pkg/radtrans/radtrans_calc.F

b55e95f1ff Oliv*0001 #include "RADTRANS_OPTIONS.h"
                0002 #ifdef ALLOW_SEAICE
                0003 #include "SEAICE_OPTIONS.h"
                0004 #endif
                0005
                0006 CBOP
                0007 C !ROUTINE: RADTRANS_CALC
                0008 C !INTERFACE: ==========================================================
                0009       SUBROUTINE RADTRANS_CALC(
                0010      I                         a, bt, bb,
                0011      O                         E0F,
                0012      I                         bi, bj, iMin, iMax, jMin, jMax,
                0013      I                         myTime, myIter, myThid )
                0014
                0015 C !DESCRIPTION:
                0016
                0017 C !USES: ===============================================================
                0018       IMPLICIT NONE
                0019 #include "SIZE.h"
                0020 #include "GRID.h"
                0021 #include "EEPARAMS.h"
                0022 #include "PARAMS.h"
0222db53b0 Oliv*0023 #ifdef ALLOW_OASIM
                0024 #include "OASIM_SIZE.h"
                0025 #include "OASIM_FIELDS.h"
                0026 #endif
b55e95f1ff Oliv*0027 #include "RADTRANS_SIZE.h"
                0028 #include "RADTRANS_PARAMS.h"
                0029 #include "RADTRANS_FIELDS.h"
                0030 #ifdef ALLOW_SEAICE
                0031 #include "SEAICE_SIZE.h"
                0032 #include "SEAICE.h"
                0033 #endif
                0034
                0035 C !INPUT PARAMETERS: ===================================================
                0036 C  myTime :: time at end of (sub)timestep
                0037 C  myThid :: thread number
                0038       _RL a(1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr, nlam)
                0039       _RL bt(1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr, nlam)
                0040       _RL bb(1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr, nlam)
                0041       _RL myTime
                0042       INTEGER bi, bj, iMin, iMax, jMin, jMax, myIter, myThid
                0043
                0044 C !OUTPUT PARAMETERS: =================================================
                0045 C  E0F :: spectral scalar irradiance at top of grid cell
                0046       _RL E0F(1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr+1, nlam)
                0047 CEOP
                0048
                0049 #ifdef ALLOW_RADTRANS
                0050
                0051 C!LOCAL VARIABLES: ====================================================
                0052 C  i,j                  :: loop indices
                0053 C  k                    :: vertical level
                0054       LOGICAL  DIAGNOSTICS_IS_ON
                0055       EXTERNAL DIAGNOSTICS_IS_ON
                0056       CHARACTER*8 diagname
                0057       INTEGER i,j,k,l,jp,klow
                0058       INTEGER iyr,imon,iday,isec,lp,wd,mydate(4)
                0059       _RL rad2deg
                0060       _RL delta
                0061       _RL solz(1-OLx:sNx+OLx, 1-OLy:sNy+OLy)
                0062       _RL rmud(1-OLx:sNx+OLx, 1-OLy:sNy+OLy)
                0063       _RL dz_k(Nr)
                0064       _RL a_k(Nr)
                0065       _RL bt_k(Nr)
                0066       _RL bb_k(Nr)
                0067       _RL Edwsf,Eswsf,Edown
                0068       _RL Edbot(Nr),Esbot(Nr),Eubot(Nr)
                0069       _RL Estop(Nr),Eutop(Nr)
                0070       _RL amp1_k(Nr), amp2_k(Nr)
                0071       _RL x_k(Nr), y_k(Nr)
                0072       _RL r1_k(Nr), r2_k(Nr)
                0073       _RL kappa1_k(Nr), kappa2_k(Nr)
0222db53b0 Oliv*0074 #ifdef ALLOW_OASIM
                0075       INTEGER loa
                0076 #endif
