source: LMDZ6/branches/blowing_snow/libf/phylmd/surf_landice_mod.F90 @ 4774

Last change on this file since 4774 was 4485, checked in by evignon, 18 months ago

premier commit pour l'ajout de la neige soufflee sur la nouvelle branche
  • 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
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 18.7 KB
Line 
1!
2MODULE surf_landice_mod
3 
4  IMPLICIT NONE
5
6CONTAINS
7!
8!****************************************************************************************
9!
10  SUBROUTINE surf_landice(itime, dtime, knon, knindex, &
11       rlon, rlat, debut, lafin, &
12       rmu0, lwdownm, albedo, pphi1, &
13       swnet, lwnet, tsurf, p1lay, &
14       cdragh, cdragm, precip_rain, precip_snow, precip_bs, temp_air, spechum, &
15       AcoefH, AcoefQ, BcoefH, BcoefQ, &
16       AcoefU, AcoefV, BcoefU, BcoefV, &
17       ps, u1, v1, gustiness, rugoro, pctsrf, &
18       snow, qsurf, qsol, qbs1, agesno, &
19       tsoil, z0m, z0h, SFRWL, alb_dir, alb_dif, evap, fluxsens, fluxlat, fluxbs, &
20       tsurf_new, dflux_s, dflux_l, &
21       alt, slope, cloudf, &
22       snowhgt, qsnow, to_ice, sissnow, &
23       alb3, runoff, &
24       flux_u1, flux_v1)
25
26    USE dimphy
27    USE geometry_mod,     ONLY : longitude,latitude
28    USE surface_data,     ONLY : type_ocean, calice, calsno, landice_opt, iflag_albcalc
29    USE fonte_neige_mod,  ONLY : fonte_neige,run_off_lic,fqcalving_global,ffonte_global,fqfonte_global,runofflic_global
30    USE cpl_mod,          ONLY : cpl_send_landice_fields
31    USE calcul_fluxs_mod
32    USE phys_local_var_mod, ONLY : zxrhoslic, zxustartlic
33    USE phys_output_var_mod, ONLY : snow_o,zfra_o
34!FC
35    USE ioipsl_getin_p_mod, ONLY : getin_p
36    USE blowing_snow_ini_mod, ONLY : zeta_bs, pbst_bs, prt_bs
37
38#ifdef CPP_INLANDSIS
39    USE surf_inlandsis_mod,  ONLY : surf_inlandsis
40#endif
41
42    USE indice_sol_mod
43
44!    INCLUDE "indicesol.h"
45    INCLUDE "dimsoil.h"
46    INCLUDE "YOMCST.h"
47    INCLUDE "clesphys.h"
48
49! Input variables
50!****************************************************************************************
51    INTEGER, INTENT(IN)                           :: itime, knon
52    INTEGER, DIMENSION(klon), INTENT(in)          :: knindex
53    REAL, INTENT(in)                              :: dtime
54    REAL, DIMENSION(klon), INTENT(IN)             :: swnet ! net shortwave radiance
55    REAL, DIMENSION(klon), INTENT(IN)             :: lwnet ! net longwave radiance
56    REAL, DIMENSION(klon), INTENT(IN)             :: tsurf
57    REAL, DIMENSION(klon), INTENT(IN)             :: p1lay
58    REAL, DIMENSION(klon), INTENT(IN)             :: cdragh, cdragm
59    REAL, DIMENSION(klon), INTENT(IN)             :: precip_rain, precip_snow, precip_bs
60    REAL, DIMENSION(klon), INTENT(IN)             :: temp_air, spechum
61    REAL, DIMENSION(klon), INTENT(IN)             :: AcoefH, AcoefQ
62    REAL, DIMENSION(klon), INTENT(IN)             :: BcoefH, BcoefQ
63    REAL, DIMENSION(klon), INTENT(IN)             :: AcoefU, AcoefV, BcoefU, BcoefV
64    REAL, DIMENSION(klon), INTENT(IN)             :: ps
65    REAL, DIMENSION(klon), INTENT(IN)             :: u1, v1, gustiness, qbs1
66    REAL, DIMENSION(klon), INTENT(IN)             :: rugoro
67    REAL, DIMENSION(klon,nbsrf), INTENT(IN)       :: pctsrf
68
69    LOGICAL,  INTENT(IN)                          :: debut   !true if first step
70    LOGICAL,  INTENT(IN)                          :: lafin   !true if last step
71    REAL, DIMENSION(klon), INTENT(IN)             :: rlon, rlat
72    REAL, DIMENSION(klon), INTENT(IN)             :: rmu0
73    REAL, DIMENSION(klon), INTENT(IN)             :: lwdownm !ylwdown
74    REAL, DIMENSION(klon), INTENT(IN)             :: albedo  !mean albedo
75    REAL, DIMENSION(klon), INTENT(IN)             :: pphi1   
76    REAL, DIMENSION(klon), INTENT(IN)             :: alt   !mean altitude of the grid box 
77    REAL, DIMENSION(klon), INTENT(IN)             :: slope   !mean slope in grid box 
78    REAL, DIMENSION(klon), INTENT(IN)             :: cloudf  !total cloud fraction
79
80! In/Output variables
81!****************************************************************************************
82    REAL, DIMENSION(klon), INTENT(INOUT)          :: snow, qsol
83    REAL, DIMENSION(klon), INTENT(INOUT)          :: agesno
84    REAL, DIMENSION(klon, nsoilmx), INTENT(INOUT) :: tsoil
85
86! Output variables
87!****************************************************************************************
88    REAL, DIMENSION(klon), INTENT(OUT)            :: qsurf
89    REAL, DIMENSION(klon), INTENT(OUT)            :: z0m, z0h
90!albedo SB >>>
91!    REAL, DIMENSION(klon), INTENT(OUT)            :: alb1  ! new albedo in visible SW interval
92!    REAL, DIMENSION(klon), INTENT(OUT)            :: alb2  ! new albedo in near IR interval
93    REAL, DIMENSION(6), INTENT(IN)                :: SFRWL
94    REAL, DIMENSION(klon,nsw), INTENT(OUT)        :: alb_dir,alb_dif
95!albedo SB <<<
96    REAL, DIMENSION(klon), INTENT(OUT)            :: evap, fluxsens, fluxlat
97    REAL, DIMENSION(klon), INTENT(OUT)            :: fluxbs
98    REAL, DIMENSION(klon), INTENT(OUT)            :: tsurf_new
99    REAL, DIMENSION(klon), INTENT(OUT)            :: dflux_s, dflux_l     
100    REAL, DIMENSION(klon), INTENT(OUT)            :: flux_u1, flux_v1
101
102    REAL, DIMENSION(klon), INTENT(OUT)           :: alb3
103    REAL, DIMENSION(klon), INTENT(OUT)           :: qsnow   !column water in snow [kg/m2]
104    REAL, DIMENSION(klon), INTENT(OUT)           :: snowhgt !Snow height (m)
105    REAL, DIMENSION(klon), INTENT(OUT)           :: to_ice
106    REAL, DIMENSION(klon), INTENT(OUT)           :: sissnow
107    REAL, DIMENSION(klon), INTENT(OUT)           :: runoff  !Land ice runoff
108 
109
110! Local variables
111!****************************************************************************************
112    REAL, DIMENSION(klon)    :: soilcap, soilflux
113    REAL, DIMENSION(klon)    :: cal, beta, dif_grnd
114    REAL, DIMENSION(klon)    :: zfra, alb_neig
115    REAL, DIMENSION(klon)    :: radsol
116    REAL, DIMENSION(klon)    :: u0, v0, u1_lay, v1_lay, ustar
117    INTEGER                  :: i,j,nt
118    REAL, DIMENSION(klon)    :: fqfonte,ffonte
119    REAL, DIMENSION(klon)    :: run_off_lic_frac
120    REAL, DIMENSION(klon)    :: emis_new                  !Emissivity
121    REAL, DIMENSION(klon)    :: swdown,lwdown
122    REAL, DIMENSION(klon)    :: precip_snow_adv, snow_adv !Snow Drift precip./advection (not used in inlandsis)
123    REAL, DIMENSION(klon)    :: erod                      !erosion of surface snow (flux, kg/m2/s like evap)
124    REAL, DIMENSION(klon)    :: zsl_height, wind_velo     !surface layer height, wind spd
125    REAL, DIMENSION(klon)    :: dens_air,  snow_cont_air  !air density; snow content air
126    REAL, DIMENSION(klon)    :: alb_soil                  !albedo of underlying ice
127    REAL, DIMENSION(klon)    :: pexner                    !Exner potential
128    REAL                     :: pref
129    REAL, DIMENSION(klon,nsoilmx) :: tsoil0               !modif
130    REAL                          :: dtis                ! subtimestep
131    LOGICAL                       :: debut_is, lafin_is  ! debut and lafin for inlandsis
132
133    CHARACTER (len = 20)                      :: modname = 'surf_landice'
134    CHARACTER (len = 80)                      :: abort_message
135
136
137    REAL,DIMENSION(klon) :: alb1,alb2
138    REAL,DIMENSION(klon) :: precip_totsnow, evap_totsnow
139    REAL, DIMENSION (klon,6) :: alb6
140    REAL                   :: rho0, rhoice, ustart0, hsalt, esalt, qsalt
141    REAL                   :: tau_dens, tau_dens0, tau_densmin, rhomax, rhohard
142    REAL, DIMENSION(klon)  :: ws1, rhos, ustart
143! End definition
144!****************************************************************************************
145!FC
146!FC
147   REAL,SAVE :: alb_vis_sno_lic
148  !$OMP THREADPRIVATE(alb_vis_sno_lic)
149   REAL,SAVE :: alb_nir_sno_lic
150  !$OMP THREADPRIVATE(alb_nir_sno_lic)
151  LOGICAL, SAVE :: firstcall = .TRUE.
152  !$OMP THREADPRIVATE(firstcall)
153
154
155!FC firtscall initializations
156!******************************************************************************************
157  IF (firstcall) THEN
158  alb_vis_sno_lic=0.77
159  CALL getin_p('alb_vis_sno_lic',alb_vis_sno_lic)
160           PRINT*, 'alb_vis_sno_lic',alb_vis_sno_lic
161  alb_nir_sno_lic=0.77
162  CALL getin_p('alb_nir_sno_lic',alb_nir_sno_lic)
163           PRINT*, 'alb_nir_sno_lic',alb_nir_sno_lic
164 
165  firstcall=.false.
166  ENDIF
167!******************************************************************************************
168
169! Initialize output variables
170    alb3(:) = 999999.
171    alb2(:) = 999999.
172    alb1(:) = 999999.
173    fluxbs(:)=0. 
174    runoff(:) = 0.
175!****************************************************************************************
176! Calculate total absorbed radiance at surface
177!
178!****************************************************************************************
179    radsol(:) = 0.0
180    radsol(1:knon) = swnet(1:knon) + lwnet(1:knon)
181
182!****************************************************************************************
183
184!****************************************************************************************
185!  landice_opt = 0 : soil_model, calcul_flux, fonte_neige, ... 
186!  landice_opt = 1  : prepare and call INterace Lmdz SISvat (INLANDSIS)
187!****************************************************************************************
188
189
190    IF (landice_opt .EQ. 1) THEN
191
192!****************************************************************************************   
193! CALL to INLANDSIS interface
194!****************************************************************************************
195#ifdef CPP_INLANDSIS
196
197        debut_is=debut
198        lafin_is=.false.
199        ! Suppose zero surface speed
200        u0(:)            = 0.0
201        v0(:)            = 0.0
202
203
204        CALL calcul_flux_wind(knon, dtime, &
205         u0, v0, u1, v1, gustiness, cdragm, &
206         AcoefU, AcoefV, BcoefU, BcoefV, &
207         p1lay, temp_air, &
208         flux_u1, flux_v1)
209
210       
211       ! Set constants and compute some input for SISVAT
212       ! = 1000 hPa
213       ! and calculate incoming flux for SW and LW interval: swdown, lwdown
214       swdown(:)        = 0.0
215       lwdown(:)        = 0.0
216       snow_cont_air(:) = 0.  ! the snow content in air is not a prognostic variable of the model     
217       alb_soil(:)      = 0.4 ! before albedo(:) but here it is the ice albedo that we have to set
218       ustar(:)         = 0.
219       pref             = 100000.       
220       DO i = 1, knon
221          swdown(i)        = swnet(i)/(1-albedo(i))
222          lwdown(i)        = lwdownm(i)
223          wind_velo(i)     = u1(i)**2 + v1(i)**2
224          wind_velo(i)     = wind_velo(i)**0.5
225          pexner(i)        = (p1lay(i)/pref)**(RD/RCPD)
226          dens_air(i)      = p1lay(i)/RD/temp_air(i)  ! dry air density
227          zsl_height(i)    = pphi1(i)/RG     
228          tsoil0(i,:)      = tsoil(i,:) 
229          ustar(i)= (cdragm(i)*(wind_velo(i)**2))**0.5   
230       END DO
231       
232
233
234        dtis=dtime
235
236          IF (lafin) THEN
237            lafin_is=.true.
238          END IF
239
240          CALL surf_inlandsis(knon, rlon, rlat, knindex, itime, dtis, debut_is, lafin_is,&
241            rmu0, swdown, lwdown, albedo, pexner, ps, p1lay, precip_rain, precip_snow,   &
242            zsl_height, wind_velo, ustar, temp_air, dens_air, spechum, tsurf,&
243            rugoro, snow_cont_air, alb_soil, alt, slope, cloudf, &
244            radsol, qsol, tsoil0, snow, zfra, snowhgt, qsnow, to_ice, sissnow,agesno,   &
245            AcoefH, AcoefQ, BcoefH, BcoefQ, cdragm, cdragh, &
246            run_off_lic, fqfonte, ffonte, evap, erod, fluxsens, fluxlat,dflux_s, dflux_l, &
247            tsurf_new, alb1, alb2, alb3, alb6, &
248            emis_new, z0m, z0h, qsurf)
249
250          debut_is=.false.
251
252
253        ! Treatment of snow melting and calving
254
255        ! for consistency with standard LMDZ, add calving to run_off_lic
256        run_off_lic(:)=run_off_lic(:) + to_ice(:)
257
258        DO i = 1, knon
259           ffonte_global(knindex(i),is_lic)    = ffonte(i)
260           fqfonte_global(knindex(i),is_lic)   = fqfonte(i)! net melting= melting - refreezing
261           fqcalving_global(knindex(i),is_lic) = to_ice(i) ! flux
262           runofflic_global(knindex(i)) = run_off_lic(i)
263        ENDDO
264        ! Here, we assume that the calving term is equal to the to_ice term
265        ! (no ice accumulation)
266
267
268#else
269       abort_message='Pb de coherence: landice_opt = 1  mais CPP_INLANDSIS = .false.'
270       CALL abort_physic(modname,abort_message,1)
271#endif
272
273
274    ELSE
275
276!****************************************************************************************
277! Soil calculations
278!
279!****************************************************************************************
280
281    ! EV: use calbeta
282    CALL calbeta(dtime, is_lic, knon, snow, qsol, beta, cal, dif_grnd)
283
284
285    ! use soil model and recalculate properly cal
286    IF (soil_model) THEN
287       CALL soil(dtime, is_lic, knon, snow, tsurf, qsol, &
288        & longitude(knindex(1:knon)), latitude(knindex(1:knon)), tsoil, soilcap, soilflux)
289       cal(1:knon) = RCPD / soilcap(1:knon)
290       radsol(1:knon)  = radsol(1:knon) + soilflux(1:knon)
291    ELSE
292       cal = RCPD * calice
293       WHERE (snow > 0.0) cal = RCPD * calsno
294    ENDIF
295
296
297!****************************************************************************************
298! Calulate fluxes
299!
300!****************************************************************************************
301!    beta(:) = 1.0
302!    dif_grnd(:) = 0.0
303
304! Suppose zero surface speed
305    u0(:)=0.0
306    v0(:)=0.0
307    u1_lay(:) = u1(:) - u0(:)
308    v1_lay(:) = v1(:) - v0(:)
309
310    CALL calcul_fluxs(knon, is_lic, dtime, &
311         tsurf, p1lay, cal, beta, cdragh, cdragh, ps, &
312         precip_rain, precip_snow, snow, qsurf,  &
313         radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, gustiness, &
314         1.,AcoefH, AcoefQ, BcoefH, BcoefQ, &
315         tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
316
317    CALL calcul_flux_wind(knon, dtime, &
318         u0, v0, u1, v1, gustiness, cdragm, &
319         AcoefU, AcoefV, BcoefU, BcoefV, &
320         p1lay, temp_air, &
321         flux_u1, flux_v1)
322
323
324!****************************************************************************************
325! Calculate albedo
326!
327!****************************************************************************************
328 
329    CALL albsno(klon,knon,dtime,agesno(:),alb_neig(:), precip_snow(:)) 
330
331
332! EV: following lines are obsolete since we set alb1 and alb2 to constant values
333! I therefore comment them
334!    alb1(1:knon) = alb_neig(1:knon)*zfra(1:knon) + &
335!         0.6 * (1.0-zfra(1:knon))
336!
337!IM: plusieurs choix/tests sur l'albedo des "glaciers continentaux"
338!       alb1(1 : knon)  = 0.6 !IM cf FH/GK
339!       alb1(1 : knon)  = 0.82
340!       alb1(1 : knon)  = 0.77 !211003 Ksta0.77
341!       alb1(1 : knon)  = 0.8 !KstaTER0.8 & LMD_ARMIP5
342!IM: KstaTER0.77 & LMD_ARMIP6   
343
344! Attantion: alb1 and alb2 are not the same!
345    alb1(1:knon)  = alb_vis_sno_lic
346    alb2(1:knon)  = alb_nir_sno_lic
347
348
349!****************************************************************************************
350! Rugosity
351!
352!****************************************************************************************
353    z0m = z0m_landice
354    z0h = z0h_landice
355    !z0m = SQRT(z0m**2+rugoro**2)
356
357
358
359  ! Simple blowing snow param
360  if (ok_bs) then
361       ustart0 = 0.211
362       rhoice = 920.0
363       rho0 = 200.0
364       rhomax=450.0
365       rhohard=400.0
366       tau_dens0=86400.0*10.  ! 10 days by default, in s
367       tau_densmin=86400.0 ! 1 days according to in situ obs by C. Amory
368       do i = 1, knon
369           ! estimation of snow density
370           ! snow density increases with snow age and
371           ! increases even faster in case of sedimentation of blowing snow
372           tau_dens=max(tau_densmin, tau_dens0*exp(-abs(precip_bs(i))/pbst_bs-abs(precip_rain(i))/prt_bs))
373           rhos(i)=rho0+(rhohard-rho0)*(1.-exp(-agesno(i)*86400.0/tau_dens))
374           ! blowing snow flux formula used in MAR
375           ws1(i)=(u1(i)**2+v1(i)**2)**0.5
376           ustar(i)=(cdragm(i)*(u1(i)**2+v1(i)**2))**0.5
377           ustart(i)=ustart0*exp(max(rhoice/rho0-rhoice/rhos(i),0.))*exp(max(0.,rhos(i)-rhomax))
378           ! we have multiplied by exp to prevent erosion when rhos>rhomax (usefull till
379           ! rhohard<450)
380           esalt=1./(3.25*max(ustar(i),0.001))
381           hsalt=0.08436*ustar(i)**1.27
382           qsalt=(max(ustar(i)**2-ustart(i)**2,0.))/(RG*hsalt)*esalt
383           !ep=qsalt*cdragm(i)*sqrt(u1(i)**2+v1(i)**2)
384           fluxbs(i)= zeta_bs*p1lay(i)/RD/temp_air(i)*ws1(i)*cdragm(i)*(qbs1(i)-qsalt)
385       enddo
386
387       ! for outputs
388       do j = 1, knon
389          i = knindex(j)
390          zxustartlic(i) = ustart(j)
391          zxrhoslic(i) = rhos(j)
392       enddo
393
394  endif
395
396
397
398!****************************************************************************************
399! Calculate surface snow amount
400!   
401!****************************************************************************************
402    IF (ok_bs) THEN
403      precip_totsnow=precip_snow+precip_bs
404      evap_totsnow=evap-fluxbs ! flux bs is positive towards the surface (snow erosion)
405    ELSE
406      precip_totsnow=precip_snow
407      evap_totsnow=evap
408    ENDIF
409
410    CALL fonte_neige(knon, is_lic, knindex, dtime, &
411         tsurf, precip_rain, precip_totsnow,  &
412         snow, qsol, tsurf_new, evap_totsnow)
413
414    WHERE (snow(1 : knon) .LT. 0.0001) agesno(1 : knon) = 0.                                         
415    zfra(1:knon) = MAX(0.0,MIN(1.0,snow(1:knon)/(snow(1:knon)+10.0))) 
416
417
418    END IF ! landice_opt
419
420
421!****************************************************************************************
422! Send run-off on land-ice to coupler if coupled ocean.
423! run_off_lic has been calculated in fonte_neige or surf_inlandsis
424! If landice_opt>=2, corresponding call is done from surf_land_orchidee
425!****************************************************************************************
426    IF (type_ocean=='couple' .AND. landice_opt .LT. 2) THEN
427       ! Compress fraction where run_off_lic is active (here all pctsrf(is_lic))
428       run_off_lic_frac(:)=0.0
429       DO j = 1, knon
430          i = knindex(j)
431          run_off_lic_frac(j) = pctsrf(i,is_lic)
432       ENDDO
433
434       CALL cpl_send_landice_fields(itime, knon, knindex, run_off_lic, run_off_lic_frac)
435    ENDIF
436
437 ! transfer runoff rate [kg/m2/s](!) to physiq for output
438    runoff(1:knon)=run_off_lic(1:knon)/dtime
439
440       snow_o=0.
441       zfra_o = 0.
442       DO j = 1, knon
443           i = knindex(j)
444           snow_o(i) = snow(j)
445           zfra_o(i) = zfra(j)
446       ENDDO
447
448
449!albedo SB >>>
450     select case(NSW)
451     case(2)
452       alb_dir(1:knon,1)=alb1(1:knon)
453       alb_dir(1:knon,2)=alb2(1:knon)
454     case(4)
455       alb_dir(1:knon,1)=alb1(1:knon)
456       alb_dir(1:knon,2)=alb2(1:knon)
457       alb_dir(1:knon,3)=alb2(1:knon)
458       alb_dir(1:knon,4)=alb2(1:knon)
459     case(6)
460       alb_dir(1:knon,1)=alb1(1:knon)
461       alb_dir(1:knon,2)=alb1(1:knon)
462       alb_dir(1:knon,3)=alb1(1:knon)
463       alb_dir(1:knon,4)=alb2(1:knon)
464       alb_dir(1:knon,5)=alb2(1:knon)
465       alb_dir(1:knon,6)=alb2(1:knon)
466
467       IF ((landice_opt .EQ. 1) .AND. (iflag_albcalc .EQ. 2)) THEN
468       alb_dir(1:knon,1)=alb6(1:knon,1)
469       alb_dir(1:knon,2)=alb6(1:knon,2)
470       alb_dir(1:knon,3)=alb6(1:knon,3)
471       alb_dir(1:knon,4)=alb6(1:knon,4)
472       alb_dir(1:knon,5)=alb6(1:knon,5)
473       alb_dir(1:knon,6)=alb6(1:knon,6)
474       ENDIF
475
476     end select
477alb_dif=alb_dir
478!albedo SB <<<
479
480
481  END SUBROUTINE surf_landice
482!
483!****************************************************************************************
484!
485END MODULE surf_landice_mod
486
487
488
Note: See TracBrowser for help on using the repository browser.