source: LMDZ5/branches/testing/libf/phylmd/rrtm/swde.F90 @ 2435

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

Merged trunk changes r1920:1997 into testing branch
  • 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: 7.0 KB
Line 
1!OPTIONS XOPT(HSFUN)
2SUBROUTINE SWDE &
3 & ( KIDIA, KFDIA, KLON,&
4 & PGG  , PREF , PRMUZ, PTO1, PW,&
5 & PRE1 , PRE2 , PTR1 , PTR2 &
6 & ) 
7
8!**** *SWDE* - DELTA-EDDINGTON IN A CLOUDY LAYER
9
10!     PURPOSE.
11!     --------
12!           COMPUTES THE REFLECTIVITY AND TRANSMISSIVITY OF A CLOUDY
13!     LAYER USING THE DELTA-EDDINGTON'S APPROXIMATION.
14
15!**   INTERFACE.
16!     ----------
17!          *SWDE* IS CALLED BY *SWR*, *SWNI*
18
19!        EXPLICIT ARGUMENTS :
20!        --------------------
21! PGG    : (KLON)             ; ASSYMETRY FACTOR
22! PREF   : (KLON)             ; REFLECTIVITY OF THE UNDERLYING LAYER
23! PRMUZ  : (KLON)             ; COSINE OF SOLAR ZENITH ANGLE
24! PTO1   : (KLON)             ; OPTICAL THICKNESS
25! PW     : (KLON)             ; SINGLE SCATTERING ALBEDO
26!     ==== OUTPUTS ===
27! PRE1   : (KLON)             ; LAYER REFLECTIVITY ASSUMING NO
28!                             ; REFLECTION FROM UNDERLYING LAYER
29! PTR1   : (KLON)             ; LAYER TRANSMISSIVITY ASSUMING NO
30!                             ; REFLECTION FROM UNDERLYING LAYER
31! PRE2   : (KLON)             ; LAYER REFLECTIVITY ASSUMING
32!                             ; REFLECTION FROM UNDERLYING LAYER
33! PTR2   : (KLON)             ; LAYER TRANSMISSIVITY ASSUMING
34!                             ; REFLECTION FROM UNDERLYING LAYER
35
36!        IMPLICIT ARGUMENTS :   NONE
37!        --------------------
38
39!     METHOD.
40!     -------
41
42!          STANDARD DELTA-EDDINGTON LAYER CALCULATIONS.
43
44!     EXTERNALS.
45!     ----------
46
47!          NONE
48
49!     REFERENCE.
50!     ----------
51
52!        SEE RADIATION'S PART OF THE MODEL'S DOCUMENTATION AND
53!        ECMWF RESEARCH DEPARTMENT DOCUMENTATION OF THE IFS
54
55!     AUTHOR.
56!     -------
57!        JEAN-JACQUES MORCRETTE  *ECMWF*
58
59!     MODIFICATIONS.
60!     --------------
61!        ORIGINAL: 88-12-15
62!                   96-05-30 Michel Deque (security in EXP())
63!        M.Hamrud      01-Oct-2003 CY28 Cleaning
64!        Modified: 03-10-10 Deborah Salmond and Marta Janiskova Optimisation
65!        Modified: 03-12-13 John Hague - MASS Vector Fns
66!        Y.Seity 06-09-09 : add modset from O.Thouron (MesoNH) under NOVLP tests
67!     ------------------------------------------------------------------
68
69!     ------------------------------------------------------------------
70
71!*       0.1   ARGUMENTS
72!              ---------
73
74USE PARKIND1  ,ONLY : JPIM     ,JPRB
75USE YOMHOOK   ,ONLY : LHOOK,   DR_HOOK
76
77USE YOERDU   , ONLY : REPLOG
78USE YOMJFH   , ONLY : N_VMASS
79!++MODIFCODE
80USE YOERAD   , ONLY : NOVLP
81!--MODIFCODE
82IMPLICIT NONE
83
84INTEGER(KIND=JPIM),INTENT(IN)    :: KLON
85INTEGER(KIND=JPIM),INTENT(IN)    :: KIDIA
86INTEGER(KIND=JPIM),INTENT(IN)    :: KFDIA
87REAL(KIND=JPRB)   ,INTENT(IN)    :: PGG(KLON)
88REAL(KIND=JPRB)   ,INTENT(IN)    :: PREF(KLON)
89REAL(KIND=JPRB)   ,INTENT(IN)    :: PRMUZ(KLON)
90REAL(KIND=JPRB)   ,INTENT(IN)    :: PTO1(KLON)
91REAL(KIND=JPRB)   ,INTENT(IN)    :: PW(KLON)
92REAL(KIND=JPRB)   ,INTENT(OUT)   :: PRE1(KLON)
93REAL(KIND=JPRB)   ,INTENT(OUT)   :: PRE2(KLON)
94REAL(KIND=JPRB)   ,INTENT(OUT)   :: PTR1(KLON)
95REAL(KIND=JPRB)   ,INTENT(OUT)   :: PTR2(KLON)
96REAL(KIND=JPRB) :: ZTMP   (4,KFDIA-KIDIA+1)
97REAL(KIND=JPRB) :: ZTMP2  (KFDIA-KIDIA+1+N_VMASS)
98REAL(KIND=JPRB) :: ZTMP3  (KFDIA-KIDIA+1+N_VMASS)
99REAL(KIND=JPRB) :: ZZARG  (KFDIA-KIDIA+1+N_VMASS)
100REAL(KIND=JPRB) :: ZZARG2 (KFDIA-KIDIA+1+N_VMASS)
101
102INTEGER(KIND=JPIM) :: JL, JLL, JLEN
103
104REAL(KIND=JPRB) :: ZA11, ZA12, ZA13, ZA21, ZA22, ZA23, ZALPHA,&
105 & ZAM2B, ZAP2B,  ZB21, ZB22, ZB23, &
106 & ZBETA, ZC1A, ZC1B, ZC2A, ZC2B, ZDENA, ZDENB, &
107 & ZDT, ZEXKM, ZEXKP, ZEXMU0, ZFF, ZGP, ZRI0A, &
108 & ZRI0B, ZRI0C, ZRI0D, ZRI1A, ZRI1B, ZRI1C, &
109 & ZRI1D, ZRK, ZRM2, ZRP, ZTOP, ZWCP, ZWM, ZX1, &
110 & ZX2, ZXM2P, ZXP2P 
111
112
113REAL(KIND=JPRB) :: MINJ, MAXJ, X, Y
114REAL(KIND=JPRB) :: ZPRMUZ,ZIDENA,ZIDENB,ZRR
115REAL(KIND=JPRB) :: ZHOOK_HANDLE
116
117!     STATEMENT DUNCTIONS
118MINJ(X,Y) = Y - 0.5_JPRB*(ABS(X-Y)-(X-Y))
119MAXJ(X,Y) = Y + 0.5_JPRB*(ABS(X-Y)+(X-Y))
120
121!     ------------------------------------------------------------------
122
123!*         1.      DELTA-EDDINGTON CALCULATIONS
124
125IF (LHOOK) CALL DR_HOOK('SWDE',0,ZHOOK_HANDLE)
126
127ZDT = 2.0_JPRB/3._JPRB
128
129  DO JL   =   KIDIA,KFDIA
130    JLL=JL-KIDIA+1
131    ZPRMUZ=1.0_JPRB/PRMUZ(JL)
132!++MODIFCODE
133  IF (NOVLP >= 5) THEN !MESONH_VERSION
134   ZGP = PGG(JL)
135   ZTOP = PTO1(JL)
136   ZWCP = PW(JL) 
137  ELSE !ECMWF VERSION
138    ZFF = PGG(JL)*PGG(JL)
139    ZGP = PGG(JL)/(1.0_JPRB+PGG(JL))
140    ZTOP = (1.0_JPRB- PW(JL) * ZFF) * PTO1(JL)
141    ZWCP = (1-ZFF)* PW(JL) /(1.0_JPRB- PW(JL) * ZFF)
142  ENDIF
143!--MODIFCODE
144    ZX1 = 1.0_JPRB-ZWCP*ZGP
145    ZWM = 1.0_JPRB-ZWCP
146    ZRM2 =  PRMUZ(JL) * PRMUZ(JL)
147    ZRK = SQRT(MAXJ(REPLOG,3._JPRB*ZWM*ZX1))
148    ZX2 = (1.0_JPRB-ZRK*ZRK*ZRM2)*(4._JPRB/3._JPRB)
149    ZRR = 1.0_JPRB/ZX2
150    ZRP=ZRK/ZX1
151    ZALPHA = ZWCP*ZRM2*(1.0_JPRB+ZGP*ZWM)*ZRR
152    ZBETA = ZWCP* PRMUZ(JL) *(1.0_JPRB+3._JPRB*ZGP*ZRM2*ZWM)*ZRR
153    ZZARG(JLL)   = -MAXJ( -200._JPRB, MINJ( ZTOP*ZPRMUZ, 200._JPRB) )
154    ZZARG2(JLL)  = MINJ( ZRK*ZTOP, 200._JPRB)
155    ZTMP(1,JLL) = ZPRMUZ
156    ZTMP(2,JLL) = ZALPHA
157    ZTMP(3,JLL) = ZBETA
158    ZTMP(4,JLL) = ZRP
159  ENDDO
160
161  IF(N_VMASS /= 0 ) THEN  !USING VECTOR MASS
162    JLEN=KFDIA-KIDIA+N_VMASS-MOD(KFDIA-KIDIA,N_VMASS)
163    IF(KFDIA-KIDIA+1 /= JLEN) THEN
164      ZZARG  (KFDIA-KIDIA+2:JLEN)=1.0_JPRB
165      ZZARG2 (KFDIA-KIDIA+2:JLEN)=1.0_JPRB
166    ENDIF
167! Commente par MPL le 21.11.08
168!   CALL VEXP(ZTMP2,ZZARG, JLEN)
169!   CALL VEXP(ZTMP3,ZZARG2,JLEN)
170  ELSE
171    DO JL   =   KIDIA,KFDIA
172      JLL=JL-KIDIA+1
173      ZTMP2(JLL) = EXP(ZZARG(JLL))
174      ZTMP3(JLL) = EXP(ZZARG2(JLL))
175    ENDDO
176  ENDIF
177
178  DO JL   =   KIDIA,KFDIA
179    JLL=JL-KIDIA+1
180    ZEXMU0 = ZTMP2(JLL)
181    ZEXKP  = ZTMP3(JLL)
182    ZPRMUZ = ZTMP(1,JLL)
183    ZALPHA = ZTMP(2,JLL)
184    ZBETA  = ZTMP(3,JLL)
185    ZRP    = ZTMP(4,JLL)
186    ZEXKM = 1.0_JPRB/ZEXKP
187    ZXP2P = 1.0_JPRB+ZDT*ZRP
188    ZXM2P = 1.0_JPRB-ZDT*ZRP
189    ZAP2B = ZALPHA+ZDT*ZBETA
190    ZAM2B = ZALPHA-ZDT*ZBETA
191
192!*         1.2     WITHOUT REFLECTION FROM THE UNDERLYING LAYER
193
194    ZA11 = ZXP2P
195    ZA12 = ZXM2P
196    ZA13 = ZAP2B
197    ZA22 = ZXP2P*ZEXKP
198    ZA21 = ZXM2P*ZEXKM
199    ZA23 = ZAM2B*ZEXMU0
200    ZDENA = ZA11 * ZA22 - ZA21 * ZA12
201    ZIDENA=1.0_JPRB/ZDENA
202    ZC1A = (ZA22*ZA13-ZA12*ZA23)*ZIDENA
203    ZC2A = (ZA11*ZA23-ZA21*ZA13)*ZIDENA
204    ZRI0A = ZC1A+ZC2A-ZALPHA
205    ZRI1A = ZRP*(ZC1A-ZC2A)-ZBETA
206    PRE1(JL) = (ZRI0A-ZDT*ZRI1A)*ZPRMUZ
207    ZRI0B = ZC1A*ZEXKM+ZC2A*ZEXKP-ZALPHA*ZEXMU0
208    ZRI1B = ZRP*(ZC1A*ZEXKM-ZC2A*ZEXKP)-ZBETA*ZEXMU0
209    PTR1(JL) = ZEXMU0+(ZRI0B+ZDT*ZRI1B)*ZPRMUZ
210
211!*         1.3     WITH REFLECTION FROM THE UNDERLYING LAYER
212
213    ZB21 = ZA21- PREF(JL) *ZXP2P*ZEXKM
214    ZB22 = ZA22- PREF(JL) *ZXM2P*ZEXKP
215    ZB23 = ZA23- PREF(JL) *ZEXMU0*(ZAP2B - PRMUZ(JL) )
216    ZDENB = ZA11 * ZB22 - ZB21 * ZA12
217    ZIDENB= 1.0_JPRB/ZDENB
218    ZC1B = (ZB22*ZA13-ZA12*ZB23)*ZIDENB
219    ZC2B = (ZA11*ZB23-ZB21*ZA13)*ZIDENB
220    ZRI0C = ZC1B+ZC2B-ZALPHA
221    ZRI1C = ZRP*(ZC1B-ZC2B)-ZBETA
222    PRE2(JL) = (ZRI0C-ZDT*ZRI1C) * ZPRMUZ
223    ZRI0D = ZC1B*ZEXKM + ZC2B*ZEXKP - ZALPHA*ZEXMU0
224    ZRI1D = ZRP * (ZC1B*ZEXKM - ZC2B*ZEXKP) - ZBETA*ZEXMU0
225    PTR2(JL) = ZEXMU0 + (ZRI0D + ZDT*ZRI1D) * ZPRMUZ
226  ENDDO
227
228IF (LHOOK) CALL DR_HOOK('SWDE',1,ZHOOK_HANDLE)
229END SUBROUTINE SWDE
230
Note: See TracBrowser for help on using the repository browser.