pkg/cheapaml/cheapaml.F

cf5b5345a0 Jean*0001 #include "CHEAPAML_OPTIONS.h"
                0002
89f7c61169 Jean*0003       SUBROUTINE CHEAPAML(
                0004      I                           myTime, myIter, myThid )
                0005
                0006 C     ==================================================================
                0007 C     SUBROUTINE cheapaml
                0008 C     ==================================================================
                0009 C
                0010 C     o Get the surface fluxes used to force ocean model
                0011 C
                0012 C       Output:
                0013 C       ------
                0014 C       ustress, vstress - wind stress
                0015 C       Qnet             - net heat flux
                0016 C       EmPmR            - net freshwater flux
                0017 C       Tair  - mean air temperature (K)  at height ht (m)
                0018 C       Qair - Specific humidity kg/kg
                0019 C       Cheaptracer - passive tracer
                0020 C       ---------
                0021 C
                0022 C       Input:
                0023 C       ------
b4dc6cd434 Jean*0024 C       uWind, vWind  - mean wind speed (m/s)
89f7c61169 Jean*0025 C       Tr - Relaxation profile for Tair on boundaries (C)
                0026 C       qr - Relaxation profile for specific humidity (kg/kg)
                0027 C       CheaptracerR - Relaxation profile for passive tracer
                0028 C     ==================================================================
                0029 C     SUBROUTINE cheapaml
                0030 C     ==================================================================
                0031
                0032       IMPLICIT NONE
                0033
                0034 C     == global variables ==
cf5b5345a0 Jean*0035 #include "SIZE.h"
ced0783fba Jean*0036 #include "EEPARAMS.h"
c7cc66b68a Jean*0037 #include "EESUPPORT.h"
cf5b5345a0 Jean*0038 #include "PARAMS.h"
                0039 #include "DYNVARS.h"
                0040 #include "GRID.h"
                0041 #include "FFIELDS.h"
                0042 #include "CHEAPAML.h"
ced0783fba Jean*0043
89f7c61169 Jean*0044 C     == routine arguments ==
cf5b5345a0 Jean*0045       _RL     myTime
89f7c61169 Jean*0046       INTEGER myIter
b4dc6cd434 Jean*0047       INTEGER myThid
cf5b5345a0 Jean*0048
                0049 C     == Local variables ==
5251e2c855 Jean*0050       INTEGER bi,bj
6e2c553d69 Jean*0051       INTEGER i,j, nt, startAB
b4dc6cd434 Jean*0052       INTEGER iMin, iMax
                0053       INTEGER jMin, jMax
f7e8d2fb69 Jean*0054       LOGICAL writeDbug
                0055       CHARACTER*10 sufx
5251e2c855 Jean*0056       LOGICAL xIsPeriodic, yIsPeriodic
89f7c61169 Jean*0057
                0058 C tendencies of atmospheric temperature, current and past
