pkg/aim_v23/phy_radiat.F

d676f916b2 Jean*0001 #include "AIM_OPTIONS.h"
                0002
670d1925ca Jean*0003 C--  File phy_radiat.F:
                0004 C--   Contents
                0005 C--   o SOL_OZ
                0006 C--   o RADSW
                0007 C--   o RADLW
                0008 C--   o RADSET
                0009
                0010 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0011 CBOP
                0012 C !ROUTINE: SOL_OZ (SOLC,TYEAR)
                0013 C !INTERFACE:
d676f916b2 Jean*0014       SUBROUTINE SOL_OZ (SOLC, TYEAR, SLAT, CLAT,
                0015      O                   FSOL, OZONE, OZUPP, ZENIT, STRATZ,
                0016      I                   bi, bj, myThid)
670d1925ca Jean*0017 C !DESCRIPTION: \bv
                0018 C   Purpose: Compute the flux of incoming solar radiation
                0019 C            and a climatological ozone profile for SW absorption
                0020 C   Input:   SOLC   = solar constant (area averaged)
                0021 C            TYEAR  = time as fraction of year (0-1, 0 = 1jan.h00)
                0022 C            SLAT   = sin(lat)
                0023 C            CLAT   = cos(lat)
                0024 C   Output:  FSOL   = flux of incoming solar radiation
                0025 C            OZONE  = flux absorbed by ozone (lower stratos.)
                0026 C            OZUPP  = flux absorbed by ozone (upper stratos.)
                0027 C            ZENIT  = function of solar zenith angle
                0028 C            STRATZ = ?
                0029 C   Updated common blocks: RADZON
                0030 C *==========================================================*
                0031 C \ev
                0032
                0033 C !USES:
d676f916b2 Jean*0034       IMPLICIT NONE
                0035
                0036 C     Resolution parameters
                0037 C-- size for MITgcm & Physics package :
                0038 #include "AIM_SIZE.h"
                0039
                0040 #include "EEPARAMS.h"
670d1925ca Jean*0041 #include "PARAMS.h"
7f98c35e47 Davi*0042 #include "AIM_PARAMS.h"
d676f916b2 Jean*0043
                0044 C     Constants + functions of sigma and latitude
                0045 #include "com_physcon.h"
                0046 C     Radiation constants
                0047 #include "com_radcon.h"
                0048
670d1925ca Jean*0049 C !INPUT/OUTPUT PARAMETERS:
d676f916b2 Jean*0050       INTEGER  bi, bj, myThid
                0051       _RL SOLC, TYEAR
                0052       _RL SLAT(NGP), CLAT(NGP)
                0053       _RL FSOL(NGP), OZONE(NGP), OZUPP(NGP), ZENIT(NGP), STRATZ(NGP)
670d1925ca Jean*0054 CEOP
d676f916b2 Jean*0055
                0056 #ifdef ALLOW_AIM
670d1925ca Jean*0057 C !LOCAL VARIABLES:
d676f916b2 Jean*0058       INTEGER J, NZEN
                0059       _RL ALPHA, CSR1, CSR2, COZ1, COZ2
                0060       _RL AZEN, RZEN, CZEN, SZEN, AST, FS0, FLAT2
7f98c35e47 Davi*0061 #ifdef ALLOW_INSOLATION
670d1925ca Jean*0062       _RL TanDelcl, cosH, HourAng, TanLat
                0063       _RL largeTan
                0064       largeTan = 1. _d 16
7f98c35e47 Davi*0065 #endif
670d1925ca Jean*0066
d676f916b2 Jean*0067 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0068
                0069 C     ALPHA = year phase ( 0 - 2pi, 0 = winter solstice = 22dec.h00 )
                0070       ALPHA = 4. _d 0*ASIN(1. _d 0)*(TYEAR+10. _d 0/365. _d 0)
                0071
                0072       CSR1=-0.796 _d 0*COS(ALPHA)
                0073       CSR2= 0.147 _d 0*COS(2. _d 0*ALPHA)-0.477 _d 0
                0074       COZ1= 1.0 _d 0 * COS(ALPHA)
                0075       COZ2= 1.8 _d 0
                0076 C
                0077       AZEN=1.0
                0078       NZEN=2
                0079
7f98c35e47 Davi*0080 #ifdef ALLOW_INSOLATION
                0081       SZEN = - SIN(  OBLIQ * PI/180. _d 0) * COS(ALPHA)
                0082       RZEN = ASIN( SZEN )
                0083       CZEN =  COS( RZEN )
                0084       IF ( SZEN .EQ. 1. _d 0 ) THEN
670d1925ca Jean*0085          TanDelcl = largeTan
7f98c35e47 Davi*0086       ELSEIF ( SZEN .EQ. -1. _d 0 ) THEN
670d1925ca Jean*0087          TanDelcl =-largeTan
7f98c35e47 Davi*0088       ELSE
                0089          TanDelcl = SZEN / CZEN
                0090       ENDIF
                0091 #else
d676f916b2 Jean*0092       RZEN=-COS(ALPHA)*23.45 _d 0*ASIN(1. _d 0)/90. _d 0
                0093       CZEN=COS(RZEN)
                0094       SZEN=SIN(RZEN)
