C $Header: /home/ubuntu/mnt/e9_copy/MITgcm_contrib/lab_sea_test/my_min_max_abs.F,v 1.1 2004/07/12 01:00:20 dimitri Exp $
C $Name:  $

#include "CPP_OPTIONS.h"

C Simple versions of Martins modified min/max/abs functions

      _RS function mymax_R4( a, b )
      implicit none
      _RS a, b
      _RS mymin_R4
      external mymin_R4
      mymax_R4 = -mymin_R4(-a, -b)
      return
      end

      _RS function mymax_R8( a, b )
      implicit none
      _RL a, b
      _RL mymin_R8
      external mymin_R8
      mymax_R8 = -mymin_R8(-a, -b)
      print*,'###max',a,b,mymax_R8
      return
      end

      _RS function mymin_R4( a, b )

      implicit none

      _RS a, b

      _RS myabs_R4
      external myabs_R4

Cml      mymin_R4 = .5*(a+b)
      mymin_R4 = .5*( a+b - myabs_R4(a-b) )
Cml      mymin_R4 = MIN(a,b)

Cml      if ( a .lt. b ) then
Cml         mymin_R4 = a
Cml      else if ( a .gt. b ) then
Cml         mymin_R4 = b
Cml      else
Cml         mymin_R4 = .5*(a+b)
Cml      end if

      return
      end

      _RL function mymin_R8( a, b )

      implicit none

      _RL a, b

      _RL myabs_R8
      external myabs_R8

Cml      mymin_R8 = .5*(a+b)
      mymin_R8 = .5*( a+b - myabs_R8(a-b) )
Cml      mymin_R8 = MIN(a,b)

Cml      if ( a .lt. b ) then
Cml         mymin_R8 = a
Cml      else if ( a .gt. b ) then
Cml         mymin_R8 = b
Cml      else
Cml         mymin_R8 = .5*(a+b)
Cml      end if
      print*,'###min',a,b,mymin_R8

      return
      end

      _RS function myabs_R4( x )
      
      implicit none
C     === Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
C     input parameter
      _RS x
c     local variable
      _RS sf, rsf
      _RS smoothAbsFuncRange

c     set default value. THhs will go in ini_parms
      smoothAbsFuncRange = 0.

      if ( smoothAbsFuncRange .lt. 0.0 ) then
c     limit of mymin(a,b) = .5*(a+b)
         myabs_R4 = 0.
      else
         if ( smoothAbsFuncRange .ne. 0.0 ) then
            sf  = 10.0/smoothAbsFuncRange
            rsf = 1./sf
         else
c     limit of mymin(a,b) = min(a,b)
            sf  = 0.
            rsf = 0.
         end if
c
         if ( x .gt. smoothAbsFuncRange ) then
            myabs_R4 = x
         else if ( x .lt. -smoothAbsFuncRange ) then
            myabs_R4 = -x
         else
            myabs_R4 = log(.5*(exp(x*sf)+exp(-x*sf)))*rsf
         end if
      end if

      return
      end 

      _RL function myabs_R8( x )
      
      implicit none
C     === Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
C     input parameter
      _RL x
c     local variable
      _RL sf, rsf
      _RL smoothAbsFuncRange

c     set default value. THhs will go in ini_parms
      smoothAbsFuncRange = 1.

      if ( smoothAbsFuncRange .lt. 0.0 ) then
c     limit of mymin(a,b) = .5*(a+b)
         myabs_R8 = 0.
      else
         if ( smoothAbsFuncRange .ne. 0.0 ) then
            sf  = 10.0D0/smoothAbsFuncRange
            rsf = 1.D0/sf
         else
c     limit of mymin(a,b) = min(a,b)
            sf  = 0.D0
            rsf = 0.D0
         end if
c     
         if ( x .ge. smoothAbsFuncRange ) then
            myabs_R8 = x
         else if ( x .le. -smoothAbsFuncRange ) then
            myabs_R8 = -x
         else
            myabs_R8 = log(.5*(exp(x*sf)+exp(-x*sf)))*rsf
         end if
      end if

      return
      end