source: LMDZ6/branches/contrails/libf/phylmd/rrtm/srtm_cldprop.F90 @ 5452

Last change on this file since 5452 was 1990, checked in by Laurent Fairhead, 11 years ago

Corrections à la version r1989 pour permettre la compilation avec RRTM
Inclusion de la licence CeCILL_V2 pour RRTM

Changes to revision r1989 to enable RRTM code compilation
RRTM part put under CeCILL_V2 licence
  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
File size: 12.2 KB
Line 
1SUBROUTINE SRTM_CLDPROP &
2 & ( KLEV, K_ICLDATM, K_INFLAG, K_ICEFLAG, K_LIQFLAG, K_NSTR,&
3 & P_CLDFRAC, P_CLDDAT1, P_CLDDAT2, P_CLDDAT3, P_CLDDAT4, P_CLDDATMOM,&
4 & P_TAUCLDORIG, P_TAUCLOUD, P_SSACLOUD, P_XMOM &
5 & ) 
6
7!  path:      $Source: /storm/rc1/cvsroot/rc/rrtm_sw/src/cldprop_sw.f,v $
8!  author:    $Author: jdelamer $
9!  revision:  $Revision: 2.6 $
10!  created:   $Date: 2002/04/04 18:29:47 $
11
12!  PURPOSE:  COMPUTE THE CLOUD OPTICAL DEPTH(S) FOR EACH CLOUDY
13!  LAYER.  NOTE:  ONLY INFLAG = 0 AND INFLAG=2/LIQFLAG=1,
14!  ICEFLAG=3
15!  (HU & STAMNES, Q. FU) ARE IMPLEMENTED.
16
17USE PARKIND1  ,ONLY : JPIM     ,JPRB
18USE YOMHOOK   ,ONLY : LHOOK,   DR_HOOK
19
20USE PARSRTM ,  ONLY : JPLAY, JPBAND, JPB1, JPB2
21USE YOESRTOP,  ONLY : EXTLIQ1, SSALIQ1, ASYLIQ1 &
22 & , EXTICE3, SSAICE3, ASYICE3, FDLICE3 &
23 & , FDELTA , EXTCOICE, SSACOICE, GICE, FORWICE &
24 & , EXTCOLIQ, SSACOLIQ, GLIQ, FORWLIQ 
25
26!     ------------------------------------------------------------------
27
28IMPLICIT NONE
29
30!-- real arguments
31
32INTEGER(KIND=JPIM),INTENT(IN)    :: KLEV
33INTEGER(KIND=JPIM),INTENT(OUT)   :: K_ICLDATM
34INTEGER(KIND=JPIM),INTENT(IN)    :: K_INFLAG
35INTEGER(KIND=JPIM),INTENT(IN)    :: K_ICEFLAG
36INTEGER(KIND=JPIM),INTENT(IN)    :: K_LIQFLAG
37INTEGER(KIND=JPIM),INTENT(IN)    :: K_NSTR
38REAL(KIND=JPRB)   ,INTENT(IN)    :: P_CLDFRAC(JPLAY)
39REAL(KIND=JPRB)   ,INTENT(IN)    :: P_CLDDAT1(JPLAY)
40REAL(KIND=JPRB)   ,INTENT(IN)    :: P_CLDDAT2(JPLAY)
41REAL(KIND=JPRB)   ,INTENT(IN)    :: P_CLDDAT3(JPLAY)
42REAL(KIND=JPRB)   ,INTENT(IN)    :: P_CLDDAT4(JPLAY)
43REAL(KIND=JPRB)   ,INTENT(IN)    :: P_CLDDATMOM(0:16,JPLAY)
44REAL(KIND=JPRB)   ,INTENT(INOUT) :: P_TAUCLDORIG(JPLAY,JPBAND)
45REAL(KIND=JPRB)   ,INTENT(INOUT) :: P_TAUCLOUD(JPLAY,JPBAND)
46REAL(KIND=JPRB)   ,INTENT(OUT)   :: P_SSACLOUD(JPLAY,JPBAND)
47REAL(KIND=JPRB)   ,INTENT(OUT)   :: P_XMOM(0:16,JPLAY,JPBAND)
48!-- integer arguments
49
50!-- real locals
51REAL(KIND=JPRB) :: Z_EPS
52REAL(KIND=JPRB) :: Z_TAUCLDORIG_A, Z_FFP, Z_FFP1, Z_FFPSSA, Z_SSACLOUD_A, Z_TAUCLOUD_A
53REAL(KIND=JPRB) :: Z_CWP, Z_FICE, Z_RADICE, Z_FACTOR, Z_FINT, Z_FLIQ, Z_RADLIQ
54REAL(KIND=JPRB) :: Z_TAUICEORIG, Z_SCATICE, Z_SSAICE, Z_TAUICE &
55 & , Z_TAULIQORIG, Z_SCATLIQ, Z_SSALIQ, Z_TAULIQ 
56
57!-- integer locals
58INTEGER(KIND=JPIM) :: I_NCBANDS, I_NLAYERS
59INTEGER(KIND=JPIM) :: IB, IB1, IB2, I_LAY, ISTR , INDEX
60INTEGER(KIND=JPIM) :: I_NDBUG
61REAL(KIND=JPRB) :: ZHOOK_HANDLE
62 
63!     INCLUDE 'param.f'
64
65!      COMMON /CONTROL/   IAER, NSTR, IOUT, ISTART, IEND, ICLD
66
67!      COMMON /PROFILE/   NLAYERS,PAVEL(MXLAY),TAVEL(MXLAY),
68!     &                   PZ(0:MXLAY),TZ(0:MXLAY)
69
70!      COMMON /CLOUDIN/   INFLAG,CLDDAT1(MXLAY),CLDDAT2(MXLAY),
71!     &                   ICEFLAG,LIQFLAG,CLDDAT3(MXLAY),CLDDAT4(MXLAY),
72!     &                   CLDDATMOM(0:16,MXLAY)
73
74!      COMMON /CLOUDDAT/  NCBANDS,CLDFRAC(MXLAY),
75!     &     TAUCLOUD(MXLAY,NBANDS),SSACLOUD(MXLAY,NBANDS),
76!     &     XMOM(0:16,MXLAY,NBANDS),TAUCLDORIG(MXLAY,NBANDS)
77
78!      COMMON /HVERSN/    HVRRTM,HVRRTR,HVRATM,HVRSET,HVRTAU,
79!     *                   HVDUM1(4),HVRUTL,HVREXT,
80!     *                   HVRD1M,HVRR1M,HVREPK,HVRLPK,HVRAER,HVRBKA,
81!     *                   HVRBKB,HVRCLD,HVRDIS,HVRLAM,HVRPAR
82
83!      CHARACTER*15 HVRRTM,HVRRTR,HVRATM,HVRSET,HVRTAU,
84!     *            HVDUM1,HVRUTL,HVREXT,
85!     *            HVRD1M,HVRR1M,HVREPK,HVRLPK,HVRAER,HVRBKA,
86!     *            HVRBKB,HVRCLD,HVRDIS,HVRLAM,HVRPAR
87
88!      DIMENSION ABSICE0(2), ABSICE1(2,5), ABSICE2(40,16)
89!      DIMENSION ABSLIQ1(58,16)
90
91!      DIMENSION ABSCOICE(NBANDS), ABSCOLIQ(NBANDS)
92!      DIMENSION EXTCOLIQ(NBANDS),SSACOLIQ(NBANDS),GLIQ(NBANDS)
93!      DIMENSION EXTCOICE(NBANDS),SSACOICE(NBANDS),GICE(NBANDS)
94!      DIMENSION FORWLIQ(NBANDS),FORWICE(NBANDS)
95
96! Added for SW
97!      DIMENSION EXTLIQ1(58,NBANDS), SSALIQ1(58,NBANDS),
98!     &     ASYLIQ1(58,NBANDS)
99!      DIMENSION EXTICE3(27,NBANDS), SSAICE3(27,NBANDS),
100!     &     ASYICE3(27,NBANDS),FDLICE3(28,NBANDS)
101!      DIMENSION FDELTA(NBANDS)
102
103!DATA EPS /1.E-6/
104
105! HVRCLD = '$Revision: 2.6 $'
106
107IF (LHOOK) CALL DR_HOOK('SRTM_CLDPROP',0,ZHOOK_HANDLE)
108Z_EPS = 1.E-06_JPRB
109I_NDBUG = 3
110
111K_ICLDATM = 0
112I_NCBANDS = 29
113I_NLAYERS = KLEV
114IB1     = JPB1
115IB2     = JPB2
116
117IF (I_NDBUG <= 2) THEN
118  print *,'cldprop before loop K_INFLAG, K_ICEFLAG, K_LIQFLAG:',K_INFLAG,K_ICEFLAG,K_LIQFLAG,IB1,IB2
119ENDIF
120
121DO I_LAY = 1, I_NLAYERS
122
123  IF (P_CLDFRAC(I_LAY) >= Z_EPS) THEN
124    K_ICLDATM = 1
125
126    IF (I_NDBUG <= 2) THEN
127      print 9101,I_LAY,K_ICLDATM,P_CLDFRAC(I_LAY),P_CLDDAT1(I_LAY),P_CLDDAT2(I_LAY),P_CLDDAT3(I_LAY)&
128       & ,P_CLDDAT4(I_LAY),(P_CLDDATMOM(ISTR,I_LAY),ISTR=0,K_NSTR) 
129      9101  format(1x,'Cld :',2I3,f7.4,7E12.5)
130    ENDIF
131
132!   Ice clouds and water clouds combined.
133    IF (K_INFLAG == 0) THEN
134      Z_TAUCLDORIG_A = P_CLDDAT1(I_LAY)
135      Z_FFP = P_CLDDATMOM(K_NSTR,I_LAY)
136      Z_FFP1 = 1.0 - Z_FFP
137      Z_FFPSSA = 1.0 - Z_FFP * P_CLDDAT2(I_LAY)
138      Z_SSACLOUD_A = Z_FFP1*P_CLDDAT2(I_LAY)/Z_FFPSSA
139      Z_TAUCLOUD_A = Z_FFPSSA*Z_TAUCLDORIG_A
140
141!       DO IB = 16,NBANDS
142      DO IB = IB1 , IB2
143        P_TAUCLDORIG(I_LAY,IB) = Z_TAUCLDORIG_A
144        P_SSACLOUD(I_LAY,IB) = Z_SSACLOUD_A
145        P_TAUCLOUD(I_LAY,IB) = Z_TAUCLOUD_A
146
147        DO ISTR = 0,K_NSTR
148          P_XMOM(ISTR,I_LAY,IB) = (P_CLDDATMOM(ISTR,I_LAY) - Z_FFP)/ &
149           & (Z_FFP1) 
150        ENDDO
151      ENDDO
152
153!   Separate treatement of ice clouds and water clouds.
154    ELSEIF(K_INFLAG == 2) THEN       
155      Z_CWP = P_CLDDAT1(I_LAY)
156      Z_FICE = P_CLDDAT2(I_LAY)
157      Z_RADICE = P_CLDDAT3(I_LAY)
158
159      IF (I_NDBUG <= 1) THEN
160        print 9102,I_LAY,Z_CWP,Z_FICE,Z_RADICE
161        9102     format(1x,'A',I3,3E13.6)
162      ENDIF
163
164!   Calculation of absorption coefficients due to ice clouds.
165      IF (Z_FICE == 0.0) THEN
166!         NCBANDS = 29
167!         DO IB = 16, NCBANDS
168        DO IB = IB1 , IB2
169          EXTCOICE(IB) = 0.0_JPRB
170          SSACOICE(IB) = 1.0_JPRB
171          GICE(IB)     = 1.0_JPRB
172          FORWICE(IB)  = 0.0_JPRB
173
174          IF (I_NDBUG <= 1) THEN
175            print 9103,I_LAY,Z_FICE,Z_CWP,Z_RADICE,IB,EXTCOICE(IB),SSACOICE(IB),GICE(IB),FORWICE(IB)
176            9103         format(1x,'B',I3,F6.3,2E13.6,I3,4E12.5)
177          ENDIF
178
179        ENDDO
180
181      ELSEIF (K_ICEFLAG == 3) THEN
182        IF (Z_RADICE < 10.0 .OR. Z_RADICE > 140.0) STOP 'ICE EFFECTIVE SIZE OUT OF BOUNDS'
183!         NCBANDS = 29
184        Z_FACTOR = (Z_RADICE - 5._JPRB)/5._JPRB
185        INDEX = INT(Z_FACTOR)
186        IF (INDEX == 27) INDEX = 26
187        Z_FINT = Z_FACTOR - REAL(INDEX)
188
189!         DO IB=16,NCBANDS
190        DO IB = IB1 , IB2
191          EXTCOICE(IB) = Z_FICE * (EXTICE3(INDEX,IB) + Z_FINT * &
192           & (EXTICE3(INDEX+1,IB) - EXTICE3(INDEX,IB))) 
193          SSACOICE(IB) = SSAICE3(INDEX,IB) + Z_FINT * &
194           & (SSAICE3(INDEX+1,IB) - SSAICE3(INDEX,IB)) 
195          GICE(IB) = ASYICE3(INDEX,IB) + Z_FINT * &
196           & (ASYICE3(INDEX+1,IB) - ASYICE3(INDEX,IB)) 
197          FDELTA(IB) = FDLICE3(INDEX,IB) + Z_FINT * &
198           & (FDLICE3(INDEX+1,IB) - FDLICE3(INDEX,IB)) 
199          if (fdelta(ib) < 0.0) STOP 'FDELTA LESS THAN 0.0'
200          if (fdelta(ib) > 1.0) STOP 'FDELTA GT THAN 1.0'                     
201          FORWICE(IB) = FDELTA(IB) + 0.5 / SSACOICE(IB)
202! See Fu 1996 p. 2067
203          IF (FORWICE(IB) > GICE(IB)) FORWICE(IB) = GICE(IB)
204! Check to ensure all calculated quantities are within physical limits.
205          if (extcoice(ib) < 0.0_JPRB) STOP 'ICE EXTINCTION LESS THAN 0.0'
206          if (ssacoice(ib) > 1.0_JPRB) STOP 'ICE SSA GRTR THAN 1.0'
207          if (ssacoice(ib) < 0.0_JPRB) STOP 'ICE SSA LESS THAN 0.0'
208          if (gice(ib) > 1.0_JPRB) STOP 'ICE ASYM GRTR THAN 1.0'
209          if (gice(ib) < 0.0_JPRB) STOP 'ICE ASYM LESS THAN 0.0'
210
211          IF (I_NDBUG <= 1) THEN
212            print 9104,I_LAY,Z_FICE,Z_CWP,Z_RADICE,IB,EXTCOICE(IB),SSACOICE(IB),GICE(IB),FORWICE(IB),FDELTA(IB)
213            9104         format(1x,'C',I3,F5.3,2E13.6,I3,5E12.5)
214          ENDIF
215
216        ENDDO
217      ENDIF
218      print *,'end of ice computations for I_LAY=',I_LAY
219                 
220!  Calculation of absorption coefficients due to water clouds.
221      Z_FLIQ = 1. - Z_FICE
222      IF (Z_FLIQ == 0.0) THEN
223!         NCBANDS = 29
224!         DO IB = 16, NCBANDS
225        DO IB = IB1 , IB2
226          EXTCOLIQ(IB) = 0.0
227          SSACOLIQ(IB) = 1.0
228          GLIQ(IB) = 1.0
229          FORWLIQ(IB) = 0.0
230
231          IF (I_NDBUG <= 1) THEN
232            print 9105,I_LAY,Z_FLIQ,Z_CWP,IB,EXTCOLIQ(IB),SSACOLIQ(IB),GLIQ(IB),FORWLIQ(IB)
233            9105         format(1x,'D',I3,F5.3,1E13.6,I3,4E12.5)
234          ENDIF
235
236        ENDDO
237
238      ELSEIF (K_LIQFLAG == 1) THEN
239        Z_RADLIQ = P_CLDDAT4(I_LAY)
240        IF (Z_RADLIQ < 1.5 .OR. Z_RADLIQ > 60.) STOP 'LIQUID EFFECTIVE RADIUS OUT OF BOUNDS'
241        INDEX = INT(Z_RADLIQ - 1.5)
242        IF (INDEX == 0) INDEX = 1
243        IF (INDEX == 58) INDEX = 57
244        Z_FINT = Z_RADLIQ - 1.5 - REAL(INDEX)
245!         NCBANDS = 29
246
247!         DO IB = 16, NCBANDS
248        DO IB = IB1 , IB2
249          EXTCOLIQ(IB) = Z_FLIQ * (EXTLIQ1(INDEX,IB) + Z_FINT * &
250           & (EXTLIQ1(INDEX+1,IB) - (EXTLIQ1(INDEX,IB)))) 
251          SSACOLIQ(IB) = SSALIQ1(INDEX,IB) + Z_FINT * &
252           & (SSALIQ1(INDEX+1,IB) - SSALIQ1(INDEX,IB)) 
253          GLIQ(IB) = ASYLIQ1(INDEX,IB) + Z_FINT * &
254           & (ASYLIQ1(INDEX+1,IB) - ASYLIQ1(INDEX,IB)) 
255          FORWLIQ(IB) = GLIQ(IB)**K_NSTR
256! Check to ensure all calculated quantities are within physical limits.
257          if (extcoliq(ib) < 0.0_JPRB) STOP 'LIQUID EXTINCTION LESS THAN 0.0'
258          if (ssacoliq(ib) > 1.0_JPRB) STOP 'LIQUID SSA GRTR THAN 1.0'
259          if (ssacoliq(ib) < 0.0_JPRB) STOP 'LIQUID SSA LESS THAN 0.0'
260          if (gliq(ib) > 1.0_JPRB) STOP 'LIQUID ASYM GRTR THAN 1.0'
261          if (gliq(ib) < 0.0_JPRB) STOP 'LIQUID ASYM LESS THAN 0.0'
262
263          IF (I_NDBUG <= 1) THEN
264            print 9106,I_LAY,Z_FLIQ,Z_CWP,Z_RADLIQ,IB,EXTCOLIQ(IB),SSACOLIQ(IB),GLIQ(IB),FORWLIQ(IB)
265            9106         format(1x,'E',I3,F5.3,2E13.6,I3,5E12.5)
266          ENDIF
267
268        ENDDO
269      ENDIF
270
271      IF (I_NDBUG <= 1) THEN
272        print *,'end of liquid water computations for I_LAY=',I_LAY
273      ENDIF
274
275!       DO IB = 16, NCBANDS
276      DO IB = IB1 , IB2
277        Z_TAULIQORIG = Z_CWP * EXTCOLIQ(IB)
278        Z_TAUICEORIG = Z_CWP * EXTCOICE(IB)
279        P_TAUCLDORIG(I_LAY,IB) = Z_TAULIQORIG + Z_TAUICEORIG
280
281        IF (I_NDBUG <= 1) THEN
282          print 9107,IB,Z_TAULIQORIG,Z_TAUICEORIG,P_TAUCLDORIG(I_LAY,IB),Z_CWP &
283           & ,EXTCOLIQ(IB),EXTCOICE(IB),SSACOLIQ(IB),SSACOICE(IB) &
284           & ,FORWLIQ(IB),FORWICE(IB) 
285          9107       format(1x,'F',I3,10E12.5)
286        ENDIF
287
288        Z_SSALIQ = SSACOLIQ(IB) * (1. - FORWLIQ(IB)) / &
289         & (1. - FORWLIQ(IB) * SSACOLIQ(IB)) 
290        Z_TAULIQ =  (1. - FORWLIQ(IB) * SSACOLIQ(IB)) * &
291         & Z_TAULIQORIG 
292        Z_SSAICE = SSACOICE(IB) * (1. - FORWICE(IB)) / &
293         & (1. - FORWICE(IB) * SSACOICE(IB)) 
294        Z_TAUICE =  (1. - FORWICE(IB) * SSACOICE(IB)) * &
295         & Z_TAUICEORIG 
296        Z_SCATLIQ = Z_SSALIQ * Z_TAULIQ
297        Z_SCATICE = Z_SSAICE * Z_TAUICE
298        P_TAUCLOUD(I_LAY,IB) = Z_TAULIQ + Z_TAUICE
299        P_SSACLOUD(I_LAY,IB) = (Z_SCATLIQ + Z_SCATICE) / &
300         & P_TAUCLOUD(I_LAY,IB) 
301        P_XMOM(0,I_LAY,IB) = 1.0
302
303        IF (I_NDBUG <= 1) THEN
304          print 9108,IB,Z_TAULIQORIG,Z_TAUICEORIG,Z_SSALIQ,Z_TAULIQ,Z_SCATLIQ,Z_SSAICE,Z_TAUICE,Z_SCATICE
305          9108       format(1x,'G',I3,8E13.6)
306        ENDIF
307
308        DO ISTR = 1, K_NSTR
309!This commented code is the standard method for delta-m scaling. In accordance
310!  with the 1996 Fu paper, equation A.3, the moments for ice were calculated
311!  as in the uncommented code.
312!                     XMOM(ISTR,LAY,IB) = (SCATLIQ *  &
313!     &                    (GLIQ(IB)**ISTR - FORWLIQ(IB)) / &
314!     &                    (1. - FORWLIQ(IB)) &
315!     &                    + SCATICE * &
316!     &                    (GICE(IB)**ISTR - FORWICE(IB)) /  &
317!     &                    (1. - FORWICE(IB)))/(SCATLIQ + SCATICE)
318
319          P_XMOM(ISTR,I_LAY,IB) = (1.0/(Z_SCATLIQ+Z_SCATICE))* &
320           & (Z_SCATLIQ*(GLIQ(IB)**ISTR - FORWLIQ(IB)) / &
321           & (1. - FORWLIQ(IB)) &
322           & + Z_SCATICE * &
323           & ((gice(ib)-forwice(ib))/(1.0-forwice(ib)))**ISTR) 
324        ENDDO
325      ENDDO
326
327    ENDIF
328
329  ENDIF
330
331ENDDO
332
333IF (I_NDBUG <= 1) THEN
334  print *,'about to leave SRTM_CLDPROP'
335ENDIF
336
337!-----------------------------------------------------------------------
338IF (LHOOK) CALL DR_HOOK('SRTM_CLDPROP',1,ZHOOK_HANDLE)
339END SUBROUTINE SRTM_CLDPROP
340
Note: See TracBrowser for help on using the repository browser.