7f98c35e47 Davi*0095 #endif
d676f916b2 Jean*0096
                0097       AST=0.025 _d 0
                0098       FS0=10. _d 0
                0099 C     FS0=16.-8.*COS(ALPHA)
                0100
                0101       DO J=1,NGP
                0102
                0103         FLAT2 = 1.5 _d 0*SLAT(J)**2 - 0.5 _d 0
                0104
7f98c35e47 Davi*0105 #ifndef ALLOW_INSOLATION
d676f916b2 Jean*0106 C       solar radiation at the top
                0107         FSOL(J) = SOLC*
                0108      &     MAX( 0. _d 0, 1. _d 0+CSR1*SLAT(J)+CSR2*FLAT2 )
7f98c35e47 Davi*0109 #else
                0110         IF ( CLAT(J) .EQ. 0. _d 0 ) THEN
670d1925ca Jean*0111            TanLat = SIGN(1. _d 0, SLAT(J) ) * largeTan
7f98c35e47 Davi*0112         ELSE
670d1925ca Jean*0113            TanLat = SLAT(J)/CLAT(J)
7f98c35e47 Davi*0114         ENDIF
                0115         cosH     = - TanLat * TanDelcl
                0116         cosH     = MAX(MIN(cosH,1. _d 0), -1. _d 0)
                0117         HourAng  =  ACOS( cosH )
                0118         FSOL(J)  = 4. _d 0 / PI * SOLC *
                0119      &   (SLAT(J)*SZEN*HourAng+CLAT(J)*CZEN*SIN(HourAng))
                0120 #endif
d676f916b2 Jean*0121
                0122 C       ozone depth in upper and lower stratosphere
                0123         OZUPP(J) = EPSSW*(1.-FLAT2)
                0124         OZONE(J) = EPSSW*(1.+COZ1*SLAT(J)+COZ2*FLAT2)
                0125
                0126 C       zenith angle correction to (downward) absorptivity
                0127         ZENIT(J) = 1. + AZEN*
                0128      &    (1. _d 0-(CLAT(J)*CZEN+SLAT(J)*SZEN))**NZEN
                0129
                0130 C       ozone absorption in upper and lower stratosphere
                0131         OZUPP(J)=FSOL(J)*OZUPP(J)*ZENIT(J)
                0132         OZONE(J)=FSOL(J)*OZONE(J)*ZENIT(J)
                0133         STRATZ(J)=AST*FSOL(J)*CLAT(J)**3
                0134      &           +MAX( FS0-FSOL(J), 0. _d 0 )
                0135
                0136       ENDDO
                0137
7f98c35e47 Davi*0138 #ifdef ALLOW_DIAGNOSTICS
670d1925ca Jean*0139       IF ( useDiagnostics ) THEN
7f98c35e47 Davi*0140         CALL DIAGNOSTICS_FILL( FSOL,
                0141      &                        'FSOL    ', 1, 1, 3,bi,bj, myThid )
670d1925ca Jean*0142       ENDIF
7f98c35e47 Davi*0143 #endif
                0144
d676f916b2 Jean*0145 #endif /* ALLOW_AIM */
                0146       RETURN
                0147       END
                0148
670d1925ca Jean*0149 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
d676f916b2 Jean*0150
670d1925ca Jean*0151 CBOP
                0152 C !ROUTINE: RADSW (PSA,QA,RH,ALB,
                0153 C    &             ICLTOP,CLOUDC,FTOP,FSFC,DFABS)
                0154 C !INTERFACE:
d676f916b2 Jean*0155       SUBROUTINE RADSW (PSA,dpFac,QA,RH,ALB,
                0156      I                  FSOL, OZONE, OZUPP, ZENIT, STRATZ,
                0157      O                  TAU2, STRATC,
7f98c35e47 Davi*0158      O                  ICLTOP,CLOUDC,FTOP,FSFC,UPSWG,DFABS,
0d5086b5bf Jean*0159      I                  absCO2, kGrd,bi,bj,myThid)
670d1925ca Jean*0160 C !DESCRIPTION: \bv
                0161 C   Purpose: Compute the absorption of shortwave radiation and
                0162 C            initialize arrays for longwave-radiation routines
                0163 C   Input:   PSA    = norm. surface pressure [p/p0]           (2-dim)
                0164 C            dpFac  = cell delta_P fraction                   (3-dim)
                0165 C            QA     = specific humidity [g/kg]                (3-dim)
                0166 C            RH     = relative humidity                       (3-dim)
                0167 C            ALB    = surface albedo                          (2-dim)
                0168 C   Output:  ICLTOP = cloud top level                             (2-dim)
                0169 C            CLOUDC = total cloud cover                           (2-dim)
                0170 C            FTOP   = net downw. flux of sw rad. at the atm. top  (2-dim)
                0171 C            FSFC   = net downw. flux of sw rad. at the surface   (2-dim)
                0172 C            UPSWG  = upward flux of sw rad. at the surface       (2-dim)
                0173 C            DFABS  = flux of sw rad. absorbed by each atm. layer (3-dim)
0d5086b5bf Jean*0174 C  Input:    absCO2 = LW absorbtion in CO2 band (uniform value)
                0175 C            kGrd   = Ground level index                      (2-dim)
670d1925ca Jean*0176 C            bi,bj  = tile index
                0177 C            myThid = Thread number for this instance of the routine
                0178
                0179 C !USES:
