source: LMDZ6/branches/contrails/libf/phylmd/inlandsis/sisvat_sno_albedo.f90 @ 5445

Last change on this file since 5445 was 5246, checked in by abarral, 2 months ago

Convert fixed-form to free-form sources .F -> .{f,F}90
(WIP: some .F remain, will be handled in subsequent commits)

File size: 30.9 KB
Line 
1subroutine SnOptP(jjtime)
2
3  ! +------------------------------------------------------------------------+
4  ! | MAR/SISVAT   SnOptP                                    12-08-2019  MAR |
5  ! |   SubRoutine SnOptP computes the Snow Pack optical Properties          |
6  ! +------------------------------------------------------------------------+
7  ! |                                                                        |
8  ! |   PARAMETERS:  knonv: Total Number of columns =                        |
9  ! |   ^^^^^^^^^^        = Total Number of continental     Grid Boxes       |
10  ! |                     X       Number of Mosaic Cell per Grid Box         |
11  ! |                                                                        |
12  ! |   INPUT:   isnoSV   = total Nb of Ice/Snow Layers                      |
13  ! |   ^^^^^    ispiSV   = 0,...,nsno: Uppermost Superimposed Ice Layer     |
14  ! |                                                                        |
15  ! |                                                                        |
16  ! |   INPUT:   G1snSV   : Dendricity (<0) or Sphericity (>0) of Snow Layer |
17  ! |   ^^^^^    G2snSV   : Sphericity (>0) or Size            of Snow Layer |
18  ! |            agsnSV   : Snow       Age                             [day] |
19  ! |            ro__SV   : Snow/Soil  Volumic Mass                  [kg/m3] |
20  ! |            eta_SV   : Water      Content                       [m3/m3] |
21  ! |            rusnSV   : Surficial  Water   Thickness   [kg/m2] .OR. [mm] |
22  ! |            SWS_SV   : Surficial  Water   Status                        |
23  ! |            dzsnSV   : Snow       Layer   Thickness                 [m] |
24  ! |                                                                        |
25  ! |            albssv   : Soil       Albedo                            [-] |
26  ! |            zzsnsv   : Snow       Pack    Thickness                 [m] |
27  ! |                                                                        |
28  ! |   OUTPUT:  albisv   : Snow/Ice/Water/Soil Integrated Albedo        [-] |
29  ! |   ^^^^^^   sEX_sv   : Verticaly Integrated Extinction Coefficient      |
30  ! |            DOPsnSV   : Snow Optical diameter [m]                       |
31  ! |                                                                        |
32  ! |   Internal Variables:                                                  |
33  ! |   ^^^^^^^^^^^^^^^^^^                                                   |
34  ! |            SnOpSV   : Snow Grain optical Size                      [m] |
35  ! |            EX1_sv   : Integrated Snow Extinction (0.3--0.8micr.m)      |
36  ! |            EX2_sv   : Integrated Snow Extinction (0.8--1.5micr.m)      |
37  ! |            EX3_sv   : Integrated Snow Extinction (1.5--2.8micr.m)      |
38  ! |                                                                        |
39  ! |   METHODE:    Calcul de la taille optique des grains ? partir de       |
40  ! |   ^^^^^^^    -leur type decrit par les deux variables descriptives     |
41  ! |                    continues sur la plage -99/+99 passees en appel.    |
42  ! |              -la taille optique (1/10mm) des etoiles,                  |
43  ! |                                          des grains fins et            |
44  ! |                                          des jeunes faces planes       |
45  ! |                                                                        |
46  ! |   METHOD:     Computation of the optical diameter of the grains        |
47  ! |   ^^^^^^      described with the CROCUS formalism G1snSV / G2snSV      |
48  ! |                                                                        |
49  ! |   REFERENCE: Brun et al.      1989, J. Glaciol 35 pp. 333--342         |
50  ! |   ^^^^^^^^^  Brun et al.      1992, J. Glaciol 38 pp.  13-- 22         |
51  ! |              Eric Martin Sept.1996                                     |
52  ! |                                                                        |
53  ! |                                                                        |
54  ! +------------------------------------------------------------------------+
55
56
57
58
59  ! +--Global Variables
60  ! +  ================
61
62
63  use VARphy
64  use VAR_SV
65  use VARdSV
66  use VARxSV
67  use VARySV
68  use VARtSV
69  USE surface_data, only: iflag_albcalc,correc_alb
70
71  IMPLICIT NONE
72
73
74  ! + -- INPUT
75  integer :: jjtime
76
77  ! +--Internal Variables
78  ! +  ==================
79
80  real :: coalb1(knonv)                      ! weighted Coalbedo, Vis.
81  real :: coalb2(knonv)                      ! weighted Coalbedo, nIR 1
82  real :: coalb3(knonv)                      ! weighted Coalbedo, nIR 2
83  real :: coalbm                             ! weighted Coalbedo, mean
84  real :: sExt_1(knonv)                      ! Extinction Coeff., Vis.
85  real :: sExt_2(knonv)                      ! Extinction Coeff., nIR 1
86  real :: sExt_3(knonv)                      ! Extinction Coeff., nIR 2
87  real :: SnOpSV(knonv,      nsno)           ! Snow Grain optical Size
88  ! #AG real      agesno
89
90  integer :: isn   ,ikl   ,isn1, i
91  real :: sbeta1,sbeta2,sbeta3,sbeta4,sbeta5
92  real :: AgeMax,AlbMin,HSnoSV,HIceSV,doptmx,SignG1,Sph_OK
93  real :: dalbed,dalbeS,dalbeW
94  real :: bsegal,czemax,csegal,csza
95  real :: RoFrez,DiffRo,SignRo,SnowOK,OpSqrt
96  real :: albSn1,albIc1,a_SnI1,a_SII1
97  real :: albSn2,albIc2,a_SnI2,a_SII2
98  real :: albSn3,albIc3,a_SnI3,a_SII3
99  real :: albSno,albIce,albSnI,albSII,albWIc,albmax
100  real :: doptic,Snow_H,SIce_H,SnownH,SIcenH
101  real :: exarg1,exarg2,exarg3,sign_0,sExt_0
102  real :: albedo_old,albCor
103  real :: ro_ave,dz_ave,minalb
104  real :: l1min,l1max,l2min,l2max,l3min,l3max
105  real :: l6min(6), l6max(6), albSn6(6), a_SII6(6)
106  real :: lmintmp,lmaxtmp,albtmp
107
108  ! +--Local   DATA
109  ! +  ============
110
111  ! +--For the computation of the solar irradiance extinction in snow
112  ! +  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
113  data      sbeta1/0.0192/,sbeta2/0.4000/,sbeta3/0.1098/
114  data      sbeta4/1.0000/
115  data      sbeta5/2.00e1/
116
117  ! +--Snow Age Maximum (Taiga, e.g. Col de Porte)
118  ! +  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
119  data      AgeMax  /60.0/
120  ! +...          AgeMax:  Snow Age Maximum                              [day]
121
122  data      AlbMin  /0.94/
123  ! +...          AlbMin:  Albedo   Minimum / visible (0.3--0.8 micrometers)
124
125  data      HSnoSV  /0.01/
126  ! +...          HSnoSV:  Snow     Thickness through witch
127  ! +                      Albedo is interpolated to Ice  Albedo
128  data      HIceSV  /0.10/
129  ! +...          HIceSV:  Snow/Ice Thickness through witch
130  ! +                      Albedo is interpolated to Soil Albedo
131  data      doptmx  /2.3e-3/
132  ! +...          doptmx:  Maximum  optical Diameter    (pi * R**2)        [m]
133  ! +
134  data      czeMAX  /0.173648178/  ! 80.deg (Segal et al., 1991 JAS)
135
136  data      bsegal  /4.00       /  !
137  data      albmax  /0.99       /  ! Albedo max
138
139  ! +-- wavelength limits [m] for each broad band
140
141  data      l1min/400.0e-9/,l1max/800.0e-9/
142  data      l2min/800.0e-9/,l2max/1500.0e-9/
143  data      l3min/1500.0e-9/,l3max/2800.0e-9/
144
145  data      l6min/185.0e-9,250.0e-9,400.0e-9, &
146        690.0e-9,1190.0e-9,2380.0e-9/
147  data      l6max/250.0e-9,400.0e-9, &
148        690.0e-9,1190.0e-9,2380.0e-9,4000.0e-9/
149
150
151  ! +--Snow Grain optical Size
152  ! +  =======================
153
154    DO ikl=1,knonv
155     DO   isn=1,max(1,isnoSV(ikl))
156
157      G2snSV(ikl,isn) =  max(epsi,G2snSV(ikl,isn))
158  ! +...    Avoid non physical Values
159
160      SignG1          = sign(unun,G1snSV(ikl,isn))
161      Sph_OK          =  max(zero,SignG1)
162
163      SnOpSV(ikl,isn) =   1.e-4 * &
164  ! +...    SI:           (from 1/10 mm to m)
165
166
167  ! +--Contribution of Non Dendritic Snow
168  ! +  ----------------------------------
169
170            (    Sph_OK *(      G2snSV(ikl,isn)*G1snSV(ikl,isn)/G1_dSV &
171            +max(demi*G2snSV(ikl,isn),DFcdSV) &
172            *(unun-G1snSV(ikl,isn)                /G1_dSV)) &
173
174
175  ! +--Contribution of     Dendritic Snow
176  ! +  ----------------------------------
177
178            +(1.-Sph_OK)*(     -G1snSV(ikl,isn)*DDcdSV         /G1_dSV &
179            +(unun+G1snSV(ikl,isn)                /G1_dSV) &
180            *    (G2snSV(ikl,isn)*DScdSV         /G1_dSV &
181            +(unun-G2snSV(ikl,isn)                /G1_dSV) &
182            *DFcdSV                 )))
183      SnOpSV(ikl,isn) =  max(zero,SnOpSV(ikl,isn))
184
185  ! + --For outputs (Etienne)
186  ! + ------------------------
187      DOPsnSV(ikl,isn)=SnOpSV(ikl,isn)
188    END DO
189  END DO
190
191
192
193
194  ! +--Snow/Ice Albedo
195  ! +  ===============
196
197
198
199  ! +--Uppermost effective Snow Layer
200  ! +  ------------------------------
201
202    DO ikl=1,knonv
203
204      isn   =  max(iun,isnoSV(ikl))
205
206      SignRo = sign(unun, rocdSV - ro__SV(ikl,isn))
207      SnowOK =  max(zero,SignRo) ! Ice Density Threshold
208
209      OpSqrt = sqrt(SnOpSV(ikl,isn))
210
211  !CA    +--Correction of snow albedo for Antarctica/Greenland
212  !CA       --------------------------------------------------
213
214
215      albCor = correc_alb
216  ! #GL     albCor = 1.01
217  ! #AC    albCor = 1.01
218
219
220      IF (iflag_albcalc .GE. 1) THEN  ! Albedo calculation according to Kokhanovsky and Zege 2004
221
222      dalbed = 0.0
223      doptic=SnOpSV(ikl,isn)
224      csza=coszSV(ikl)
225
226      CALL albedo_kokhanovsky(l1min,l1max,csza,doptic,albSn1)
227      CALL albedo_kokhanovsky(l2min,l2max,csza,doptic,albSn2)
228      CALL albedo_kokhanovsky(l3min,l3max,csza,doptic,albSn3)
229
230      DO i=1,6
231         lmintmp=l6min(i)
232         lmaxtmp=l6max(i)
233         CALL albedo_kokhanovsky(lmintmp,lmaxtmp,csza,doptic,albtmp)
234         albSn6(i)=albtmp
235      ENDDO
236
237
238      ELSE ! Default calculation in SISVAT
239
240  !    Zenith Angle Correction (Segal et al.,         1991, JAS 48, p.1025)
241  !--------------------------- (Wiscombe & Warren, dec1980, JAS   , p.2723)
242  !                            (Warren,               1982,  RG   , p.  81)
243                         ! -------------------------------------------------
244
245      dalbed = 0.0
246
247      csegal = max(czemax                   ,coszSV(ikl))
248  ! #cz     dalbeS =   ((bsegal+unun)/(unun+2.0*bsegal*csegal)
249  ! #cz.                -       unun                          )*0.32
250  ! #cz.              /  bsegal
251  ! #cz     dalbeS = max(dalbeS,zero)
252  ! #cz     dalbed =     dalbeS      *       min(1,isnoSV(ikl))
253
254      dalbeW =(0.64 - csegal  )*0.0625  ! Warren 1982, RevGeo, fig.12b
255                                        ! ! 0.0625 = 5% * 1/0.8,   p.81
256                                        ! ! 0.64   = cos(50)
257      dalbed =     dalbeW      *       min(1,isnoSV(ikl))
258  !-------------------------------------------------------------------------
259
260      albSn1 =  0.96-1.580*OpSqrt
261      albSn1 =  max(albSn1,AlbMin)
262
263      albSn1 =  max(albSn1,zero)
264      albSn1 =  min(albSn1*albCor,unun)
265
266      albSn2 =  0.95-15.40*OpSqrt
267      albSn2 =  max(albSn2,zero)
268      albSn2 =  min(albSn2*albCor,unun)
269
270      doptic =  min(SnOpSV(ikl,isn),doptmx)
271      albSn3 =  346.3*doptic -32.31*OpSqrt +0.88
272      albSn3 =  max(albSn3,zero)
273      albSn3 =  min(albSn3*albCor,unun)
274
275      albSn6(1:3)=albSn1
276      albSn6(4:6)=albSn2
277
278      ! !snow albedo corection if wetsnow
279  ! #GL     albSn1 =  albSn1*max(0.9,(1.-1.5*eta_SV(ikl,isn)))
280  ! #GL     albSn2 =  albSn2*max(0.9,(1.-1.5*eta_SV(ikl,isn)))
281  ! #GL     albSn3 =  albSn3*max(0.9,(1.-1.5*eta_SV(ikl,isn)))
282
283      ENDIF
284
285
286      albSno =  So1dSV*albSn1 &
287            +  So2dSV*albSn2 &
288            +  So3dSV*albSn3
289
290  !XF
291      minalb =  (aI2dSV + (aI3dSV -aI2dSV) &
292            *  (ro__SV(ikl,isn)-ro_Ice)/(roSdSV-ro_Ice))
293      minalb =  min(aI3dSV,max(aI2dSV,minalb)) ! pure/firn albedo
294
295      SnowOK =  SnowOK*max(zero,sign(unun,     albSno-minalb))
296      albSn1 =  SnowOK*albSn1+(1.0-SnowOK)*max(albSno,minalb)
297      albSn2 =  SnowOK*albSn2+(1.0-SnowOK)*max(albSno,minalb)
298      albSn3 =  SnowOK*albSn3+(1.0-SnowOK)*max(albSno,minalb)
299      albSn6(:) =  SnowOK*albSn6(:)+(1.0-SnowOK)*max(albSno,minalb)
300
301
302  ! +           ro < roSdSV |          min al > aI3dSV
303  ! +  roSdSV < ro < rocdSV | aI2dSV < min al < aI3dSV (fct of density)
304
305
306  ! +--Snow/Ice Pack Thickness
307  ! +  -----------------------
308
309      isn    =    max(min(isnoSV(ikl) ,ispiSV(ikl)),0)
310      Snow_H =  zzsnsv(ikl,isnoSV(ikl))-zzsnsv(ikl,isn)
311      SIce_H =  zzsnsv(ikl,isnoSV(ikl))
312      SnownH =  Snow_H  /  HSnoSV
313      SnownH =  min(unun,  SnownH)
314      SIcenH =  SIce_H  / (HIceSV)
315      SIcenH =  min(unun,  SIcenH)
316
317  ! +       The value of SnownH is set to 1 in case of ice lenses above
318  ! +       1m of dry snow (ro<600kg/m3) for using CROCUS albedo
319
320      ! ro_ave =  0.
321      ! dz_ave =  0.
322      ! SnowOK =  1.
323   ! do isn    =  isnoSV(ikl),1,-1
324   !    ro_ave =  ro_ave + ro__SV(ikl,isn) * dzsnSV(ikl,isn) * SnowOK
325   !    dz_ave =  dz_ave +                   dzsnSV(ikl,isn) * SnowOK
326   !    SnowOK =  max(zero,sign(unun,1.-dz_ave))
327   ! enddo
328
329   !    ro_ave =  ro_ave / max(dz_ave,epsi)
330   !    SnowOK =  max(zero,sign(unun,600.-ro_ave))
331   !    SnownH =  SnowOK + SnownH * (1. - SnowOK)
332
333  ! +--Integrated Snow/Ice Albedo: Case of Water on Bare Ice
334  ! +  -----------------------------------------------------
335
336      isn    =  max(min(isnoSV(ikl) ,ispiSV(ikl)),0)
337
338      albWIc =  aI1dSV-(aI1dSV-aI2dSV) &
339            *  exp(-(rusnSV(ikl) & !
340            *  (1. -SWS_SV(ikl) & ! 0 <=> freezing
341            *  (1  -min(1,iabs(isn-isnoSV(ikl))))) & ! 1 <=> isn=isnoSV
342            /   ru_dSV)**0.50)                        !
343      ! albWIc = max(aI1dSV,min(aI2dSV,albWIc+slopSV(ikl)*
344  !    .             min(5.,max(1.,dx/10000.))))
345
346      SignRo = sign(unun,ro_Ice-5.-ro__SV(ikl,isn))    ! RoSN<920kg/m3
347      SnowOK =  max(zero,SignRo)
348
349      albWIc = (1. - SnowOK) * albWIc + SnowOK &
350            * (aI2dSV + (aI3dSV -aI2dSV) &
351            * (ro__SV(ikl,isn)-ro_Ice)/(roSdSV-ro_Ice))
352
353  ! +  rocdSV < ro < ro_ice | aI2dSV< al <aI3dSV (fct of density)
354  ! +           ro > ro_ice | aI1dSV< al <aI2dSV (fct of superficial water content
355
356
357  ! +--Integrated Snow/Ice      Albedo
358  ! +  -------------------------------
359
360      a_SII1      =     albWIc      +(albSn1-albWIc)     *SnownH
361      a_SII1      = min(a_SII1       ,albSn1)
362
363      a_SII2      =     albWIc      +(albSn2-albWIc)     *SnownH
364      a_SII2      = min(a_SII2       ,albSn2)
365
366      a_SII3      =     albWIc      +(albSn3-albWIc)     *SnownH
367      a_SII3      = min(a_SII3       ,albSn3)
368
369      DO i=1,6
370      a_SII6(i)   = albWIc      +(albSn6(i)-albWIc)     *SnownH
371      a_SII6(i)   = min(a_SII6(i)       ,albSn6(i))
372      ENDDO
373
374  !c #AG     agesno =      min(agsnSV(ikl,isn)          ,AgeMax)
375  !c #AG     a_SII1      =     a_SII1      -0.175*agesno/AgeMax
376  ! +...                                   Impurities: Col de Porte Parameter.
377
378
379
380  ! +--Elsewhere    Integrated Snow/Ice Albedo
381  ! +  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
382  ! #cp     ELSE
383        albSII =     So1dSV*a_SII1 &
384              + So2dSV*a_SII2 &
385              + So3dSV*a_SII3
386  ! #cp     END IF
387
388
389  ! +--Integrated Snow/Ice/Soil Albedo
390  ! +  -------------------------------
391
392        alb1sv(ikl) =     albssv(ikl) +(a_SII1-albssv(ikl))*SIcenH
393        alb1sv(ikl) = min(alb1sv(ikl)  ,a_SII1)
394
395        alb2sv(ikl) =     albssv(ikl) +(a_SII2-albssv(ikl))*SIcenH
396        alb2sv(ikl) = min(alb2sv(ikl)  ,a_SII2)
397
398        alb3sv(ikl) =     albssv(ikl) +(a_SII3-albssv(ikl))*SIcenH
399        alb3sv(ikl) = min(alb3sv(ikl)  ,a_SII3)
400
401        albisv(ikl) =     albssv(ikl) +(albSII-albssv(ikl))*SIcenH
402        albisv(ikl) = min(albisv(ikl)  ,albSII)
403
404        DO i=1,6
405        alb6sv(ikl,i) = albssv(ikl) +(a_SII6(i)-albssv(ikl))*SIcenH
406        alb6sv(ikl,i) = min(alb6sv(ikl,i)  ,a_SII6(i))
407        ENDDO
408
409
410  ! +--Integrated Snow/Ice/Soil Albedo: Clouds Correction! Greuell & all., 1994
411  ! +  --------------------------------------------------! Glob.&t Planet.Change
412                                                   ! ! (9):91-114
413        alb1sv(ikl) = alb1sv(ikl) + 0.05 *(cld_SV(ikl)-0.5)*SIcenH &
414              + dalbed      *    (1.-cld_SV(ikl))
415        alb2sv(ikl) = alb2sv(ikl) + 0.05 *(cld_SV(ikl)-0.5)*SIcenH &
416              + dalbed      *    (1.-cld_SV(ikl))
417        alb3sv(ikl) = alb3sv(ikl) + 0.05 *(cld_SV(ikl)-0.5)*SIcenH &
418              + dalbed      *    (1.-cld_SV(ikl))
419        alb6sv(ikl,:) =alb6sv(ikl,:)+0.05 *(cld_SV(ikl)-0.5)*SIcenH &
420              + dalbed      *    (1.-cld_SV(ikl))
421        albisv(ikl) = albisv(ikl) + 0.05 *(cld_SV(ikl)-0.5)*SIcenH &
422              + dalbed      *    (1.-cld_SV(ikl))
423
424  ! +--Integrated Snow/Ice/Soil Albedo: Minimum snow albedo = aI1dSV
425  ! +  -------------------------------------------------------------
426
427        albedo_old  = albisv(ikl)
428        albisv(ikl) = max(albisv(ikl),aI1dSV   * SIcenH &
429              + albssv(ikl) *(1.0        - SIcenH))
430        alb1sv(ikl) = alb1sv(ikl) - 1.0/3.0 & ! 33 %
431              * (albedo_old-albisv(ikl)) / So1dSV
432        alb2sv(ikl) = alb2sv(ikl) - 1.0/3.0 & ! 33 %
433              * (albedo_old-albisv(ikl)) / So2dSV
434        alb3sv(ikl) = alb3sv(ikl) - 1.0/3.0 & ! 33 %
435              * (albedo_old-albisv(ikl)) / So3dSV
436        alb6sv(ikl,1:3) = alb6sv(ikl,1:3) - 1.0/6.0 & ! 16 %
437              * (albedo_old-albisv(ikl)) / (So1dSV/3)
438        alb6sv(ikl,4:6) = alb6sv(ikl,4:6) - 1.0/6.0 & ! 16 %
439              * (albedo_old-albisv(ikl)) / (So2dSV/3)
440
441
442  ! +--Integrated Snow/Ice/Soil Albedo: Maximum albedo = 95%
443  ! +  -----------------------------------------------------
444
445        albedo_old  = albisv(ikl)
446        albisv(ikl) = min(albisv(ikl),0.95)
447        alb1sv(ikl) = alb1sv(ikl) - 1.0/3.0 & ! 33 %
448              * (albedo_old-albisv(ikl)) / So1dSV
449        alb2sv(ikl) = alb2sv(ikl) - 1.0/3.0 & ! 33 %
450              * (albedo_old-albisv(ikl)) / So2dSV
451        alb3sv(ikl) = alb3sv(ikl) - 1.0/3.0 & ! 33 %
452              * (albedo_old-albisv(ikl)) / So3dSV
453        alb6sv(ikl,1:3) = alb6sv(ikl,1:3) - 1.0/6.0 & ! 16 %
454              * (albedo_old-albisv(ikl)) / (So1dSV/3)
455        alb6sv(ikl,4:6) = alb6sv(ikl,4:6) - 1.0/6.0 & ! 16 %
456              * (albedo_old-albisv(ikl)) / (So2dSV/3)
457
458
459    !Sea Ice/snow permanent-interractive prescription from Nemo
460    !AO_CK 20/02/2020
461
462    ! ! No check if coupling update since MAR and NEMO albedo are too different
463    ! and since MAR albedo is computed on properties that are not in NEMO
464    ! ! prescription for each time step with NEMO values
465
466  ! #AO      if (LSmask(ikl) .eq. 0 .and. coupling_ao .eq. .true.) then
467  ! #AO       if (AOmask(ikl) .eq. 0) then
468  ! #AO       albisv(ikl) =  (1.-AOmask(ikl))* albAOsisv(ikl)
469  ! #AO.                    +(AOmask(ikl)*albisv(ikl))
470  ! #AO       alb1sv(ikl) =   (1.-AOmask(ikl))* albAOsisv(ikl)
471  ! #AO.                    +(AOmask(ikl)*alb1sv(ikl))
472  ! #AO       alb2sv(ikl) =   (1.-AOmask(ikl))* albAOsisv(ikl)
473  ! #AO.                    +(AOmask(ikl)*alb2sv(ikl))
474  ! #AO       alb3sv(ikl) =   (1.-AOmask(ikl))* albAOsisv(ikl)
475  ! #AO.                    +(AOmask(ikl)*alb3sv(ikl))
476  ! #AO       endif
477  ! #AO      endif
478
479
480        alb1sv(ikl) = min(max(zero,alb1sv(ikl)),albmax)
481        alb2sv(ikl) = min(max(zero,alb2sv(ikl)),albmax)
482        alb3sv(ikl) = min(max(zero,alb3sv(ikl)),albmax)
483
484        DO i=1,6
485            alb6sv(ikl,i) = min(max(zero,alb6sv(ikl,i)),albmax)
486        ENDDO
487    END DO
488
489
490  ! +--Extinction Coefficient: Exponential Factor
491  ! +  ==========================================
492
493    DO ikl=1,knonv
494      sExt_1(ikl)        = 1.
495      sExt_2(ikl)        = 1.
496      sExt_3(ikl)        = 1.
497      sEX_sv(ikl,nsno+1) = 1.
498
499      coalb1(ikl) = (1.          -alb1sv(ikl))*So1dSV
500      coalb2(ikl) = (1.          -alb2sv(ikl))*So2dSV
501      coalb3(ikl) = (1.          -alb3sv(ikl))*So3dSV
502      coalbm      =  coalb1(ikl) +coalb2(ikl) +coalb3(ikl)
503      coalb1(ikl) =  coalb1(ikl)              /coalbm
504      coalb2(ikl) =  coalb2(ikl)              /coalbm
505      coalb3(ikl) =  coalb3(ikl)              /coalbm
506    END DO
507
508  !XF
509
510    DO   isn=nsno,1,-1
511      DO ikl=1,knonv
512        sEX_sv(ikl,isn) = 1.0
513       ! !sEX_sv(ikl,isn) = 0.95 ! if MAR is too warm in summer
514      END DO
515    END DO
516
517    DO ikl=1,knonv
518     DO isn=max(1,isnoSV(ikl)),1,-1
519
520      SignRo = sign(unun, rocdSV - ro__SV(ikl,isn))
521      SnowOK =  max(zero,SignRo) ! Ice Density Threshold
522
523      RoFrez =  1.e-3      * ro__SV(ikl,isn) * (1.0-eta_SV(ikl,isn))
524
525      OpSqrt = sqrt(max(epsi,SnOpSV(ikl,isn)))
526      exarg1 =      SnowOK  *1.e2 *max(sbeta1*RoFrez/OpSqrt,sbeta2) &
527            +(1.0-SnowOK)           *sbeta5
528      exarg2 =      SnowOK  *1.e2 *max(sbeta3*RoFrez/OpSqrt,sbeta4) &
529            +(1.0-SnowOK)           *sbeta5
530      exarg3 =      SnowOK  *1.e2     *sbeta5 &
531            +(1.0-SnowOK)           *sbeta5
532
533
534  ! +--Integrated Extinction of Solar Irradiance (Normalized Value)
535  ! +  ============================================================
536
537      sExt_1(ikl) = sExt_1(ikl) &
538            * exp(min(0.0,-exarg1 *dzsnSV(ikl,isn)))
539      sign_0      =              sign(unun,eps9   -sExt_1(ikl))
540      sExt_0      =               max(zero,sign_0)*sExt_1(ikl)
541      sExt_1(ikl) = sExt_1(ikl)                   -sExt_0
542
543      sExt_2(ikl) = sExt_2(ikl) &
544            * exp(min(0.0,-exarg2 *dzsnSV(ikl,isn)))
545      sign_0      =              sign(unun,eps9   -sExt_2(ikl))
546      sExt_0      =               max(zero,sign_0)*sExt_2(ikl)
547      sExt_2(ikl) = sExt_2(ikl)                   -sExt_0
548
549      sExt_3(ikl) = sExt_3(ikl) &
550            * exp(min(0.0,-exarg3 *dzsnSV(ikl,isn)))
551      sign_0      =              sign(unun,eps9   -sExt_3(ikl))
552      sExt_0      =               max(zero,sign_0)*sExt_3(ikl)
553      sExt_3(ikl) = sExt_3(ikl)                   -sExt_0
554
555      sEX_sv(ikl,isn) = coalb1(ikl) * sExt_1(ikl) &
556            + coalb2(ikl) * sExt_2(ikl) &
557            + coalb3(ikl) * sExt_3(ikl)
558    END DO
559  END DO
560
561  DO   isn=0,-nsol,-1
562    DO ikl=1,knonv
563      sEX_sv(ikl,isn) = 0.0
564    END DO
565  END DO
566
567
568  return
569
570
571end subroutine snoptp
572
573
574!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
575SUBROUTINE albedo_kokhanovsky(lambdamin,lambdamax, &
576        cossza,dopt,albint)
577  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
578  ! Authors: Hajar El Habchi El Fenniri, Etienne Vignon, Cecile Agosta
579       ! Ghislain Picard
580  ! Routine that calculates the  integrated snow spectral albedo between
581  ! lambdamin and lambdamax following Kokhanisky and Zege 2004,
582  ! Scattering optics of snow, Applied Optics, Vol 43, No7
583  ! Code inspired from the snowoptics package of Ghislain Picard:
584  ! https://github.com/ghislainp/snowoptics
585  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
586
587
588  USE VARphy
589
590  IMPLICIT NONE
591
592  ! Inputs
593  !--------
594  REAL :: lambdamin   ! minimum wavelength for integration [m]
595  REAL :: lambdamax   ! maximum wavelength for integration [m]
596  REAL :: cossza     ! solar zenith angle cosinus
597  REAL :: dopt        ! optical diameter [m]
598
599  !Outputs
600  !-------
601  REAL :: albint
602
603  ! Local Variables
604  !-----------------
605
606  REAL :: ropt,cosalb,norm,Pas
607  REAL :: SSA,alpha,gamm,R,cos30,alb30
608  INTEGER :: i
609
610
611  REAL :: B_amp                       ! amplification factor
612  PARAMETER (B_amp=1.6)
613
614  REAL :: g_asy                       ! asymetry factor
615  PARAMETER (g_asy=0.845)
616
617  INTEGER :: nlambda                  ! length of wavelength vector
618  PARAMETER(nlambda=200)
619
620  REAL :: lmin
621  PARAMETER(lmin=185.0E-9)
622
623  REAL :: lmax
624  PARAMETER(lmax=4000.0E-9)
625
626  REAL :: albmax
627  PARAMETER(albmax=0.99)
628
629  REAL :: wavelengths(nlambda)
630  REAL :: ni(nlambda)
631
632  DATA wavelengths / 1.85000000e-07, 2.04170854e-07, &
633        2.23341709e-07, 2.42512563e-07, &
634        2.61683417e-07, 2.80854271e-07, 3.00025126e-07, 3.19195980e-07, &
635        3.38366834e-07, 3.57537688e-07, 3.76708543e-07, 3.95879397e-07, &
636        4.15050251e-07, 4.34221106e-07, 4.53391960e-07, 4.72562814e-07, &
637        4.91733668e-07, 5.10904523e-07, 5.30075377e-07, 5.49246231e-07, &
638        5.68417085e-07, 5.87587940e-07, 6.06758794e-07, 6.25929648e-07, &
639        6.45100503e-07, 6.64271357e-07, 6.83442211e-07, 7.02613065e-07, &
640        7.21783920e-07, 7.40954774e-07, 7.60125628e-07, 7.79296482e-07, &
641        7.98467337e-07, 8.17638191e-07, 8.36809045e-07, 8.55979899e-07, &
642        8.75150754e-07, 8.94321608e-07, 9.13492462e-07, 9.32663317e-07, &
643        9.51834171e-07, 9.71005025e-07, 9.90175879e-07, 1.00934673e-06, &
644        1.02851759e-06, 1.04768844e-06, 1.06685930e-06, 1.08603015e-06, &
645        1.10520101e-06, 1.12437186e-06, 1.14354271e-06, 1.16271357e-06, &
646        1.18188442e-06, 1.20105528e-06, 1.22022613e-06, 1.23939698e-06, &
647        1.25856784e-06, 1.27773869e-06, 1.29690955e-06, 1.31608040e-06, &
648        1.33525126e-06, 1.35442211e-06, 1.37359296e-06, 1.39276382e-06, &
649        1.41193467e-06, 1.43110553e-06, 1.45027638e-06, 1.46944724e-06, &
650        1.48861809e-06, 1.50778894e-06, 1.52695980e-06, 1.54613065e-06, &
651        1.56530151e-06, 1.58447236e-06, 1.60364322e-06, 1.62281407e-06, &
652        1.64198492e-06, 1.66115578e-06, 1.68032663e-06, 1.69949749e-06, &
653        1.71866834e-06, 1.73783920e-06, 1.75701005e-06, 1.77618090e-06, &
654        1.79535176e-06, 1.81452261e-06, 1.83369347e-06, 1.85286432e-06, &
655        1.87203518e-06, 1.89120603e-06, 1.91037688e-06, 1.92954774e-06, &
656        1.94871859e-06, 1.96788945e-06, 1.98706030e-06, 2.00623116e-06, &
657        2.02540201e-06, 2.04457286e-06, 2.06374372e-06, 2.08291457e-06, &
658        2.10208543e-06, 2.12125628e-06, 2.14042714e-06, 2.15959799e-06, &
659        2.17876884e-06, 2.19793970e-06, 2.21711055e-06, 2.23628141e-06, &
660        2.25545226e-06, 2.27462312e-06, 2.29379397e-06, 2.31296482e-06, &
661        2.33213568e-06, 2.35130653e-06, 2.37047739e-06, 2.38964824e-06, &
662        2.40881910e-06, 2.42798995e-06, 2.44716080e-06, 2.46633166e-06, &
663        2.48550251e-06, 2.50467337e-06, 2.52384422e-06, 2.54301508e-06, &
664        2.56218593e-06, 2.58135678e-06, 2.60052764e-06, 2.61969849e-06, &
665        2.63886935e-06, 2.65804020e-06, 2.67721106e-06, 2.69638191e-06, &
666        2.71555276e-06, 2.73472362e-06, 2.75389447e-06, 2.77306533e-06, &
667        2.79223618e-06, 2.81140704e-06, 2.83057789e-06, 2.84974874e-06, &
668        2.86891960e-06, 2.88809045e-06, 2.90726131e-06, 2.92643216e-06, &
669        2.94560302e-06, 2.96477387e-06, 2.98394472e-06, 3.00311558e-06, &
670        3.02228643e-06, 3.04145729e-06, 3.06062814e-06, 3.07979899e-06, &
671        3.09896985e-06, 3.11814070e-06, 3.13731156e-06, 3.15648241e-06, &
672        3.17565327e-06, 3.19482412e-06, 3.21399497e-06, 3.23316583e-06, &
673        3.25233668e-06, 3.27150754e-06, 3.29067839e-06, 3.30984925e-06, &
674        3.32902010e-06, 3.34819095e-06, 3.36736181e-06, 3.38653266e-06, &
675        3.40570352e-06, 3.42487437e-06, 3.44404523e-06, 3.46321608e-06, &
676        3.48238693e-06, 3.50155779e-06, 3.52072864e-06, 3.53989950e-06, &
677        3.55907035e-06, 3.57824121e-06, 3.59741206e-06, 3.61658291e-06, &
678        3.63575377e-06, 3.65492462e-06, 3.67409548e-06, 3.69326633e-06, &
679        3.71243719e-06, 3.73160804e-06, 3.75077889e-06, 3.76994975e-06, &
680        3.78912060e-06, 3.80829146e-06, 3.82746231e-06, 3.84663317e-06, &
681        3.86580402e-06, 3.88497487e-06, 3.90414573e-06, 3.92331658e-06, &
682        3.94248744e-06, 3.96165829e-06, 3.98082915e-06, 4.00000000e-06/
683
684
685  DATA ni /7.74508407e-10, 7.74508407e-10, &
686        7.74508407e-10, 7.74508407e-10, &
687        7.74508407e-10, 7.74508407e-10, 7.74508407e-10, 7.74508407e-10, &
688        6.98381122e-10, 6.23170274e-10, 5.97655992e-10, 5.84106004e-10, &
689        5.44327597e-10, 5.71923510e-10, 6.59723827e-10, 8.05183870e-10, &
690        1.03110161e-09, 1.36680386e-09, 1.85161253e-09, 2.56487751e-09, &
691        3.56462855e-09, 4.89450926e-09, 6.49252022e-09, 9.62029335e-09, &
692        1.32335015e-08, 1.75502184e-08, 2.19240625e-08, 3.03304156e-08, &
693        4.07715972e-08, 5.00414911e-08, 7.09722331e-08, 1.00773751e-07, &
694        1.31427409e-07, 1.42289041e-07, 1.49066787e-07, 2.01558515e-07, &
695        2.99106105e-07, 4.03902086e-07, 4.54292169e-07, 5.21906983e-07, &
696        6.27643362e-07, 9.43955678e-07, 1.33464494e-06, 1.97278315e-06, &
697        2.31801329e-06, 2.20584676e-06, 1.85568138e-06, 1.73395193e-06, &
698        1.73101406e-06, 1.91333936e-06, 2.26413019e-06, 3.23959718e-06, &
699        4.94316963e-06, 6.89378896e-06, 1.02237444e-05, 1.21439656e-05, &
700        1.31567585e-05, 1.33448288e-05, 1.32000000e-05, 1.31608040e-05, &
701        1.33048369e-05, 1.40321464e-05, 1.51526244e-05, 1.80342858e-05, &
702        3.82875736e-05, 1.07325259e-04, 2.11961637e-04, 3.82008054e-04, &
703        5.30897470e-04, 5.29244735e-04, 4.90876605e-04, 4.33905427e-04, &
704        3.77795349e-04, 3.17633099e-04, 2.81078564e-04, 2.57579485e-04, &
705        2.42203100e-04, 2.23789060e-04, 2.04306870e-04, 1.87909255e-04, &
706        1.73117146e-04, 1.61533186e-04, 1.53420328e-04, 1.47578033e-04, &
707        1.42334776e-04, 1.35691466e-04, 1.30495414e-04, 1.36065123e-04, &
708        1.70928821e-04, 2.66389730e-04, 4.80957955e-04, 8.25041961e-04, &
709        1.21654792e-03, 1.50232875e-03, 1.62316078e-03, 1.61649750e-03, &
710        1.53736801e-03, 1.42343711e-03, 1.24459117e-03, 1.02388611e-03, &
711        7.89112523e-04, 5.97204264e-04, 4.57152413e-04, 3.62341259e-04, &
712        2.99128332e-04, 2.57035569e-04, 2.26992203e-04, 2.07110355e-04, &
713        2.05835688e-04, 2.25108810e-04, 2.64262893e-04, 3.23594011e-04, &
714        3.93061117e-04, 4.62789970e-04, 5.19664416e-04, 5.59739628e-04, &
715        5.93476084e-04, 6.22797885e-04, 6.57484833e-04, 6.92849600e-04, &
716        7.26584901e-04, 7.56604648e-04, 7.68009488e-04, 7.65579073e-04, &
717        7.50526164e-04, 7.39809972e-04, 7.55622847e-04, 8.05099514e-04, &
718        9.67279246e-04, 1.16281559e-03, 1.42570247e-03, 2.04986585e-03, &
719        2.93971170e-03, 4.49827711e-03, 7.32537532e-03, 1.18889734e-02, &
720        1.85851805e-02, 2.86242532e-02, 4.34131035e-02, 6.37828307e-02, &
721        9.24145850e-02, 1.35547945e-01, 1.94143245e-01, 2.54542814e-01, &
722        3.02282024e-01, 3.42214181e-01, 3.85475065e-01, 4.38000000e-01, &
723        4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, &
724        4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, &
725        4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, &
726        4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, &
727        4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, &
728        4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, &
729        4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, &
730        4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, &
731        4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, &
732        4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, &
733        4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, &
734        4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, &
735        4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01/
736
737
738  Pas     =(lmax-lmin)/nlambda
739  ropt    = dopt/2.0
740  SSA     = 3.0/(rhoIce*ropt)
741  cos30   = cos(30.0/180.0*pi)
742
743
744  albint=0.
745  norm=0.
746
747  DO i = 1,nlambda
748      gamm = 4.0 * pi * ni(i) / wavelengths(i)
749      cosalb = 2.0 / (SSA * rhoice) * B_amp * gamm
750      alpha = 16. / 3 * cosalb / (1.0 - g_asy)
751      R = exp(-(alpha**0.5) * 3.0 / 7.0 * (1.0 + 2.0 * cossza))
752      alb30 = exp(-(alpha**0.5)* 3.0 / 7.0 * (1.0 + 20 * cos30))
753
754      IF ((wavelengths(i).GE.lambdamin).AND. &
755            (wavelengths(i).LT.lambdamax)) THEN
756      albint=albint+R*Pas  ! rectangle integration -> can be improved with trapezintegration
757      norm=norm+Pas
758      ENDIF
759
760  END DO
761
762  albint=max(0.,min(albint/max(norm,1E-10),albmax))
763
764END SUBROUTINE albedo_kokhanovsky
765
Note: See TracBrowser for help on using the repository browser.