source: LMDZ6/branches/contrails/libf/phylmd/ecrad/ifsrrtm/srtm_taumol28.F90 @ 5431

Last change on this file since 5431 was 4773, checked in by idelkadi, 12 months ago
  • Update of Ecrad in LMDZ The same organization of the Ecrad offline version is retained in order to facilitate the updating of Ecrad in LMDZ and the comparison between online and offline results. version 1.6.1 of Ecrad (​https://github.com/lguez/ecrad.git)
  • Implementation of the double call of Ecrad in LMDZ


File size: 8.0 KB
Line 
1SUBROUTINE SRTM_TAUMOL28 &
2 & ( KIDIA   , KFDIA    , KLEV,&
3 & P_FAC00   , P_FAC01  , P_FAC10   , P_FAC11,&
4 & K_JP      , K_JT     , K_JT1     , P_ONEMINUS,&
5 & P_COLMOL  , P_COLO2  , P_COLO3,&
6 & K_LAYTROP,&
7 & P_SFLUXZEN, P_TAUG   , P_TAUR    , PRMU0   &
8 & ) 
9
10!     Written by Eli J. Mlawer, Atmospheric & Environmental Research.
11
12!     BAND 28:  38000-50000 cm-1 (low - O3,O2; high - O3,O2)
13
14! Modifications
15!        M.Hamrud      01-Oct-2003 CY28 Cleaning
16
17!     JJMorcrette 2003-02-24 adapted to ECMWF environment
18!        D.Salmond  31-Oct-2007 Vector version in the style of RRTM from Meteo France & NEC
19!     JJMorcrette 20010610 Flexible configuration for number of g-points
20
21USE PARKIND1 , ONLY : JPIM, JPRB
22USE YOMHOOK  , ONLY : LHOOK, DR_HOOK, JPHOOK
23USE PARSRTM  , ONLY : JPG
24USE YOESRTM  , ONLY : NG28
25USE YOESRTA28, ONLY : ABSA, ABSB, SFLUXREFC, RAYL, LAYREFFR, STRRAT 
26USE YOESRTWN , ONLY : NSPA, NSPB
27
28IMPLICIT NONE
29
30!-- Output
31INTEGER(KIND=JPIM),INTENT(IN)    :: KIDIA, KFDIA
32INTEGER(KIND=JPIM),INTENT(IN)    :: KLEV
33REAL(KIND=JPRB)   ,INTENT(IN)    :: P_FAC00(KIDIA:KFDIA,KLEV)
34REAL(KIND=JPRB)   ,INTENT(IN)    :: P_FAC01(KIDIA:KFDIA,KLEV)
35REAL(KIND=JPRB)   ,INTENT(IN)    :: P_FAC10(KIDIA:KFDIA,KLEV)
36REAL(KIND=JPRB)   ,INTENT(IN)    :: P_FAC11(KIDIA:KFDIA,KLEV)
37INTEGER(KIND=JPIM),INTENT(IN)    :: K_JP(KIDIA:KFDIA,KLEV)
38INTEGER(KIND=JPIM),INTENT(IN)    :: K_JT(KIDIA:KFDIA,KLEV)
39INTEGER(KIND=JPIM),INTENT(IN)    :: K_JT1(KIDIA:KFDIA,KLEV)
40REAL(KIND=JPRB)   ,INTENT(IN)    :: P_ONEMINUS(KIDIA:KFDIA)
41REAL(KIND=JPRB)   ,INTENT(IN)    :: P_COLMOL(KIDIA:KFDIA,KLEV)
42REAL(KIND=JPRB)   ,INTENT(IN)    :: P_COLO2(KIDIA:KFDIA,KLEV)
43REAL(KIND=JPRB)   ,INTENT(IN)    :: P_COLO3(KIDIA:KFDIA,KLEV)
44INTEGER(KIND=JPIM),INTENT(IN)    :: K_LAYTROP(KIDIA:KFDIA)
45
46REAL(KIND=JPRB)   ,INTENT(OUT)   :: P_SFLUXZEN(KIDIA:KFDIA,JPG)
47REAL(KIND=JPRB)   ,INTENT(OUT)   :: P_TAUG(KIDIA:KFDIA,KLEV,JPG)
48REAL(KIND=JPRB)   ,INTENT(OUT)   :: P_TAUR(KIDIA:KFDIA,KLEV,JPG)
49REAL(KIND=JPRB)   ,INTENT(IN)    :: PRMU0(KIDIA:KFDIA)
50!- from INTFAC     
51!- from INTIND
52!- from PRECISE             
53!- from PROFDATA             
54!- from SELF             
55INTEGER(KIND=JPIM) :: IG, IND0, IND1, JS, I_LAY, I_LAYSOLFR(KIDIA:KFDIA), I_NLAYERS, IPLON
56
57REAL(KIND=JPRB) :: Z_FAC000, Z_FAC001, Z_FAC010, Z_FAC011, Z_FAC100, Z_FAC101,&
58 & Z_FAC110, Z_FAC111, Z_FS, Z_SPECCOMB, Z_SPECMULT, Z_SPECPARM, &
59 & Z_TAURAY 
60REAL(KIND=JPHOOK) :: ZHOOK_HANDLE
61
62IF (LHOOK) CALL DR_HOOK('SRTM_TAUMOL28',0,ZHOOK_HANDLE)
63
64I_NLAYERS = KLEV
65
66!     Compute the optical depth by interpolating in ln(pressure),
67!     temperature, and appropriate species.  Below LAYTROP, the water
68!     vapor self-continuum is interpolated (in temperature) separately. 
69
70DO I_LAY = 1, I_NLAYERS
71  DO IPLON = KIDIA, KFDIA
72    IF (PRMU0(IPLON) > 0.0_JPRB) THEN
73      IF (I_LAY <= K_LAYTROP(IPLON)) THEN
74        Z_SPECCOMB = P_COLO3(IPLON,I_LAY) + STRRAT*P_COLO2(IPLON,I_LAY)
75        Z_SPECPARM = P_COLO3(IPLON,I_LAY)/Z_SPECCOMB
76        IF (Z_SPECPARM >= P_ONEMINUS(IPLON)) Z_SPECPARM = P_ONEMINUS(IPLON)
77        Z_SPECMULT = 8.*(Z_SPECPARM)
78        JS = 1 + INT(Z_SPECMULT)
79        Z_FS = MOD(Z_SPECMULT, 1.0_JPRB )
80        ! Z_FAC000 = (1. - Z_FS) * P_FAC00(I_LAY)
81        ! Z_FAC010 = (1. - Z_FS) * P_FAC10(I_LAY)
82        ! Z_FAC100 = Z_FS * P_FAC00(I_LAY)
83        ! Z_FAC110 = Z_FS * P_FAC10(I_LAY)
84        ! Z_FAC001 = (1. - Z_FS) * P_FAC01(I_LAY)
85        ! Z_FAC011 = (1. - Z_FS) * P_FAC11(I_LAY)
86        ! Z_FAC101 = Z_FS * P_FAC01(I_LAY)
87        ! Z_FAC111 = Z_FS * P_FAC11(I_LAY)
88        IND0 = ((K_JP(IPLON,I_LAY)-1)*5+(K_JT(IPLON,I_LAY)-1))*NSPA(28) + JS
89        IND1 = (K_JP(IPLON,I_LAY)*5+(K_JT1(IPLON,I_LAY)-1))*NSPA(28) + JS
90        Z_TAURAY = P_COLMOL(IPLON,I_LAY) * RAYL
91
92        !  DO IG = 1, NG(28)
93!CDIR UNROLL=NG28
94        DO IG = 1 , NG28
95          P_TAUG(IPLON,I_LAY,IG) = Z_SPECCOMB * &
96           !    & (Z_FAC000 * ABSA(IND0,IG) + &
97           !    & Z_FAC100 * ABSA(IND0+1,IG) + &
98           !    & Z_FAC010 * ABSA(IND0+9,IG) + &
99           !    & Z_FAC110 * ABSA(IND0+10,IG) + &
100           !    & Z_FAC001 * ABSA(IND1,IG) + &
101           !    & Z_FAC101 * ABSA(IND1+1,IG) + &
102           !    & Z_FAC011 * ABSA(IND1+9,IG) + &
103           !    & Z_FAC111 * ABSA(IND1+10,IG))   
104           & (&
105           & (1. - Z_FS) * ( ABSA(IND0,IG) * P_FAC00(IPLON,I_LAY) + &
106           &                 ABSA(IND0+9,IG) * P_FAC10(IPLON,I_LAY) + &
107           &                 ABSA(IND1,IG) * P_FAC01(IPLON,I_LAY) + &
108           &                 ABSA(IND1+9,IG) * P_FAC11(IPLON,I_LAY) ) + &
109           & Z_FS        * ( ABSA(IND0+1,IG) * P_FAC00(IPLON,I_LAY) + &
110           &                 ABSA(IND0+10,IG) * P_FAC10(IPLON,I_LAY) + &
111           &                 ABSA(IND1+1,IG) * P_FAC01(IPLON,I_LAY) + &
112           &                 ABSA(IND1+10,IG) * P_FAC11(IPLON,I_LAY) ) &
113           & )
114          !     &           + TAURAY
115          !    SSA(LAY,IG) = TAURAY/TAUG(LAY,IG)
116          P_TAUR(IPLON,I_LAY,IG) = Z_TAURAY
117        ENDDO
118      ENDIF
119    ENDIF
120  ENDDO
121ENDDO
122
123I_LAYSOLFR(:) = I_NLAYERS
124
125DO I_LAY = 1, I_NLAYERS
126  DO IPLON = KIDIA, KFDIA
127    IF (PRMU0(IPLON) > 0.0_JPRB) THEN
128      IF (I_LAY >= K_LAYTROP(IPLON)+1) THEN
129        IF (K_JP(IPLON,I_LAY-1) < LAYREFFR .AND. K_JP(IPLON,I_LAY) >= LAYREFFR) &
130         & I_LAYSOLFR(IPLON) = I_LAY 
131        Z_SPECCOMB = P_COLO3(IPLON,I_LAY) + STRRAT*P_COLO2(IPLON,I_LAY)
132        Z_SPECPARM = P_COLO3(IPLON,I_LAY)/Z_SPECCOMB
133        IF (Z_SPECPARM >= P_ONEMINUS(IPLON)) Z_SPECPARM = P_ONEMINUS(IPLON)
134        Z_SPECMULT = 4.*(Z_SPECPARM)
135        JS = 1 + INT(Z_SPECMULT)
136        Z_FS = MOD(Z_SPECMULT, 1.0_JPRB )
137        ! Z_FAC000 = (1. - Z_FS) * P_FAC00(I_LAY)
138        ! Z_FAC010 = (1. - Z_FS) * P_FAC10(I_LAY)
139        ! Z_FAC100 = Z_FS * P_FAC00(I_LAY)
140        ! Z_FAC110 = Z_FS * P_FAC10(I_LAY)
141        ! Z_FAC001 = (1. - Z_FS) * P_FAC01(I_LAY)
142        ! Z_FAC011 = (1. - Z_FS) * P_FAC11(I_LAY)
143        ! Z_FAC101 = Z_FS * P_FAC01(I_LAY)
144        ! Z_FAC111 = Z_FS * P_FAC11(I_LAY)
145        IND0 = ((K_JP(IPLON,I_LAY)-13)*5+(K_JT(IPLON,I_LAY)-1))*NSPB(28) + JS
146        IND1 = ((K_JP(IPLON,I_LAY)-12)*5+(K_JT1(IPLON,I_LAY)-1))*NSPB(28) + JS
147        Z_TAURAY = P_COLMOL(IPLON,I_LAY) * RAYL
148
149        !  DO IG = 1, NG(28)
150!CDIR UNROLL=NG28
151        DO IG = 1 , NG28
152          P_TAUG(IPLON,I_LAY,IG) = Z_SPECCOMB * &
153           !    & (Z_FAC000 * ABSB(IND0,IG) + &
154           !    & Z_FAC100 * ABSB(IND0+1,IG) + &
155           !    & Z_FAC010 * ABSB(IND0+5,IG) + &
156           !    & Z_FAC110 * ABSB(IND0+6,IG) + &
157           !    & Z_FAC001 * ABSB(IND1,IG) + &
158           !    & Z_FAC101 * ABSB(IND1+1,IG) + &
159           !    & Z_FAC011 * ABSB(IND1+5,IG) + &
160           !    & Z_FAC111 * ABSB(IND1+6,IG))   
161           & (&
162           & (1. - Z_FS) * ( ABSB(IND0,IG) * P_FAC00(IPLON,I_LAY) + &
163           &                 ABSB(IND0+5,IG) * P_FAC10(IPLON,I_LAY) + &
164           &                 ABSB(IND1,IG) * P_FAC01(IPLON,I_LAY) + &
165           &                 ABSB(IND1+5,IG) * P_FAC11(IPLON,I_LAY) ) + &
166           & Z_FS        * ( ABSB(IND0+1,IG) * P_FAC00(IPLON,I_LAY) + &
167           &                 ABSB(IND0+6,IG) * P_FAC10(IPLON,I_LAY) + &
168           &                 ABSB(IND1+1,IG) * P_FAC01(IPLON,I_LAY) + &
169           &                 ABSB(IND1+6,IG) * P_FAC11(IPLON,I_LAY) ) &
170           & )
171          !     &           + TAURAY
172          !    SSA(LAY,IG) = TAURAY/TAUG(LAY,IG)
173          IF (I_LAY == I_LAYSOLFR(IPLON)) P_SFLUXZEN(IPLON,IG) = SFLUXREFC(IG,JS) &
174           & + Z_FS * (SFLUXREFC(IG,JS+1) - SFLUXREFC(IG,JS)) 
175! The following actually improves this band by setting the solar
176! spectrum at each g point equal to what would be computed if
177! molecular oxygen was set to zero. But it is worse overall due to a
178! compensating error with the previous band 27.
179!          IF (I_LAY == I_LAYSOLFR) P_SFLUXZEN(IPLON,IG) = SFLUXREFC(IG,5)
180          P_TAUR(IPLON,I_LAY,IG) = Z_TAURAY
181        ENDDO
182      ENDIF
183    ENDIF
184  ENDDO
185ENDDO
186
187!-----------------------------------------------------------------------
188IF (LHOOK) CALL DR_HOOK('SRTM_TAUMOL28',1,ZHOOK_HANDLE)
189
190END SUBROUTINE SRTM_TAUMOL28
Note: See TracBrowser for help on using the repository browser.