d676f916b2 Jean*0180       IMPLICIT NONE
                0181
                0182 C     Resolution parameters
                0183
                0184 C-- size for MITgcm & Physics package :
                0185 #include "AIM_SIZE.h"
                0186
                0187 #include "EEPARAMS.h"
670d1925ca Jean*0188 #include "PARAMS.h"
d676f916b2 Jean*0189
                0190 C     Constants + functions of sigma and latitude
                0191 #include "com_physcon.h"
                0192 C     Radiation parameters
                0193 #include "com_radcon.h"
                0194
670d1925ca Jean*0195 C !INPUT/OUTPUT PARAMETERS:
b3097ed02d Jean*0196       _RL PSA(NGP),dpFac(NGP,NLEV),QA(NGP,NLEV),RH(NGP,NLEV)
                0197       _RL ALB(NGP,0:3)
d676f916b2 Jean*0198       INTEGER  ICLTOP(NGP)
b3097ed02d Jean*0199       _RL CLOUDC(NGP), FTOP(NGP), FSFC(NGP,0:3), DFABS(NGP,NLEV)
7f98c35e47 Davi*0200       _RL UPSWG(NGP)
d676f916b2 Jean*0201
                0202       _RL FSOL(NGP), OZONE(NGP), OZUPP(NGP), ZENIT(NGP), STRATZ(NGP)
                0203       _RL TAU2(NGP,NLEV,NBAND),STRATC(NGP)
                0204 c     _RL FLUX(NGP,4)
                0205
0d5086b5bf Jean*0206       _RL absCO2
d676f916b2 Jean*0207       INTEGER  kGrd(NGP)
                0208       INTEGER  bi, bj, myThid
670d1925ca Jean*0209 CEOP
d676f916b2 Jean*0210
                0211 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0212
                0213 #ifdef ALLOW_AIM
670d1925ca Jean*0214 C !LOCAL VARIABLES:
d676f916b2 Jean*0215 c     _RL  QCLOUD(NGP), ACLOUD(NGP), PSAZ(NGP),
                0216       _RL  QCLOUD(NGP), ACLOUD(NGP),
                0217      &     ALBTOP(NGP,NLEV), FREFL(NGP,NLEV), FLUX(NGP,2)
bd70561229 Davi*0218 #ifdef ALLOW_CLOUD_3D
670d1925ca Jean*0219 C_DE       CLDCLW = Local cloud cover                           (3-dim)
bd70561229 Davi*0220       _RL CLDCLW(NGP,NLEV), ACLDLW(NGP,NLEV)
                0221 #endif
d676f916b2 Jean*0222
670d1925ca Jean*0223 C-  jmc: define "FLUX" as a local variable & remove Equivalences:
d676f916b2 Jean*0224 c     EQUIVALENCE (ALBTOP(1,1),TAU2(1,1,3))
                0225 c     EQUIVALENCE ( FREFL(1,1),TAU2(1,1,4))
                0226       INTEGER NL1(NGP)
                0227       INTEGER K, J
e749d70ece Jean*0228       LOGICAL makeClouds
d676f916b2 Jean*0229       _RL FBAND1, FBAND2, RRCL, RQCL, DQACL, QACL3
                0230       _RL ABS1, DELTAP
670d1925ca Jean*0231
d676f916b2 Jean*0232 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0233
                0234       FBAND2=0.05 _d 0
                0235       FBAND1=1.-FBAND2
                0236
                0237       DO J=1,NGP
                0238         NL1(J)=kGrd(J)-1
                0239       ENDDO
                0240
                0241 C--   1.  Cloud cover:
                0242 C         defined as a linear fun. of the maximum relative humidity
                0243 C         in all tropospheric layers above PBL:
                0244 C         CLOUDC =  0 for RHmax < RHCL1, = 1 for RHmax > RHCL2.
                0245 C         This value is reduced by a factor (Qbase/QACL) if the
                0246 C         cloud-base absolute humidity Qbase < QACL.
                0247 C
e749d70ece Jean*0248       makeClouds = ICLTOP(1) .GE. 0
d676f916b2 Jean*0249       RRCL=1./(RHCL2-RHCL1)
                0250       RQCL=1./QACL2
                0251 C
                0252       DO J=1,NGP
                0253         CLOUDC(J)=0.
                0254         QCLOUD(J)=0.
e749d70ece Jean*0255         ICLTOP(J)=NLEV+1
d676f916b2 Jean*0256         FREFL(J,1)=0.
                0257       ENDDO
                0258
                0259       DO K=1,NLEV
                0260        DO J=1,NGP
                0261          ALBTOP(J,K)=0.
bd70561229 Davi*0262 #ifdef ALLOW_CLOUD_3D
                0263          CLDCLW(J,K)=0.
                0264 #endif
d676f916b2 Jean*0265        ENDDO
                0266       ENDDO
                0267 C
e749d70ece Jean*0268       IF ( makeClouds ) THEN
                0269 C-    skipp this part for clear-sky diagnostics
                0270
                0271        DQACL=(QACL2-QACL1)/(0.5 _d 0 - SIG(2))
                0272        DO J=1,NGP
                0273         ICLTOP(J)= kGrd(J)
                0274         DO K=NL1(J),2,-1
d676f916b2 Jean*0275          QACL3=MIN(QACL2,QACL1+DQACL*(SIG(K)-SIG(2)))
                0276          IF (RH(J,K).GT.RHCL1.AND.QA(J,K).GT.QACL1) THEN
                0277             CLOUDC(J)=MAX(CLOUDC(J),RH(J,K)-RHCL1)
                0278             IF (QA(J,K).GT.QACL3) ICLTOP(J)=K
