C $Header: /u/gcmpack/MITgcm/model/src/update_masks_etc.F,v 1.11 2017/04/04 23:22:38 jmc Exp $
C $Name: $
#include "PACKAGES_CONFIG.h"
#include "CPP_OPTIONS.h"
#ifdef ALLOW_AUTODIFF
# include "AUTODIFF_OPTIONS.h"
#endif
C-- File update_masks_etc.F:
C-- Contents
C-- o S/R UPDATE_MASKS_ETC
C-- o FCT SMOOTHMIN_RS( a, b )
C-- o FCT SMOOTHMIN_RL( a, b )
C-- o FCT SMOOTHABS_RS( x )
C-- o FCT SMOOTHABS_RL( x )
Cml o S/R LIMIT_HFACC_TO_ONE
Cml o S/R ADLIMIT_HFACC_TO_ONE
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
CBOP
C !ROUTINE: UPDATE_MASKS_ETC
C !INTERFACE:
SUBROUTINE UPDATE_MASKS_ETC( myThid )
C !DESCRIPTION: \bv
C *==========================================================*
C | SUBROUTINE UPDATE_MASKS_ETC
C | o Re-initialise masks and topography factors after a new
C | hFacC has been calculated by the minimizer
C *==========================================================*
C | These arrays are used throughout the code and describe
C | the topography of the domain through masks (0s and 1s)
C | and fractional height factors (0
C | distinguish between the lopped-cell and full-step
C | topographic representations.
C *==========================================================*
C | code taken from ini_masks_etc.F
C *==========================================================*
C \ev
C !USES:
IMPLICIT NONE
C === Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "SURFACE.h"
Cml we need optimcycle for storing the new hFaC(C/W/S) and depth
#ifdef ALLOW_AUTODIFF
# include "optim.h"
#endif
C !INPUT/OUTPUT PARAMETERS:
C == Routine arguments ==
C myThid - Number of this instance of INI_MASKS_ETC
INTEGER myThid
#ifdef ALLOW_DEPTH_CONTROL
C !FUNCTIONS:
_RS SMOOTHMIN_RS
EXTERNAL
C !LOCAL VARIABLES:
C == Local variables ==
C bi,bj :: Loop counters
C I,J,K
C tmpfld :: Temporary array used to compute & write Total Depth
INTEGER bi, bj
INTEGER I, J, K
_RS tmpfld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
CHARACTER*(MAX_LEN_MBUF) suff
Cml(
INTEGER Im1, Jm1
_RL hFacCtmp, hFacCtmp2
_RL hFacMnSz
Cml)
CEOP
C- Calculate lopping factor hFacC : over-estimate the part inside of the domain
C taking into account the lower_R Boundary (Bathymetrie / Top of Atmos)
DO bj=myByLo(myThid), myByHi(myThid)
DO bi=myBxLo(myThid), myBxHi(myThid)
DO K=1, Nr
hFacMnSz=max( hFacMin, min(hFacMinDr*recip_drF(k),1. _d 0) )
DO J=1-OLy,sNy+OLy
DO I=1-OLx,sNx+OLx
C o Non-dimensional distance between grid bound. and domain lower_R bound.
#ifdef ALLOW_DEPTH_CONTROL
hFacCtmp = (rF(K)-xx_r_low(I,J,bi,bj))*recip_drF(K)
#else
hFacCtmp = (rF(K)-R_low(I,J,bi,bj))*recip_drF(K)
#endif /* ALLOW_DEPTH_CONTROL */
Cml IF ( hFacCtmp .LE. 0. _d 0 ) THEN
CmlC IF ( hFacCtmp .LT. 0.5*hfacMnSz ) THEN
Cml hFacCtmp2 = 0. _d 0
Cml ELSE
Cml hFacCtmp2 = hFacCtmp + hFacMnSz*(
Cml & EXP(-hFacCtmp/hFacMnSz)-EXP(-1./hFacMnSz) )
Cml ENDIF
Cml CALL limit_hfacc_to_one( hFacCtmp2 )
Cml hFacC(I,J,K,bi,bj) = hFacCtmp2
IF ( hFacCtmp .LE. 0. _d 0 ) THEN
C IF ( hFacCtmp .LT. 0.5*hfacMnSz ) THEN
hFacC(I,J,K,bi,bj) = 0. _d 0
ELSEIF ( hFacCtmp .GT. 1. _d 0 ) THEN
hFacC(I,J,K,bi,bj) = 1. _d 0
ELSE
hFacC(I,J,K,bi,bj) = hFacCtmp + hFacMnSz*(
& EXP(-hFacCtmp/hFacMnSz)-EXP(-1./hFacMnSz) )
ENDIF
Cml print '(A,3I5,F20.16)', 'ml-hfac:', I,J,K,hFacC(I,J,K,bi,bj)
CmlC o Select between, closed, open or partial (0,1,0-1)
Cml hFacCtmp=min( max( hFacCtmp, 0. _d 0) , 1. _d 0)
CmlC o Impose minimum fraction and/or size (dimensional)
Cml IF (hFacCtmp.LT.hFacMnSz) THEN
Cml IF (hFacCtmp.LT.hFacMnSz*0.5) THEN
Cml hFacC(I,J,K,bi,bj)=0.
Cml ELSE
Cml hFacC(I,J,K,bi,bj)=hFacMnSz
Cml ENDIF
Cml ELSE
Cml hFacC(I,J,K,bi,bj)=hFacCtmp
Cml ENDIF
Cml ENDIF
Cml print '(A,F15.4,F20.16)', 'ml-hfac:', R_low(i,j,bi,bj),hFacC(I,J,K,bi,bj)
ENDDO
ENDDO
ENDDO
C - end bi,bj loops.
ENDDO
ENDDO
C _EXCH_XYZ_RS(hFacC,myThid)
C- Re-calculate lower-R Boundary position, taking into account hFacC
DO bj=myByLo(myThid), myByHi(myThid)
DO bi=myBxLo(myThid), myBxHi(myThid)
DO J=1-OLy,sNy+OLy
DO I=1-OLx,sNx+OLx
R_low(i,j,bi,bj) = rF(1)
ENDDO
ENDDO
DO K=Nr,1,-1
DO J=1-OLy,sNy+OLy
DO I=1-OLx,sNx+OLx
R_low(I,J,bi,bj) = R_low(I,J,bi,bj)
& - drF(K)*hFacC(I,J,K,bi,bj)
ENDDO
ENDDO
ENDDO
C - end bi,bj loops.
ENDDO
ENDDO
Cml DO bj=myByLo(myThid), myByHi(myThid)
Cml DO bi=myBxLo(myThid), myBxHi(myThid)
CmlC- Re-calculate Reference surface position, taking into account hFacC
Cml DO J=1-OLy,sNy+OLy
Cml DO I=1-OLx,sNx+OLx
Cml Ro_surf(I,J,bi,bj) = R_low(I,J,bi,bj)
Cml DO K=Nr,1,-1
Cml Ro_surf(I,J,bi,bj) = Ro_surf(I,J,bi,bj)
Cml & + drF(k)*hFacC(I,J,K,bi,bj)
Cml ENDDO
Cml ENDDO
Cml ENDDO
CmlC - end bi,bj loops.
Cml ENDDO
Cml ENDDO
IF ( plotLevel.GE.debLevC ) THEN
_BARRIER
CALL PLOT_FIELD_XYRS( R_low,
& 'Model R_low (update_masks_etc)', 1, myThid )
CML I assume that Ro_surf is not changed anywhere else in the code
CML and since it is not changed in this routine, we do not need to
CML print it again.
CML CALL PLOT_FIELD_XYRS( Ro_surf,
CML & 'Model Ro_surf (update_masks_etc)', 1, myThid )
ENDIF
C Calculate quantities derived from XY depth map
DO bj = myByLo(myThid), myByHi(myThid)
DO bi = myBxLo(myThid), myBxHi(myThid)
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
C Total fluid column thickness (r_unit) :
tmpfld(i,j,bi,bj) = Ro_surf(i,j,bi,bj) - R_low(i,j,bi,bj)
C Inverse of fluid column thickness (1/r_unit)
IF ( tmpfld(i,j,bi,bj) .LE. 0. ) THEN
recip_Rcol(i,j,bi,bj) = 0.
ELSE
recip_Rcol(i,j,bi,bj) = 1. _d 0 / tmpfld(i,j,bi,bj)
ENDIF
ENDDO
ENDDO
ENDDO
ENDDO
C _EXCH_XY_RS( recip_Rcol, myThid )
C hFacW and hFacS (at U and V points)
CML This will be the crucial part of the code, because here the minimum
CML function MIN is involved which does not have a continuous derivative
CML for MIN(x,y) at y=x.
CML The thin walls representation has been moved into this loop, that is
CML before the call to EXCH_UV_XVY_RS, because TAMC will prefer it this
CML way. On the other hand, this might cause difficulties in some
CML configurations.
DO bj=myByLo(myThid), myByHi(myThid)
DO bi=myBxLo(myThid), myBxHi(myThid)
DO K=1, Nr
DO J=1-OLy,sNy+OLy
DO I=1-OLx,sNx+OLx
Im1=MAX(I-1,1-OLx)
Jm1=MAX(J-1,1-OLy)
IF (DYG(I,J,bi,bj).EQ.0.) THEN
C thin walls representation of non-periodic
C boundaries such as happen on the lat-lon grid at the N/S poles.
C We should really supply a flag for doing this.
hFacW(I,J,K,bi,bj)=0.
ELSE
hFacW(I,J,K,bi,bj)=maskW(I,J,K,bi,bj)*
#ifdef USE_SMOOTH_MIN
& SMOOTHMIN_RS(hFacC(I,J,K,bi,bj),hFacC(Im1,J,K,bi,bj))
#else
& MIN(hFacC(I,J,K,bi,bj),hFacC(Im1,J,K,bi,bj))
#endif /* USE_SMOOTH_MIN */
ENDIF
IF (DXG(I,J,bi,bj).EQ.0.) THEN
hFacS(I,J,K,bi,bj)=0.
ELSE
hFacS(I,J,K,bi,bj)=maskS(I,J,K,bi,bj)*
#ifdef USE_SMOOTH_MIN
& SMOOTHMIN_RS(hFacC(I,J,K,bi,bj),hFacC(I,Jm1,K,bi,bj))
#else
& MIN(hFacC(I,J,K,bi,bj),hFacC(I,Jm1,K,bi,bj))
#endif /* USE_SMOOTH_MIN */
ENDIF
ENDDO
ENDDO
ENDDO
ENDDO
ENDDO
#if ( defined (ALLOW_AUTODIFF)
defined (ALLOW_AUTODIFF_MONITOR)
defined (ALLOW_DEPTH_CONTROL) )
C Include call to a dummy routine. Its adjoint will be called at the proper
C place in the adjoint code. The adjoint routine will print out adjoint
C values if requested. The location of the call is important, it has to be
C after the adjoint of the exchanges (DO_GTERM_BLOCKING_EXCHANGES).
Cml CALL DUMMY_IN_HFAC( 'W', 0, myThid )
Cml CALL DUMMY_IN_HFAC( 'S', 0, myThid )
#endif
CALL EXCH_UV_XYZ_RS(hFacW,hFacS,.FALSE.,myThid)
#if ( defined (ALLOW_AUTODIFF)
defined (ALLOW_AUTODIFF_MONITOR)
defined (ALLOW_DEPTH_CONTROL) )
C Include call to a dummy routine. Its adjoint will be called at the proper
C place in the adjoint code. The adjoint routine will print out adjoint
C values if requested. The location of the call is important, it has to be
C after the adjoint of the exchanges (DO_GTERM_BLOCKING_EXCHANGES).
Cml CALL DUMMY_IN_HFAC( 'W', 1, myThid )
Cml CALL DUMMY_IN_HFAC( 'S', 1, myThid )
#endif
C- Write to disk: Total Column Thickness & hFac(C,W,S):
WRITE(suff,'(I10.10)') optimcycle
CALL WRITE_FLD_XY_RS( 'Depth.',suff,tmpfld,optimcycle,myThid)
CALL WRITE_FLD_XYZ_RS( 'hFacC.',suff,hFacC,optimcycle,myThid)
CALL WRITE_FLD_XYZ_RS( 'hFacW.',suff,hFacW,optimcycle,myThid)
CALL WRITE_FLD_XYZ_RS( 'hFacS.',suff,hFacS,optimcycle,myThid)
IF ( plotLevel.GE.debLevC ) THEN
_BARRIER
C-- Write to monitor file (standard output)
CALL PLOT_FIELD_XYZRS( hFacC,'hFacC (update_masks_etc)',
& Nr, 1, myThid )
CALL PLOT_FIELD_XYZRS( hFacW,'hFacW (update_masks_etc)',
& Nr, 1, myThid )
CALL PLOT_FIELD_XYZRS( hFacS,'hFacS (update_masks_etc)',
& Nr, 1, myThid )
ENDIF
C Masks and reciprocals of hFac[CWS]
Cml The masks should stay constant, so they are not recomputed at this time
Cml implicitly implying that no cell that is wet in the begin will ever dry
Cml up! This is a strong constraint and should be implementent as a hard
Cml inequality contraint when performing optimization (m1qn3 cannot do that)
Cml Also, I am assuming here that the new hFac(s) never become zero during
Cml optimization!
DO bj = myByLo(myThid), myByHi(myThid)
DO bi = myBxLo(myThid), myBxHi(myThid)
DO K=1,Nr
DO J=1-OLy,sNy+OLy
DO I=1-OLx,sNx+OLx
IF (hFacC(I,J,K,bi,bj) .NE. 0. ) THEN
Cml IF (maskC(I,J,K,bi,bj) .NE. 0. ) THEN
recip_hFacC(I,J,K,bi,bj) = 1. _d 0 / hFacC(I,J,K,bi,bj)
Cml maskC(I,J,K,bi,bj) = 1.
ELSE
recip_hFacC(I,J,K,bi,bj) = 0.
Cml maskC(I,J,K,bi,bj) = 0.
ENDIF
IF (hFacW(I,J,K,bi,bj) .NE. 0. ) THEN
Cml IF (maskW(I,J,K,bi,bj) .NE. 0. ) THEN
recip_hFacW(I,J,K,bi,bj) = 1. _d 0 / hFacw(I,J,K,bi,bj)
Cml maskW(I,J,K,bi,bj) = 1.
ELSE
recip_hFacW(I,J,K,bi,bj) = 0.
Cml maskW(I,J,K,bi,bj) = 0.
ENDIF
IF (hFacS(I,J,K,bi,bj) .NE. 0. ) THEN
Cml IF (maskS(I,J,K,bi,bj) .NE. 0. ) THEN
recip_hFacS(I,J,K,bi,bj) = 1. _d 0 / hFacS(I,J,K,bi,bj)
Cml maskS(I,J,K,bi,bj) = 1.
ELSE
recip_hFacS(I,J,K,bi,bj) = 0.
Cml maskS(I,J,K,bi,bj) = 0.
ENDIF
ENDDO
ENDDO
ENDDO
CmlCml(
Cml ENDDO
Cml ENDDO
Cml _EXCH_XYZ_RS(recip_hFacC , myThid )
Cml _EXCH_XYZ_RS(recip_hFacW , myThid )
Cml _EXCH_XYZ_RS(recip_hFacS , myThid )
Cml _EXCH_XYZ_RS(maskC , myThid )
Cml _EXCH_XYZ_RS(maskW , myThid )
Cml _EXCH_XYZ_RS(maskS , myThid )
Cml DO bj = myByLo(myThid), myByHi(myThid)
Cml DO bi = myBxLo(myThid), myBxHi(myThid)
CmlCml)
#ifdef NONLIN_FRSURF
C-- Save initial geometrical hFac factor into h0Fac (fixed in time):
C Note: In case 1 pkg modifies hFac (from packages_init_fixed, called
C later in sequence of calls) this pkg would need also to update h0Fac.
DO k=1,Nr
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
h0FacC(i,j,k,bi,bj) = _hFacC(i,j,k,bi,bj)
h0FacW(i,j,k,bi,bj) = _hFacW(i,j,k,bi,bj)
h0FacS(i,j,k,bi,bj) = _hFacS(i,j,k,bi,bj)
ENDDO
ENDDO
ENDDO
#endif /* NONLIN_FRSURF */
C - end bi,bj loops.
ENDDO
ENDDO
#endif /* ALLOW_DEPTH_CONTROL */
RETURN
END
#ifdef USE_SMOOTH_MIN
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
_RS FUNCTION SMOOTHMIN_RS( a, b )
IMPLICIT NONE
_RS a, b
_RS SMOOTHABS_RS
EXTERNAL
Cml smoothMin_R4 = .5*(a+b)
SMOOTHMIN_RS = .5*( a+b - SMOOTHABS_RS(a-b) )
CML smoothMin_R4 = MIN(a,b)
RETURN
END
_RL FUNCTION SMOOTHMIN_RL( a, b )
IMPLICIT NONE
_RL a, b
_RL SMOOTHABS_RL
EXTERNAL
Cml smoothMin_R8 = .5*(a+b)
SMOOTHMIN_RL = .5*( a+b - SMOOTHABS_RL(a-b) )
Cml smoothMin_R8 = MIN(a,b)
RETURN
END
_RS FUNCTION SMOOTHABS_RS( 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
IF ( smoothAbsFuncRange .LT. 0.0 ) THEN
c limit of smoothMin(a,b) = .5*(a+b)
SMOOTHABS_RS = 0.
ELSE
IF ( smoothAbsFuncRange .NE. 0.0 ) THEN
sf = 10.0/smoothAbsFuncRange
rsf = 1./sf
ELSE
c limit of smoothMin(a,b) = min(a,b)
sf = 0.
rsf = 0.
ENDIF
c
IF ( x .GT. smoothAbsFuncRange ) THEN
SMOOTHABS_RS = x
ELSEIF ( x .LT. -smoothAbsFuncRange ) THEN
SMOOTHABS_RS = -x
ELSE
SMOOTHABS_RS = log(.5*(exp(x*sf)+exp(-x*sf)))*rsf
ENDIF
ENDIF
RETURN
END
_RL FUNCTION SMOOTHABS_RL( 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
IF ( smoothAbsFuncRange .LT. 0.0 ) THEN
c limit of smoothMin(a,b) = .5*(a+b)
SMOOTHABS_RL = 0.
ELSE
IF ( smoothAbsFuncRange .NE. 0.0 ) THEN
sf = 10.0D0/smoothAbsFuncRange
rsf = 1.D0/sf
ELSE
c limit of smoothMin(a,b) = min(a,b)
sf = 0.D0
rsf = 0.D0
ENDIF
c
IF ( x .GE. smoothAbsFuncRange ) THEN
SMOOTHABS_RL = x
ELSEIF ( x .LE. -smoothAbsFuncRange ) THEN
SMOOTHABS_RL = -x
ELSE
SMOOTHABS_RL = log(.5*(exp(x*sf)+exp(-x*sf)))*rsf
ENDIF
ENDIF
RETURN
END
#endif /* USE_SMOOTH_MIN */
Cml#ifdef ALLOW_DEPTH_CONTROL
Cmlcadj SUBROUTINE limit_hfacc_to_one INPUT = 1
Cmlcadj SUBROUTINE limit_hfacc_to_one OUTPUT = 1
Cmlcadj SUBROUTINE limit_hfacc_to_one ACTIVE = 1
Cmlcadj SUBROUTINE limit_hfacc_to_one DEPEND = 1
Cmlcadj SUBROUTINE limit_hfacc_to_one REQUIRED
Cmlcadj SUBROUTINE limit_hfacc_to_one ADNAME = adlimit_hfacc_to_one
Cml#endif /* ALLOW_DEPTH_CONTROL */
Cml SUBROUTINE LIMIT_HFACC_TO_ONE( hf )
Cml
Cml _RL hf
Cml
Cml IF ( hf .GT. 1. _d 0 ) THEN
Cml hf = 1. _d 0
Cml ENDIF
Cml
Cml RETURN
Cml END
Cml
Cml SUBROUTINE ADLIMIT_HFACC_TO_ONE( hf, adhf )
Cml
Cml _RL hf, adhf
Cml
Cml RETURN
Cml END
#ifdef ALLOW_DEPTH_CONTROL
cadj SUBROUTINE dummy_in_hfac INPUT = 1, 2, 3
cadj SUBROUTINE dummy_in_hfac OUTPUT =
cadj SUBROUTINE dummy_in_hfac ACTIVE =
cadj SUBROUTINE dummy_in_hfac DEPEND = 1, 2, 3
cadj SUBROUTINE dummy_in_hfac REQUIRED
cadj SUBROUTINE dummy_in_hfac INFLUENCED
cadj SUBROUTINE dummy_in_hfac ADNAME = addummy_in_hfac
cadj SUBROUTINE dummy_in_hfac FTLNAME = g_dummy_in_hfac
#endif /* ALLOW_DEPTH_CONTROL */