model/src/ini_forcing.F

212a8d049e Ed H*0001 #include "PACKAGES_CONFIG.h"
1dbaea09ee Chri*0002 #include "CPP_OPTIONS.h"
924557e60a Chri*0003
9366854e02 Chri*0004 CBOP
                0005 C     !ROUTINE: INI_FORCING
                0006 C     !INTERFACE:
924557e60a Chri*0007       SUBROUTINE INI_FORCING( myThid )
                0008
9366854e02 Chri*0009 C     !DESCRIPTION: \bv
                0010 C     *==========================================================*
d197c88195 Jean*0011 C     | SUBROUTINE INI_FORCING
                0012 C     | o Set model initial forcing fields.
9366854e02 Chri*0013 C     *==========================================================*
                0014 C     \ev
                0015
                0016 C     !USES:
                0017       IMPLICIT NONE
924557e60a Chri*0018 C     === Global variables ===
                0019 #include "SIZE.h"
                0020 #include "EEPARAMS.h"
                0021 #include "PARAMS.h"
                0022 #include "GRID.h"
a5e98ae15f Jean*0023 #include "SURFACE.h"
924557e60a Chri*0024 #include "FFIELDS.h"
                0025
9366854e02 Chri*0026 C     !INPUT/OUTPUT PARAMETERS:
924557e60a Chri*0027 C     == Routine arguments ==
983485b08a Jean*0028 C     myThid :: my Thread Id number
924557e60a Chri*0029       INTEGER myThid
                0030
9366854e02 Chri*0031 C     !LOCAL VARIABLES:
924557e60a Chri*0032 C     == Local variables ==
983485b08a Jean*0033 C     bi,bj  :: Tile indices
                0034 C     i, j   :: Loop counters
924557e60a Chri*0035       INTEGER bi, bj
d197c88195 Jean*0036       INTEGER  i, j
00c7090dc0 Mart*0037 #ifdef SHORTWAVE_HEATING
                0038       INTEGER  k, km
                0039       _RL SWFracK(Nr+1), swfac
                0040 #endif
9366854e02 Chri*0041 CEOP
924557e60a Chri*0042
d197c88195 Jean*0043 C-    Initialise all arrays in common blocks
d067a44b7e Jean*0044 C     <-- moved to new S/R INI_FFIELDS
d197c88195 Jean*0045
28659cf1dc Patr*0046       DO bj = myByLo(myThid), myByHi(myThid)
                0047        DO bi = myBxLo(myThid), myBxHi(myThid)
4eb48150b4 Jean*0048         DO j=1-OLy,sNy+OLy
                0049          DO i=1-OLx,sNx+OLx
d197c88195 Jean*0050           IF ( doThetaClimRelax .AND.
                0051      &         ABS(yC(i,j,bi,bj)).LE.latBandClimRelax ) THEN
                0052            lambdaThetaClimRelax(i,j,bi,bj) = 1. _d 0/tauThetaClimRelax
28659cf1dc Patr*0053           ELSE
d197c88195 Jean*0054            lambdaThetaClimRelax(i,j,bi,bj) = 0. _d 0
28659cf1dc Patr*0055           ENDIF
                0056           IF ( doSaltClimRelax .AND.
d197c88195 Jean*0057      &         ABS(yC(i,j,bi,bj)).LE.latBandClimRelax ) THEN
                0058            lambdaSaltClimRelax(i,j,bi,bj) = 1. _d 0/tauSaltClimRelax
28659cf1dc Patr*0059           ELSE
d197c88195 Jean*0060            lambdaSaltClimRelax(i,j,bi,bj) = 0. _d 0
28659cf1dc Patr*0061           ENDIF
                0062          ENDDO
                0063         ENDDO
                0064        ENDDO
                0065       ENDDO
d197c88195 Jean*0066
                0067 C-    every-one waits before master thread loads from file
3365bdc872 Jean*0068 C     this is done within IO routines => no longer needed
                0069 c     _BARRIER
