tutorial_global_oce_latlon/code/ptracers_forcing_surf.F

4307c62f79 Jean*0001 #include "PTRACERS_OPTIONS.h"
                0002
                0003 CBOP
                0004 C !ROUTINE: PTRACERS_FORCING_SURF
                0005
                0006 C !INTERFACE: ==========================================================
                0007       SUBROUTINE PTRACERS_FORCING_SURF(
3f4989e08f Jean*0008      I                            relaxForcingS,
4307c62f79 Jean*0009      I                            bi, bj, iMin, iMax, jMin, jMax,
                0010      I                            myTime,myIter,myThid )
                0011
                0012 C !DESCRIPTION:
                0013 C     Precomputes surface forcing term for pkg/ptracers.
                0014 C     Precomputation is needed because of non-local KPP transport term,
                0015 C     routine KPP_TRANSPORT_PTR.
                0016
                0017 C !USES: ===============================================================
                0018       IMPLICIT NONE
                0019 #include "SIZE.h"
                0020 #include "EEPARAMS.h"
                0021 #include "PARAMS.h"
                0022 #include "GRID.h"
                0023 #include "SURFACE.h"
3f4989e08f Jean*0024 #include "DYNVARS.h"
4307c62f79 Jean*0025 #include "FFIELDS.h"
                0026 #include "PTRACERS_SIZE.h"
22e0cff85e Jean*0027 #include "PTRACERS_PARAMS.h"
d6215f7b79 Jean*0028 #include "PTRACERS_START.h"
22e0cff85e Jean*0029 #include "PTRACERS_FIELDS.h"
4307c62f79 Jean*0030
                0031 C !INPUT PARAMETERS: ===================================================
3f4989e08f Jean*0032 C  relaxForcingS        :: Salt forcing due to surface relaxation
4307c62f79 Jean*0033 C  bi,bj                :: tile indices
                0034 C  myTime               :: model time
                0035 C  myIter               :: time-step number
                0036 C  myThid               :: thread number