bd70561229 Davi*0279 #ifdef ALLOW_CLOUD_3D
                0280             CLDCLW(J,K)=MAX(0. _d 0,RH(J,K)-RHCL1)
                0281             CLDCLW(J,K)=MIN(1. _d 0,CLDCLW(J,K)*RRCL)
                0282 #endif
d676f916b2 Jean*0283          ENDIF
e749d70ece Jean*0284         ENDDO
d676f916b2 Jean*0285        ENDDO
                0286
e749d70ece Jean*0287        DO J=1,NGP
                0288         IF (kGrd(J).NE.0)
                0289      &  QCLOUD(J)= MAX( QA(J,kGrd(J)), QA(J,NL1(J)) )
d676f916b2 Jean*0290         CLOUDC(J)=MIN(1. _d 0,CLOUDC(J)*RRCL)
                0291         IF (CLOUDC(J).GT.0.0) THEN
                0292           CLOUDC(J)=CLOUDC(J)*MIN(1. _d 0,QCLOUD(J)*RQCL)
bd70561229 Davi*0293 #ifdef ALLOW_CLOUD_3D
                0294           DO K=NL1(J),2,-1
                0295              CLDCLW(J,K)=CLDCLW(J,K)*MIN(1. _d 0,QCLOUD(J)*RQCL)
                0296           ENDDO
                0297 #endif
d676f916b2 Jean*0298           ALBTOP(J,ICLTOP(J))=ALBCL*CLOUDC(J)
                0299         ELSE
                0300           ICLTOP(J)=NLEV+1
                0301         ENDIF
e749d70ece Jean*0302        ENDDO
                0303
                0304 C-    end if makeClouds
                0305       ENDIF
d676f916b2 Jean*0306
                0307 C
                0308 C--   2. Shortwave transmissivity:
                0309 C        function of layer mass, ozone (in the statosphere),
                0310 C        abs. humidity and cloud cover (in the troposphere)
                0311
                0312       DO J=1,NGP
                0313 c_FM    PSAZ(J)=PSA(J)*ZENIT(J)
                0314         ACLOUD(J)=CLOUDC(J)*(ABSCL1+ABSCL2*QCLOUD(J))
                0315       ENDDO
                0316
                0317       DO J=1,NGP
                0318 c_FM    DELTAP=PSAZ(J)*DSIG(1)
                0319         DELTAP=ZENIT(J)*DSIG(1)*dpFac(J,1)
                0320         TAU2(J,1,1)=EXP(-DELTAP*ABSDRY)
                0321       ENDDO
                0322 C
                0323       DO J=1,NGP
                0324        DO K=2,NL1(J)
                0325 c_FM     ABS1=ABSDRY+ABSAER*SIG(K)**2
                0326 c_FM     DELTAP=PSAZ(J)*DSIG(K)
                0327          ABS1=ABSDRY+ABSAER*(SIG(K)/PSA(J))**2
                0328          DELTAP=ZENIT(J)*DSIG(K)*dpFac(J,K)
                0329          IF (K.EQ.ICLTOP(J)) THEN
                0330            TAU2(J,K,1)=EXP(-DELTAP*
                0331      &                 (ABS1+ABSWV1*QA(J,K)+2.*ACLOUD(J)))
                0332          ELSE IF (K.GT.ICLTOP(J)) THEN
                0333            TAU2(J,K,1)=EXP(-DELTAP*
                0334      &                 (ABS1+ABSWV1*QA(J,K)+ACLOUD(J)))
                0335          ELSE
                0336            TAU2(J,K,1)=EXP(-DELTAP*(ABS1+ABSWV1*QA(J,K)))
                0337          ENDIF
                0338        ENDDO
                0339       ENDDO
                0340
                0341 c_FM  ABS1=ABSDRY+ABSAER*SIG(NLEV)**2
                0342       DO J=1,NGP
                0343        K = kGrd(J)
                0344        ABS1=ABSDRY+ABSAER*(SIG(K)/PSA(J))**2
                0345 c_FM    DELTAP=PSAZ(J)*DSIG(NLEV)
                0346         DELTAP=ZENIT(J)*DSIG(K)*dpFac(J,K)
                0347         TAU2(J,K,1)=EXP(-DELTAP*(ABS1+ABSWV1*QA(J,K)))
                0348       ENDDO
                0349
                0350       DO J=1,NGP
                0351        DO K=2,kGrd(J)
                0352          DELTAP=ZENIT(J)*DSIG(K)*dpFac(J,K)
                0353          TAU2(J,K,2)=EXP(-DELTAP*ABSWV2*QA(J,K))
                0354        ENDDO
                0355       ENDDO
                0356 C
                0357 C---  3. Shortwave downward flux
                0358 C
                0359 C     3.1  Absorption in the stratosphere
                0360
                0361 C     3.1.1 Initialization of fluxes (subtracting
                0362 C           ozone absorption in the upper stratosphere)
                0363
                0364       DO J=1,NGP
                0365         FTOP(J)  =FSOL(J)
                0366         FLUX(J,1)=FSOL(J)*FBAND1-OZUPP(J)
                0367         FLUX(J,2)=FSOL(J)*FBAND2
                0368         STRATC(J)=STRATZ(J)*PSA(J)
                0369       ENDDO
                0370
                0371 C     3.1.2 Ozone and dry-air absorption
                0372 C           in the lower (modelled) stratosphere
                0373
                0374       DO J=1,NGP
                0375         DFABS(J,1)=FLUX(J,1)
                0376         FLUX (J,1)=TAU2(J,1,1)*(FLUX(J,1)-OZONE(J)*PSA(J))
                0377         DFABS(J,1)=DFABS(J,1)-FLUX(J,1)
                0378       ENDDO
                0379
                0380 C     3.3  Absorption and reflection in the troposphere
                0381 C
                0382       DO J=1,NGP
                0383        DO K=2,kGrd(J)
                0384          FREFL(J,K)=FLUX(J,1)*ALBTOP(J,K)
                0385          FLUX (J,1)=FLUX(J,1)-FREFL(J,K)
                0386          DFABS(J,K)=FLUX(J,1)
                0387          FLUX (J,1)=TAU2(J,K,1)*FLUX(J,1)
                0388          DFABS(J,K)=DFABS(J,K)-FLUX(J,1)
                0389        ENDDO
                0390       ENDDO
                0391
                0392       DO J=1,NGP
                0393        DO K=2,kGrd(J)
                0394          DFABS(J,K)=DFABS(J,K)+FLUX(J,2)
                0395          FLUX (J,2)=TAU2(J,K,2)*FLUX(J,2)
                0396          DFABS(J,K)=DFABS(J,K)-FLUX(J,2)
                0397        ENDDO
                0398       ENDDO
                0399
                0400 C
                0401 C---  4. Shortwave upward flux
                0402 C
                0403 C     4.1  Absorption and reflection at the surface
                0404 C
                0405       DO J=1,NGP
