source: LMDZ6/branches/cirrus/libf/phylmd/ecrad/ifsrrtm/srtm_taumol16.F90

Last change on this file was 4773, checked in by idelkadi, 11 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: 7.7 KB
Line 
1SUBROUTINE SRTM_TAUMOL16 &
2 & ( KIDIA   , KFDIA     , KLEV,&
3 & P_FAC00   , P_FAC01   , P_FAC10   , P_FAC11,&
4 & K_JP      , K_JT      , K_JT1     , P_ONEMINUS,&
5 & P_COLH2O  , P_COLCH4  , P_COLMOL,&
6 & K_LAYTROP , P_SELFFAC , P_SELFFRAC, K_INDSELF , P_FORFAC  , P_FORFRAC, K_INDFOR,&
7 & P_SFLUXZEN, P_TAUG    , P_TAUR    , PRMU0     &
8 & ) 
9
10!     Written by Eli J. Mlawer, Atmospheric & Environmental Research.
11
12!     BAND 16:  2600-3250 cm-1 (low - H2O,CH4; high - CH4)
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 : NG16
25USE YOESRTA16, ONLY : ABSA, ABSB, FORREFC, SELFREFC, SFLUXREFC, RAYL, LAYREFFR, STRRAT1 
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_COLH2O(KIDIA:KFDIA,KLEV)
42REAL(KIND=JPRB)   ,INTENT(IN)    :: P_COLCH4(KIDIA:KFDIA,KLEV)
43REAL(KIND=JPRB)   ,INTENT(IN)    :: P_COLMOL(KIDIA:KFDIA,KLEV)
44INTEGER(KIND=JPIM),INTENT(IN)    :: K_LAYTROP(KIDIA:KFDIA)
45REAL(KIND=JPRB)   ,INTENT(IN)    :: P_SELFFAC(KIDIA:KFDIA,KLEV)
46REAL(KIND=JPRB)   ,INTENT(IN)    :: P_SELFFRAC(KIDIA:KFDIA,KLEV)
47INTEGER(KIND=JPIM),INTENT(IN)    :: K_INDSELF(KIDIA:KFDIA,KLEV)
48REAL(KIND=JPRB)   ,INTENT(IN)    :: P_FORFAC(KIDIA:KFDIA,KLEV)
49REAL(KIND=JPRB)   ,INTENT(IN)    :: P_FORFRAC(KIDIA:KFDIA,KLEV)
50INTEGER(KIND=JPIM),INTENT(IN)    :: K_INDFOR(KIDIA:KFDIA,KLEV)
51
52REAL(KIND=JPRB)   ,INTENT(OUT)   :: P_SFLUXZEN(KIDIA:KFDIA,JPG)
53REAL(KIND=JPRB)   ,INTENT(OUT)   :: P_TAUG(KIDIA:KFDIA,KLEV,JPG)
54REAL(KIND=JPRB)   ,INTENT(OUT)   :: P_TAUR(KIDIA:KFDIA,KLEV,JPG)
55REAL(KIND=JPRB)   ,INTENT(IN)    :: PRMU0(KIDIA:KFDIA)
56!- from INTFAC     
57!- from INTIND
58!- from PRECISE             
59!- from PROFDATA             
60!- from SELF             
61INTEGER(KIND=JPIM) :: IG, IND0, IND1, INDS, INDF, JS, I_LAY, I_LAYSOLFR(KIDIA:KFDIA), I_NLAYERS, IPLON
62
63REAL(KIND=JPRB) :: Z_FAC000, Z_FAC001, Z_FAC010, Z_FAC011, Z_FAC100, Z_FAC101,&
64 & Z_FAC110, Z_FAC111, Z_FS, Z_SPECCOMB, Z_SPECMULT, Z_SPECPARM, &
65 & Z_TAURAY 
66REAL(KIND=JPHOOK) :: ZHOOK_HANDLE
67
68IF (LHOOK) CALL DR_HOOK('SRTM_TAUMOL16',0,ZHOOK_HANDLE)
69
70I_NLAYERS = KLEV
71
72!     Compute the optical depth by interpolating in ln(pressure),
73!     temperature, and appropriate species.  Below LAYTROP, the water
74!     vapor self-continuum is interpolated (in temperature) separately. 
75
76DO I_LAY = 1, I_NLAYERS
77  DO IPLON = KIDIA, KFDIA
78    IF (PRMU0(IPLON) > 0.0_JPRB) THEN
79      IF (I_LAY <= K_LAYTROP(IPLON)) THEN
80        Z_SPECCOMB = P_COLH2O(IPLON,I_LAY) + STRRAT1*P_COLCH4(IPLON,I_LAY)
81        Z_SPECPARM = P_COLH2O(IPLON,I_LAY)/Z_SPECCOMB
82        IF (Z_SPECPARM >= P_ONEMINUS(IPLON)) Z_SPECPARM = P_ONEMINUS(IPLON)
83        Z_SPECMULT = 8._JPRB*(Z_SPECPARM)
84        JS = 1 + INT(Z_SPECMULT)
85        Z_FS = MOD(Z_SPECMULT, 1.0_JPRB )
86        ! Z_FAC000 = (1.0_JPRB - Z_FS) * P_FAC00(I_LAY)
87        ! Z_FAC010 = (1.0_JPRB - Z_FS) * P_FAC10(I_LAY)
88        ! Z_FAC100 = Z_FS * P_FAC00(I_LAY)
89        ! Z_FAC110 = Z_FS * P_FAC10(I_LAY)
90        ! Z_FAC001 = (1.0_JPRB - Z_FS) * P_FAC01(I_LAY)
91        ! Z_FAC011 = (1.0_JPRB - Z_FS) * P_FAC11(I_LAY)
92        ! Z_FAC101 = Z_FS * P_FAC01(I_LAY)
93        ! Z_FAC111 = Z_FS * P_FAC11(I_LAY)
94        IND0 = ((K_JP(IPLON,I_LAY)-1)*5+(K_JT(IPLON,I_LAY)-1))*NSPA(16) + JS
95        IND1 = (K_JP(IPLON,I_LAY)*5+(K_JT1(IPLON,I_LAY)-1))*NSPA(16) + JS
96        INDS = K_INDSELF(IPLON,I_LAY)
97        INDF = K_INDFOR(IPLON,I_LAY)
98        Z_TAURAY = P_COLMOL(IPLON,I_LAY) * RAYL
99
100        !  print 9001,LAY,IND0,IND1,INDS,INDF,FAC000,FAC010,FAC100,FAC110,FAC001,FAC011,FAC101,FAC111 &
101        !   &,TAURAY,SELFFAC(LAY),SELFFRAC(LAY),FORFAC(LAY),FORFRAC(LAY)
1029001    format(1x,'T16 ',5I4,13E12.3)
103
104        !  DO IG = 1, NG(16)
105!CDIR UNROLL=NG16
106        DO IG = 1, NG16
107          P_TAUG(IPLON,I_LAY,IG) = Z_SPECCOMB * &
108           !    & (Z_FAC000 * ABSA(IND0   ,IG) + &
109           !    & Z_FAC100 * ABSA(IND0 +1,IG) + &
110           !    & Z_FAC010 * ABSA(IND0 +9,IG) + &
111           !    & Z_FAC110 * ABSA(IND0+10,IG) + &
112           !    & Z_FAC001 * ABSA(IND1   ,IG) + &
113           !    & Z_FAC101 * ABSA(IND1 +1,IG) + &
114           !    & Z_FAC011 * ABSA(IND1 +9,IG) + &
115           !    & Z_FAC111 * ABSA(IND1+10,IG)) + &
116           & (&
117           & (1. - Z_FS) * ( ABSA(IND0,IG) * P_FAC00(IPLON,I_LAY) + &
118           &                 ABSA(IND0+9,IG) * P_FAC10(IPLON,I_LAY) + &
119           &                 ABSA(IND1,IG) * P_FAC01(IPLON,I_LAY) + &
120           &                 ABSA(IND1+9,IG) * P_FAC11(IPLON,I_LAY) ) + &
121           & Z_FS        * ( ABSA(IND0+1,IG) * P_FAC00(IPLON,I_LAY) + &
122           &                 ABSA(IND0+10,IG) * P_FAC10(IPLON,I_LAY) + &
123           &                 ABSA(IND1+1,IG) * P_FAC01(IPLON,I_LAY) + &
124           &                 ABSA(IND1+10,IG) * P_FAC11(IPLON,I_LAY) ) &
125           & ) + &
126           & P_COLH2O(IPLON,I_LAY) * &
127           & (P_SELFFAC(IPLON,I_LAY) * (SELFREFC(INDS,IG) + &
128           & P_SELFFRAC(IPLON,I_LAY) * &
129           & (SELFREFC(INDS+1,IG) - SELFREFC(INDS,IG))) + &
130           & P_FORFAC(IPLON,I_LAY) * (FORREFC(INDF,IG) + &
131           & P_FORFRAC(IPLON,I_LAY) * &
132           & (FORREFC(INDF+1,IG) - FORREFC(INDF,IG))))   
133          !     &           + TAURAY
134          !    SSA(LAY,IG) = TAURAY/TAUG(LAY,IG)
135          P_TAUR(IPLON,I_LAY,IG) = Z_TAURAY
136          !    print 9002,LAY,IG,ABSA(IND0,IG),ABSA(IND0+1,IG),ABSA(IND0+9,IG),ABSA(IND0+10,IG) &
137          !      &, ABSA(IND1,IG),ABSA(IND1+1,IG),ABSA(IND1+9,IG),ABSA(IND1+10,IG) &
138          !      &, SELFREFC(INDS+1,IG),SELFREFC(INDS,IG),FORREFC(INDF+1,IG),FORREFC(INDF,IG)
1399002      format(1x,'U16 ',2I3,12E12.3)
140        ENDDO
141      ENDIF
142    ENDIF
143  ENDDO
144ENDDO
145
146I_LAYSOLFR(:) = I_NLAYERS
147
148DO I_LAY = 1, I_NLAYERS
149  DO IPLON = KIDIA, KFDIA
150    IF (PRMU0(IPLON) > 0.0_JPRB) THEN
151      IF (I_LAY >= K_LAYTROP(IPLON)+1) THEN
152        IF (K_JP(IPLON,I_LAY-1) < LAYREFFR .AND. K_JP(IPLON,I_LAY) >= LAYREFFR) &
153         & I_LAYSOLFR(IPLON) = I_LAY 
154        IND0 = ((K_JP(IPLON,I_LAY)-13)*5+(K_JT(IPLON,I_LAY)-1))*NSPB(16) + 1
155        IND1 = ((K_JP(IPLON,I_LAY)-12)*5+(K_JT1(IPLON,I_LAY)-1))*NSPB(16) + 1
156        Z_TAURAY = P_COLMOL(IPLON,I_LAY) * RAYL
157
158        !  DO IG = 1, NG(16)
159!CDIR UNROLL=NG16
160        DO IG = 1, NG16
161          P_TAUG(IPLON,I_LAY,IG) = P_COLCH4(IPLON,I_LAY) * &
162           & (P_FAC00(IPLON,I_LAY) * ABSB(IND0  ,IG) + &
163           & P_FAC10(IPLON,I_LAY) * ABSB(IND0+1,IG) + &
164           & P_FAC01(IPLON,I_LAY) * ABSB(IND1  ,IG) + &
165           & P_FAC11(IPLON,I_LAY) * ABSB(IND1+1,IG))   
166          !     &           + TAURAY
167          !    SSA(LAY,IG) = TAURAY/TAUG(LAY,IG)
168          IF (I_LAY == I_LAYSOLFR(IPLON)) P_SFLUXZEN(IPLON,IG) = SFLUXREFC(IG)
169          P_TAUR(IPLON,I_LAY,IG) = Z_TAURAY 
170        ENDDO
171      ENDIF
172    ENDIF
173  ENDDO
174ENDDO
175
176!DO LAY=1,NLAYERS
177!  print 9003,LAY,(TAUG(LAY,IG),IG=1,NG16)
1789003 format(1x,'O16 ',I3,16E13.5)
179!END DO
180
181!----------------------------------------------------------------------
182IF (LHOOK) CALL DR_HOOK('SRTM_TAUMOL16',1,ZHOOK_HANDLE)
183
184END SUBROUTINE SRTM_TAUMOL16
Note: See TracBrowser for help on using the repository browser.