d197c88195 Jean*0070
ab42872a05 Alis*0071       IF ( zonalWindFile .NE. ' '  ) THEN
                0072        CALL READ_FLD_XY_RS( zonalWindFile, ' ', fu, 0, myThid )
                0073       ENDIF
                0074       IF ( meridWindFile .NE. ' '  ) THEN
                0075        CALL READ_FLD_XY_RS( meridWindFile, ' ', fv, 0, myThid )
                0076       ENDIF
                0077       IF ( surfQFile .NE. ' '  ) THEN
                0078        CALL READ_FLD_XY_RS( surfQFile, ' ', Qnet, 0, myThid )
2d2cc93d4f Jean*0079       ELSEIF ( surfQnetFile .NE. ' '  ) THEN
                0080        CALL READ_FLD_XY_RS( surfQnetFile, ' ', Qnet, 0, myThid )
ab42872a05 Alis*0081       ENDIF
                0082       IF ( EmPmRfile .NE. ' '  ) THEN
                0083        CALL READ_FLD_XY_RS( EmPmRfile, ' ', EmPmR, 0, myThid )
62fd6ae4e5 Jean*0084 c      IF ( convertEmP2rUnit.EQ.mass2rUnit ) THEN
b5f408f39d Jean*0085 C-     EmPmR is now (after c59h) expressed in kg/m2/s (fresh water mass flux)
                0086         DO bj = myByLo(myThid), myByHi(myThid)
                0087          DO bi = myBxLo(myThid), myBxHi(myThid)
4eb48150b4 Jean*0088           DO j=1-OLy,sNy+OLy
                0089            DO i=1-OLx,sNx+OLx
b5f408f39d Jean*0090             EmPmR(i,j,bi,bj) = EmPmR(i,j,bi,bj)*rhoConstFresh
                0091            ENDDO
                0092           ENDDO
                0093          ENDDO
                0094         ENDDO
62fd6ae4e5 Jean*0095 c      ENDIF
ab42872a05 Alis*0096       ENDIF
1e273d1bf5 Jean*0097       IF ( saltFluxFile .NE. ' '  ) THEN
                0098        CALL READ_FLD_XY_RS( saltFluxFile, ' ', saltFlux, 0, myThid )
                0099       ENDIF
ab42872a05 Alis*0100       IF ( thetaClimFile .NE. ' '  ) THEN
                0101        CALL READ_FLD_XY_RS( thetaClimFile, ' ', SST, 0, myThid )
                0102       ENDIF
                0103       IF ( saltClimFile .NE. ' '  ) THEN
                0104        CALL READ_FLD_XY_RS( saltClimFile, ' ', SSS, 0, myThid )
                0105       ENDIF
28659cf1dc Patr*0106       IF ( lambdaThetaFile .NE. ' '  ) THEN
d197c88195 Jean*0107        CALL READ_FLD_XY_RS( lambdaThetaFile, ' ',
28659cf1dc Patr*0108      &  lambdaThetaClimRelax, 0, myThid )
                0109       ENDIF
                0110       IF ( lambdaSaltFile .NE. ' '  ) THEN
d197c88195 Jean*0111        CALL READ_FLD_XY_RS( lambdaSaltFile, ' ',
28659cf1dc Patr*0112      &  lambdaSaltClimRelax, 0, myThid )
                0113       ENDIF
310851b9c0 Patr*0114 #ifdef SHORTWAVE_HEATING
00c7090dc0 Mart*0115       IF ( surfQswFile .NE. ' ' ) THEN
310851b9c0 Patr*0116        CALL READ_FLD_XY_RS( surfQswFile, ' ', Qsw, 0, myThid )
2d2cc93d4f Jean*0117        IF ( surfQFile .NE. ' '  ) THEN
                0118 C-     Qnet is now (after c54) the net Heat Flux (including SW)
d197c88195 Jean*0119         DO bj = myByLo(myThid), myByHi(myThid)
                0120          DO bi = myBxLo(myThid), myBxHi(myThid)
