Back to home page

MITgcm

 
 

    


File indexing completed on 2025-02-02 06:10:51 UTC

view on githubraw file Latest commit 701e10a9 on 2025-02-01 19:15:20 UTC
6d54cf9ca1 Ed H*0001 #include "BULK_FORCE_OPTIONS.h"
7753507405 Curt*0002 
dd80d278b6 Jean*0003 CBOP
                0004 C     !ROUTINE: BULKF_FORCING
                0005 C     !INTERFACE:
679d149d01 Jean*0006       SUBROUTINE BULKF_FORCING(
dd80d278b6 Jean*0007      I                          myTime, myIter, myThid )
7753507405 Curt*0008 
dd80d278b6 Jean*0009 C     !DESCRIPTION: \bv
                0010 C     *==========================================================*
                0011 C     | SUBROUTINE BULKF_FORCING
                0012 C     *==========================================================*
                0013 C     \ev
7753507405 Curt*0014 
dd80d278b6 Jean*0015 C     o Get the surface fluxes used to force ocean model
                0016 C       Output:
                0017 C       ------
                0018 C       ustress, vstress :: wind stress
                0019 C       qnet             :: net heat flux
                0020 C       empmr            :: freshwater flux
                0021 C       ---------
                0022 C
                0023 C       Input:
                0024 C       ------
                0025 C       uwind, vwind  :: mean wind speed (m/s)     at height hu (m)
                0026 C       Tair  :: mean air temperature (K)  at height ht (m)
                0027 C       Qair  :: mean air humidity (kg/kg) at height hq (m)
                0028 C       theta(k=1) :: sea surface temperature (K)
                0029 C       rain   :: precipitation
                0030 C       runoff :: river(ice) runoff
                0031 C
                0032 C     ==================================================================
                0033 C     SUBROUTINE bulkf_forcing
                0034 C     ==================================================================
7753507405 Curt*0035 
dd80d278b6 Jean*0036 C     !USES:
                0037       IMPLICIT NONE
                0038 C     == global variables ==
7753507405 Curt*0039 #include "EEPARAMS.h"
                0040 #include "SIZE.h"
                0041 #include "PARAMS.h"
                0042 #include "DYNVARS.h"
                0043 #include "GRID.h"
                0044 #include "FFIELDS.h"
6a1d3c464b Jean*0045 #include "BULKF_PARAMS.h"
7753507405 Curt*0046 #include "BULKF.h"
f4245d1665 Curt*0047 #include "BULKF_INT.h"
dd80d278b6 Jean*0048 #include "BULKF_TAVE.h"
7753507405 Curt*0049 
dd80d278b6 Jean*0050 C     !INPUT/OUTPUT PARAMETERS:
                0051 C     == routine arguments ==
6a1d3c464b Jean*0052       _RL     myTime
dd80d278b6 Jean*0053       INTEGER myIter
                0054       INTEGER myThid
                0055 CEOP
7753507405 Curt*0056 
6a1d3c464b Jean*0057 #ifdef ALLOW_BULK_FORCE
7753507405 Curt*0058 C     == Local variables ==
679d149d01 Jean*0059       INTEGER bi,bj
e5b783de15 Jean*0060       INTEGER i,j
                0061       INTEGER ks, iceornot
7753507405 Curt*0062 
679d149d01 Jean*0063       _RL     df0dT, hfl, evp, dEvdT
                0064 #ifdef ALLOW_FORMULA_AIM
                0065       _RL     SHF(1), EVPloc(1), SLRU(1)
                0066       _RL     dEvp(1), sFlx(0:2)
                0067 #endif
7753507405 Curt*0068 
e5b783de15 Jean*0069 C-    surface level index:
                0070       ks = 1
701e10a905 Mart*0071       IF ( usingPcoords ) ks = Nr
7753507405 Curt*0072 
e5b783de15 Jean*0073 C-    Compute surface fluxes over ice-free ocean only:
                0074       iceornot = 0
7753507405 Curt*0075 
6a1d3c464b Jean*0076       DO bj=myByLo(myThid),myByHi(myThid)
                0077        DO bi=myBxLo(myThid),myBxHi(myThid)
e5b783de15 Jean*0078 
70964a532e Jean*0079          DO j = 1-OLy,sNy+OLy
                0080           DO i = 1-OLx,sNx+OLx