b55e95f1ff Oliv*0077
                0078       _RL Ed(1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr, nlam)
                0079       _RL Es(1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr, nlam)
                0080       _RL Eu(1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr, nlam)
                0081       _RL Ef(1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr, nlam)
                0082
                0083 #ifdef ALLOW_DIAGNOSTICS
                0084       _RL Rirr(1-OLx:sNx+OLx, 1-OLy:sNy+OLy, nlam)
                0085 #ifdef RADTRANS_DIAG_SOLUTION
                0086       _RL Est(1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr, nlam)
                0087       _RL Eub(1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr, nlam)
                0088       _RL amp1(1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr, nlam)
                0089       _RL amp2(1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr, nlam)
                0090       _RL x3d(1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr, nlam)
                0091       _RL y3d(1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr, nlam)
                0092       _RL r1(1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr, nlam)
                0093       _RL r2(1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr, nlam)
                0094       _RL kap1(1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr, nlam)
                0095       _RL kap2(1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr, nlam)
                0096 #endif
                0097 #endif
                0098
                0099       rad2deg = 180 _d 0 / PI
                0100
                0101 C ======================================================================
                0102
                0103 #ifdef ALLOW_SEAICE
                0104       IF (RT_useSEAICE) THEN
                0105        DO j=jMin,jMax
                0106         DO i=iMin,iMax
                0107          RT_iceFrac(i,j,bi,bj) = AREA(i,j,bi,bj)
                0108         ENDDO
                0109        ENDDO
                0110       ENDIF
                0111 #endif
                0112
                0113 C ======================================================================
                0114 C--   compute solar zenith angles
                0115
0222db53b0 Oliv*0116 C     OASIM computes its own
                0117       IF (.NOT.(useOASIM .AND. RT_useOASIMrmud)) THEN
                0118
                0119        IF (RT_useMeanCosSolz) THEN
                0120 C      from day-time average of cosine of zenith angle
b55e95f1ff Oliv*0121
                0122         CALL RADTRANS_DECLINATION_SPENCER(delta, myTime, myIter, myThid)
                0123         CALL RADTRANS_SOLZ_DAYTIME(solz, delta, bi, bj,
                0124      &                             iMin, iMax, jMin, jMax, myThid)
                0125
0222db53b0 Oliv*0126        ELSE
                0127 C      compute zenith angle at local noon
b55e95f1ff Oliv*0128
e8b7a8749d Oliv*0129         IF (RT_useNoonSolz) THEN
                0130          isec = 12*3600
                0131         ELSE
                0132          isec = -1
                0133         ENDIF
1c7c72345b Oliv*0134 #ifdef ALLOW_SUN
b55e95f1ff Oliv*0135         CALL SUN_SFCSOLZ(
                0136      O                    solz,
                0137      I                    isec, bi, bj, iMin, iMax, jMin, jMax,
                0138      I                    myTime, myIter, myThid )
1c7c72345b Oliv*0139 #endif
b55e95f1ff Oliv*0140
0222db53b0 Oliv*0141        ENDIF
b55e95f1ff Oliv*0142
0222db53b0 Oliv*0143 c      compute 1/cos(zenith) for direct light below surface
                0144        CALL RADTRANS_RMUD_BELOW(rmud,solz,iMin,iMax,jMin,jMax,myThid)
                0145
                0146 C     not useOASIM
                0147       ENDIF
b55e95f1ff Oliv*0148
                0149 C ======================================================================
                0150       DO j=jMin,jMax
                0151        DO i=iMin,iMax
                0152 C ----------------------------------------------------------------------
                0153         DO k=1,Nr
                0154          dz_k(k) = drF(k)*HFacC(i,j,k,bi,bj)
                0155         ENDDO
                0156
                0157         klow = MIN(RT_kmax, kLowC(i,j,bi,bj))
                0158
                0159         DO l = 1,nlam
                0160
0222db53b0 Oliv*0161 #ifdef ALLOW_OASIM
                0162          IF (useOASIM) THEN
                0163
                0164 C use light computed in the oasim package
                0165           DO loa = 1,nlt
                0166            IF (nint(RT_wbRefWLs(l)) .EQ. oasim_lam(loa)) EXIT
                0167           ENDDO
                0168           IF (loa .gt. nlt) THEN
                0169             print*,'RADTRANS_CALC: waveband not found in OASIM:',
                0170      &             RT_wbRefWLs(l)
                0171             STOP 'RADTRANS_CALC: waveband not found in OASIM'
                0172           ENDIF
                0173           Edwsf = oasim_edbelow(i,j,bi,bj,loa)*RT_oasimWgt(l)
                0174           Eswsf = oasim_esbelow(i,j,bi,bj,loa)*RT_oasimWgt(l)
                0175
                0176           IF (RT_useOASIMrmud) THEN
                0177            rmud(i,j) = OASIM_rmud(i,j,bi,bj)
                0178           ENDIF
                0179
                0180          ELSE
                0181 #else
                0182          IF (.TRUE.) THEN
                0183 #endif /* ALLOW_OASIM */
                0184