4e530425d3 Jean*0121           DO j=1-OLy,sNy+OLy
                0122            DO i=1-OLx,sNx+OLx
2d2cc93d4f Jean*0123             Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj) + Qsw(i,j,bi,bj)
4e530425d3 Jean*0124            ENDDO
2d2cc93d4f Jean*0125           ENDDO
                0126          ENDDO
4e530425d3 Jean*0127         ENDDO
2d2cc93d4f Jean*0128        ENDIF
310851b9c0 Patr*0129       ENDIF
00c7090dc0 Mart*0130 C--   initialisation for the case of no shortwave penetration
                0131 C--   (shortwave radiation only heats the top layer)
                0132       DO bj = myByLo(myThid), myByHi(myThid)
                0133        DO bi = myBxLo(myThid), myBxHi(myThid)
                0134         DO k = 1,Nr+1
                0135          swfac = 0. _d 0
                0136          IF ( usingZCoords ) THEN
                0137           IF ( k .EQ. 1 ) swfac = 1. _d 0
                0138          ELSE
                0139           IF ( k .EQ. Nr+1 ) swfac = 1. _d 0
                0140          ENDIF
                0141          DO j=1-OLy,sNy+OLy
                0142           DO i=1-OLx,sNx+OLx
                0143            SWFrac3D(i,j,k,bi,bj) = swfac
                0144           ENDDO
                0145          ENDDO
                0146         ENDDO
                0147        ENDDO
                0148       ENDDO
                0149       IF ( selectPenetratingSW .GT. 0 ) THEN
                0150 C     For now SWFrac3D is held fixed throughout the simulation, so it is
                0151 C     initialised here instead of somewhere in S/R INITIALISE_VARIA
                0152        DO k=1,Nr+1
                0153         IF ( usingZCoords ) THEN
                0154          SWFracK(k) = rF(k) - rF(1)
                0155         ELSE
                0156 C     this is the oceanic pressure coordinate case
                0157          SWFracK(k) = ( rF(Nr+1) - rF(k) )
                0158      &              * recip_rhoConst*recip_gravity
                0159 C     valid also with Mass-Coordinate (but not with variable gravity):
                0160 c        SWFracK(k) = ( phiRef(2*k-1) - phiRef(2*Nr+1) )*recip_gravity
                0161         ENDIF
                0162        ENDDO
                0163        CALL SWFRAC(
                0164      I             Nr+1, oneRL,
                0165      U             SWFracK,
                0166      I             zeroRL, 0, myThid )
                0167        DO bj = myByLo(myThid), myByHi(myThid)
                0168         DO bi = myBxLo(myThid), myBxHi(myThid)
                0169          DO k = 1,Nr+1
                0170           swfac = 1. _d 0
                0171 C     Here, km is the index for the mask physically below the interface k.
                0172           IF ( usingZCoords ) THEN
                0173            km = MIN(k,Nr)
                0174            IF ( k .EQ. Nr+1 ) swfac = 0. _d 0
                0175           ELSE
                0176            km = MAX(k-1,1)
                0177            IF ( k .EQ. 1 ) swfac = 0. _d 0
                0178           ENDIF
                0179           DO j=1-OLy,sNy+OLy
                0180            DO i=1-OLx,sNx+OLx
                0181             SWFrac3D(i,j,k,bi,bj) = SWFracK(k)*swfac
                0182      &           *maskC(i,j,km,bi,bj)
                0183            ENDDO
                0184           ENDDO
                0185          ENDDO
                0186         ENDDO
                0187        ENDDO
                0188 C     endif selectPenetratingSW
                0189       ENDIF
                0190 #endif /* SHORTWAVE_HEATING */
395b093796 Mart*0191 #ifdef ATMOSPHERIC_LOADING
                0192       IF ( pLoadFile .NE. ' '  ) THEN