b3097ed02d Jean*0406 C      for each surface type:
                0407         FSFC(J,1)=FLUX(J,1)*(1.-ALB(J,1))+FLUX(J,2)
                0408         FSFC(J,2)=FLUX(J,1)*(1.-ALB(J,2))+FLUX(J,2)
                0409         FSFC(J,3)=FLUX(J,1)*(1.-ALB(J,3))+FLUX(J,2)
                0410 C      weighted average according to land/sea/sea-ice fraction:
                0411         FSFC(J,0)=FLUX(J,1)+FLUX(J,2)
                0412         FLUX(J,1)=FLUX(J,1)*ALB(J,0)
                0413         FSFC(J,0)=FSFC(J,0)-FLUX(J,1)
d676f916b2 Jean*0414       ENDDO
7f98c35e47 Davi*0415
                0416 C--   Store upward shortwave flux at the surface for diagnostics purpose
                0417       DO J=1,NGP
                0418         UPSWG(J)=FLUX(J,1)
                0419       ENDDO
d676f916b2 Jean*0420 C
                0421 C     4.2  Absorption of upward flux
                0422 C
78542777fd Davi*0423       DO K=NLEV,1,-1
                0424        DO J=1,NGP
                0425         IF ( K .LE. kGrd(J) ) THEN
d676f916b2 Jean*0426          DFABS(J,K)=DFABS(J,K)+FLUX(J,1)
                0427          FLUX (J,1)=TAU2(J,K,1)*FLUX(J,1)
                0428          DFABS(J,K)=DFABS(J,K)-FLUX(J,1)
                0429          FLUX (J,1)=FLUX(J,1)+FREFL(J,K)
78542777fd Davi*0430         ELSE
                0431          DFABS(J,K)= 0. _d 0
                0432         ENDIF
d676f916b2 Jean*0433        ENDDO
                0434       ENDDO
                0435 C
                0436 C     4.3  Net solar radiation = incoming - outgoing
                0437 C
                0438       DO J=1,NGP
                0439         FTOP(J)=FTOP(J)-FLUX(J,1)
                0440       ENDDO
                0441
                0442 C
                0443 C---  5.  Initialization of longwave radiation model
                0444 C
                0445 C     5.1  Longwave transmissivity:
                0446 C          function of layer mass, abs. humidity and cloud cover.
                0447
bd70561229 Davi*0448 #ifdef ALLOW_CLOUD_3D
                0449       DO K=2,NLEV
                0450         DO J=1,NGP
                0451           ACLDLW(J,K)=CLDCLW(J,K)*(ABLCL1+ABLCL2*QCLOUD(J))
                0452         ENDDO
                0453       ENDDO
                0454 #else
d676f916b2 Jean*0455       DO J=1,NGP
                0456         ACLOUD(J)=CLOUDC(J)*(ABLCL1+ABLCL2*QCLOUD(J))
                0457       ENDDO
bd70561229 Davi*0458 #endif
d676f916b2 Jean*0459
                0460       DO J=1,NGP
                0461 c_FM    DELTAP=PSA(J)*DSIG(1)
                0462         DELTAP=DSIG(1)*dpFac(J,1)
                0463         TAU2(J,1,1)=EXP(-DELTAP*ABLWIN)
0d5086b5bf Jean*0464         TAU2(J,1,2)=EXP(-DELTAP*absCO2)
d676f916b2 Jean*0465         TAU2(J,1,3)=1.
                0466         TAU2(J,1,4)=1.
                0467       ENDDO
                0468
                0469       DO K=2,NLEV
                0470        DO J=1,NGP
                0471 c_FM     DELTAP=PSA(J)*DSIG(K)
                0472          DELTAP=DSIG(K)*dpFac(J,K)
                0473          IF ( K.GE.ICLTOP(J).AND.K.NE.kGrd(J) ) THEN