e5b783de15 Jean*0081            IF ( maskC(i,j,ks,bi,bj).NE.0. _d 0 ) THEN
f4245d1665 Curt*0082 
679d149d01 Jean*0083 #ifdef ALLOW_FORMULA_AIM
                0084              IF ( useFluxFormula_AIM ) THEN
                0085                CALL BULKF_FORMULA_AIM(
701e10a905 Mart*0086      I            gcmSST(i,j,bi,bj), flwdwn(i,j,bi,bj),
679d149d01 Jean*0087      I            thAir(i,j,bi,bj), Tair(i,j,bi,bj),
                0088      I            Qair(i,j,bi,bj), wspeed(i,j,bi,bj),
                0089      O            SHF, EVPloc, SLRU,
                0090      O            dEvp, sFlx,
                0091      I            iceornot, myThid )
                0092 
                0093                   flwup(i,j,bi,bj)= ocean_emissivity*SLRU(1)
                0094 C-    reverse sign (AIM convention -> BULKF convention):
                0095                   fsh(i,j,bi,bj) = -SHF(1)
                0096                   flh(i,j,bi,bj) = -Lvap*EVPloc(1)
                0097 C-    Convert from g/m2/s to m/s
f664a6d8bb Jean*0098                   evap(i,j,bi,bj) = EVPloc(1) * 1. _d -3 / rhoFW
679d149d01 Jean*0099                   dEvdT = dEvp(1) * 1. _d -3
                0100                   df0dT = sFlx(2)
e5b783de15 Jean*0101 
                0102              ELSEIF ( blk_nIter.EQ.0 ) THEN
679d149d01 Jean*0103 #else  /* ALLOW_FORMULA_AIM */
e5b783de15 Jean*0104              IF ( blk_nIter.EQ.0 ) THEN
679d149d01 Jean*0105 #endif /* ALLOW_FORMULA_AIM */
548c63e38c Jean*0106                CALL BULKF_FORMULA_LANL(
e5b783de15 Jean*0107      I            uwind(i,j,bi,bj),vwind(i,j,bi,bj),wspeed(i,j,bi,bj),
                0108      I            Tair(i,j,bi,bj), Qair(i,j,bi,bj),
701e10a905 Mart*0109      I            cloud(i,j,bi,bj), gcmSST(i,j,bi,bj),
679d149d01 Jean*0110      O            flwup(i,j,bi,bj), flh(i,j,bi,bj),
                0111      O            fsh(i,j,bi,bj), df0dT,
                0112      O            ustress(i,j,bi,bj), vstress(i,j,bi,bj),
                0113      O            evp, savssq(i,j,bi,bj), dEvdT,
                0114      I            iceornot, myThid )
f4245d1665 Curt*0115 C               Note that the LANL flux conventions are opposite
                0116 C               of what they are in the model.
                0117 
e5b783de15 Jean*0118 C-             Convert from kg/m2/s to m/s
                0119                evap(i,j,bi,bj) = evp/rhoFW
                0120 
                0121              ELSE
f664a6d8bb Jean*0122                CALL BULKF_FORMULA_LAY(
                0123      I            uwind(i,j,bi,bj), vwind(i,j,bi,bj),
                0124      I            wspeed(i,j,bi,bj), Tair(i,j,bi,bj),
701e10a905 Mart*0125      I            Qair(i,j,bi,bj), gcmSST(i,j,bi,bj),
f664a6d8bb Jean*0126      O            flwup(i,j,bi,bj), flh(i,j,bi,bj),
                0127      O            fsh(i,j,bi,bj), df0dT,
                0128      O            ustress(i,j,bi,bj), vstress(i,j,bi,bj),
                0129      O            evp, savssq(i,j,bi,bj), dEvdT,
                0130      I            iceornot, i,j,bi,bj,myThid )
                0131 
548c63e38c Jean*0132 C-             Convert from kg/m2/s to m/s
f664a6d8bb Jean*0133                evap(i,j,bi,bj) = evp/rhoFW
548c63e38c Jean*0134 
679d149d01 Jean*0135              ENDIF
                0136 
e5b783de15 Jean*0137 C- use down long wave data
                0138              flwupnet(i,j,bi,bj)=flwup(i,j,bi,bj)-flwdwn(i,j,bi,bj)
                0139 C- using down solar, need to have water albedo -- .1
                0140              fswnet(i,j,bi,bj) = solar(i,j,bi,bj)
                0141      &                         *( 1. _d 0 - ocean_albedo )
