source: lmdz_wrf/trunk/WRFV3/phys/module_ra_sw.F @ 915

Last change on this file since 915 was 1, checked in by lfita, 10 years ago
  • -- --- Opening of the WRF+LMDZ coupling repository --- -- -

WRF: version v3.3
LMDZ: version v1818

More details in:

File size: 17.8 KB
Line 
1!WRF:MODEL_LAYER:PHYSICS
2!
3MODULE module_ra_sw
4
5      REAL,PRIVATE,SAVE :: CSSCA
6
7CONTAINS
8
9!------------------------------------------------------------------
10   SUBROUTINE SWRAD(dt,RTHRATEN,GSW,XLAT,XLONG,ALBEDO,            &
11                    rho_phy,T3D,QV3D,QC3D,QR3D,                   &
12                    QI3D,QS3D,QG3D,P3D,pi3D,dz8w,GMT,             &
13                    R,CP,G,JULDAY,                                &
14                    XTIME,DECLIN,SOLCON,                          &
15                    F_QV,F_QC,F_QR,F_QI,F_QS,F_QG,                &
16                    pm2_5_dry,pm2_5_water,pm2_5_dry_ec,           &
17                    RADFRQ,ICLOUD,DEGRAD,warm_rain,               &
18                    ids,ide, jds,jde, kds,kde,                    &
19                    ims,ime, jms,jme, kms,kme,                    &
20                    its,ite, jts,jte, kts,kte                     &
21                    )
22!------------------------------------------------------------------
23   IMPLICIT NONE
24!------------------------------------------------------------------
25   INTEGER,    INTENT(IN   ) ::        ids,ide, jds,jde, kds,kde, &
26                                       ims,ime, jms,jme, kms,kme, &
27                                       its,ite, jts,jte, kts,kte
28
29   LOGICAL,    INTENT(IN   ) ::        warm_rain
30   INTEGER,    INTENT(IN   ) ::        icloud
31
32   REAL, INTENT(IN    )      ::        RADFRQ,DEGRAD,             &
33                                       XTIME,DECLIN,SOLCON
34!
35   REAL, DIMENSION( ims:ime, kms:kme, jms:jme ),                  &
36         INTENT(IN    ) ::                                   P3D, &
37                                                            pi3D, &
38                                                         rho_phy, &
39                                                            dz8w, &
40                                                             T3D
41   REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), OPTIONAL ,       &
42         INTENT(IN    ) ::                             pm2_5_dry, &
43                                                     pm2_5_water, &
44                                                    pm2_5_dry_ec
45
46
47   REAL, DIMENSION( ims:ime, kms:kme, jms:jme ),                  &
48         INTENT(INOUT)  ::                              RTHRATEN
49!
50   REAL, DIMENSION( ims:ime, jms:jme ),                           &
51         INTENT(IN   )  ::                                  XLAT, &
52                                                           XLONG, &
53                                                          ALBEDO
54!
55   REAL, DIMENSION( ims:ime, jms:jme ),                           &
56         INTENT(INOUT)  ::                                   GSW
57!
58   REAL, INTENT(IN   )   ::                        GMT,R,CP,G,dt
59!
60   INTEGER, INTENT(IN  ) ::                               JULDAY 
61
62
63
64!
65! Optional
66!
67   REAL, DIMENSION( ims:ime, kms:kme, jms:jme ),                  &
68         OPTIONAL,                                                &
69         INTENT(IN    ) ::                                        &
70                                                            QV3D, &
71                                                            QC3D, &
72                                                            QR3D, &
73                                                            QI3D, &
74                                                            QS3D, &
75                                                            QG3D
76
77   LOGICAL, OPTIONAL, INTENT(IN )      ::        F_QV,F_QC,F_QR,F_QI,F_QS,F_QG
78 
79! LOCAL VARS
80 
81   REAL, DIMENSION( kts:kte ) ::                                  &
82                                                          TTEN1D, &
83                                                          RHO01D, &
84                                                             P1D, &
85                                                              DZ, &
86                                                             T1D, &
87                                                            QV1D, &
88                                                            QC1D, &
89                                                            QR1D, &
90                                                            QI1D, &
91                                                            QS1D, &
92                                                            QG1D
93!
94   REAL::      XLAT0,XLONG0,ALB0,GSW0
95
96!
97   INTEGER :: i,j,K,NK
98   LOGICAL :: predicate , do_topo_shading
99   real :: aer_dry1(kts:kte),aer_water1(kts:kte)
100
101!------------------------------------------------------------------
102
103   j_loop: DO J=jts,jte
104   i_loop: DO I=its,ite
105
106! reverse vars
107         DO K=kts,kte
108            QV1D(K)=0.
109            QC1D(K)=0.
110            QR1D(K)=0.
111            QI1D(K)=0.
112            QS1D(K)=0.
113            QG1D(K)=0.
114         ENDDO
115
116         DO K=kts,kte
117            NK=kme-1-K+kms
118            TTEN1D(K)=0.
119
120            T1D(K)=T3D(I,NK,J)
121            P1D(K)=P3D(I,NK,J)
122            RHO01D(K)=rho_phy(I,NK,J)
123            DZ(K)=dz8w(I,NK,J)
124         ENDDO
125
126         IF( PRESENT(pm2_5_dry) .AND. PRESENT(pm2_5_water) )THEN
127            DO K=kts,kte
128               NK=kme-1-K+kms
129               aer_dry1(k)   = pm2_5_dry(i,nk,j)
130               aer_water1(k) = pm2_5_water(i,nk,j)
131            ENDDO
132         ELSE
133            DO K=kts,kte
134               aer_dry1(k)   = 0.
135               aer_water1(k) = 0.
136            ENDDO
137         ENDIF
138
139         IF (PRESENT(F_QV) .AND. PRESENT(QV3D)) THEN
140            IF (F_QV) THEN
141               DO K=kts,kte
142                  NK=kme-1-K+kms
143                  QV1D(K)=QV3D(I,NK,J)
144                  QV1D(K)=max(0.,QV1D(K))
145               ENDDO
146            ENDIF
147         ENDIF
148
149         IF (PRESENT(F_QC) .AND. PRESENT(QC3D)) THEN
150            IF (F_QC) THEN
151               DO K=kts,kte
152                  NK=kme-1-K+kms
153                  QC1D(K)=QC3D(I,NK,J)
154                  QC1D(K)=max(0.,QC1D(K))
155               ENDDO
156            ENDIF
157         ENDIF
158
159         IF (PRESENT(F_QR) .AND. PRESENT(QR3D)) THEN
160            IF (F_QR) THEN
161               DO K=kts,kte
162                  NK=kme-1-K+kms
163                  QR1D(K)=QR3D(I,NK,J)
164                  QR1D(K)=max(0.,QR1D(K))
165               ENDDO
166            ENDIF
167         ENDIF
168
169!
170         IF ( PRESENT( F_QI ) ) THEN
171            predicate = F_QI
172         ELSE
173            predicate = .FALSE.
174         ENDIF
175
176         IF ( predicate .AND. PRESENT( QI3D ) ) THEN
177            DO K=kts,kte
178               NK=kme-1-K+kms
179               QI1D(K)=QI3D(I,NK,J)
180               QI1D(K)=max(0.,QI1D(K))
181            ENDDO
182         ELSE
183            IF (.not. warm_rain) THEN
184               DO K=kts,kte
185               IF(T1D(K) .lt. 273.15) THEN
186                  QI1D(K)=QC1D(K)
187                  QC1D(K)=0.
188                  QS1D(K)=QR1D(K)
189                  QR1D(K)=0.
190               ENDIF
191               ENDDO
192            ENDIF
193         ENDIF
194
195         IF (PRESENT(F_QS) .AND. PRESENT(QS3D)) THEN
196            IF (F_QS) THEN
197               DO K=kts,kte         
198                  NK=kme-1-K+kms
199                  QS1D(K)=QS3D(I,NK,J)
200                  QS1D(K)=max(0.,QS1D(K))
201               ENDDO
202            ENDIF
203         ENDIF
204
205         IF (PRESENT(F_QG) .AND. PRESENT(QG3D)) THEN
206            IF (F_QG) THEN
207               DO K=kts,kte         
208                  NK=kme-1-K+kms
209                  QG1D(K)=QG3D(I,NK,J)
210                  QG1D(K)=max(0.,QG1D(K))
211               ENDDO
212            ENDIF
213         ENDIF
214
215         XLAT0=XLAT(I,J)
216         XLONG0=XLONG(I,J)
217         ALB0=ALBEDO(I,J)
218! slope code removed - factor now done in surface driver
219           CALL SWPARA(TTEN1D,GSW0,XLAT0,XLONG0,ALB0,              &
220                       T1D,QV1D,QC1D,QR1D,QI1D,QS1D,QG1D,P1D,      &
221                       XTIME,GMT,RHO01D,DZ,                        &
222                       R,CP,G,DECLIN,SOLCON,                       &
223                       RADFRQ,ICLOUD,DEGRAD,aer_dry1,aer_water1,   &
224                       kts,kte      )
225         GSW(I,J)=GSW0
226         DO K=kts,kte         
227            NK=kme-1-K+kms
228            RTHRATEN(I,K,J)=RTHRATEN(I,K,J)+TTEN1D(NK)/pi3D(I,K,J)
229         ENDDO
230!
231   ENDDO i_loop
232   ENDDO j_loop                                         
233
234   END SUBROUTINE SWRAD
235
236!------------------------------------------------------------------
237   SUBROUTINE SWPARA(TTEN,GSW,XLAT,XLONG,ALBEDO,                  &
238                     T,QV,QC,QR,QI,QS,QG,P,                       &
239                     XTIME, GMT, RHO0, DZ,                        &
240                     R,CP,G,DECLIN,SOLCON,                        &
241                     RADFRQ,ICLOUD,DEGRAD,aer_dry1,aer_water1,    &
242                     kts,kte,slope_rad,shadow,slp_azi,slope       )
243!------------------------------------------------------------------
244!     TO CALCULATE SHORT-WAVE ABSORPTION AND SCATTERING IN CLEAR
245!     AIR AND REFLECTION AND ABSORPTION IN CLOUD LAYERS (STEPHENS,
246!     1984)
247!     CHANGES:
248!       REDUCE EFFECTS OF ICE CLOUDS AND PRECIP ON LIQUID WATER PATH
249!       ADD EFFECT OF GRAUPEL
250!------------------------------------------------------------------
251
252  IMPLICIT NONE
253
254  INTEGER, INTENT(IN ) ::                 kts,kte
255!
256  REAL, DIMENSION( kts:kte ), INTENT(IN   )  ::                   &
257                                                            RHO0, &
258                                                               T, &
259                                                               P, &
260                                                              DZ, &
261                                                              QV, &
262                                                              QC, &
263                                                              QR, &
264                                                              QI, &
265                                                              QS, &
266                                                              QG
267
268   REAL, DIMENSION( kts:kte ), INTENT(INOUT)::              TTEN
269!
270   REAL, INTENT(IN  )   ::               XTIME,GMT,R,CP,G,DECLIN, &
271                                        SOLCON,XLAT,XLONG,ALBEDO, &
272                                                  RADFRQ, DEGRAD
273!
274   INTEGER, INTENT(IN) :: icloud
275   REAL, INTENT(INOUT)  ::                                   GSW
276! For slope-dependent radiation
277
278   INTEGER, OPTIONAL, INTENT(IN) :: slope_rad,shadow
279   REAL, OPTIONAL,    INTENT(IN) :: slp_azi,slope
280
281! LOCAL VARS
282!
283   REAL, DIMENSION( kts:kte+1 ) ::                         SDOWN
284
285   REAL, DIMENSION( kts:kte )   ::                          XLWP, &
286                                                            XATP, &
287                                                            XWVP, &
288                                             aer_dry1,aer_water1, &
289                                                              RO
290!
291   REAL, DIMENSION( 4, 5 ) ::                             ALBTAB, &
292                                                          ABSTAB
293
294   REAL, DIMENSION( 4    ) ::                             XMUVAL
295
296   REAL :: beta
297
298!------------------------------------------------------------------
299
300      DATA ALBTAB/0.,0.,0.,0., &
301           69.,58.,40.,15.,    &
302           90.,80.,70.,60.,    &
303           94.,90.,82.,78.,    &
304           96.,92.,85.,80./
305
306      DATA ABSTAB/0.,0.,0.,0., &
307           0.,2.5,4.,5.,       &
308           0.,2.6,7.,10.,      &
309           0.,3.3,10.,14.,     &
310           0.,3.7,10.,15./
311
312      DATA XMUVAL/0.,0.2,0.5,1.0/
313
314      REAL :: bext340, absc, alba, alw, csza,dabsa,dsca,dabs
315      REAL :: bexth2o, dscld, hrang,ff,oldalb,oldabs,oldabc
316      REAL :: soltop, totabs, tloctm, ugcm, uv,xabs,xabsa,wv
317      REAL :: wgm, xalb, xi, xsca, xt24,xmu,xabsc,trans0,yj
318      REAL :: xxlat,ww
319      INTEGER :: iil,ii,jjl,ju,k,iu
320
321! For slope-dependent radiation
322
323   REAL :: diffuse_frac, corr_fac, csza_slp
324
325
326      GSW=0.0
327      bext340=5.E-6
328      bexth2o=5.E-6
329      SOLTOP=SOLCON
330      XT24=MOD(XTIME+RADFRQ*0.5,1440.)
331      TLOCTM=GMT+XT24/60.+XLONG/15.
332      HRANG=15.*(TLOCTM-12.)*DEGRAD
333      XXLAT=XLAT*DEGRAD
334      CSZA=SIN(XXLAT)*SIN(DECLIN)+COS(XXLAT)*COS(DECLIN)*COS(HRANG)
335
336!     RETURN IF NIGHT       
337      IF(CSZA.LE.1.E-9)GOTO 7
338!
339      DO K=kts, kte
340
341! P in the unit of 10mb
342         RO(K)=P(K)/(R*T(K))
343         XWVP(K)=RO(K)*QV(K)*DZ(K)*1000.
344! KG/M**2
345          XATP(K)=RO(K)*DZ(K)
346      ENDDO
347!
348!     G/M**2
349!     REDUCE WEIGHT OF LIQUID AND ICE IN SHORT-WAVE SCHEME
350!     ADD GRAUPEL EFFECT (ASSUMED SAME AS RAIN)
351!
352      IF (ICLOUD.EQ.0)THEN
353         DO K=kts, kte
354            XLWP(K)=0.
355         ENDDO
356      ELSE
357         DO K=kts, kte
358            XLWP(K)=RO(K)*1000.*DZ(K)*(QC(K)+0.1*QI(K)+0.05* &
359                    QR(K)+0.02*QS(K)+0.05*QG(K))
360         ENDDO
361      ENDIF
362!
363      XMU=CSZA
364      SDOWN(1)=SOLTOP*XMU
365!     SET WW (G/M**2) LIQUID WATER PATH INTEGRATED DOWN
366!     SET UV (G/M**2) WATER VAPOR PATH INTEGRATED DOWN
367      WW=0.
368      UV=0.
369      OLDALB=0.
370      OLDABC=0.
371      TOTABS=0.
372!     CONTRIBUTIONS DUE TO CLEAR AIR AND CLOUD
373      DSCA=0.
374      DABS=0.
375      DSCLD=0.
376!
377! CONTRIBUTION DUE TO AEROSOLS (FOR CHEMISTRY)
378      DABSA=0.
379!
380      DO 200 K=kts,kte
381         WW=WW+XLWP(K)
382         UV=UV+XWVP(K)
383!     WGM IS WW/COS(THETA) (G/M**2)
384!     UGCM IS UV/COS(THETA) (G/CM**2)
385         WGM=WW/XMU
386         UGCM=UV*0.0001/XMU
387!
388         OLDABS=TOTABS
389!     WATER VAPOR ABSORPTION AS IN LACIS AND HANSEN (1974)
390         TOTABS=2.9*UGCM/((1.+141.5*UGCM)**0.635+5.925*UGCM)
391!     APPROXIMATE RAYLEIGH + AEROSOL SCATTERING
392!        XSCA=1.E-5*XATP(K)/XMU
393!          XSCA=(1.E-5*XATP(K)+aer_dry1(K)*bext340+aer_water1(K)*bexth2o)/XMU
394         beta=0.4*(1.0-XMU)+0.1
395!     CSSCA - CLEAR-SKY SCATTERING SET FROM NAMELIST SWRAD_SCAT
396         XSCA=(cssca*XATP(K)+beta*aer_dry1(K)*bext340*DZ(K) &
397              +beta*aer_water1(K)*bexth2o*DZ(K))/XMU   
398
399!     LAYER VAPOR ABSORPTION DONE FIRST
400         XABS=(TOTABS-OLDABS)*(SDOWN(1)-DSCLD-DSCA-DABSA)/SDOWN(K)
401!rs   AEROSOL ABSORB (would be elemental carbon). So far XABSA = 0.
402         XABSA=0.
403         IF(XABS.LT.0.)XABS=0.
404!
405         ALW=ALOG10(WGM+1.)
406         IF(ALW.GT.3.999)ALW=3.999
407!
408         DO II=1,3
409            IF(XMU.GT.XMUVAL(II))THEN
410              IIL=II
411              IU=II+1
412              XI=(XMU-XMUVAL(II))/(XMUVAL(II+1)-XMUVAL(II))+FLOAT(IIL)
413            ENDIF
414         ENDDO
415!
416         JJL=IFIX(ALW)+1
417         JU=JJL+1
418         YJ=ALW+1.
419!     CLOUD ALBEDO
420         ALBA=(ALBTAB(IU,JU)*(XI-IIL)*(YJ-JJL)   &
421              +ALBTAB(IIL,JU)*(IU-XI)*(YJ-JJL)   &
422              +ALBTAB(IU,JJL)*(XI-IIL)*(JU-YJ)   &
423              +ALBTAB(IIL,JJL)*(IU-XI)*(JU-YJ))  &
424             /((IU-IIL)*(JU-JJL))
425!     CLOUD ABSORPTION
426         ABSC=(ABSTAB(IU,JU)*(XI-IIL)*(YJ-JJL)   &
427              +ABSTAB(IIL,JU)*(IU-XI)*(YJ-JJL)   &
428              +ABSTAB(IU,JJL)*(XI-IIL)*(JU-YJ)   &
429              +ABSTAB(IIL,JJL)*(IU-XI)*(JU-YJ))  &
430             /((IU-IIL)*(JU-JJL))
431!     LAYER ALBEDO AND ABSORPTION
432         XALB=(ALBA-OLDALB)*(SDOWN(1)-DSCA-DABS)/SDOWN(K)
433         XABSC=(ABSC-OLDABC)*(SDOWN(1)-DSCA-DABS)/SDOWN(K)
434         IF(XALB.LT.0.)XALB=0.
435         IF(XABSC.LT.0.)XABSC=0.
436         DSCLD=DSCLD+(XALB+XABSC)*SDOWN(K)*0.01
437         DSCA=DSCA+XSCA*SDOWN(K)
438         DABS=DABS+XABS*SDOWN(K)
439         DABSA=DABSA+XABSA*SDOWN(K)
440         OLDALB=ALBA
441         OLDABC=ABSC
442!     LAYER TRANSMISSIVITY
443         TRANS0=100.-XALB-XABSC-XABS*100.-XSCA*100.
444         IF(TRANS0.LT.1.)THEN
445           FF=99./(XALB+XABSC+XABS*100.+XSCA*100.)
446           XALB=XALB*FF
447           XABSC=XABSC*FF
448           XABS=XABS*FF
449           XSCA=XSCA*FF
450           TRANS0=1.
451         ENDIF
452         SDOWN(K+1)=AMAX1(1.E-9,SDOWN(K)*TRANS0*0.01)
453         TTEN(K)=SDOWN(K)*(XABSC+XABS*100.+XABSA*100.)*0.01/( &
454                 RO(K)*CP*DZ(K))
455  200   CONTINUE
456!
457        GSW=(1.-ALBEDO)*SDOWN(kte+1)
458
459    IF (PRESENT(slope_rad)) THEN
460! Slope-dependent solar radiation part
461
462      if (slope_rad.eq.1) then
463
464!  Parameterize diffuse fraction of global solar radiation as a function of the ratio between TOA radiation and surface global radiation
465
466        diffuse_frac = min(1.,1/(max(0.1,2.1-2.8*log(log(SDOWN(kts)/max(SDOWN(kte+1),1.e-3))))))
467        if ((slope.eq.0).or.(diffuse_frac.eq.1).or.(csza.lt.1.e-2)) then  ! no topographic effects when all radiation is diffuse or the sun is too close to the horizon
468        corr_fac = 1
469        goto 140
470        endif
471
472! cosine of zenith angle over sloping topography
473
474        csza_slp = ((SIN(XXLAT)*COS(HRANG))*                                          &
475                    (-cos(slp_azi)*sin(slope))-SIN(HRANG)*(sin(slp_azi)*sin(slope))+  &
476                    (COS(XXLAT)*COS(HRANG))*cos(slope))*                              &
477                   COS(DECLIN)+(COS(XXLAT)*(cos(slp_azi)*sin(slope))+                 &
478                   SIN(XXLAT)*cos(slope))*SIN(DECLIN)
479        IF(csza_slp.LE.1.E-4) csza_slp = 0
480
481! Topographic shading
482
483        if (shadow.eq.1) csza_slp = 0
484
485! Correction factor for sloping topography; the diffuse fraction of solar radiation is assumed to be unaffected by the slope
486        corr_fac = diffuse_frac + (1-diffuse_frac)*csza_slp/csza
487
488 140    continue   
489
490        GSW=(1.-ALBEDO)*SDOWN(kte+1)*corr_fac
491       
492      endif
493    ENDIF
494
495    7 CONTINUE
496!
497   END SUBROUTINE SWPARA
498
499!====================================================================
500   SUBROUTINE swinit(swrad_scat,                                    &
501                     allowed_to_read ,                              &
502                     ids, ide, jds, jde, kds, kde,                  &
503                     ims, ime, jms, jme, kms, kme,                  &
504                     its, ite, jts, jte, kts, kte                   )
505!--------------------------------------------------------------------
506   IMPLICIT NONE
507!--------------------------------------------------------------------
508   LOGICAL , INTENT(IN)           :: allowed_to_read
509   INTEGER , INTENT(IN)           :: ids, ide, jds, jde, kds, kde,  &
510                                     ims, ime, jms, jme, kms, kme,  &
511                                     its, ite, jts, jte, kts, kte
512
513   REAL , INTENT(IN)              :: swrad_scat
514
515!     CSSCA - CLEAR-SKY SCATTERING SET FROM NAMELIST SWRAD_SCAT
516   cssca = swrad_scat * 1.e-5
517
518   END SUBROUTINE swinit
519
520END MODULE module_ra_sw
Note: See TracBrowser for help on using the repository browser.