pkg/cheapaml/cheapaml.F

cf5b5345a0 Jean*0001 #include "CHEAPAML_OPTIONS.h"
                0002
89f7c61169 Jean*0003       SUBROUTINE CHEAPAML(
d861cd501f Jean*0004      I                     myTime, myIter, myThid )
89f7c61169 Jean*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
d861cd501f Jean*0071         _RL uRelWind(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0072         _RL vRelWind(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
8fd83faf35 Jean*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
d861cd501f Jean*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
d54b0079d9 Brun*0307
                0308          IF (useRelativeWind) THEN
                0309            DO j = 1-OLy, sNy+OLy
d861cd501f Jean*0310             DO i = 1-OLx, sNx+OLx
                0311               uRelWind(i,j,bi,bj) = uWind(i,j,bi,bj)-uVel(i,j,1,bi,bj)
                0312               vRelWind(i,j,bi,bj) = vWind(i,j,bi,bj)-vVel(i,j,1,bi,bj)
                0313             ENDDO
d54b0079d9 Brun*0314            ENDDO
                0315          ELSE
d861cd501f Jean*0316            DO j = 1-OLy, sNy+OLy
                0317             DO i = 1-OLx, sNx+OLx
                0318               uRelWind(i,j,bi,bj) = uWind(i,j,bi,bj)
                0319               vRelWind(i,j,bi,bj) = vWind(i,j,bi,bj)
                0320             ENDDO
d54b0079d9 Brun*0321            ENDDO
08ad22b4de Jean*0322          ENDIF
d861cd501f Jean*0323          DO j = jMin,jMax
                0324           DO i = iMin,iMax
                0325             windSq(i,j) = ( uRelWind( i ,j,bi,bj)*uRelWind( i ,j,bi,bj)
                0326      &                    + uRelWind(i+1,j,bi,bj)*uRelWind(i+1,j,bi,bj)
                0327      &                    + vRelWind(i, j ,bi,bj)*vRelWind(i, j ,bi,bj)
                0328      &                    + vRelWind(i,j+1,bi,bj)*vRelWind(i,j+1,bi,bj)
                0329      &                    )*halfRL
                0330 #ifdef INCONSISTENT_WIND_LOCATION
                0331             windSq(i,j) = uRelWind(i,j,bi,bj)*uRelWind(i,j,bi,bj)
                0332      &                  + vRelWind(i,j,bi,bj)*vRelWind(i,j,bi,bj)
                0333 #endif
                0334           ENDDO
                0335          ENDDO
d54b0079d9 Brun*0336
8fd83faf35 Jean*0337          IF ( useThSIce ) THEN
                0338            CALL CHEAPAML_SEAICE(
                0339      I                    solar(1-OLx,1-OLy,bi,bj),
                0340      I                    cheapdlongwave(1-OLx,1-OLy,bi,bj),
                0341      I                    uWind(1-OLx,1-OLy,bi,bj),
                0342      I                    vWind(1-OLx,1-OLy,bi,bj), lath,
                0343      O                    fsha, flha, evp, xolw, ssqt, q100, cdq,
                0344      O                    Tsurf, iceFrac, Qsw(1-OLx,1-OLy,bi,bj),
                0345      I                    bi, bj, myTime, myIter, myThid )
                0346            DO j = jMin,jMax
                0347             DO i = iMin,iMax
                0348               CALL CHEAPAML_COARE3_FLUX(
                0349      I                      i, j, bi, bj, 0,
                0350      I                      theta(1-OLx,1-OLy,1,bi,bj), windSq,
                0351      O                      fsha0, flha0, evp_0, xolw0,
                0352      O                      ssqt0, q100_0, cdq_0,
                0353      O                      surfDrag(i,j,bi,bj),
                0354      O                      dumArg(1), dumArg(2), dumArg(3), dumArg(4),
                0355      I                      myIter, myThid )
                0356               Qnet(i,j,bi,bj) = (
                0357      &                           -solar(i,j,bi,bj)
                0358      &                           +xolw0 - cheapdlongwave(i,j,bi,bj)
                0359      &                           +fsha0
                0360      &                           +flha0
                0361      &                          )*maskC(i,j,1,bi,bj)
                0362               EmPmR(i,j,bi,bj) = evp_0
                0363               icFrac  = iceFrac(i,j)
                0364               opFrac = 1. _d 0 - icFrac
                0365 C-     Qsw (from FFIELDS.h) has opposite (and annoying) sign convention:
                0366               Qsw(i,j,bi,bj) = - ( icFrac*Qsw(i,j,bi,bj)
                0367      &                           + opFrac*solar(i,j,bi,bj) )
                0368               fsha(i,j) = icFrac*fsha(i,j) + opFrac*fsha0
                0369               flha(i,j) = icFrac*flha(i,j) + opFrac*flha0
                0370               evp(i,j)  = icFrac*evp(i,j)  + opFrac*evp_0
                0371               xolw(i,j) = icFrac*xolw(i,j) + opFrac*xolw0
                0372               ssqt(i,j) = icFrac*ssqt(i,j) + opFrac*ssqt0
                0373               q100(i,j) = icFrac*q100(i,j) + opFrac*q100_0
                0374               cdq(i,j)  = icFrac*cdq(i,j)  + opFrac*cdq_0
                0375             ENDDO
                0376            ENDDO
                0377          ELSE
                0378            DO j = jMin,jMax
                0379             DO i = iMin,iMax
                0380              IF (FluxFormula.EQ.'LANL') THEN
d861cd501f Jean*0381               CALL CHEAPAML_LANL_FLUX(
                0382      I                      i, j, bi, bj, 0,
                0383      I                      theta(1-OLx,1-OLy,1,bi,bj), windSq,
                0384      O                      fsha(i,j), flha(i,j), evp(i,j), xolw(i,j),
                0385      O                      ssqt(i,j), q100(i,j), cdq(i,j),
                0386      O                      surfDrag(i,j,bi,bj),
                0387 c    O                      dumArg(1), dumArg(2), dumArg(3), dumArg(4),
                0388      I                      myIter, myThid )
8fd83faf35 Jean*0389              ELSEIF (FluxFormula.EQ.'COARE3') THEN
                0390               CALL CHEAPAML_COARE3_FLUX(
                0391      I                      i, j, bi, bj, 0,
                0392      I                      theta(1-OLx,1-OLy,1,bi,bj), windSq,
                0393      O                      fsha(i,j), flha(i,j), evp(i,j), xolw(i,j),
                0394      O                      ssqt(i,j), q100(i,j), cdq(i,j),
                0395      O                      surfDrag(i,j,bi,bj),
                0396      O                      dumArg(1), dumArg(2), dumArg(3), dumArg(4),
                0397      I                      myIter, myThid )
                0398              ENDIF
                0399              IF (useFreshWaterFlux) THEN
                0400               EmPmR(i,j,bi,bj) = evp(i,j)
                0401              ENDIF
                0402             ENDDO
                0403            ENDDO
                0404          ENDIF
                0405
                0406          DO j = jMin,jMax
                0407           DO i = iMin,iMax
c7cc66b68a Jean*0408
89f7c61169 Jean*0409 C atmospheric upwelled long wave
8fd83faf35 Jean*0410            ttt = Tair(i,j,bi,bj)-gamma_blk*(CheapHgrid(i,j,bi,bj)-zt)
                0411 c          xalwu = stefan*(ttt+celsius2K)**4*0.5 _d 0
                0412            xalwu = stefan*(0.5*Tair(i,j,bi,bj)+0.5*ttt+celsius2K)**4
d861cd501f Jean*0413      &                   *halfRL
89f7c61169 Jean*0414 C atmospheric downwelled long wave
d861cd501f Jean*0415            xalwd = stefan*(Tair(i,j,bi,bj)+celsius2K)**4*halfRL
89f7c61169 Jean*0416 C total flux at upper atmospheric layer interface
8fd83faf35 Jean*0417            xflu = ( -solar(i,j,bi,bj) + xalwu + flha(i,j)
                0418      &            )*xef*maskC(i,j,1,bi,bj)
89f7c61169 Jean*0419 C lower flux calculation.
8fd83faf35 Jean*0420            xfld = ( -solar(i,j,bi,bj) - xalwd + xolw(i,j)
                0421      &              + fsha(i,j) + flha(i,j)
                0422      &            )*xef*maskC(i,j,1,bi,bj)
0c111b8b6e Nico*0423
b4dc6cd434 Jean*0424            IF (useDLongWave) THEN
8fd83faf35 Jean*0425              xlwnet = xolw(i,j)-cheapdlongwave(i,j,bi,bj)
b4dc6cd434 Jean*0426            ELSE
0c111b8b6e Nico*0427 C net long wave (see Josey et al. JGR 1997)
                0428 C coef lambda replaced by 0.5+lat/230
                0429 C convert spec humidity in water vapor pressure (mbar) using coef 1000/0.622=1607.7
b4dc6cd434 Jean*0430              xlwnet = 0.98 _d 0*stefan*(theta(i,j,1,bi,bj)+celsius2K)**4
cf6b9ab292 Brun*0431      &          *(0.39 _d 0 - 0.05 _d 0*SQRT(ABS(qair(i,j,bi,bj))
                0432      &          *     1607.7 _d 0))
8fd83faf35 Jean*0433      &        *( oneRL - (halfRL+ABS(yG(i,j,bi,bj))/230. _d 0)
                0434      &                   *cheapclouds(i,j,bi,bj)**2 )
                0435      &        + 4.0*0.98 _d 0*stefan*(theta(i,j,1,bi,bj)+celsius2K)**3
                0436      &          *(theta(i,j,1,bi,bj)-Tair(i,j,bi,bj))
89f7c61169 Jean*0437
b4dc6cd434 Jean*0438 c            xlwnet = xolw-stefan*(theta(i,j,1,bi,bj)+celsius2K)**4.
89f7c61169 Jean*0439 c     &       *(0.65+11.22*qair(i,j,bi,bj) + 0.25*cheapclouds(i,j,bi,bj)
                0440 c     &       -8.23*qair(i,j,bi,bj)*cheapclouds(i,j,bi,bj))
b4dc6cd434 Jean*0441            ENDIF
2616d73cb2 Nico*0442 C clouds
b4dc6cd434 Jean*0443 c          ttt2=Tair(i,j,bi,bj)-1.5*gamma_blk*CheapHgrid(i,j,bi,bj)
                0444 c          Fclouds = stefan*ttt2**4*(0.4*cheapclouds(i,j,bi,bj)+1-0.4)/2
                0445 c          ttt2=Tair(i,j,bi,bj)-3*gamma_blk*CheapHgrid(i,j,bi,bj)+celsius2K
                0446 c          Fclouds = 0.3*stefan*ttt2**4 + 0.22*xolw*cheapclouds(i,j,bi,bj)
89f7c61169 Jean*0447 C add flux divergences into atmospheric temperature tendency
b4dc6cd434 Jean*0448            gTair(i,j,bi,bj)= (xfld-xflu)/CheapHgrid(i,j,bi,bj)
8fd83faf35 Jean*0449            IF ( .NOT.useThSIce ) THEN
                0450             Qnet(i,j,bi,bj) = (
b4dc6cd434 Jean*0451      &                         -solar(i,j,bi,bj)
                0452 c    &                         -xalwd
                0453 c    &                         -Fclouds
                0454 c    &                         +xolw
                0455      &                         +xlwnet
8fd83faf35 Jean*0456      &                         +fsha(i,j)
                0457      &                         +flha(i,j)
b4dc6cd434 Jean*0458      &                        )*maskC(i,j,1,bi,bj)
8fd83faf35 Jean*0459              Qsw(i,j,bi,bj) = -solar(i,j,bi,bj)
                0460            ENDIF
2616d73cb2 Nico*0461
89f7c61169 Jean*0462 C need to precip?
b4dc6cd434 Jean*0463            IF (useFreshWaterFlux) THEN
cf6b9ab292 Brun*0464              q=q100(i,j)
89f7c61169 Jean*0465 C compute saturation specific humidity at atmospheric
                0466 C layer top
                0467 C first, what is the pressure there?
                0468 C ts is surface atmospheric temperature
b4dc6cd434 Jean*0469             ts=Tair(i,j,bi,bj)+gamma_blk*zt+celsius2K
                0470             pt=p0*(1-gamma_blk*CheapHgrid(i,j,bi,bj)/ts)
cf6b9ab292 Brun*0471      &          **(gravity/gamma_blk/gasR)
89f7c61169 Jean*0472
d861cd501f Jean*0473             IF (.NOT.usePrecip) THEN
89f7c61169 Jean*0474 C factor to compute rainfall from specific humidity
cf6b9ab292 Brun*0475               dm=100.*(p0-pt)*recip_gravity
0c111b8b6e Nico*0476 C     Large scale precip
cf6b9ab292 Brun*0477               precip = 0.
                0478               IF ( wWind(i,j,bi,bj).GT.0. .AND.
                0479      &             q.GT.ssqt(i,j)*0.7 _d 0 ) THEN
                0480                 precip = precip
8fd83faf35 Jean*0481      &               + ( (q-ssqt(i,j)*0.7 _d 0)*dm/cheap_pr2 )
                0482      &                 *(wWind(i,j,bi,bj)/0.75 _d -5)**2
cf6b9ab292 Brun*0483               ENDIF
2616d73cb2 Nico*0484
0c111b8b6e Nico*0485 C     Convective precip
cf6b9ab292 Brun*0486               IF (q.GT.0.0214 _d 0 .AND. q.GT.ssqt(i,j)*0.9 _d 0) THEN
                0487                 precip = precip + ((q-ssqt(i,j)*0.9 _d 0)*dm/cheap_pr1)
                0488               ENDIF
0c111b8b6e Nico*0489
cf6b9ab292 Brun*0490               cheapPrecip(i,j,bi,bj) = precip*1200/CheapHgrid(i,j,bi,bj)
                0491             ENDIF
08ad22b4de Jean*0492
d861cd501f Jean*0493             entrain = cdq(i,j)*q*0.25 _d 0
0c111b8b6e Nico*0494
b4dc6cd434 Jean*0495 c           gqair(i,j,bi,bj)=(evp-precip-entrain)/CheapHgrid(i,j,bi,bj)
8fd83faf35 Jean*0496             gqair(i,j,bi,bj) = (evp(i,j)-entrain)/CheapHgrid(i,j,bi,bj)
                0497      &                        /rhoa*maskC(i,j,1,bi,bj)
                0498             EmPmR(i,j,bi,bj) = ( EmPmR(i,j,bi,bj)
                0499      &                          -cheapPrecip(i,j,bi,bj)
b4dc6cd434 Jean*0500      &                         )*maskC(i,j,1,bi,bj)
                0501            ENDIF
ced0783fba Jean*0502
89f7c61169 Jean*0503           ENDDO
                0504          ENDDO
                0505
                0506 C it is not necessary to use the Adams2d subroutine as
                0507 C the forcing is always computed at the current time step.
5251e2c855 Jean*0508 C note: full oceanic time step deltaT is used below
                0509          CALL CHEAPAML_TIMESTEP(
b4dc6cd434 Jean*0510      I           bi, bj, iMin, iMax, jMin, jMax, deltaT,
5251e2c855 Jean*0511      I           gTair,
                0512      U           Tair,
                0513      I           0, myIter, myThid )
89f7c61169 Jean*0514 C       do implicit time stepping over land
                0515          DO j=1-OLy,sNy+OLy
                0516           DO i=1-OLx,sNx+OLx
d861cd501f Jean*0517             dtemp = tr(i,j,bi,bj)-Tair(i,j,bi,bj)
                0518             Tair(i,j,bi,bj) = Tair(i,j,bi,bj) + dtemp*xrelf(i,j,bi,bj)
89f7c61169 Jean*0519           ENDDO
                0520          ENDDO
51132e5783 Nico*0521
89f7c61169 Jean*0522 C do water
                0523         IF (useFreshWaterFlux) THEN
5251e2c855 Jean*0524          CALL CHEAPAML_TIMESTEP(
b4dc6cd434 Jean*0525      I           bi, bj, iMin, iMax, jMin, jMax, deltaT,
5251e2c855 Jean*0526      I           gqair,
                0527      U           qair,
                0528      I           0, myIter, myThid )
89f7c61169 Jean*0529 C     do implicit time stepping over land and or buffer
                0530          DO j=1-OLy,sNy+OLy
                0531           DO i=1-OLx,sNx+OLx
d861cd501f Jean*0532             dq = qr(i,j,bi,bj)-qair(i,j,bi,bj)
                0533             qair(i,j,bi,bj) = qair(i,j,bi,bj) + dq*xrelf(i,j,bi,bj)
                0534             IF (qair(i,j,bi,bj).LT.zeroRL) qair(i,j,bi,bj) = 0. _d 0
89f7c61169 Jean*0535           ENDDO
                0536          ENDDO
                0537         ENDIF
                0538
                0539 C do tracer
                0540         IF (useCheapTracer) THEN
                0541 C     do implicit time stepping over land and or buffer
                0542          DO j=1-OLy,sNy+OLy
                0543           DO i=1-OLx,sNx+OLx
d861cd501f Jean*0544             dtr = CheaptracerR(i,j,bi,bj)-Cheaptracer(i,j,bi,bj)
89f7c61169 Jean*0545             Cheaptracer(i,j,bi,bj) = Cheaptracer(i,j,bi,bj)
                0546      &                             + dtr*xrelf(i,j,bi,bj)
                0547           ENDDO
                0548          ENDDO
                0549         ENDIF
ced0783fba Jean*0550
8fd83faf35 Jean*0551 #ifdef ALLOW_DIAGNOSTICS
                0552         IF ( useDiagnostics ) THEN
                0553          CALL DIAGNOSTICS_FILL( fsha,'CH_SH   ',0,1,2,bi,bj,myThid )
                0554          CALL DIAGNOSTICS_FILL( flha,'CH_LH   ',0,1,2,bi,bj,myThid )
                0555          CALL DIAGNOSTICS_FILL( q100,'CH_q100 ',0,1,2,bi,bj,myThid )
                0556          CALL DIAGNOSTICS_FILL( ssqt,'CH_ssqt ',0,1,2,bi,bj,myThid )
                0557         ENDIF
                0558 #endif /* ALLOW_DIAGNOSTICS */
                0559
89f7c61169 Jean*0560 C close bi,bj loops
                0561        ENDDO
                0562       ENDDO
                0563
                0564 C update edges
b4dc6cd434 Jean*0565        _EXCH_XY_RL(Tair,myThid)
                0566        _EXCH_XY_RS(Qnet,myThid)
89f7c61169 Jean*0567       IF (useFreshWaterFlux) THEN
b4dc6cd434 Jean*0568        _EXCH_XY_RL(qair,myThid)
                0569        _EXCH_XY_RS(EmPmR,myThid)
89f7c61169 Jean*0570       ENDIF
                0571       IF (useCheapTracer) THEN
b4dc6cd434 Jean*0572        _EXCH_XY_RL(Cheaptracer,myThid)
                0573       ENDIF
d861cd501f Jean*0574       IF ( .NOT.useStressOption ) THEN
                0575        _EXCH_XY_RL( surfDrag, myThid )
89f7c61169 Jean*0576       ENDIF
                0577
                0578 C reset edges to open boundary profiles
d25d6ad15e Jean*0579 c     IF ( .NOT.(cheapamlXperiodic.AND.cheapamlYperiodic) ) THEN
5251e2c855 Jean*0580       IF ( notUsingXPeriodicity.OR.notUsingYPeriodicity ) THEN
                0581         xIsPeriodic = .NOT.notUsingXPeriodicity
                0582         yIsPeriodic = .NOT.notUsingYPeriodicity
                0583         DO bj=myByLo(myThid),myByHi(myThid)
                0584          DO bi=myBxLo(myThid),myBxHi(myThid)
                0585            CALL CHEAPAML_COPY_EDGES(
                0586 c    I                   cheapamlXperiodic, cheapamlYperiodic,
                0587      I                   xIsPeriodic, yIsPeriodic,
                0588      I                   Tr(1-OLx,1-OLy,bi,bj),
                0589      U                   Tair(1-OLx,1-OLy,bi,bj),
                0590      I                   bi, bj, myIter, myThid )
                0591           IF (useFreshWaterFlux) THEN
                0592            CALL CHEAPAML_COPY_EDGES(
                0593 c    I                   cheapamlXperiodic, cheapamlYperiodic,
                0594      I                   xIsPeriodic, yIsPeriodic,
                0595      I                   qr(1-OLx,1-OLy,bi,bj),
                0596      U                   qair(1-OLx,1-OLy,bi,bj),
                0597      I                   bi, bj, myIter, myThid )
                0598           ENDIF
                0599           IF (useCheapTracer) THEN
                0600            CALL CHEAPAML_COPY_EDGES(
                0601 c    I                   cheapamlXperiodic, cheapamlYperiodic,
                0602      I                   xIsPeriodic, yIsPeriodic,
                0603      I                   CheaptracerR(1-OLx,1-OLy,bi,bj),
                0604      U                   Cheaptracer(1-OLx,1-OLy,bi,bj),
                0605      I                   bi, bj, myIter, myThid )
                0606           ENDIF
89f7c61169 Jean*0607          ENDDO
51132e5783 Nico*0608         ENDDO
89f7c61169 Jean*0609       ENDIF
c7cc66b68a Jean*0610
8fd83faf35 Jean*0611 c     CALL PLOT_FIELD_XYRS( gTair, 'S/R CHEAPAML gTair',1,myThid)
                0612 c     CALL PLOT_FIELD_XYRS( Tair, 'S/R CHEAPAML Tair',1,myThid)
                0613 c     CALL PLOT_FIELD_XYRS( Qnet, 'S/R CHEAPAML Qnet',1,myThid)
cf5b5345a0 Jean*0614
                0615       DO bj=myByLo(myThid),myByHi(myThid)
                0616        DO bi=myBxLo(myThid),myBxHi(myThid)
b4dc6cd434 Jean*0617
d861cd501f Jean*0618         IF ( .NOT.useStressOption ) THEN
d54b0079d9 Brun*0619           DO j = 1-OLy+1,sNy+OLy
d861cd501f Jean*0620            DO i = 1-OLx+1,sNx+OLx
                0621             fu(i,j,bi,bj) = maskW(i,j,1,bi,bj)*halfRL
b4dc6cd434 Jean*0622      &          *( surfDrag(i-1,j,bi,bj) + surfDrag(i,j,bi,bj) )
d861cd501f Jean*0623      &          * uRelWind(i,j,bi,bj)
                0624             fv(i,j,bi,bj) = maskS(i,j,1,bi,bj)*halfRL
b4dc6cd434 Jean*0625      &          *( surfDrag(i,j-1,bi,bj) + surfDrag(i,j,bi,bj) )
d861cd501f Jean*0626      &          * vRelWind(i,j,bi,bj)
                0627            ENDDO
b4dc6cd434 Jean*0628           ENDDO
                0629 #ifdef INCONSISTENT_WIND_LOCATION
d861cd501f Jean*0630           DO j = 1-OLy,sNy+OLy
                0631            DO i = 1-OLx+1,sNx+OLx
                0632             fu(i,j,bi,bj) = maskW(i,j,1,bi,bj)*halfRL
                0633      &          *( surfDrag(i-1,j,bi,bj)*uRelWind(i-1,j,bi,bj)
                0634      &           + surfDrag(i,j,bi,bj) * uRelWind(i,j,bi,bj) )
                0635            ENDDO
b4dc6cd434 Jean*0636           ENDDO
d54b0079d9 Brun*0637           DO j = 1-OLy+1,sNy+OLy
d861cd501f Jean*0638            DO i = 1-OLx,sNx+OLx
                0639             fv(i,j,bi,bj) = maskS(i,j,1,bi,bj)*halfRL
                0640      &          *( surfDrag(i,j-1,bi,bj)*vRelWind(i,j-1,bi,bj)
                0641      &           + surfDrag(i,j,bi,bj) * vRelWind(i,j,bi,bj) )
                0642            ENDDO
d54b0079d9 Brun*0643           ENDDO
                0644 #endif /* INCONSISTENT_WIND_LOCATION */
                0645
d861cd501f Jean*0646         ELSE ! useStressOption
89f7c61169 Jean*0647 Cswd move wind stresses to u and v points
                0648          DO j = 1-OLy,sNy+OLy
                0649           DO i = 1-OLx+1,sNx+OLx
d861cd501f Jean*0650             fu(i,j,bi,bj) = maskW(i,j,1,bi,bj)*halfRL
b4dc6cd434 Jean*0651      &          *( ustress(i,j,bi,bj) + ustress(i-1,j,bi,bj) )
89f7c61169 Jean*0652           ENDDO
                0653          ENDDO
                0654          DO j = 1-OLy+1,sNy+OLy
                0655           DO i = 1-OLx,sNx+OLx
d861cd501f Jean*0656             fv(i,j,bi,bj) = maskS(i,j,1,bi,bj)*halfRL
b4dc6cd434 Jean*0657      &          *( vstress(i,j,bi,bj) + vstress(i,j-1,bi,bj) )
89f7c61169 Jean*0658           ENDDO
                0659          ENDDO
b4dc6cd434 Jean*0660         ENDIF
cf5b5345a0 Jean*0661
                0662 C--   end bi,bj loops
                0663        ENDDO
                0664       ENDDO
2616d73cb2 Nico*0665
4fa4901be6 Nico*0666 #ifdef ALLOW_DIAGNOSTICS
bf944b1865 Jean*0667 C- note: with thSIce, CH_QNET and CH_EmP correspond to the fluxes over the
                0668 C     open-ocean fraction of the grid-cell (similar to EXFqnet and EXFempmr).
                0669 C     Use instead diagnostics SIflxAtm & SIfrwAtm to get the grid-cell
                0670 C      averaged (open-ocean fraction + ice-covered fraction).
b4dc6cd434 Jean*0671       IF ( useDiagnostics ) THEN
                0672        CALL DIAGNOSTICS_FILL(uWind,  'CH_Uwind',0,1,0,1,1,myThid)
                0673        CALL DIAGNOSTICS_FILL(vWind,  'CH_Vwind',0,1,0,1,1,myThid)
                0674        CALL DIAGNOSTICS_FILL_RS(Qnet,'CH_QNET ',0,1,0,1,1,myThid)
                0675        IF (useFreshWaterFlux) THEN
                0676         CALL DIAGNOSTICS_FILL_RS( EmPmR, 'CH_EmP  ', 0,1,0,1,1,myThid)
                0677         CALL DIAGNOSTICS_FILL(cheapPrecip,'CH_Prec ',0,1,0,1,1,myThid)
                0678        ENDIF
                0679        IF (useCheapTracer) THEN
89f7c61169 Jean*0680         CALL DIAGNOSTICS_FILL(Cheaptracer,'CH_Trace',0,1,0,1,1,myThid)
b4dc6cd434 Jean*0681        ENDIF
51132e5783 Nico*0682       ENDIF
4fa4901be6 Nico*0683 #endif /* ALLOW_DIAGNOSTICS */
c7cc66b68a Jean*0684
89f7c61169 Jean*0685 c     DO bj=myByLo(myThid),myByHi(myThid)
                0686 c      DO bi=myBxLo(myThid),myBxHi(myThid)
                0687 c        DO j = 1-OLy,sNy+OLy
                0688 c         DO i = 1-OLx+1,sNx+OLx
                0689 c           fu(i,j,bi,bj) = 0.0
                0690 c           fv(i,j,bi,bj) = 0.0
                0691 c           Qnet(i,j,bi,bj) = 0.0
                0692 c           EmPmR(i,j,bi,bj) = 0.0
                0693 c         ENDDO
                0694 c        ENDDO
                0695 c      ENDDO
                0696 c     ENDDO
cf5b5345a0 Jean*0697
                0698       RETURN
                0699       END