679d149d01 Jean*0142            ElSE
                0143              ustress(i,j,bi,bj) = 0. _d 0
                0144              vstress(i,j,bi,bj) = 0. _d 0
                0145              fsh(i,j,bi,bj)     = 0. _d 0
                0146              flh(i,j,bi,bj)     = 0. _d 0
                0147              flwup(i,j,bi,bj)   = 0. _d 0
                0148              evap(i,j,bi,bj)    = 0. _d 0
                0149              fswnet(i,j,bi,bj)  = 0. _d 0
                0150              savssq(i,j,bi,bj)  = 0. _d 0
                0151            ENDIF
                0152           ENDDO
                0153          ENDDO
                0154 
                0155          IF ( calcWindStress ) THEN
e5b783de15 Jean*0156 C-  move wind stresses to u and v points
70964a532e Jean*0157            DO j = 1-OLy,sNy+OLy
                0158             DO i = 1-OLx+1,sNx+OLx
6a1d3c464b Jean*0159               fu(i,j,bi,bj) = maskW(i,j,1,bi,bj)
                0160      &          *(ustress(i,j,bi,bj)+ustress(i-1,j,bi,bj))*0.5 _d 0
                0161             ENDDO
                0162            ENDDO
70964a532e Jean*0163            DO j = 1-OLy+1,sNy+OLy
                0164             DO i = 1-OLx,sNx+OLx
6a1d3c464b Jean*0165               fv(i,j,bi,bj) = maskS(i,j,1,bi,bj)
                0166      &          *(vstress(i,j,bi,bj)+vstress(i,j-1,bi,bj))*0.5 _d 0
                0167             ENDDO
                0168            ENDDO
679d149d01 Jean*0169          ENDIF
6a1d3c464b Jean*0170 
e5b783de15 Jean*0171 C-    Add all contributions.
70964a532e Jean*0172          DO j = 1-OLy,sNy+OLy
                0173           DO i = 1-OLx,sNx+OLx
e5b783de15 Jean*0174             IF ( maskC(i,j,ks,bi,bj).NE.0. _d 0 ) THEN
                0175 C-       Net downward surface heat flux :
7753507405 Curt*0176               hfl = 0. _d 0
f4245d1665 Curt*0177               hfl = hfl + fsh(i,j,bi,bj)
                0178               hfl = hfl + flh(i,j,bi,bj)
6a1d3c464b Jean*0179               hfl = hfl - flwupnet(i,j,bi,bj)
                0180               hfl = hfl + fswnet(i,j,bi,bj)
e5b783de15 Jean*0181 C- Heat flux:
6a1d3c464b Jean*0182               Qnet(i,j,bi,bj) = -hfl
e96c64fcd5 Jean*0183               Qsw (i,j,bi,bj) = -fswnet(i,j,bi,bj)
7753507405 Curt*0184 #ifdef COUPLE_MODEL
                0185               dFdT(i,j,bi,bj) = df0dT
                0186 #endif
e5b783de15 Jean*0187 C- Fresh-water flux from Precipitation and Evaporation.
6a1d3c464b Jean*0188               EmPmR(i,j,bi,bj) = (evap(i,j,bi,bj)-rain(i,j,bi,bj)
a5003302cb Jean*0189      &                           - runoff(i,j,bi,bj))*rhoConstFresh
e5b783de15 Jean*0190 C---- cheating: now done in S/R BULKF_FLUX_ADJUST, over ice-free ocean only
7753507405 Curt*0191 c            Qnet(i,j,bi,bj) = Qnetch(i,j,bi,bj)
                0192 c            EmPmR(i,j,bi,bj) = EmPch(i,j,bi,bj)
e5b783de15 Jean*0193 C----
679d149d01 Jean*0194             ELSE
6a1d3c464b Jean*0195               Qnet(i,j,bi,bj) = 0. _d 0
e96c64fcd5 Jean*0196               Qsw (i,j,bi,bj) = 0. _d 0
6a1d3c464b Jean*0197               EmPmR(i,j,bi,bj)= 0. _d 0
7753507405 Curt*0198 #ifdef COUPLE_MODEL
6a1d3c464b Jean*0199               dFdT(i,j,bi,bj) = 0. _d 0
7753507405 Curt*0200 #endif
679d149d01 Jean*0201             ENDIF
6a1d3c464b Jean*0202           ENDDO
                0203          ENDDO