b55e95f1ff Oliv*0185 C use read-in light
                0186           Edwsf = RT_Ed_sfc(i,j,bi,bj,l)
                0187           Eswsf = RT_Es_sfc(i,j,bi,bj,l)
                0188
0222db53b0 Oliv*0189 C endif useOASIM
                0190          ENDIF
                0191
                0192          IF (myiter .ge. 0) THEN
                0193           Edwsf = Edwsf*(1.0 _d 0 - RT_iceFrac(i,j,bi,bj))
                0194           Eswsf = Eswsf*(1.0 _d 0 - RT_iceFrac(i,j,bi,bj))
                0195          ENDIF
b55e95f1ff Oliv*0196
                0197          DO k=1,Nr
                0198           a_k(k) = a(i,j,k,l)
                0199           bt_k(k) = bt(i,j,k,l)
                0200           bb_k(k) = bb(i,j,k,l)
                0201          ENDDO
                0202
                0203          CALL RADTRANS_SOLVE(
                0204      I                dz_k,rmud(i,j),Edwsf,Eswsf,a_k,bt_k,bb_k,klow,
                0205      O                Edbot,Esbot,Eubot,Estop,Eutop,
                0206      O                amp1_k,amp2_k, x_k, y_k,
                0207      O                r1_k,r2_k,kappa1_k,kappa2_k,
                0208      I                myThid)
                0209
                0210          Ed(i,j,1,l) = Edwsf
                0211          Es(i,j,1,l) = Eswsf
                0212          DO k=1,Nr-1
                0213           Ed(i,j,k+1,l) = Edbot(k)
                0214           Es(i,j,k+1,l) = Esbot(k)
                0215          ENDDO
                0216          DO k=1,Nr
                0217           Eu(i,j,k,l) = Eutop(k)
                0218          ENDDO
                0219
                0220          DO k=1,Nr
                0221 C convert to scalar irradiance in quanta
                0222           E0F(i,j,k,l) = rmud(i,j)*Ed(i,j,k,l) + RT_rmus*Es(i,j,k,l)
                0223      &                 + RT_rmuu*Eu(i,j,k,l)
                0224          ENDDO
                0225          E0F(i,j,Nr+1,l) = rmud(i,j)*Edbot(Nr) + RT_rmus*Esbot(Nr)
                0226      &                   + RT_rmuu*Eubot(Nr)
                0227
                0228 #ifdef ALLOW_DIAGNOSTICS
                0229          Edown = Edwsf + Eswsf
                0230          IF (Edown.GT.0 _d 0)THEN
                0231           Rirr(i,j,l) = Eutop(1)/Edown
                0232          ELSE
                0233           Rirr(i,j,l) = 0 _d 0
                0234          ENDIF
                0235 #ifdef RADTRANS_DIAG_SOLUTION
                0236          Eub(i,j,1,l) = 0.0 _d 0
                0237 #endif
                0238          DO k=1,Nr-1
                0239 #ifdef RADTRANS_DIAG_SOLUTION
                0240           Eub(i,j,k+1,l) = Eubot(k)
                0241 #endif
                0242          ENDDO
                0243          DO k=1,Nr
                0244 #ifdef RADTRANS_DIAG_SOLUTION
                0245           Est(i,j,k,l) = Estop(k)
                0246           amp1(i,j,k,l) = amp1_k(k)
                0247           amp2(i,j,k,l) = amp2_k(k)
                0248           x3d(i,j,k,l) = x_k(k)
                0249           y3d(i,j,k,l) = y_k(k)
                0250           r1(i,j,k,l) = r1_k(k)
                0251           r2(i,j,k,l) = r2_k(k)
                0252           kap1(i,j,k,l) = kappa1_k(k)
                0253           kap2(i,j,k,l) = kappa2_k(k)
                0254 #endif
                0255          ENDDO
                0256 #endif
                0257
                0258 C       l
                0259         ENDDO
                0260
                0261 C      i,j
                0262        ENDDO
                0263       ENDDO
                0264 C ======================================================================
                0265
                0266 #ifdef ALLOW_DIAGNOSTICS
                0267       IF (useDIAGNOSTICS .AND. myIter .GE.0) THEN
                0268        CALL DIAGNOSTICS_FILL(rmud,'rmud    ',1,1,2,bi,bj,myThid)
                0269        DO l=1,nlam
                0270         WRITE(diagname, '(A,I3.3)') 'Rirr', l
                0271         CALL DIAGNOSTICS_FILL(
                0272      &             Rirr(1-OLx,1-OLy,l),diagname,0,1,2,bi,bj,myThid)
                0273         WRITE(diagname, '(A,I3.3)') 'Ed', l
                0274         CALL DIAGNOSTICS_FILL(
                0275      &             Ed(1-OLx,1-OLy,1,l),diagname,0,Nr,2,bi,bj,myThid)
                0276         WRITE(diagname, '(A,I3.3)') 'Es', l
                0277         CALL DIAGNOSTICS_FILL(
                0278      &             Es(1-OLx,1-OLy,1,l),diagname,0,Nr,2,bi,bj,myThid)
                0279         WRITE(diagname, '(A,I3.3)') 'Eu', l
                0280         CALL DIAGNOSTICS_FILL(
                0281      &             Eu(1-OLx,1-OLy,1,l),diagname,0,Nr,2,bi,bj,myThid)
                0282         WRITE(diagname, '(A,I3.3)') 'E0F', l
                0283         CALL DIAGNOSTICS_FILL(
                0284      &             E0F(1-OLx,1-OLy,1,l),diagname,0,Nr,2,bi,bj,myThid)
                0285 #ifdef RADTRANS_DIAG_SOLUTION
                0286         WRITE(diagname, '(A,I3.3)') 'Estop', l
                0287         CALL DIAGNOSTICS_FILL(
                0288      &             Est(1-OLx,1-OLy,1,l),diagname,0,Nr,2,bi,bj,myThid)
                0289         WRITE(diagname, '(A,I3.3)') 'Eubot', l
                0290         CALL DIAGNOSTICS_FILL(
                0291      &             Eub(1-OLx,1-OLy,1,l),diagname,0,Nr,2,bi,bj,myThid)
                0292         WRITE(diagname, '(A,I3.3)') 'amp1_', l
                0293         CALL DIAGNOSTICS_FILL(
                0294      &             amp1(1-OLx,1-OLy,1,l),diagname,0,Nr,2,bi,bj,myThid)
                0295         WRITE(diagname, '(A,I3.3)') 'amp2_', l
                0296         CALL DIAGNOSTICS_FILL(
                0297      &             amp2(1-OLx,1-OLy,1,l),diagname,0,Nr,2,bi,bj,myThid)
                0298         WRITE(diagname, '(A,I3.3)') 'x_', l
                0299         CALL DIAGNOSTICS_FILL(
                0300      &             x3d(1-OLx,1-OLy,1,l),diagname,0,Nr,2,bi,bj,myThid)
                0301         WRITE(diagname, '(A,I3.3)') 'y_', l
                0302         CALL DIAGNOSTICS_FILL(
                0303      &             y3d(1-OLx,1-OLy,1,l),diagname,0,Nr,2,bi,bj,myThid)
                0304         WRITE(diagname, '(A,I3.3)') 'r1_', l
                0305         CALL DIAGNOSTICS_FILL(
                0306      &             r1(1-OLx,1-OLy,1,l),diagname,0,Nr,2,bi,bj,myThid)
                0307         WRITE(diagname, '(A,I3.3)') 'r2_', l
                0308         CALL DIAGNOSTICS_FILL(
                0309      &             r2(1-OLx,1-OLy,1,l),diagname,0,Nr,2,bi,bj,myThid)
                0310         WRITE(diagname, '(A,I3.3)') 'att1_', l
                0311         CALL DIAGNOSTICS_FILL(
                0312      &             kap1(1-OLx,1-OLy,1,l),diagname,0,Nr,2,bi,bj,myThid)
                0313         WRITE(diagname, '(A,I3.3)') 'att2_', l
                0314         CALL DIAGNOSTICS_FILL(
                0315      &             kap2(1-OLx,1-OLy,1,l),diagname,0,Nr,2,bi,bj,myThid)
                0316 #endif
                0317        ENDDO
                0318       ENDIF
                0319 #endif
                0320
                0321 #endif /* ALLOW_RADTRANS */
                0322
                0323       RETURN
                0324       END