bd70561229 Davi*0474 #ifdef ALLOW_CLOUD_3D
                0475            TAU2(J,K,1)=EXP(-DELTAP*(ABLWIN+ACLDLW(J,K)))
                0476 #else
d676f916b2 Jean*0477            TAU2(J,K,1)=EXP(-DELTAP*(ABLWIN+ACLOUD(J)))
bd70561229 Davi*0478 #endif
d676f916b2 Jean*0479          ELSE
                0480            TAU2(J,K,1)=EXP(-DELTAP*ABLWIN)
                0481          ENDIF
0d5086b5bf Jean*0482          TAU2(J,K,2)=EXP(-DELTAP*absCO2)
d676f916b2 Jean*0483          TAU2(J,K,3)=EXP(-DELTAP*ABLWV1*QA(J,K))
                0484          TAU2(J,K,4)=EXP(-DELTAP*ABLWV2*QA(J,K))
                0485        ENDDO
                0486       ENDDO
bd70561229 Davi*0487
                0488 #if (defined ALLOW_CLOUD_3D && defined ALLOW_DIAGNOSTICS)
670d1925ca Jean*0489       IF ( useDiagnostics ) THEN
bd70561229 Davi*0490         CALL DIAGNOSTICS_FILL( CLDCLW,
                0491      &                        'CLDCLW  ',-1, Nr, 3,bi,bj, myThid )
670d1925ca Jean*0492       ENDIF
bd70561229 Davi*0493 #endif
                0494
d676f916b2 Jean*0495 #endif /* ALLOW_AIM */
                0496
                0497       RETURN
                0498       END
                0499
670d1925ca Jean*0500 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
d676f916b2 Jean*0501
670d1925ca Jean*0502 CBOP
                0503 C !ROUTINE: RADLW (IMODE,TA,TS,ST4S,
                0504 C    &             FTOP,FSFC,DFABS)
                0505 C !INTERFACE:
d676f916b2 Jean*0506       SUBROUTINE RADLW (IMODE,TA,TS,ST4S,
                0507      I                  OZUPP, STRATC, TAU2,
                0508      &                  FLUX, ST4A,
                0509      &                  FTOP,FSFC,DFABS,
                0510      I                  kGrd,bi,bj,myThid)
                0511
670d1925ca Jean*0512 C !DESCRIPTION: \bv
                0513 C   Purpose: Compute the absorption of longwave radiation
                0514 C   Input:   IMODE  = index for operation mode
                0515 C                     -1 : downward flux only
                0516 C                      0 : downward + upward flux
                0517 C                     +1 : upward flux only
                0518 C            TA     = absolute temperature (3-dim)
                0519 C            TS     = surface temperature                  [if IMODE=0,1]
                0520 C            ST4S   = surface blackbody emission             [if IMODE=1]
                0521 C            FSFC   = FSFC  output from RADLW(-1,... )       [if IMODE=1]
                0522 C            DFABS  = DFABS output from RADLW(-1,... )       [if IMODE=1]
                0523 C   Output:  ST4S   = surface blackbody emission             [if IMODE=0]
                0524 C            FTOP   = outgoing flux of lw rad. at the top  [if IMODE=0,1]
                0525 C            FSFC   = downward flux of lw rad. at the sfc. [if IMODE= -1]
                0526 C                     net upw. flux of lw rad. at the sfc. [if IMODE=0,1]
                0527 C            DFABS  = flux of lw rad. absorbed by each atm. layer (3-dim)
                0528 C  Input:    kGrd   = Ground level index                      (2-dim)
                0529 C            bi,bj  = tile index
                0530 C            myThid = Thread number for this instance of the routine
                0531
                0532 C !USES:
d676f916b2 Jean*0533       IMPLICIT NONE
                0534
                0535 C     Resolution parameters
                0536 C-- size for MITgcm & Physics package :
                0537 #include "AIM_SIZE.h"
                0538 #include "EEPARAMS.h"
                0539
                0540 C     Number of radiation bands with tau < 1
                0541 c     INTEGER NBAND
                0542 c     PARAMETER ( NBAND=4 )
                0543 C     Constants + functions of sigma and latitude
                0544 #include "com_physcon.h"
                0545 C     Radiation parameters
                0546 #include "com_radcon.h"
                0547
670d1925ca Jean*0548 C !INPUT/OUTPUT PARAMETERS:
d676f916b2 Jean*0549       INTEGER IMODE
                0550       _RL TA(NGP,NLEV), TS(NGP), ST4S(NGP)
                0551       _RL FTOP(NGP), FSFC(NGP), DFABS(NGP,NLEV)
                0552       _RL OZUPP(NGP), STRATC(NGP)
                0553       _RL TAU2(NGP,NLEV,NBAND), FLUX(NGP,NBAND), ST4A(NGP,NLEV,2)
                0554
                0555       INTEGER kGrd(NGP)
                0556       INTEGER bi,bj,myThid
670d1925ca Jean*0557 CEOP
d676f916b2 Jean*0558
                0559 #ifdef ALLOW_AIM