7753507405 Curt*0204 
70964a532e Jean*0205          IF ( temp_EvPrRn .NE. UNSET_RL ) THEN
                0206 C--   Account for energy content of Precip + RunOff & Evap. Assumes:
                0207 C     1) Rain has same temp as Air
                0208 C     2) Snow has no heat capacity (consistent with seaice & thsice pkgs)
                0209 C     3) No distinction between sea-water Cp and fresh-water Cp
                0210 C     4) Run-Off comes at the temp of surface water (with same Cp)
                0211 C     5) Evap is released to the Atmos @ surf-temp (=SST); should be using
                0212 C        the water-vapor heat capacity here and consistently in Bulk-Formulae;
                0213 C        Could also be put directly into Latent Heat flux.
                0214 c         IF ( SnowFile .NE. ' ' ) THEN
                0215 C--   Melt snow (if provided) into the ocean and account for rain-temp
                0216 c          DO j = 1-OLy,sNy+OLy
                0217 c           DO i = 1-OLx,sNx+OLx
                0218 c             Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj)
                0219 c    &              + Lfresh*snowPrecip(i,j,bi,bj)*rhoConstFresh
                0220 c    &              - HeatCapacity_Cp
                0221 c    &               *( Tair(i,j,bi,bj) - Tf0kel - temp_EvPrRn )
                0222 c    &               *( rain(i,j,bi,bj)- snowPrecip(i,j,bi,bj) )
                0223 c    &               *rhoConstFresh
                0224 c           ENDDO
                0225 c          ENDDO
                0226 c         ELSE
                0227 C--   Make snow (according to Air Temp) and melt it in the ocean
                0228 C     note: here we just use the same criteria as over seaice but would be
                0229 C           better to consider a higher altitude air temp, e.g., 850.mb
                0230            DO j = 1-OLy,sNy+OLy
                0231             DO i = 1-OLx,sNx+OLx
                0232              IF ( Tair(i,j,bi,bj).LE.Tf0kel ) THEN
                0233                Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj)
                0234      &              + Lfresh*rain(i,j,bi,bj)*rhoConstFresh
                0235               ELSE
                0236 C--   Account for rain-temp
                0237                Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj)
                0238      &              - HeatCapacity_Cp
                0239      &               *( Tair(i,j,bi,bj) - Tf0kel - temp_EvPrRn )
                0240      &               *rain(i,j,bi,bj)*rhoConstFresh
                0241               ENDIF
                0242             ENDDO
                0243            ENDDO
                0244 c         ENDIF
                0245 C--   Account for energy content of Evap and RunOff:
                0246           DO j = 1-OLy,sNy+OLy
                0247             DO i = 1-OLx,sNx+OLx
                0248 c             Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj)
                0249 c    &              + ( theta(i,j,ks,bi,bj) - temp_EvPrRn )
                0250 c    &               *( evap(i,j,bi,bj)*cpwv
                0251 c    &                - runoff(i,j,bi,bj)*HeatCapacity_Cp
                0252 c    &                )*rhoConstFresh
                0253               Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj)
                0254      &              + ( theta(i,j,ks,bi,bj) - temp_EvPrRn )
                0255      &               *( evap(i,j,bi,bj) - runoff(i,j,bi,bj) )
                0256      &               *HeatCapacity_Cp*rhoConstFresh
                0257               Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj)*maskC(i,j,ks,bi,bj)
                0258             ENDDO
                0259           ENDDO
                0260          ENDIF
                0261 
679d149d01 Jean*0262          IF ( blk_taveFreq.GT.0. _d 0 )
70964a532e Jean*0263      &     CALL BULKF_AVE( bi, bj, myThid )
7753507405 Curt*0264 
6a1d3c464b Jean*0265 C--   end bi,bj loops
                0266        ENDDO
                0267       ENDDO
7753507405 Curt*0268 
6a1d3c464b Jean*0269 C--   Update the tile edges.
                0270 C jmc: Is it necessary ?
12ffad7671 Jean*0271 c     _EXCH_XY_RS(Qnet,   myThid)
                0272 c     _EXCH_XY_RS(EmPmR,   myThid)
6a1d3c464b Jean*0273 c     CALL EXCH_UV_XY_RS(fu, fv, .TRUE., myThid)
7753507405 Curt*0274 
                0275 #endif  /*ALLOW_BULK_FORCE*/
                0276 
6a1d3c464b Jean*0277       RETURN
                0278       END