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