670d1925ca Jean*0560 C !LOCAL VARIABLES:
d676f916b2 Jean*0561       INTEGER K, J, JB
                0562 c     INTEGER J0, Jl, I2
                0563       INTEGER NL1(NGP)
                0564       _RL REFSFC, BRAD, EMIS
670d1925ca Jean*0565
d676f916b2 Jean*0566 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0567
                0568       DO J=1,NGP
                0569         NL1(J)=kGrd(J)-1
                0570       ENDDO
                0571
                0572       REFSFC=1.-EMISFC
                0573
                0574       IF (IMODE.EQ.1) GO TO 410
                0575
                0576 C---  1. Blackbody emission from atmospheric full and half levels.
                0577 C        Temperature is interpolated as a linear function of ln sigma.
                0578 C        At the lower boundary, the emission is linearly extrapolated;
                0579 C        at the upper boundary, the atmosphere is assumed isothermal.
                0580
                0581       DO K=1,NLEV
                0582        DO J=1,NGP
                0583          ST4A(J,K,1)=TA(J,K)*TA(J,K)
                0584          ST4A(J,K,1)=SBC*ST4A(J,K,1)*ST4A(J,K,1)
                0585        ENDDO
                0586       ENDDO
                0587 C
                0588       DO K=1,NLEV-1
                0589        DO J=1,NGP
                0590          ST4A(J,K,2)=TA(J,K)+WVI(K,2)*(TA(J,K+1)-TA(J,K))
                0591          ST4A(J,K,2)=ST4A(J,K,2)*ST4A(J,K,2)
                0592          ST4A(J,K,2)=SBC*ST4A(J,K,2)*ST4A(J,K,2)
                0593        ENDDO
                0594       ENDDO
                0595 C
                0596       DO J=1,NGP
                0597 c       ST4A(J,NLEV,2)=ST4A(J,NLEV,1)
                0598         K=kGrd(J)
                0599         ST4A(J,K,2)=2.*ST4A(J,K,1)-ST4A(J,NL1(J),2)
                0600       ENDDO
                0601
                0602 C---  2. Initialization
                0603 C---     (including the stratospheric correction term)
                0604
                0605       DO J=1,NGP
                0606         FTOP(J)   = 0.
                0607         FSFC(J)   = STRATC(J)
                0608         DFABS(J,1)=-STRATC(J)
                0609       ENDDO
                0610
                0611       DO K=2,NLEV
                0612        DO J=1,NGP
                0613          DFABS(J,K)=0.
                0614        ENDDO
                0615       ENDDO
                0616
                0617 C---  3. Emission ad absorption of longwave downward flux.
                0618 C        Downward emission is an average of the emission from the full level
                0619 C        and the half-level below, weighted according to the transmissivity
                0620 C        of the layer.
                0621
                0622 C     3.1  Stratosphere
                0623
                0624       K=1
                0625       DO JB=1,2
                0626        DO J=1,NGP
                0627          BRAD=ST4A(J,K,2)+TAU2(J,K,JB)*(ST4A(J,K,1)-ST4A(J,K,2))
                0628          EMIS=FBAND(NINT(TA(J,K)),JB)*(1.-TAU2(J,K,JB))
                0629          FLUX(J,JB)=EMIS*BRAD
                0630          DFABS(J,K)=DFABS(J,K)-FLUX(J,JB)
                0631        ENDDO
                0632       ENDDO
                0633
                0634       DO JB=3,NBAND
                0635        DO J=1,NGP
                0636          FLUX(J,JB)=0.
                0637        ENDDO
                0638       ENDDO
                0639
                0640 C     3.2  Troposphere
                0641
                0642       DO JB=1,NBAND
                0643        DO J=1,NGP
                0644         DO K=2,kGrd(J)
                0645           BRAD=ST4A(J,K,2)+TAU2(J,K,JB)*(ST4A(J,K,1)-ST4A(J,K,2))
                0646           EMIS=FBAND(NINT(TA(J,K)),JB)*(1.-TAU2(J,K,JB))
                0647           DFABS(J,K)=DFABS(J,K)+FLUX(J,JB)
                0648           FLUX(J,JB)=TAU2(J,K,JB)*FLUX(J,JB)+EMIS*BRAD
                0649           DFABS(J,K)=DFABS(J,K)-FLUX(J,JB)
                0650         ENDDO
                0651        ENDDO
                0652       ENDDO
                0653
                0654       DO JB=1,NBAND
                0655        DO J=1,NGP
                0656          FSFC(J)=FSFC(J)+EMISFC*FLUX(J,JB)
                0657        ENDDO
                0658       ENDDO
                0659
                0660       IF (IMODE.EQ.-1) RETURN
                0661
                0662 C---  4. Emission ad absorption of longwave upward flux
                0663 C        Upward emission is an average of the emission from the full level
                0664 C        and the half-level above, weighted according to the transmissivity
                0665 C        of the layer (for the top layer, full-level emission is used).
                0666 C        Surface lw emission in "band 0" goes directly into FTOP.
                0667
                0668 C     4.1  Surface
                0669
                0670       DO J=1,NGP
                0671         ST4S(J)=TS(J)*TS(J)
                0672         ST4S(J)=SBC*ST4S(J)*ST4S(J)