3f4989e08f Jean*0037       _RL relaxForcingS(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
4307c62f79 Jean*0038       INTEGER bi, bj, iMin, iMax, jMin, jMax
                0039       _RL myTime
                0040       INTEGER myIter
                0041       INTEGER myThid
                0042
                0043 #ifdef ALLOW_PTRACERS
                0044
                0045 C !LOCAL VARIABLES: ====================================================
                0046 C  i,j                  :: loop indices
                0047 C  iTrc                 :: tracer index
                0048 C  ks                   :: surface level index
                0049       INTEGER i, j
                0050       INTEGER iTrc, ks
3f4989e08f Jean*0051       _RL add2EmP(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0052       _RL epsil, cutoff, tmpVar
4307c62f79 Jean*0053 CEOP
                0054
                0055       IF ( usingPCoords ) THEN
                0056         ks = Nr
                0057       ELSE
                0058         ks = 1
                0059       ENDIF
                0060
                0061 C Example of how to add forcing at the surface
                0062       DO iTrc=1,PTRACERS_numInUse
d6215f7b79 Jean*0063 c       IF ( PTRACERS_StepFwd(iTrc) ) THEN
4307c62f79 Jean*0064           DO j = jMin, jMax
                0065            DO i = iMin, iMax
d6215f7b79 Jean*0066 C-    Case of age-tracer: at the surface, add 10.days relaxation towards zero:
22e0cff85e Jean*0067              surfaceForcingPTr(i,j,bi,bj,iTrc) =
62fd6ae4e5 Jean*0068      &        + 1. _d 0 / (10. _d 0 * 86400. _d 0)
                0069      &                  * ( 0. _d 0 - pTracer(i,j,ks,bi,bj,iTrc) )
97f6d163ab Davi*0070      &                  * drF(ks) * _hFacC(i,j,ks,bi,bj)
4307c62f79 Jean*0071            ENDDO
                0072           ENDDO
d6215f7b79 Jean*0073 c       ENDIF
4307c62f79 Jean*0074       ENDDO
                0075
3f4989e08f Jean*0076 C--   Option to convert Salt-relaxation into additional EmP contribution
                0077       IF ( PTRACERS_addSrelax2EmP ) THEN
                0078 C-    here we assume that salt_EvPrRn = 0
                0079 C     set cutoff value to prevent too large additional EmP:
                0080 C       current limit is set to 0.1 CFL
                0081         epsil = 1. _d -10
                0082         cutoff = 0.1 _d 0 *drF(ks)/PTRACERS_dTLev(ks)
                0083         IF ( ( (nonlinFreeSurf.GT.0 .OR. usingPCoords)
                0084      &         .AND. useRealFreshWaterFlux )
                0085      &     .OR.convertFW2Salt .EQ. -1. ) THEN
                0086          DO j = jMin, jMax
                0087           DO i = iMin, iMax
                0088             tmpVar = MAX( salt(i,j,ks,bi,bj), epsil )
                0089             add2EmP(i,j) = relaxForcingS(i,j)/tmpVar
                0090             add2EmP(i,j) = rUnit2mass
                0091      &                  *MAX( -cutoff, MIN( add2EmP(i,j), cutoff ) )
                0092           ENDDO
                0093          ENDDO
                0094         ELSE
                0095          DO j = jMin, jMax
                0096           DO i = iMin, iMax
                0097             add2EmP(i,j) = relaxForcingS(i,j)/convertFW2Salt
                0098             add2EmP(i,j) = rUnit2mass
                0099      &                  *MAX( -cutoff, MIN( add2EmP(i,j), cutoff ) )
                0100           ENDDO
                0101          ENDDO
                0102         ENDIF
                0103 #ifdef ALLOW_DIAGNOSTICS
                0104         IF ( useDiagnostics ) THEN
                0105          CALL DIAGNOSTICS_FILL(add2EmP,'Add2EmP ',0,1,2,bi,bj,myThid)
                0106         ENDIF
                0107 #endif /* ALLOW_DIAGNOSTICS */
                0108       ELSE
                0109         DO j = jMin, jMax
                0110           DO i = iMin, iMax
                0111             add2EmP(i,j) = 0. _d 0
                0112           ENDDO
                0113         ENDDO
                0114       ENDIF
                0115 C-- end of "addEmP" setting
                0116
4307c62f79 Jean*0117 #ifdef EXACT_CONSERV
                0118       IF ( (nonlinFreeSurf.GT.0 .OR. usingPCoords)
                0119      &     .AND. useRealFreshWaterFlux ) THEN
                0120
                0121        DO iTrc=1,PTRACERS_numInUse
                0122
                0123 c-  NonLin_FrSurf and RealFreshWaterFlux : PmEpR effectively changes
                0124 c   the water column height ; temp., salt, (tracer) flux associated
                0125 c   with this input/output of water is added here to the surface tendency.
                0126 c
d6215f7b79 Jean*0127          IF ( PTRACERS_StepFwd(iTrc) .AND.
                0128      &        PTRACERS_EvPrRn(iTrc).NE.UNSET_RL ) THEN
4307c62f79 Jean*0129           DO j = jMin, jMax
                0130            DO i = iMin, iMax
22e0cff85e Jean*0131              surfaceForcingPTr(i,j,bi,bj,iTrc) =
                0132      &          surfaceForcingPTr(i,j,bi,bj,iTrc)
3f4989e08f Jean*0133      &        + ( PmEpR(i,j,bi,bj) - add2EmP(i,j) )
4307c62f79 Jean*0134      &          *( PTRACERS_EvPrRn(iTrc) - pTracer(i,j,ks,bi,bj,iTrc) )
62fd6ae4e5 Jean*0135      &          *mass2rUnit
4307c62f79 Jean*0136            ENDDO
                0137           ENDDO
                0138          ENDIF
                0139
                0140        ENDDO
a09a74749d Davi*0141
                0142 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0143       ELSE
                0144 #else /* EXACT_CONSERV */
                0145       IF (.TRUE.) THEN
4307c62f79 Jean*0146 #endif /* EXACT_CONSERV */
                0147
a09a74749d Davi*0148 C--   EmPmR does not really affect the water column height (for tracer budget)
                0149 C     and is converted to a salt tendency.
                0150
                0151        IF (convertFW2Salt .EQ. -1.) THEN
                0152 C-    use local surface tracer field to calculate forcing term:
                0153
                0154         DO iTrc=1,PTRACERS_numInUse
                0155
d6215f7b79 Jean*0156          IF ( PTRACERS_StepFwd(iTrc) .AND.
                0157      &        PTRACERS_EvPrRn(iTrc).NE.UNSET_RL ) THEN
22e0cff85e Jean*0158 C        account for Rain/Evap tracer content (PTRACERS_EvPrRn) using
a09a74749d Davi*0159 C        local surface tracer
                0160           DO j = jMin, jMax
                0161            DO i = iMin, iMax
22e0cff85e Jean*0162             surfaceForcingPTr(i,j,bi,bj,iTrc) =
                0163      &          surfaceForcingPTr(i,j,bi,bj,iTrc)
3f4989e08f Jean*0164      &        + ( EmPmR(i,j,bi,bj) + add2EmP(i,j) )
22e0cff85e Jean*0165      &          *( pTracer(i,j,ks,bi,bj,iTrc) - PTRACERS_EvPrRn(iTrc) )
a09a74749d Davi*0166      &          *mass2rUnit
                0167            ENDDO
                0168           ENDDO
                0169          ENDIF
                0170
                0171         ENDDO
                0172
                0173        ELSE
                0174 C-    use uniform tracer value to calculate forcing term:
                0175
                0176         DO iTrc=1,PTRACERS_numInUse
                0177
d6215f7b79 Jean*0178          IF ( PTRACERS_StepFwd(iTrc) .AND.
                0179      &        PTRACERS_EvPrRn(iTrc).NE.UNSET_RL ) THEN
22e0cff85e Jean*0180 C     account for Rain/Evap tracer content (PTRACERS_EvPrRn) assuming uniform
a09a74749d Davi*0181 C     surface tracer (=PTRACERS_ref)
                0182           DO j = jMin, jMax
                0183            DO i = iMin, iMax
22e0cff85e Jean*0184             surfaceForcingPTr(i,j,bi,bj,iTrc) =
                0185      &          surfaceForcingPTr(i,j,bi,bj,iTrc)
3f4989e08f Jean*0186      &        + ( EmPmR(i,j,bi,bj) + add2EmP(i,j) )
22e0cff85e Jean*0187      &            *( PTRACERS_ref(ks,iTrc) - PTRACERS_EvPrRn(iTrc) )
a09a74749d Davi*0188      &            *mass2rUnit
                0189            ENDDO
                0190           ENDDO
                0191          ENDIF
                0192
                0193         ENDDO
                0194
                0195 C-    end local-surface-tracer / uniform-value distinction
                0196        ENDIF
                0197
                0198       ENDIF
                0199
                0200 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0201
4307c62f79 Jean*0202 #endif /* ALLOW_PTRACERS */
                0203
                0204       RETURN
                0205       END