6e2c553d69 Jean*0059         _RL gTair(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0060         _RL gqair(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0061         _RL gCheaptracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
89f7c61169 Jean*0062 C zonal and meridional transports
                0063         _RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0064         _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
cf5b5345a0 Jean*0065 C       AML timestep
ced0783fba Jean*0066         _RL deltaTTracer,deltaTm,ts,xalwu
8fd83faf35 Jean*0067         _RL dm,pt,xalwd,xlwnet
2616d73cb2 Nico*0068         _RL dtemp,xflu,xfld,dq,dtr
c7cc66b68a Jean*0069 c       _RL Fclouds, ttt2
8fd83faf35 Jean*0070         _RL q,precip,ttt,entrain
                0071         _RL uRelWind(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0072         _RL vRelWind(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0073         _RL windSq  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0074         _RL fsha(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0075         _RL flha(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0076         _RL evp (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0077         _RL xolw(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0078         _RL ssqt(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0079         _RL q100(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0080         _RL cdq (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0081 C     surfDrag   :: surface drag coeff (for wind stress)
                0082         _RL surfDrag(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0083         _RL dumArg(6)
                0084         _RL fsha0, flha0, evp_0, xolw0, ssqt0, q100_0, cdq_0
                0085         _RL Tsurf  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0086         _RL iceFrac(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0087         _RL icFrac, opFrac
2616d73cb2 Nico*0088 C temp var
fe3b1ad426 Jean*0089         _RL tmpFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
2616d73cb2 Nico*0090
89f7c61169 Jean*0091 C useful values
                0092 C inverse of time step
ced0783fba Jean*0093         deltaTm=1. _d 0/deltaT
cf5b5345a0 Jean*0094
ced0783fba Jean*0095 C atmospheric timestep
cf5b5345a0 Jean*0096         deltaTtracer = deltaT/FLOAT(cheapaml_ntim)
                0097
f7e8d2fb69 Jean*0098 c       writeDbug = debugLevel.GE.debLevC .AND.
                0099 c     &             DIFFERENT_MULTIPLE(diagFreq, myTime, deltaTClock)
                0100         writeDbug = debugLevel.GE.debLevC .AND. diagFreq.GT.0.
                0101
b4dc6cd434 Jean*0102 #ifdef ALLOW_DIAGNOSTICS
                0103 C--   fill-in diagnostics for cheapAML state variables:
                0104       IF ( useDiagnostics ) THEN
                0105        CALL DIAGNOSTICS_FILL( Tair, 'CH_TAIR ', 0,1,0,1,1, myThid )
                0106        IF (useFreshWaterFlux)
                0107      & CALL DIAGNOSTICS_FILL( Qair, 'CH_QAIR ', 0,1,0,1,1, myThid )
                0108       ENDIF
                0109 #endif /* ALLOW_DIAGNOSTICS */
                0110
cf5b5345a0 Jean*0111       DO bj=myByLo(myThid),myByHi(myThid)
                0112        DO bi=myBxLo(myThid),myBxHi(myThid)
89f7c61169 Jean*0113 C initialize net heat flux and fresh water flux arrays
                0114          DO j = 1-OLy,sNy+OLy
                0115           DO i = 1-OLx,sNx+OLx
                0116             Qnet(i,j,bi,bj) = 0. _d 0
                0117             EmPmR(i,j,bi,bj)= 0. _d 0
ced0783fba Jean*0118           ENDDO
                0119          ENDDO
89f7c61169 Jean*0120        ENDDO
                0121       ENDDO
ced0783fba Jean*0122
89f7c61169 Jean*0123 C this is a reprogramming to speed up cheapaml
                0124 C the short atmospheric time step is applied to
                0125 C advection and diffusion only.  diabatic forcing is computed
                0126 C once and used for the entire oceanic time step.
c7cc66b68a Jean*0127
89f7c61169 Jean*0128 C cycle through atmospheric advective/diffusive
                0129 C surface temperature evolution
cf5b5345a0 Jean*0130
89f7c61169 Jean*0131       DO nt=1,cheapaml_ntim
cf5b5345a0 Jean*0132
89f7c61169 Jean*0133         DO bj=myByLo(myThid),myByHi(myThid)
                0134          DO bi=myBxLo(myThid),myBxHi(myThid)
                0135
                0136 C compute advective and diffusive flux divergence
                0137           DO j=1-OLy,sNy+OLy
                0138            DO i=1-OLx,sNx+OLx
                0139              gTair(i,j,bi,bj)=0. _d 0
58426debb4 Jean*0140              uTrans(i,j) = uWind(i,j,bi,bj)*dyG(i,j,bi,bj)
                0141              vTrans(i,j) = vWind(i,j,bi,bj)*dxG(i,j,bi,bj)
89f7c61169 Jean*0142            ENDDO
                0143           ENDDO
2b5bd8961b Jean*0144           CALL CHEAPAML_CALC_RHS(
89f7c61169 Jean*0145      I           bi, bj, 1-OLx,sNx+OLx, 1-OLy,sNy+OLy,
                0146      I           uTrans, vTrans,
b4dc6cd434 Jean*0147      I           uWind, vWind,
5251e2c855 Jean*0148      I           cheapaml_kdiff, Tair,
                0149      I           deltaTtracer, zu, useFluxLimit,
                0150      I           cheapamlXperiodic, cheapamlYperiodic,
b4dc6cd434 Jean*0151      O           wWind,
5251e2c855 Jean*0152      U           gTair,
cf5b5345a0 Jean*0153      I           myTime, myIter, myThid )
451f8a1a2d Jean*0154          IF  (.NOT.useFluxLimit ) THEN
                0155            startAB = cheapTairStartAB + nt - 1
                0156            CALL ADAMS_BASHFORTH2(
fe3b1ad426 Jean*0157      I           bi, bj, 1, 1,
9d0e3cbad3 Jean*0158      U           gTair(1-OLx,1-OLy,bi,bj),
                0159      U           gTairm(1-OLx,1-OLy,bi,bj), tmpFld,
89f7c61169 Jean*0160      I           startAB, myIter, myThid )
451f8a1a2d Jean*0161          ENDIF
5251e2c855 Jean*0162          CALL CHEAPAML_TIMESTEP(
89f7c61169 Jean*0163      I           bi, bj, 1-OLx,sNx+OLx, 1-OLy,sNy+OLy, deltaTtracer,
5251e2c855 Jean*0164      I           gTair,
                0165      U           Tair,
                0166      I           nt, myIter, myThid )
89f7c61169 Jean*0167 C close bi,bj loops
                0168          ENDDO
                0169         ENDDO
                0170 C update edges
b4dc6cd434 Jean*0171         _EXCH_XY_RL(Tair,myThid)
89f7c61169 Jean*0172
                0173        IF (useFreshWaterFlux) THEN
                0174 C do water
                0175         DO bj=myByLo(myThid),myByHi(myThid)
                0176          DO bi=myBxLo(myThid),myBxHi(myThid)
                0177           DO j=1-OLy,sNy+OLy
                0178            DO i=1-OLx,sNx+OLx
                0179              gqair(i,j,bi,bj)=0. _d 0
58426debb4 Jean*0180              uTrans(i,j) = uWind(i,j,bi,bj)*dyG(i,j,bi,bj)
                0181              vTrans(i,j) = vWind(i,j,bi,bj)*dxG(i,j,bi,bj)
89f7c61169 Jean*0182            ENDDO
                0183           ENDDO
2b5bd8961b Jean*0184           CALL CHEAPAML_CALC_RHS(
89f7c61169 Jean*0185      I           bi, bj, 1-OLx,sNx+OLx, 1-OLy,sNy+OLy,
                0186      I           uTrans, vTrans,
b4dc6cd434 Jean*0187      I           uWind, vWind,
5251e2c855 Jean*0188      I           cheapaml_kdiff, qair,
                0189      I           deltaTtracer, zu, useFluxLimit,
                0190      I           cheapamlXperiodic, cheapamlYperiodic,
b4dc6cd434 Jean*0191      O           wWind,
5251e2c855 Jean*0192      U           gqair,
ced0783fba Jean*0193      I           myTime, myIter, myThid )
451f8a1a2d Jean*0194           IF  (.NOT.useFluxLimit ) THEN
                0195            startAB = cheapTairStartAB + nt - 1
                0196            CALL ADAMS_BASHFORTH2(
fe3b1ad426 Jean*0197      I           bi, bj, 1, 1,
9d0e3cbad3 Jean*0198      U           gqair(1-OLx,1-OLy,bi,bj),
                0199      U           gqairm(1-OLx,1-OLy,bi,bj), tmpFld,
89f7c61169 Jean*0200      I           startAB, myIter, myThid )
451f8a1a2d Jean*0201           ENDIF
5251e2c855 Jean*0202           CALL CHEAPAML_TIMESTEP(
89f7c61169 Jean*0203      I           bi, bj, 1-OLx,sNx+OLx, 1-OLy,sNy+OLy, deltaTtracer,
5251e2c855 Jean*0204      I           gqair,
                0205      U           qair,
                0206      I           nt, myIter, myThid )
89f7c61169 Jean*0207 C close bi, bj loops
                0208          ENDDO
                0209         ENDDO
                0210 C update edges
b4dc6cd434 Jean*0211         _EXCH_XY_RL(qair,myThid)
89f7c61169 Jean*0212        ENDIF         ! if use freshwater
51132e5783 Nico*0213
89f7c61169 Jean*0214        IF (useCheapTracer) THEN
                0215 C     do tracer
                0216         DO bj=myByLo(myThid),myByHi(myThid)
                0217          DO bi=myBxLo(myThid),myBxHi(myThid)
                0218           DO j=1-OLy,sNy+OLy
                0219            DO i=1-OLx,sNx+OLx
                0220              gCheaptracer(i,j,bi,bj)=0. _d 0
58426debb4 Jean*0221              uTrans(i,j) = uWind(i,j,bi,bj)*dyG(i,j,bi,bj)
                0222              vTrans(i,j) = vWind(i,j,bi,bj)*dxG(i,j,bi,bj)
89f7c61169 Jean*0223            ENDDO
                0224           ENDDO
2b5bd8961b Jean*0225           CALL CHEAPAML_CALC_RHS(
89f7c61169 Jean*0226      I           bi, bj, 1-OLx,sNx+OLx, 1-OLy,sNy+OLy,
                0227      I           uTrans, vTrans,
b4dc6cd434 Jean*0228      I           uWind, vWind,
5251e2c855 Jean*0229      I           cheapaml_kdiff, Cheaptracer,
                0230      I           deltaTtracer, zu, useFluxLimit,
                0231      I           cheapamlXperiodic, cheapamlYperiodic,
b4dc6cd434 Jean*0232      O           wWind,
5251e2c855 Jean*0233      U           gCheaptracer,
51132e5783 Nico*0234      I           myTime, myIter, myThid )
451f8a1a2d Jean*0235           IF  (.NOT.useFluxLimit ) THEN
                0236            startAB = cheapTracStartAB + nt - 1
                0237            CALL ADAMS_BASHFORTH2(
fe3b1ad426 Jean*0238      I           bi, bj, 1, 1,
9d0e3cbad3 Jean*0239      U           gCheaptracer(1-OLx,1-OLy,bi,bj),
                0240      U           gCheaptracerm(1-OLx,1-OLy,bi,bj), tmpFld,
fe3b1ad426 Jean*0241      I           startAB, myIter, myThid )
451f8a1a2d Jean*0242           ENDIF
5251e2c855 Jean*0243           CALL CHEAPAML_TIMESTEP(
89f7c61169 Jean*0244      I           bi, bj, 1-OLx,sNx+OLx, 1-OLy,sNy+OLy, deltaTtracer,
5251e2c855 Jean*0245      I           gCheaptracer,
                0246      U           Cheaptracer,
                0247      I           nt, myIter, myThid )
89f7c61169 Jean*0248 C     close bi, bj loops
                0249          ENDDO
                0250         ENDDO
                0251 C     update edges
b4dc6cd434 Jean*0252         _EXCH_XY_RL(Cheaptracer,myThid)
89f7c61169 Jean*0253        ENDIF                   ! if use tracer
c7cc66b68a Jean*0254
89f7c61169 Jean*0255 C reset boundaries to open boundary profile
d25d6ad15e Jean*0256        IF ( .NOT.(cheapamlXperiodic.AND.cheapamlYperiodic) ) THEN
4fa4901be6 Nico*0257         DO bj=myByLo(myThid),myByHi(myThid)
89f7c61169 Jean*0258          DO bi=myBxLo(myThid),myBxHi(myThid)
5251e2c855 Jean*0259            CALL CHEAPAML_COPY_EDGES(
                0260      I                   cheapamlXperiodic, cheapamlYperiodic,
                0261      I                   Tr(1-OLx,1-OLy,bi,bj),
                0262      U                   Tair(1-OLx,1-OLy,bi,bj),
                0263      I                   bi, bj, myIter, myThid )
                0264           IF (useFreshWaterFlux) THEN
                0265            CALL CHEAPAML_COPY_EDGES(
                0266      I                   cheapamlXperiodic, cheapamlYperiodic,
                0267      I                   qr(1-OLx,1-OLy,bi,bj),
                0268      U                   qair(1-OLx,1-OLy,bi,bj),
                0269      I                   bi, bj, myIter, myThid )
                0270           ENDIF
                0271           IF (useCheapTracer) THEN
                0272            CALL CHEAPAML_COPY_EDGES(
                0273      I                   cheapamlXperiodic, cheapamlYperiodic,
                0274      I                   CheaptracerR(1-OLx,1-OLy,bi,bj),
                0275      U                   Cheaptracer(1-OLx,1-OLy,bi,bj),
                0276      I                   bi, bj, myIter, myThid )
                0277           ENDIF
89f7c61169 Jean*0278          ENDDO
                0279         ENDDO
                0280        ENDIF
c7cc66b68a Jean*0281
89f7c61169 Jean*0282 C--   end loop on nt (short time-step loop)
                0283       ENDDO
f7e8d2fb69 Jean*0284       IF ( writeDbug ) THEN
                0285        WRITE(sufx,'(I10.10)') myIter
                0286        CALL WRITE_FLD_XY_RL('tAir_afAdv.', sufx, Tair, myIter, myThid )
                0287        IF (useFreshWaterFlux)
                0288      & CALL WRITE_FLD_XY_RL('qAir_afAdv.', sufx, qair, myIter, myThid )
                0289       ENDIF
ced0783fba Jean*0290
89f7c61169 Jean*0291 C cycling on short atmospheric time step is now done
cf5b5345a0 Jean*0292
89f7c61169 Jean*0293 C     now continue with diabatic forcing
b4dc6cd434 Jean*0294       iMin = 1
                0295       iMax = sNx
                0296       jMin = 1
                0297       jMax = sNy
                0298
51132e5783 Nico*0299       DO bj=myByLo(myThid),myByHi(myThid)
89f7c61169 Jean*0300        DO bi=myBxLo(myThid),myBxHi(myThid)
b4dc6cd434 Jean*0301
d54b0079d9 Brun*0302            DO j = 1-OLy, sNy+OLy
                0303              DO i = 1-OLx, sNx+OLx
                0304                surfDrag(i,j,bi,bj) = 0.
                0305              ENDDO
                0306            ENDDO
                0307
                0308          IF (useRelativeWind) THEN
                0309            DO j = 1-OLy, sNy+OLy
                0310              DO i = 1-OLx, sNx+OLx
                0311                uRelWind(i,j) = uWind(i,j,bi,bj)-uVel(i,j,1,bi,bj)
                0312                vRelWind(i,j) = vWind(i,j,bi,bj)-vVel(i,j,1,bi,bj)
                0313              ENDDO
                0314            ENDDO
                0315
                0316            DO j = jMin,jMax
                0317              DO i = iMin,iMax
                0318                windSq(i,j) = ( uRelWind( i ,j)*uRelWind( i ,j)
                0319      &            + uRelWind(i+1,j)*uRelWind(i+1,j)
                0320      &            + vRelWind(i, j )*vRelWind(i, j )
                0321      &            + vRelWind(i,j+1)*vRelWind(i,j+1)
                0322      &            )*0.5 _d 0
b4dc6cd434 Jean*0323 #ifdef INCONSISTENT_WIND_LOCATION
d54b0079d9 Brun*0324                windSq(i,j) = uRelWind(i,j)*uRelWind(i,j)
                0325      &            + vRelWind(i,j)*vRelWind(i,j)
b4dc6cd434 Jean*0326 #endif
d54b0079d9 Brun*0327              ENDDO
                0328            ENDDO
                0329          ELSE
                0330            DO j = jMin,jMax
                0331              DO i = iMin,iMax
                0332                windSq(i,j) = ( uWind( i ,j,bi,bj)*uWind( i ,j,bi,bj)
                0333      &            + uWind(i+1,j,bi,bj)*uWind(i+1,j,bi,bj)
                0334      &            + vWind(i, j ,bi,bj)*vWind(i, j,bi,bj )
                0335      &            + vWind(i,j+1,bi,bj)*vWind(i,j+1,bi,bj)
                0336      &            )*0.5 _d 0
                0337 #ifdef INCONSISTENT_WIND_LOCATION
                0338                windSq(i,j) = uWind(i,j,bi,bj)*uWind(i,j,bi,bj)
                0339      &            + vWind(i,j,bi,bj)*vWind(i,j,bi,bj)
                0340 #endif
                0341              ENDDO
                0342            ENDDO
                0343
08ad22b4de Jean*0344          ENDIF
d54b0079d9 Brun*0345
8fd83faf35 Jean*0346          IF ( useThSIce ) THEN
                0347            CALL CHEAPAML_SEAICE(
                0348      I                    solar(1-OLx,1-OLy,bi,bj),
                0349      I                    cheapdlongwave(1-OLx,1-OLy,bi,bj),
                0350      I                    uWind(1-OLx,1-OLy,bi,bj),
                0351      I                    vWind(1-OLx,1-OLy,bi,bj), lath,
                0352      O                    fsha, flha, evp, xolw, ssqt, q100, cdq,
                0353      O                    Tsurf, iceFrac, Qsw(1-OLx,1-OLy,bi,bj),
                0354      I                    bi, bj, myTime, myIter, myThid )
                0355            DO j = jMin,jMax
                0356             DO i = iMin,iMax
                0357               CALL CHEAPAML_COARE3_FLUX(
                0358      I                      i, j, bi, bj, 0,
                0359      I                      theta(1-OLx,1-OLy,1,bi,bj), windSq,
                0360      O                      fsha0, flha0, evp_0, xolw0,
                0361      O                      ssqt0, q100_0, cdq_0,
                0362      O                      surfDrag(i,j,bi,bj),
                0363      O                      dumArg(1), dumArg(2), dumArg(3), dumArg(4),
                0364      I                      myIter, myThid )
                0365               Qnet(i,j,bi,bj) = (
                0366      &                           -solar(i,j,bi,bj)
                0367      &                           +xolw0 - cheapdlongwave(i,j,bi,bj)
                0368      &                           +fsha0
                0369      &                           +flha0
                0370      &                          )*maskC(i,j,1,bi,bj)
                0371               EmPmR(i,j,bi,bj) = evp_0
                0372               icFrac  = iceFrac(i,j)
                0373               opFrac = 1. _d 0 - icFrac
                0374 C-     Qsw (from FFIELDS.h) has opposite (and annoying) sign convention:
                0375               Qsw(i,j,bi,bj) = - ( icFrac*Qsw(i,j,bi,bj)
                0376      &                           + opFrac*solar(i,j,bi,bj) )
                0377               fsha(i,j) = icFrac*fsha(i,j) + opFrac*fsha0
                0378               flha(i,j) = icFrac*flha(i,j) + opFrac*flha0
                0379               evp(i,j)  = icFrac*evp(i,j)  + opFrac*evp_0
                0380               xolw(i,j) = icFrac*xolw(i,j) + opFrac*xolw0
                0381               ssqt(i,j) = icFrac*ssqt(i,j) + opFrac*ssqt0
                0382               q100(i,j) = icFrac*q100(i,j) + opFrac*q100_0
                0383               cdq(i,j)  = icFrac*cdq(i,j)  + opFrac*cdq_0
                0384             ENDDO
                0385            ENDDO
                0386          ELSE
                0387            DO j = jMin,jMax
                0388             DO i = iMin,iMax
                0389              IF (FluxFormula.EQ.'LANL') THEN
                0390               CALL cheapaml_LANL_flux(
                0391      I                      i, j, bi, bj,
                0392      O                      fsha(i,j), flha(i,j), evp(i,j),
                0393      O                      xolw(i,j), ssqt(i,j), q100(i,j) )
                0394              ELSEIF (FluxFormula.EQ.'COARE3') THEN
                0395               CALL CHEAPAML_COARE3_FLUX(
                0396      I                      i, j, bi, bj, 0,
                0397      I                      theta(1-OLx,1-OLy,1,bi,bj), windSq,
                0398      O                      fsha(i,j), flha(i,j), evp(i,j), xolw(i,j),
                0399      O                      ssqt(i,j), q100(i,j), cdq(i,j),
                0400      O                      surfDrag(i,j,bi,bj),
                0401      O                      dumArg(1), dumArg(2), dumArg(3), dumArg(4),
                0402      I                      myIter, myThid )
                0403              ENDIF
                0404              IF (useFreshWaterFlux) THEN
                0405               EmPmR(i,j,bi,bj) = evp(i,j)
                0406              ENDIF
                0407             ENDDO
                0408            ENDDO
                0409          ENDIF
                0410
                0411          DO j = jMin,jMax
                0412           DO i = iMin,iMax
c7cc66b68a Jean*0413
89f7c61169 Jean*0414 C atmospheric upwelled long wave
8fd83faf35 Jean*0415            ttt = Tair(i,j,bi,bj)-gamma_blk*(CheapHgrid(i,j,bi,bj)-zt)
                0416 c          xalwu = stefan*(ttt+celsius2K)**4*0.5 _d 0
                0417            xalwu = stefan*(0.5*Tair(i,j,bi,bj)+0.5*ttt+celsius2K)**4
                0418      &                   *0.5 _d 0
89f7c61169 Jean*0419 C atmospheric downwelled long wave
8fd83faf35 Jean*0420            xalwd = stefan*(Tair(i,j,bi,bj)+celsius2K)**4*0.5 _d 0
89f7c61169 Jean*0421 C total flux at upper atmospheric layer interface
8fd83faf35 Jean*0422            xflu = ( -solar(i,j,bi,bj) + xalwu + flha(i,j)
                0423      &            )*xef*maskC(i,j,1,bi,bj)
89f7c61169 Jean*0424 C lower flux calculation.
8fd83faf35 Jean*0425            xfld = ( -solar(i,j,bi,bj) - xalwd + xolw(i,j)
                0426      &              + fsha(i,j) + flha(i,j)
                0427      &            )*xef*maskC(i,j,1,bi,bj)
0c111b8b6e Nico*0428
b4dc6cd434 Jean*0429            IF (useDLongWave) THEN
8fd83faf35 Jean*0430              xlwnet = xolw(i,j)-cheapdlongwave(i,j,bi,bj)
b4dc6cd434 Jean*0431            ELSE
0c111b8b6e Nico*0432 C net long wave (see Josey et al. JGR 1997)
                0433 C coef lambda replaced by 0.5+lat/230
                0434 C convert spec humidity in water vapor pressure (mbar) using coef 1000/0.622=1607.7
b4dc6cd434 Jean*0435              xlwnet = 0.98 _d 0*stefan*(theta(i,j,1,bi,bj)+celsius2K)**4
cf6b9ab292 Brun*0436      &          *(0.39 _d 0 - 0.05 _d 0*SQRT(ABS(qair(i,j,bi,bj))
                0437      &          *     1607.7 _d 0))
8fd83faf35 Jean*0438      &        *( oneRL - (halfRL+ABS(yG(i,j,bi,bj))/230. _d 0)
                0439      &                   *cheapclouds(i,j,bi,bj)**2 )
                0440      &        + 4.0*0.98 _d 0*stefan*(theta(i,j,1,bi,bj)+celsius2K)**3
                0441      &          *(theta(i,j,1,bi,bj)-Tair(i,j,bi,bj))
89f7c61169 Jean*0442
b4dc6cd434 Jean*0443 c            xlwnet = xolw-stefan*(theta(i,j,1,bi,bj)+celsius2K)**4.
89f7c61169 Jean*0444 c     &       *(0.65+11.22*qair(i,j,bi,bj) + 0.25*cheapclouds(i,j,bi,bj)
                0445 c     &       -8.23*qair(i,j,bi,bj)*cheapclouds(i,j,bi,bj))
b4dc6cd434 Jean*0446            ENDIF
2616d73cb2 Nico*0447 C clouds
b4dc6cd434 Jean*0448 c          ttt2=Tair(i,j,bi,bj)-1.5*gamma_blk*CheapHgrid(i,j,bi,bj)
                0449 c          Fclouds = stefan*ttt2**4*(0.4*cheapclouds(i,j,bi,bj)+1-0.4)/2
                0450 c          ttt2=Tair(i,j,bi,bj)-3*gamma_blk*CheapHgrid(i,j,bi,bj)+celsius2K
                0451 c          Fclouds = 0.3*stefan*ttt2**4 + 0.22*xolw*cheapclouds(i,j,bi,bj)
89f7c61169 Jean*0452 C add flux divergences into atmospheric temperature tendency
b4dc6cd434 Jean*0453            gTair(i,j,bi,bj)= (xfld-xflu)/CheapHgrid(i,j,bi,bj)
8fd83faf35 Jean*0454            IF ( .NOT.useThSIce ) THEN
                0455             Qnet(i,j,bi,bj) = (
b4dc6cd434 Jean*0456      &                         -solar(i,j,bi,bj)
                0457 c    &                         -xalwd
                0458 c    &                         -Fclouds
                0459 c    &                         +xolw
                0460      &                         +xlwnet
8fd83faf35 Jean*0461      &                         +fsha(i,j)
                0462      &                         +flha(i,j)
b4dc6cd434 Jean*0463      &                        )*maskC(i,j,1,bi,bj)
8fd83faf35 Jean*0464              Qsw(i,j,bi,bj) = -solar(i,j,bi,bj)
                0465            ENDIF
2616d73cb2 Nico*0466
89f7c61169 Jean*0467 C need to precip?
b4dc6cd434 Jean*0468            IF (useFreshWaterFlux) THEN
cf6b9ab292 Brun*0469              q=q100(i,j)
89f7c61169 Jean*0470 C compute saturation specific humidity at atmospheric
                0471 C layer top
                0472 C first, what is the pressure there?
                0473 C ts is surface atmospheric temperature
b4dc6cd434 Jean*0474             ts=Tair(i,j,bi,bj)+gamma_blk*zt+celsius2K
                0475             pt=p0*(1-gamma_blk*CheapHgrid(i,j,bi,bj)/ts)
cf6b9ab292 Brun*0476      &          **(gravity/gamma_blk/gasR)
89f7c61169 Jean*0477
cf6b9ab292 Brun*0478              IF (.NOT.usePrecip) THEN
89f7c61169 Jean*0479 C factor to compute rainfall from specific humidity
cf6b9ab292 Brun*0480               dm=100.*(p0-pt)*recip_gravity
0c111b8b6e Nico*0481 C     Large scale precip
cf6b9ab292 Brun*0482               precip = 0.
                0483               IF ( wWind(i,j,bi,bj).GT.0. .AND.
                0484      &             q.GT.ssqt(i,j)*0.7 _d 0 ) THEN
                0485                 precip = precip
8fd83faf35 Jean*0486      &               + ( (q-ssqt(i,j)*0.7 _d 0)*dm/cheap_pr2 )
                0487      &                 *(wWind(i,j,bi,bj)/0.75 _d -5)**2
cf6b9ab292 Brun*0488               ENDIF
2616d73cb2 Nico*0489
0c111b8b6e Nico*0490 C     Convective precip
cf6b9ab292 Brun*0491               IF (q.GT.0.0214 _d 0 .AND. q.GT.ssqt(i,j)*0.9 _d 0) THEN
                0492                 precip = precip + ((q-ssqt(i,j)*0.9 _d 0)*dm/cheap_pr1)
                0493               ENDIF
0c111b8b6e Nico*0494
cf6b9ab292 Brun*0495               cheapPrecip(i,j,bi,bj) = precip*1200/CheapHgrid(i,j,bi,bj)
                0496             ENDIF
08ad22b4de Jean*0497
cf6b9ab292 Brun*0498               entrain = cdq(i,j)*q*0.25
0c111b8b6e Nico*0499
b4dc6cd434 Jean*0500 c           gqair(i,j,bi,bj)=(evp-precip-entrain)/CheapHgrid(i,j,bi,bj)
8fd83faf35 Jean*0501             gqair(i,j,bi,bj) = (evp(i,j)-entrain)/CheapHgrid(i,j,bi,bj)
                0502      &                        /rhoa*maskC(i,j,1,bi,bj)
                0503             EmPmR(i,j,bi,bj) = ( EmPmR(i,j,bi,bj)
                0504      &                          -cheapPrecip(i,j,bi,bj)
b4dc6cd434 Jean*0505      &                         )*maskC(i,j,1,bi,bj)
                0506            ENDIF
ced0783fba Jean*0507
89f7c61169 Jean*0508           ENDDO
                0509          ENDDO
                0510
                0511 C it is not necessary to use the Adams2d subroutine as
                0512 C the forcing is always computed at the current time step.
5251e2c855 Jean*0513 C note: full oceanic time step deltaT is used below
                0514          CALL CHEAPAML_TIMESTEP(
b4dc6cd434 Jean*0515      I           bi, bj, iMin, iMax, jMin, jMax, deltaT,
5251e2c855 Jean*0516      I           gTair,
                0517      U           Tair,
                0518      I           0, myIter, myThid )
89f7c61169 Jean*0519 C       do implicit time stepping over land
                0520          DO j=1-OLy,sNy+OLy
                0521           DO i=1-OLx,sNx+OLx
                0522             dtemp=tr(i,j,bi,bj)-Tair(i,j,bi,bj)
                0523             Tair(i,j,bi,bj)=Tair(i,j,bi,bj)+dtemp*xrelf(i,j,bi,bj)
                0524           ENDDO
                0525          ENDDO
51132e5783 Nico*0526
89f7c61169 Jean*0527 C do water
                0528         IF (useFreshWaterFlux) THEN
5251e2c855 Jean*0529          CALL CHEAPAML_TIMESTEP(
b4dc6cd434 Jean*0530      I           bi, bj, iMin, iMax, jMin, jMax, deltaT,
5251e2c855 Jean*0531      I           gqair,
                0532      U           qair,
                0533      I           0, myIter, myThid )
89f7c61169 Jean*0534 C     do implicit time stepping over land and or buffer
                0535          DO j=1-OLy,sNy+OLy
                0536           DO i=1-OLx,sNx+OLx
                0537             dq=qr(i,j,bi,bj)-qair(i,j,bi,bj)
                0538             qair(i,j,bi,bj)=qair(i,j,bi,bj)+dq*xrelf(i,j,bi,bj)
                0539             IF (qair(i,j,bi,bj).LT.0.0) qair(i,j,bi,bj) = 0.0 _d 0
                0540           ENDDO
                0541          ENDDO
                0542         ENDIF
                0543
                0544 C do tracer
                0545         IF (useCheapTracer) THEN
                0546 C     do implicit time stepping over land and or buffer
                0547          DO j=1-OLy,sNy+OLy
                0548           DO i=1-OLx,sNx+OLx
                0549             dtr=CheaptracerR(i,j,bi,bj)-Cheaptracer(i,j,bi,bj)
                0550             Cheaptracer(i,j,bi,bj) = Cheaptracer(i,j,bi,bj)
                0551      &                             + dtr*xrelf(i,j,bi,bj)
                0552           ENDDO
                0553          ENDDO
                0554         ENDIF
ced0783fba Jean*0555
8fd83faf35 Jean*0556 #ifdef ALLOW_DIAGNOSTICS
                0557         IF ( useDiagnostics ) THEN
                0558          CALL DIAGNOSTICS_FILL( fsha,'CH_SH   ',0,1,2,bi,bj,myThid )
                0559          CALL DIAGNOSTICS_FILL( flha,'CH_LH   ',0,1,2,bi,bj,myThid )
                0560          CALL DIAGNOSTICS_FILL( q100,'CH_q100 ',0,1,2,bi,bj,myThid )
                0561          CALL DIAGNOSTICS_FILL( ssqt,'CH_ssqt ',0,1,2,bi,bj,myThid )
                0562         ENDIF
                0563 #endif /* ALLOW_DIAGNOSTICS */
                0564
89f7c61169 Jean*0565 C close bi,bj loops
                0566        ENDDO
                0567       ENDDO
                0568
                0569 C update edges
b4dc6cd434 Jean*0570        _EXCH_XY_RL(Tair,myThid)
                0571        _EXCH_XY_RS(Qnet,myThid)
89f7c61169 Jean*0572       IF (useFreshWaterFlux) THEN
b4dc6cd434 Jean*0573        _EXCH_XY_RL(qair,myThid)
                0574        _EXCH_XY_RS(EmPmR,myThid)
89f7c61169 Jean*0575       ENDIF
                0576       IF (useCheapTracer) THEN
b4dc6cd434 Jean*0577        _EXCH_XY_RL(Cheaptracer,myThid)
                0578       ENDIF
                0579       IF (.NOT.useStressOption) THEN
                0580        IF ( FluxFormula.EQ.'COARE3' ) THEN
                0581         _EXCH_XY_RL( surfDrag, myThid )
                0582        ELSE
                0583         CALL EXCH_UV_AGRID_3D_RL( ustress, vstress, .TRUE., 1, myThid )
                0584        ENDIF
89f7c61169 Jean*0585       ENDIF
                0586
                0587 C reset edges to open boundary profiles
d25d6ad15e Jean*0588 c     IF ( .NOT.(cheapamlXperiodic.AND.cheapamlYperiodic) ) THEN
5251e2c855 Jean*0589       IF ( notUsingXPeriodicity.OR.notUsingYPeriodicity ) THEN
                0590         xIsPeriodic = .NOT.notUsingXPeriodicity
                0591         yIsPeriodic = .NOT.notUsingYPeriodicity
                0592         DO bj=myByLo(myThid),myByHi(myThid)
                0593          DO bi=myBxLo(myThid),myBxHi(myThid)
                0594            CALL CHEAPAML_COPY_EDGES(
                0595 c    I                   cheapamlXperiodic, cheapamlYperiodic,
                0596      I                   xIsPeriodic, yIsPeriodic,
                0597      I                   Tr(1-OLx,1-OLy,bi,bj),
                0598      U                   Tair(1-OLx,1-OLy,bi,bj),
                0599      I                   bi, bj, myIter, myThid )
                0600           IF (useFreshWaterFlux) THEN
                0601            CALL CHEAPAML_COPY_EDGES(
                0602 c    I                   cheapamlXperiodic, cheapamlYperiodic,
                0603      I                   xIsPeriodic, yIsPeriodic,
                0604      I                   qr(1-OLx,1-OLy,bi,bj),
                0605      U                   qair(1-OLx,1-OLy,bi,bj),
                0606      I                   bi, bj, myIter, myThid )
                0607           ENDIF
                0608           IF (useCheapTracer) THEN
                0609            CALL CHEAPAML_COPY_EDGES(
                0610 c    I                   cheapamlXperiodic, cheapamlYperiodic,
                0611      I                   xIsPeriodic, yIsPeriodic,
                0612      I                   CheaptracerR(1-OLx,1-OLy,bi,bj),
                0613      U                   Cheaptracer(1-OLx,1-OLy,bi,bj),
                0614      I                   bi, bj, myIter, myThid )
                0615           ENDIF
89f7c61169 Jean*0616          ENDDO
51132e5783 Nico*0617         ENDDO
89f7c61169 Jean*0618       ENDIF
c7cc66b68a Jean*0619
8fd83faf35 Jean*0620 c     CALL PLOT_FIELD_XYRS( gTair, 'S/R CHEAPAML gTair',1,myThid)
                0621 c     CALL PLOT_FIELD_XYRS( Tair, 'S/R CHEAPAML Tair',1,myThid)
                0622 c     CALL PLOT_FIELD_XYRS( Qnet, 'S/R CHEAPAML Qnet',1,myThid)
cf5b5345a0 Jean*0623
                0624       DO bj=myByLo(myThid),myByHi(myThid)
                0625        DO bi=myBxLo(myThid),myBxHi(myThid)
b4dc6cd434 Jean*0626
                0627         IF ( .NOT.useStressOption .AND. FluxFormula.EQ.'COARE3' ) THEN
d54b0079d9 Brun*0628           IF (useRelativeWind) THEN
08ad22b4de Jean*0629
d54b0079d9 Brun*0630           DO j = 1-OLy+1,sNy+OLy
b4dc6cd434 Jean*0631           DO i = 1-OLx+1,sNx+OLx
                0632             fu(i,j,bi,bj) = maskW(i,j,1,bi,bj)*0.5 _d 0
                0633      &          *( surfDrag(i-1,j,bi,bj) + surfDrag(i,j,bi,bj) )
                0634      &          *( uWind(i,j,bi,bj)-uVel(i,j,1,bi,bj) )
                0635             fv(i,j,bi,bj) = maskS(i,j,1,bi,bj)*0.5 _d 0
                0636      &          *( surfDrag(i,j-1,bi,bj) + surfDrag(i,j,bi,bj) )
                0637      &          *( vWind(i,j,bi,bj)-vVel(i,j,1,bi,bj) )
                0638           ENDDO
                0639          ENDDO
                0640 #ifdef INCONSISTENT_WIND_LOCATION
                0641          DO j = 1-OLy,sNy+OLy
                0642           DO i = 1-OLx+1,sNx+OLx
                0643             fu(i,j,bi,bj) = maskW(i,j,1,bi,bj)*0.5 _d 0
                0644      &          *( surfDrag(i-1,j,bi,bj)
                0645      &             *( uWind(i-1,j,bi,bj)-uVel(i-1,j,1,bi,bj) )
                0646      &           + surfDrag(i,j,bi,bj)
                0647      &             *( uWind(i,j,bi,bj) - uVel(i,j,1,bi,bj) ) )
                0648           ENDDO
                0649          ENDDO
                0650          DO j = 1-OLy+1,sNy+OLy
                0651           DO i = 1-OLx,sNx+OLx
                0652             fv(i,j,bi,bj) = maskS(i,j,1,bi,bj)*0.5 _d 0
                0653      &          *( surfDrag(i,j-1,bi,bj)
                0654      &             *( vWind(i,j-1,bi,bj)-vVel(i,j-1,1,bi,bj) )
                0655      &           + surfDrag(i,j,bi,bj)
                0656      &             *( vWind(i,j,bi,bj) - vVel(i,j,1,bi,bj) ) )
                0657           ENDDO
                0658          ENDDO
                0659 #endif /* INCONSISTENT_WIND_LOCATION */
d54b0079d9 Brun*0660        ELSE ! relative wind
                0661
                0662           DO j = 1-OLy+1,sNy+OLy
                0663           DO i = 1-OLx+1,sNx+OLx
                0664             fu(i,j,bi,bj) = maskW(i,j,1,bi,bj)*0.5 _d 0
                0665      &          *( surfDrag(i-1,j,bi,bj) + surfDrag(i,j,bi,bj) )
                0666      &          *uWind(i,j,bi,bj)
                0667             fv(i,j,bi,bj) = maskS(i,j,1,bi,bj)*0.5 _d 0
                0668      &          *( surfDrag(i,j-1,bi,bj) + surfDrag(i,j,bi,bj) )
                0669      &          *vWind(i,j,bi,bj)
                0670           ENDDO
                0671          ENDDO
                0672 #ifdef INCONSISTENT_WIND_LOCATION
                0673          DO j = 1-OLy,sNy+OLy
                0674           DO i = 1-OLx+1,sNx+OLx
                0675             fu(i,j,bi,bj) = maskW(i,j,1,bi,bj)*0.5 _d 0
                0676      &          *( surfDrag(i-1,j,bi,bj)
                0677      &             *uWind(i-1,j,bi,bj)
                0678      &           + surfDrag(i,j,bi,bj)
                0679      &             *uWind(i,j,bi,bj))
                0680           ENDDO
                0681          ENDDO
                0682          DO j = 1-OLy+1,sNy+OLy
                0683           DO i = 1-OLx,sNx+OLx
                0684             fv(i,j,bi,bj) = maskS(i,j,1,bi,bj)*0.5 _d 0
                0685      &          *( surfDrag(i,j-1,bi,bj)
                0686      &             *vWind(i,j-1,bi,bj)
                0687      &           + surfDrag(i,j,bi,bj)
                0688      &             *vWind(i,j,bi,bj))
                0689           ENDDO
                0690          ENDDO
                0691 #endif /* INCONSISTENT_WIND_LOCATION */
                0692          ENDIF !relative wind
                0693
                0694        ELSE
89f7c61169 Jean*0695 Cswd move wind stresses to u and v points
                0696          DO j = 1-OLy,sNy+OLy
                0697           DO i = 1-OLx+1,sNx+OLx
b4dc6cd434 Jean*0698             fu(i,j,bi,bj) = maskW(i,j,1,bi,bj)*0.5 _d 0
                0699      &          *( ustress(i,j,bi,bj) + ustress(i-1,j,bi,bj) )
89f7c61169 Jean*0700           ENDDO
                0701          ENDDO
                0702          DO j = 1-OLy+1,sNy+OLy
                0703           DO i = 1-OLx,sNx+OLx
b4dc6cd434 Jean*0704             fv(i,j,bi,bj) = maskS(i,j,1,bi,bj)*0.5 _d 0
                0705      &          *( vstress(i,j,bi,bj) + vstress(i,j-1,bi,bj) )
89f7c61169 Jean*0706           ENDDO
                0707          ENDDO
b4dc6cd434 Jean*0708         ENDIF
cf5b5345a0 Jean*0709
                0710 C--   end bi,bj loops
                0711        ENDDO
                0712       ENDDO
2616d73cb2 Nico*0713
4fa4901be6 Nico*0714 #ifdef ALLOW_DIAGNOSTICS
bf944b1865 Jean*0715 C- note: with thSIce, CH_QNET and CH_EmP correspond to the fluxes over the
                0716 C     open-ocean fraction of the grid-cell (similar to EXFqnet and EXFempmr).
                0717 C     Use instead diagnostics SIflxAtm & SIfrwAtm to get the grid-cell
                0718 C      averaged (open-ocean fraction + ice-covered fraction).
b4dc6cd434 Jean*0719       IF ( useDiagnostics ) THEN
                0720        CALL DIAGNOSTICS_FILL(uWind,  'CH_Uwind',0,1,0,1,1,myThid)
                0721        CALL DIAGNOSTICS_FILL(vWind,  'CH_Vwind',0,1,0,1,1,myThid)
                0722        CALL DIAGNOSTICS_FILL_RS(Qnet,'CH_QNET ',0,1,0,1,1,myThid)
                0723        IF (useFreshWaterFlux) THEN
                0724         CALL DIAGNOSTICS_FILL_RS( EmPmR, 'CH_EmP  ', 0,1,0,1,1,myThid)
                0725         CALL DIAGNOSTICS_FILL(cheapPrecip,'CH_Prec ',0,1,0,1,1,myThid)
                0726        ENDIF
                0727        IF (useCheapTracer) THEN
89f7c61169 Jean*0728         CALL DIAGNOSTICS_FILL(Cheaptracer,'CH_Trace',0,1,0,1,1,myThid)
b4dc6cd434 Jean*0729        ENDIF
51132e5783 Nico*0730       ENDIF
4fa4901be6 Nico*0731 #endif /* ALLOW_DIAGNOSTICS */
c7cc66b68a Jean*0732
89f7c61169 Jean*0733 c     DO bj=myByLo(myThid),myByHi(myThid)
                0734 c      DO bi=myBxLo(myThid),myBxHi(myThid)
                0735 c        DO j = 1-OLy,sNy+OLy
                0736 c         DO i = 1-OLx+1,sNx+OLx
                0737 c           fu(i,j,bi,bj) = 0.0
                0738 c           fv(i,j,bi,bj) = 0.0
                0739 c           Qnet(i,j,bi,bj) = 0.0
                0740 c           EmPmR(i,j,bi,bj) = 0.0
                0741 c         ENDDO
                0742 c        ENDDO
                0743 c      ENDDO
                0744 c     ENDDO
cf5b5345a0 Jean*0745
                0746       RETURN
                0747       END