e749d70ece Jean*0673         ST4S(J)=EMISFC*ST4S(J)
d676f916b2 Jean*0674       ENDDO
                0675
                0676 C     Entry point for upward-only mode (IMODE=1)
                0677  410  CONTINUE
                0678
                0679       DO J=1,NGP
                0680         FSFC(J)=ST4S(J)-FSFC(J)
                0681         FTOP(J)=FTOP(J)+FBAND(NINT(TS(J)),0)*ST4S(J)
                0682       ENDDO
                0683
                0684       DO JB=1,NBAND
                0685        DO J=1,NGP
                0686          FLUX(J,JB)=FBAND(NINT(TS(J)),JB)*ST4S(J)
                0687      &              +REFSFC*FLUX(J,JB)
                0688        ENDDO
                0689       ENDDO
                0690
                0691 C     4.2  Troposphere
                0692
                0693       DO JB=1,NBAND
                0694        DO J=1,NGP
                0695         DO K=kGrd(J),2,-1
                0696           BRAD=ST4A(J,K-1,2)+TAU2(J,K,JB)*(ST4A(J,K,1)-ST4A(J,K-1,2))
                0697           EMIS=FBAND(NINT(TA(J,K)),JB)*(1.-TAU2(J,K,JB))
                0698           DFABS(J,K)=DFABS(J,K)+FLUX(J,JB)
                0699           FLUX(J,JB)=TAU2(J,K,JB)*FLUX(J,JB)+EMIS*BRAD
                0700           DFABS(J,K)=DFABS(J,K)-FLUX(J,JB)
                0701         ENDDO
                0702        ENDDO
                0703       ENDDO
                0704
                0705 C     4.3  Stratosphere
                0706
                0707       K=1
                0708       DO JB=1,2
                0709        DO J=1,NGP
                0710          EMIS=FBAND(NINT(TA(J,K)),JB)*(1.-TAU2(J,K,JB))
                0711          DFABS(J,K)=DFABS(J,K)+FLUX(J,JB)
                0712          FLUX(J,JB)=TAU2(J,K,JB)*FLUX(J,JB)+EMIS*ST4A(J,K,1)
                0713          DFABS(J,K)=DFABS(J,K)-FLUX(J,JB)
                0714        ENDDO
                0715       ENDDO
                0716
                0717 C     4.4  Outgoing longwave radiation
                0718
                0719       DO JB=1,NBAND
                0720        DO J=1,NGP
                0721          FTOP(J)=FTOP(J)+FLUX(J,JB)
                0722        ENDDO
                0723       ENDDO
                0724
                0725       DO J=1,NGP
                0726         FTOP(J)=FTOP(J)+OZUPP(J)
                0727       ENDDO
                0728
                0729 #endif /* ALLOW_AIM */
                0730
                0731       RETURN
                0732       END
                0733
670d1925ca Jean*0734 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
d676f916b2 Jean*0735
670d1925ca Jean*0736 CBOP
                0737 C !ROUTINE: RADSET
                0738 C !INTERFACE:
d676f916b2 Jean*0739       SUBROUTINE RADSET( myThid )
                0740
670d1925ca Jean*0741 C !DESCRIPTION: \bv
                0742 C   Purpose: compute energy fractions in LW bands
                0743 C            as a function of temperature
                0744 C   Initialized common blocks: RADFIX
d676f916b2 Jean*0745
670d1925ca Jean*0746 C !USES:
d676f916b2 Jean*0747       IMPLICIT NONE
                0748
                0749 C     Resolution parameters
                0750 C-- size for MITgcm & Physics package :
                0751 #include "AIM_SIZE.h"
                0752 #include "EEPARAMS.h"
                0753
                0754 C     Radiation constants
                0755 #include "com_radcon.h"
                0756
670d1925ca Jean*0757 C !INPUT/OUTPUT PARAMETERS:
d676f916b2 Jean*0758       INTEGER myThid
670d1925ca Jean*0759 CEOP
d676f916b2 Jean*0760
                0761 #ifdef ALLOW_AIM
670d1925ca Jean*0762 C !LOCAL VARIABLES:
d676f916b2 Jean*0763       INTEGER JTEMP, JB
                0764       _RL EPS3
                0765
                0766 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0767
                0768       EPS3=0.95 _d 0
                0769
                0770       DO JTEMP=200,320
                0771         FBAND(JTEMP,0)= EPSLW
                0772         FBAND(JTEMP,2)= 0.148 _d 0 - 3.0 _d -6 *(JTEMP-247)**2
                0773         FBAND(JTEMP,3)=(0.375 _d 0 - 5.5 _d -6 *(JTEMP-282)**2)*EPS3
                0774         FBAND(JTEMP,4)= 0.314 _d 0 + 1.0 _d -5 *(JTEMP-315)**2
670d1925ca Jean*0775         FBAND(JTEMP,1)= 1. _d 0 -(FBAND(JTEMP,0)+FBAND(JTEMP,2)
d676f916b2 Jean*0776      &                           +FBAND(JTEMP,3)+FBAND(JTEMP,4))
                0777       ENDDO
                0778
                0779       DO JB=0,NBAND
                0780         DO JTEMP=lwTemp1,199
                0781           FBAND(JTEMP,JB)=FBAND(200,JB)
                0782         ENDDO
                0783         DO JTEMP=321,lwTemp2
                0784           FBAND(JTEMP,JB)=FBAND(320,JB)
                0785         ENDDO
                0786       ENDDO
                0787
                0788 #endif /* ALLOW_AIM */
                0789
                0790       RETURN
                0791       END