a8bcab80b7 Jean*0193        CALL READ_FLD_XY_RS( pLoadFile, ' ', pLoad, 0, myThid )
395b093796 Mart*0194       ENDIF
00c7090dc0 Mart*0195 #endif /* ATMOSPHERIC_LOADING */
4eb48150b4 Jean*0196 #ifdef ALLOW_ADDFLUID
                0197       IF ( addMassFile .NE. ' ' ) THEN
                0198        CALL READ_FLD_XYZ_RL( addMassFile, ' ', addMass, 0, myThid )
                0199        CALL EXCH_XYZ_RL( addMass, myThid )
                0200       ENDIF
                0201 #endif /* ALLOW_ADDFLUID */
90929f8806 Patr*0202 #ifdef ALLOW_GEOTHERMAL_FLUX
                0203       IF ( geothermalFile .NE. ' ' ) THEN
d067a44b7e Jean*0204        CALL READ_FLD_XY_RS( geothermalFile, ' ',
90929f8806 Patr*0205      &  geothermalFlux, 0, myThid )
                0206        CALL EXCH_XY_RS( geothermalFlux, myThid )
                0207 # ifdef ALLOW_MONITOR
                0208        CALL MON_PRINTSTATS_RS(
                0209      &  1,geothermalFlux,'geothermalFlux',myThid)
                0210 # endif
                0211       ENDIF
                0212 #endif /* ALLOW_GEOTHERMAL_FLUX */
7e00d7e8f9 Jean*0213 #ifdef ALLOW_BALANCE_FLUXES
                0214       IF ( selectBalanceEmPmR.EQ.2 ) THEN
                0215 C-    set default weight to 1 (i.e., same correction as selectBalanceEmPmR=1 )
                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             weight2BalanceFlx(i,j,bi,bj) = oneRS
                0221           ENDDO
                0222          ENDDO
                0223         ENDDO
                0224        ENDDO
                0225       ENDIF
                0226       IF ( wghtBalanceFile .NE. ' ' ) THEN
                0227        CALL READ_FLD_XY_RS( wghtBalanceFile, ' ',
                0228      &                      weight2BalanceFlx, 0, myThid )
                0229        CALL EXCH_XY_RS( weight2BalanceFlx, myThid )
                0230       ENDIF
                0231 #endif /* ALLOW_GEOTHERMAL_FLUX */
339a1b85b2 Patr*0232
023bfd6664 Jean*0233       CALL EXCH_UV_XY_RS( fu,fv, .TRUE., myThid )
                0234       CALL EXCH_XY_RS( Qnet , myThid )
                0235       CALL EXCH_XY_RS( EmPmR, myThid )
                0236       CALL EXCH_XY_RS( saltFlux, myThid )
                0237       CALL EXCH_XY_RS( SST  , myThid )
                0238       CALL EXCH_XY_RS( SSS  , myThid )
                0239       CALL EXCH_XY_RS( lambdaThetaClimRelax, myThid )
                0240       CALL EXCH_XY_RS( lambdaSaltClimRelax , myThid )
395b093796 Mart*0241 #ifdef SHORTWAVE_HEATING
00c7090dc0 Mart*0242       IF ( surfQswFile .NE. ' ' )
                0243      &     CALL EXCH_XY_RS( Qsw  , myThid )
395b093796 Mart*0244 #endif
                0245 #ifdef ATMOSPHERIC_LOADING
0320e25227 Mart*0246       CALL EXCH_XY_RS( pLoad  , myThid )
a8bcab80b7 Jean*0247 C     CALL PLOT_FIELD_XYRS( pLoad, 'S/R INI_FORCING pLoad',1,myThid)
395b093796 Mart*0248 #endif
0127add478 Alis*0249 C     CALL PLOT_FIELD_XYRS( fu, 'S/R INI_FORCING FU',1,myThid)
                0250 C     CALL PLOT_FIELD_XYRS( fv, 'S/R INI_FORCING FV',1,myThid)
c1dd0647a3 Chri*0251
924557e60a Chri*0252       RETURN
                0253       END