source: LMDZ6/trunk/libf/phylmd/surf_landice_mod.F90 @ 5503

Last change on this file since 5503 was 5486, checked in by evignon, 6 days ago

inclusion d'un diagnostique de la sublimation de la glace sur les landice
pour la conservation de l'eau
  • 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: 27.9 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       AcoefQBS, BcoefQBS, &
18       ps, u1, v1, gustiness, rugoro, pctsrf, &
19       snow, qsurf, qsol, qbs1, agesno, &
20       tsoil, z0m, z0h, SFRWL, alb_dir, alb_dif, evap, icesub_lic, fluxsens, fluxlat, fluxbs, &
21       tsurf_new, dflux_s, dflux_l, &
22       alt, slope, cloudf, &
23       snowhgt, qsnow, to_ice, sissnow, &
24       alb3, runoff, &
25       flux_u1, flux_v1 &
26#ifdef ISO
27         &      ,xtprecip_rain, xtprecip_snow,xtspechum,Rland_ice &
28         &      ,xtsnow,xtsol,xtevap &
29#endif               
30           &    )
31
32    USE dimphy
33    USE geometry_mod,     ONLY : longitude,latitude
34    USE surface_data,     ONLY : type_ocean, calice, calsno, landice_opt, iflag_albcalc
35    USE fonte_neige_mod,  ONLY : fonte_neige,run_off_lic,fqcalving_global,ffonte_global,fqfonte_global,runofflic_global
36    USE cpl_mod,          ONLY : cpl_send_landice_fields
37    USE calcul_fluxs_mod
38    USE phys_local_var_mod, ONLY : zxrhoslic, zxustartlic, zxqsaltlic, tempsmoothlic
39    USE phys_output_var_mod, ONLY : snow_o,zfra_o
40#ifdef ISO   
41    USE fonte_neige_mod,  ONLY : xtrun_off_lic
42    USE infotrac_phy,     ONLY : ntiso,niso
43    USE isotopes_routines_mod, ONLY: calcul_iso_surf_lic_vectall
44#ifdef ISOVERIF
45    USE isotopes_mod, ONLY: iso_eau,ridicule
46    USE isotopes_verif_mod
47#endif
48#endif
49 
50    USE clesphys_mod_h
51    USE yomcst_mod_h
52    USE ioipsl_getin_p_mod, ONLY : getin_p
53    USE lmdz_blowing_snow_ini, ONLY : c_esalt_bs, zeta_bs, pbst_bs, prt_bs, rhoice_bs, rhohard_bs
54    USE lmdz_blowing_snow_ini, ONLY : rhofresh_bs, tau_eqsalt_bs, tau_dens0_bs, tau_densmin_bs
55    USE surf_inlandsis_mod,  ONLY : surf_inlandsis
56
57    USE indice_sol_mod
58    USE lmdz_cppkeys_wrapper, ONLY: CPPKEY_INLANDSIS
59    USE dimsoil_mod_h, ONLY: nsoilmx
60
61
62
63! Input variables
64!****************************************************************************************
65    INTEGER, INTENT(IN)                           :: itime, knon
66    INTEGER, DIMENSION(klon), INTENT(in)          :: knindex
67    REAL, INTENT(in)                              :: dtime
68    REAL, DIMENSION(klon), INTENT(IN)             :: swnet ! net shortwave radiance
69    REAL, DIMENSION(klon), INTENT(IN)             :: lwnet ! net longwave radiance
70    REAL, DIMENSION(klon), INTENT(IN)             :: tsurf
71    REAL, DIMENSION(klon), INTENT(IN)             :: p1lay
72    REAL, DIMENSION(klon), INTENT(IN)             :: cdragh, cdragm
73    REAL, DIMENSION(klon), INTENT(IN)             :: precip_rain, precip_snow, precip_bs
74    REAL, DIMENSION(klon), INTENT(IN)             :: temp_air, spechum
75    REAL, DIMENSION(klon), INTENT(IN)             :: AcoefH, AcoefQ
76    REAL, DIMENSION(klon), INTENT(IN)             :: BcoefH, BcoefQ
77    REAL, DIMENSION(klon), INTENT(IN)             :: AcoefU, AcoefV, BcoefU, BcoefV
78    REAL, DIMENSION(klon), INTENT(IN)             :: AcoefQBS, BcoefQBS
79    REAL, DIMENSION(klon), INTENT(IN)             :: ps
80    REAL, DIMENSION(klon), INTENT(IN)             :: u1, v1, gustiness, qbs1
81    REAL, DIMENSION(klon), INTENT(IN)             :: rugoro
82    REAL, DIMENSION(klon,nbsrf), INTENT(IN)       :: pctsrf
83#ifdef ISO
84    REAL, DIMENSION(ntiso,klon), INTENT(IN)       :: xtprecip_rain, xtprecip_snow
85    REAL, DIMENSION(ntiso,klon), INTENT(IN)       :: xtspechum
86#endif
87
88
89    LOGICAL,  INTENT(IN)                          :: debut   !true if first step
90    LOGICAL,  INTENT(IN)                          :: lafin   !true if last step
91    REAL, DIMENSION(klon), INTENT(IN)             :: rlon, rlat
92    REAL, DIMENSION(klon), INTENT(IN)             :: rmu0
93    REAL, DIMENSION(klon), INTENT(IN)             :: lwdownm !ylwdown
94    REAL, DIMENSION(klon), INTENT(IN)             :: albedo  !mean albedo
95    REAL, DIMENSION(klon), INTENT(IN)             :: pphi1   
96    REAL, DIMENSION(klon), INTENT(IN)             :: alt   !mean altitude of the grid box 
97    REAL, DIMENSION(klon), INTENT(IN)             :: slope   !mean slope in grid box 
98    REAL, DIMENSION(klon), INTENT(IN)             :: cloudf  !total cloud fraction
99
100! In/Output variables
101!****************************************************************************************
102    REAL, DIMENSION(klon), INTENT(INOUT)          :: snow, qsol
103    REAL, DIMENSION(klon), INTENT(INOUT)          :: agesno
104    REAL, DIMENSION(klon, nsoilmx), INTENT(INOUT) :: tsoil
105#ifdef ISO
106    REAL, DIMENSION(niso,klon), INTENT(INOUT)     :: xtsnow, xtsol
107    REAL, DIMENSION(niso,klon), INTENT(INOUT)     :: Rland_ice
108#endif
109
110
111! Output variables
112!****************************************************************************************
113    REAL, DIMENSION(klon), INTENT(OUT)            :: qsurf
114    REAL, DIMENSION(klon), INTENT(OUT)            :: z0m, z0h
115!albedo SB >>>
116!    REAL, DIMENSION(klon), INTENT(OUT)            :: alb1  ! new albedo in visible SW interval
117!    REAL, DIMENSION(klon), INTENT(OUT)            :: alb2  ! new albedo in near IR interval
118    REAL, DIMENSION(6), INTENT(IN)                :: SFRWL
119    REAL, DIMENSION(klon,nsw), INTENT(OUT)        :: alb_dir,alb_dif
120!albedo SB <<<
121    REAL, DIMENSION(klon), INTENT(OUT)            :: evap, fluxsens, fluxlat, icesub_lic
122    REAL, DIMENSION(klon), INTENT(OUT)            :: fluxbs
123    REAL, DIMENSION(klon), INTENT(OUT)            :: tsurf_new
124    REAL, DIMENSION(klon), INTENT(OUT)            :: dflux_s, dflux_l     
125    REAL, DIMENSION(klon), INTENT(OUT)            :: flux_u1, flux_v1
126
127    REAL, DIMENSION(klon), INTENT(OUT)           :: alb3
128    REAL, DIMENSION(klon), INTENT(OUT)           :: qsnow   !column water in snow [kg/m2]
129    REAL, DIMENSION(klon), INTENT(OUT)           :: snowhgt !Snow height (m)
130    REAL, DIMENSION(klon), INTENT(OUT)           :: to_ice
131    REAL, DIMENSION(klon), INTENT(OUT)           :: sissnow
132    REAL, DIMENSION(klon), INTENT(OUT)           :: runoff  !Land ice runoff
133#ifdef ISO
134    REAL, DIMENSION(ntiso,klon), INTENT(OUT)     :: xtevap     
135#endif
136 
137
138! Local variables
139!****************************************************************************************
140    REAL, DIMENSION(klon)    :: soilcap, soilflux
141    REAL, DIMENSION(klon)    :: cal, beta, dif_grnd
142    REAL, DIMENSION(klon)    :: zfra, alb_neig
143    REAL, DIMENSION(klon)    :: radsol
144    REAL, DIMENSION(klon)    :: u0, v0, u1_lay, v1_lay, ustar
145    INTEGER                  :: i,j,nt
146    REAL, DIMENSION(klon)    :: fqfonte,ffonte
147    REAL, DIMENSION(klon)    :: run_off_lic_frac
148#ifdef ISO       
149    REAL, PARAMETER          :: t_coup = 273.15
150    REAL, DIMENSION(klon)    :: fqfonte_diag
151    REAL, DIMENSION(klon)    :: fq_fonte_diag
152    REAL, DIMENSION(klon)    ::  snow_evap_diag
153    REAL, DIMENSION(klon)    ::  fqcalving_diag
154    REAL max_eau_sol_diag 
155    REAL, DIMENSION(klon)    ::  runoff_diag
156    REAL, DIMENSION(klon)    ::    run_off_lic_diag
157    REAL                     ::  coeff_rel_diag
158    INTEGER                  :: ixt
159    REAL, DIMENSION(niso,klon) :: xtsnow_prec,xtsol_prec
160    REAL, DIMENSION(klon) :: snow_prec,qsol_prec
161#endif
162
163
164    REAL, DIMENSION(klon)    :: emis_new                  !Emissivity
165    REAL, DIMENSION(klon)    :: swdown,lwdown
166    REAL, DIMENSION(klon)    :: precip_snow_adv, snow_adv !Snow Drift precip./advection (not used in inlandsis)
167    REAL, DIMENSION(klon)    :: erod                      !erosion of surface snow (flux, kg/m2/s like evap)
168    REAL, DIMENSION(klon)    :: zsl_height, wind_velo     !surface layer height, wind spd
169    REAL, DIMENSION(klon)    :: dens_air,  snow_cont_air  !air density; snow content air
170    REAL, DIMENSION(klon)    :: alb_soil                  !albedo of underlying ice
171    REAL, DIMENSION(klon)    :: pexner                    !Exner potential
172    REAL                     :: pref
173    REAL, DIMENSION(klon,nsoilmx) :: tsoil0               !modif
174    REAL                          :: dtis                ! subtimestep
175    LOGICAL                       :: debut_is, lafin_is  ! debut and lafin for inlandsis
176
177    CHARACTER (len = 20)                      :: modname = 'surf_landice'
178    CHARACTER (len = 80)                      :: abort_message
179
180
181    REAL,DIMENSION(klon) :: alb1,alb2
182    REAL                 :: time_tempsmooth,coef_tempsmooth
183    REAL,DIMENSION(klon) :: precip_totsnow, evap_totsnow
184    REAL, DIMENSION (klon,6) :: alb6
185    REAL                   :: esalt
186    REAL                   :: lambdasalt,fluxsalt, csalt, nunu, aa, bb, cc
187    REAL                   :: tau_dens, maxerosion
188    REAL, DIMENSION(klon)  :: ws1, rhod, rhos, ustart0, ustart, qsalt, hsalt
189    REAL, DIMENSION(klon)  :: fluxbs_1, fluxbs_2, bsweight_fresh
190    LOGICAL, DIMENSION(klon) :: ok_remaining_freshsnow
191    REAL  :: ta1, ta2, ta3, z01, z02, z03, coefa, coefb, coefc, coefd
192
193
194! End definition
195!****************************************************************************************
196!FC
197!FC
198   REAL,SAVE :: alb_vis_sno_lic
199  !$OMP THREADPRIVATE(alb_vis_sno_lic)
200   REAL,SAVE :: alb_nir_sno_lic
201  !$OMP THREADPRIVATE(alb_nir_sno_lic)
202  LOGICAL, SAVE :: firstcall = .TRUE.
203  !$OMP THREADPRIVATE(firstcall)
204
205
206!FC firtscall initializations
207!******************************************************************************************
208#ifdef ISO
209#ifdef ISOVERIF
210!     write(*,*) 'surf_land_ice 1499'   
211  DO i=1,knon
212    IF (iso_eau > 0) THEN
213      CALL iso_verif_egalite_choix(xtsnow(iso_eau,i),snow(i), &
214    &                              'surf_land_ice 126',errmax,errmaxrel)
215    ENDIF !IF (iso_eau > 0) THEN     
216  ENDDO !DO i=1,knon 
217#endif
218#endif
219
220  IF (firstcall) THEN
221  alb_vis_sno_lic=0.77
222  CALL getin_p('alb_vis_sno_lic',alb_vis_sno_lic)
223           PRINT*, 'alb_vis_sno_lic',alb_vis_sno_lic
224  alb_nir_sno_lic=0.77
225  CALL getin_p('alb_nir_sno_lic',alb_nir_sno_lic)
226           PRINT*, 'alb_nir_sno_lic',alb_nir_sno_lic
227 
228  DO j=1,knon
229       i = knindex(j)
230       tempsmoothlic(i) = temp_air(j)
231  ENDDO
232  firstcall=.false.
233  ENDIF
234!******************************************************************************************
235
236! Initialize output variables
237    alb3(:) = 999999.
238    alb2(:) = 999999.
239    alb1(:) = 999999.
240    fluxbs(:)=0. 
241    runoff(:) = 0.
242!****************************************************************************************
243! Calculate total absorbed radiance at surface
244!
245!****************************************************************************************
246    radsol(:) = 0.0
247    radsol(1:knon) = swnet(1:knon) + lwnet(1:knon)
248
249!****************************************************************************************
250
251!****************************************************************************************
252!  landice_opt = 0 : soil_model, calcul_flux, fonte_neige, ... 
253!  landice_opt = 1  : prepare and call INterace Lmdz SISvat (INLANDSIS)
254!  landice_opt = 2  : skip surf_landice and use orchidee over all land surfaces
255!****************************************************************************************
256
257
258    IF (landice_opt .EQ. 1) THEN
259
260!****************************************************************************************   
261! CALL to INLANDSIS interface
262!****************************************************************************************
263IF (CPPKEY_INLANDSIS) THEN
264
265#ifdef ISO
266        CALL abort_physic('surf_landice 235','isotopes pas dans INLANDSIS',1)
267#endif
268
269        debut_is=debut
270        lafin_is=.false.
271        ! Suppose zero surface speed
272        u0(:)            = 0.0
273        v0(:)            = 0.0
274
275
276        CALL calcul_flux_wind(knon, dtime, &
277         u0, v0, u1, v1, gustiness, cdragm, &
278         AcoefU, AcoefV, BcoefU, BcoefV, &
279         p1lay, temp_air, &
280         flux_u1, flux_v1)
281
282       
283       ! Set constants and compute some input for SISVAT
284       ! = 1000 hPa
285       ! and calculate incoming flux for SW and LW interval: swdown, lwdown
286       swdown(:)        = 0.0
287       lwdown(:)        = 0.0
288       snow_cont_air(:) = 0.  ! the snow content in air is not a prognostic variable of the model     
289       alb_soil(:)      = 0.4 ! before albedo(:) but here it is the ice albedo that we have to set
290       ustar(:)         = 0.
291       pref             = 100000.       
292       DO i = 1, knon
293          swdown(i)        = swnet(i)/(1-albedo(i))
294          lwdown(i)        = lwdownm(i)
295          wind_velo(i)     = u1(i)**2 + v1(i)**2
296          wind_velo(i)     = wind_velo(i)**0.5
297          pexner(i)        = (p1lay(i)/pref)**(RD/RCPD)
298          dens_air(i)      = p1lay(i)/RD/temp_air(i)  ! dry air density
299          zsl_height(i)    = pphi1(i)/RG     
300          tsoil0(i,:)      = tsoil(i,:) 
301          ustar(i)= (cdragm(i)*(wind_velo(i)**2))**0.5   
302       END DO
303       
304
305
306        dtis=dtime
307
308          IF (lafin) THEN
309            lafin_is=.true.
310          END IF
311
312          CALL surf_inlandsis(knon, rlon, rlat, knindex, itime, dtis, debut_is, lafin_is,&
313            rmu0, swdown, lwdown, albedo, pexner, ps, p1lay, precip_rain, precip_snow,   &
314            zsl_height, wind_velo, ustar, temp_air, dens_air, spechum, tsurf,&
315            rugoro, snow_cont_air, alb_soil, alt, slope, cloudf, &
316            radsol, qsol, tsoil0, snow, zfra, snowhgt, qsnow, to_ice, sissnow,agesno,   &
317            AcoefH, AcoefQ, BcoefH, BcoefQ, cdragm, cdragh, &
318            run_off_lic, fqfonte, ffonte, evap, erod, fluxsens, fluxlat,dflux_s, dflux_l, &
319            tsurf_new, alb1, alb2, alb3, alb6, &
320            emis_new, z0m, z0h, qsurf)
321
322          debut_is=.false.
323
324
325        ! Treatment of snow melting and calving
326
327        ! for consistency with standard LMDZ, add calving to run_off_lic
328        run_off_lic(:)=run_off_lic(:) + to_ice(:)
329
330        DO i = 1, knon
331           ffonte_global(knindex(i),is_lic)    = ffonte(i)
332           fqfonte_global(knindex(i),is_lic)   = fqfonte(i)! net melting= melting - refreezing
333           fqcalving_global(knindex(i),is_lic) = to_ice(i) ! flux
334           runofflic_global(knindex(i)) = run_off_lic(i)
335        ENDDO
336        ! Here, we assume that the calving term is equal to the to_ice term
337        ! (no ice accumulation)
338
339
340ELSE
341       abort_message='Pb de coherence: landice_opt = 1  mais CPP_INLANDSIS = .false.'
342       CALL abort_physic(modname,abort_message,1)
343END IF
344
345
346    ELSE
347
348!****************************************************************************************
349! Soil calculations
350!
351!****************************************************************************************
352
353    ! EV: use calbeta
354    CALL calbeta(dtime, is_lic, knon, snow, qsol, beta, cal, dif_grnd)
355
356
357    ! use soil model and recalculate properly cal
358    IF (soil_model) THEN
359       CALL soil(dtime, is_lic, knon, snow, tsurf, qsol, &
360        & longitude(knindex(1:knon)), latitude(knindex(1:knon)), tsoil, soilcap, soilflux)
361       cal(1:knon) = RCPD / soilcap(1:knon)
362       radsol(1:knon)  = radsol(1:knon) + soilflux(1:knon)
363    ELSE
364       cal = RCPD * calice
365       WHERE (snow > 0.0) cal = RCPD * calsno
366    ENDIF
367
368
369!****************************************************************************************
370! Calulate fluxes
371!
372!****************************************************************************************
373
374! Suppose zero surface speed
375    u0(:)=0.0
376    v0(:)=0.0
377    u1_lay(:) = u1(:) - u0(:)
378    v1_lay(:) = v1(:) - v0(:)
379
380    CALL calcul_fluxs(knon, is_lic, dtime, &
381         tsurf, p1lay, cal, beta, cdragh, cdragh, ps, &
382         precip_rain, precip_snow, snow, qsurf,  &
383         radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, gustiness, &
384         1.,AcoefH, AcoefQ, BcoefH, BcoefQ, &
385         tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
386
387#ifdef ISO
388#ifdef ISOVERIF
389     DO i=1,knon
390       IF (iso_eau > 0) THEN
391         IF (snow(i) > ridicule) THEN
392           CALL iso_verif_egalite_choix(xtsnow(iso_eau,i),snow(i), &
393    &                                   'surf_land_ice 1151',errmax,errmaxrel)
394         ENDIF !IF ((snow(i) > ridicule)) THEN
395       ENDIF !IF (iso_eau > 0) THEN
396     ENDDO !DO i=1,knon 
397#endif
398
399    DO i=1,knon
400      snow_prec(i)=snow(i)
401      DO ixt=1,niso
402        xtsnow_prec(ixt,i)=xtsnow(ixt,i)
403      ENDDO !DO ixt=1,niso
404      ! initialisation:
405      fq_fonte_diag(i)=0.0
406      fqfonte_diag(i)=0.0
407      snow_evap_diag(i)=0.0
408    ENDDO !DO i=1,knon
409#endif         
410
411    CALL calcul_flux_wind(knon, dtime, &
412         u0, v0, u1, v1, gustiness, cdragm, &
413         AcoefU, AcoefV, BcoefU, BcoefV, &
414         p1lay, temp_air, &
415         flux_u1, flux_v1)
416
417
418!****************************************************************************************
419! Calculate albedo
420!
421!****************************************************************************************
422
423! Attantion: alb1 and alb2 are not the same!
424    alb1(1:knon)  = alb_vis_sno_lic
425    alb2(1:knon)  = alb_nir_sno_lic
426
427
428!****************************************************************************************
429! Rugosity
430!
431!****************************************************************************************
432
433if (z0m_landice .GT. 0.) then
434    z0m(1:knon) = z0m_landice
435    z0h(1:knon) = z0m_landice*ratio_z0hz0m_landice
436else
437    ! parameterization of z0=f(T) following measurements in Adelie Land by Amory et al 2017
438    coefa = 0.1658 !0.1862 !Ant
439    coefb = -50.3869 !-55.7718 !Ant
440    ta1 = 253.15 !255. Ant
441    ta2 = 273.15
442    ta3 = 273.15+3
443    z01 = exp(coefa*ta1 + coefb) !~0.2 ! ~0.25 mm
444    z02 = exp(coefa*ta2 + coefb) !~6  !~7 mm
445    z03 = z01
446    coefc = log(z03/z02)/(ta3-ta2)
447    coefd = log(z03)-coefc*ta3
448    time_tempsmooth=2.*86400.
449    coef_tempsmooth=min(1.,dtime/time_tempsmooth)
450    !coef_tempsmooth=0.
451    do j=1,knon
452      i=knindex(j)
453      tempsmoothlic(i)=temp_air(j)*coef_tempsmooth+tempsmoothlic(i)*(1.-coef_tempsmooth)
454      if (tempsmoothlic(i) .lt. ta1) then
455        z0m(j) = z01
456      else if (tempsmoothlic(i).ge.ta1 .and. tempsmoothlic(i).lt.ta2) then
457        z0m(j) = exp(coefa*tempsmoothlic(i) + coefb)
458      else if (tempsmoothlic(i).ge.ta2 .and. tempsmoothlic(i).lt.ta3) then
459        ! if st > 0, melting induce smooth surface
460        z0m(j) = exp(coefc*tempsmoothlic(i) + coefd)
461      else
462        z0m(j) = z03
463      endif
464      z0h(j)=z0m(j)
465    enddo
466
467endif   
468 
469
470!****************************************************************************************
471! Simple blowing snow param
472!****************************************************************************************
473! we proceed in 2 steps:
474! first we erode - if possible -the accumulated snow during the time step
475! then we update the density of the underlying layer and see if we can also erode
476! this layer
477
478
479   if (ok_bs) then
480       fluxbs(:)=0.
481       do j=1,knon
482          ws1(j)=(u1(j)**2+v1(j)**2)**0.5
483          ustar(j)=(cdragm(j)*(u1(j)**2+v1(j)**2))**0.5
484          rhod(j)=p1lay(j)/RD/temp_air(j)
485          ustart0(j) =(log(2.868)-log(1.625))/0.085*sqrt(cdragm(j))
486       enddo
487
488       ! 1st step: erosion of fresh snow accumulated during the time step
489       do j=1, knon
490       if (precip_snow(j) .GT. 0.) then
491           rhos(j)=rhofresh_bs
492           ! blowing snow flux formula used in MAR
493           ustart(j)=ustart0(j)*exp(max(rhoice_bs/rhofresh_bs-rhoice_bs/rhos(j),0.))*exp(max(0.,rhos(j)-rhohard_bs))
494           ! we have multiplied by exp to prevent erosion when rhos>rhohard_bs
495           ! computation of qbs at the top of the saltation layer
496           ! default formulation from MAR model (Amory et al. 2021, Gallee et al. 2001)
497           esalt=1./(c_esalt_bs*max(1.e-6,ustar(j)))
498           hsalt(j)=0.08436*(max(1.e-6,ustar(j))**1.27)
499           qsalt(j)=(max(ustar(j)**2-ustart(j)**2,0.))/(RG*hsalt(j))*esalt
500           ! calculation of erosion (flux positive towards the surface here)
501           ! consistent with implicit resolution of turbulent mixing equation
502           ! Nemoto and Nishimura 2004 show that steady-state saltation is achieved within a time tau_eqsalt_bs of about 10s
503           ! we thus prevent snowerosion (snow particle transfer from the saltation layer to the first model level)
504           ! integrated over tau_eqsalt_bs to exceed the total mass of snow particle in the saltation layer
505           ! (rho*qsalt*hsalt)
506           ! during this first step we also lower bound the erosion to the amount of fresh snow accumulated during the time step
507           maxerosion=min(precip_snow(j),hsalt(j)*qsalt(j)*rhod(j)/tau_eqsalt_bs)
508
509           fluxbs_1(j)=rhod(j)*ws1(j)*cdragh(j)*zeta_bs*(AcoefQBS(j)-qsalt(j)) &
510                   / (1.-rhod(j)*ws1(j)*cdragh(j)*zeta_bs*BcoefQBS(j)*dtime)
511           fluxbs_1(j)=max(-maxerosion,fluxbs_1(j))
512
513           if (precip_snow(j) .gt. abs(fluxbs_1(j))) then
514               ok_remaining_freshsnow(j)=.true.
515               bsweight_fresh(j)=1.
516           else
517               ok_remaining_freshsnow(j)=.false.
518               bsweight_fresh(j)=exp(-(abs(fluxbs_1(j))-precip_snow(j))/precip_snow(j))
519           endif
520       else
521           ok_remaining_freshsnow(j)=.false.
522           fluxbs_1(j)=0.
523           bsweight_fresh(j)=0.
524       endif
525       enddo
526
527
528       ! we now compute the snow age of the overlying layer (snow surface after erosion of the fresh snow accumulated during the time step)
529       ! this is done through the routine albsno
530       CALL albsno(klon,knon,dtime,agesno(:),alb_neig(:), precip_snow(:)+fluxbs_1(:))
531
532       ! 2nd step:
533       ! computation of threshold friction velocity
534       ! which depends on surface snow density
535       do j = 1, knon
536        if (ok_remaining_freshsnow(j)) then
537           fluxbs_2(j)=0.
538        else
539           ! we start eroding the underlying layer
540           ! estimation of snow density
541           ! snow density increases with snow age and
542           ! increases even faster in case of sedimentation of blowing snow or rain
543           tau_dens=max(tau_densmin_bs, tau_dens0_bs*exp(-abs(precip_bs(j))/pbst_bs - &
544                    abs(precip_rain(j))/prt_bs)*exp(-max(tsurf(j)-RTT,0.)))
545           rhos(j)=rhofresh_bs+(rhohard_bs-rhofresh_bs)*(1.-exp(-agesno(j)*86400.0/tau_dens))
546           ! blowing snow flux formula used in MAR
547           ustart(j)=ustart0(j)*exp(max(rhoice_bs/rhofresh_bs-rhoice_bs/rhos(j),0.))*exp(max(0.,rhos(j)-rhohard_bs))
548           ! we have multiplied by exp to prevent erosion when rhos>rhohard_bs
549           ! computation of qbs at the top of the saltation layer
550           ! default formulation from MAR model (Amory et al. 2021, Gallee et al. 2001)
551           esalt=1./(c_esalt_bs*max(1.e-6,ustar(j)))
552           hsalt(j)=0.08436*(max(1.e-6,ustar(j))**1.27)
553           qsalt(j)=(max(ustar(j)**2-ustart(j)**2,0.))/(RG*hsalt(j))*esalt
554           ! calculation of erosion (flux positive towards the surface here)
555           ! consistent with implicit resolution of turbulent mixing equation
556           ! Nemoto and Nishimura 2004 show that steady-state saltation is achieved within a time tau_eqsalt_bs of about 10s
557           ! we thus prevent snowerosion (snow particle transfer from the saltation layer to the first model level)
558           ! integrated over tau_eqsalt_bs to exceed the total mass of snow particle in the saltation layer
559           ! (rho*qsalt*hsalt)
560           maxerosion=hsalt(j)*qsalt(j)*rhod(j)/tau_eqsalt_bs
561           fluxbs_2(j)=rhod(j)*ws1(j)*cdragh(j)*zeta_bs*(AcoefQBS(j)-qsalt(j)) &
562                   / (1.-rhod(j)*ws1(j)*cdragh(j)*zeta_bs*BcoefQBS(j)*dtime)
563           fluxbs_2(j)=max(-maxerosion,fluxbs_2(j))
564         endif
565       enddo
566
567
568
569
570       ! final flux and outputs       
571        do j=1, knon
572              ! total flux is the erosion of fresh snow +
573              ! a fraction of the underlying snow (if all the fresh snow has been eroded)
574              ! the calculation of the fraction is quite delicate since we do not know
575              ! how much time was needed to erode the fresh snow. We assume that this time
576              ! is dt*exp(-(abs(fluxbs1)-precipsnow)/precipsnow)=dt*bsweight_fresh
577
578              fluxbs(j)=fluxbs_1(j)+fluxbs_2(j)*(1.-bsweight_fresh(j))
579              i = knindex(j)
580              zxustartlic(i) = ustart(j)
581              zxrhoslic(i) = rhos(j)
582              zxqsaltlic(i)=qsalt(j)
583        enddo
584
585
586  else ! not ok_bs
587  ! those lines are useful to calculate the snow age
588       CALL albsno(klon,knon,dtime,agesno(:),alb_neig(:), precip_snow(:))
589
590  endif ! if ok_bs
591
592
593
594!****************************************************************************************
595! Calculate snow amount
596!   
597!****************************************************************************************
598    IF (ok_bs) THEN
599      precip_totsnow(:)=precip_snow(:)+precip_bs(:)
600      evap_totsnow(:)=evap(:)-fluxbs(:) ! flux bs is positive towards the surface (snow erosion)
601    ELSE
602      precip_totsnow(:)=precip_snow(:)
603      evap_totsnow(:)=evap(:)
604    ENDIF
605   
606 
607    CALL fonte_neige(knon, is_lic, knindex, dtime, &
608         tsurf, precip_rain, precip_totsnow, &
609         snow, qsol, tsurf_new, evap_totsnow, icesub_lic &
610#ifdef ISO   
611     &  ,fq_fonte_diag,fqfonte_diag,snow_evap_diag,fqcalving_diag     &
612     &  ,max_eau_sol_diag,runoff_diag,run_off_lic_diag,coeff_rel_diag &
613#endif
614     &   )
615
616
617#ifdef ISO
618#ifdef ISOVERIF
619    DO i=1,knon 
620      IF (iso_eau > 0) THEN 
621        CALL iso_verif_egalite_choix(Rland_ice(iso_eau,i),1.0, &
622     &                               'surf_landice_mod 217',errmax,errmaxrel)
623      ENDIF !IF (iso_eau > 0) THEN
624    ENDDO !DO i=1,knon
625#endif
626
627    CALL calcul_iso_surf_lic_vectall(klon,knon, &
628     &    evap,snow_evap_diag,Tsurf_new,snow, &
629     &    fq_fonte_diag,fqfonte_diag,dtime,t_coup, &
630     &    precip_snow,xtprecip_snow,precip_rain,xtprecip_rain, snow_prec,xtsnow_prec, &
631     &    xtspechum,spechum,ps,Rland_ice, &
632     &    xtevap,xtsnow,fqcalving_diag, &
633     &    knindex,is_lic,run_off_lic_diag,coeff_rel_diag &
634     &   )
635
636!        call fonte_neige_export_xtrun_off_lic_0(knon,xtrun_off_lic_0_diag)
637
638#endif
639   
640    WHERE (snow(1 : knon) .LT. 0.0001) agesno(1 : knon) = 0.                                         
641    zfra(1:knon) = MAX(0.0,MIN(1.0,snow(1:knon)/(snow(1:knon)+10.0))) 
642
643
644    END IF ! landice_opt
645
646
647!****************************************************************************************
648! Send run-off on land-ice to coupler if coupled ocean.
649! run_off_lic has been calculated in fonte_neige or surf_inlandsis
650! If landice_opt>=2, corresponding call is done from surf_land_orchidee
651!****************************************************************************************
652    IF (type_ocean=='couple' .AND. landice_opt .LT. 2) THEN
653       ! Compress fraction where run_off_lic is active (here all pctsrf(is_lic))
654       run_off_lic_frac(:)=0.0
655       DO j = 1, knon
656          i = knindex(j)
657          run_off_lic_frac(j) = pctsrf(i,is_lic)
658       ENDDO
659
660       CALL cpl_send_landice_fields(itime, knon, knindex, run_off_lic, run_off_lic_frac)
661    ENDIF
662
663 ! transfer runoff rate [kg/m2/s](!) to physiq for output
664    runoff(1:knon)=run_off_lic(1:knon)/dtime
665
666       snow_o=0.
667       zfra_o = 0.
668       DO j = 1, knon
669           i = knindex(j)
670           snow_o(i) = snow(j)
671           zfra_o(i) = zfra(j)
672       ENDDO
673
674
675!albedo SB >>>
676     select case(NSW)
677     case(2)
678       alb_dir(1:knon,1)=alb1(1:knon)
679       alb_dir(1:knon,2)=alb2(1:knon)
680     case(4)
681       alb_dir(1:knon,1)=alb1(1:knon)
682       alb_dir(1:knon,2)=alb2(1:knon)
683       alb_dir(1:knon,3)=alb2(1:knon)
684       alb_dir(1:knon,4)=alb2(1:knon)
685     case(6)
686       alb_dir(1:knon,1)=alb1(1:knon)
687       alb_dir(1:knon,2)=alb1(1:knon)
688       alb_dir(1:knon,3)=alb1(1:knon)
689       alb_dir(1:knon,4)=alb2(1:knon)
690       alb_dir(1:knon,5)=alb2(1:knon)
691       alb_dir(1:knon,6)=alb2(1:knon)
692
693       IF ((landice_opt .EQ. 1) .AND. (iflag_albcalc .EQ. 2)) THEN
694       alb_dir(1:knon,1)=alb6(1:knon,1)
695       alb_dir(1:knon,2)=alb6(1:knon,2)
696       alb_dir(1:knon,3)=alb6(1:knon,3)
697       alb_dir(1:knon,4)=alb6(1:knon,4)
698       alb_dir(1:knon,5)=alb6(1:knon,5)
699       alb_dir(1:knon,6)=alb6(1:knon,6)
700       ENDIF
701
702     end select
703alb_dif=alb_dir
704!albedo SB <<<
705
706
707  END SUBROUTINE surf_landice
708!
709!****************************************************************************************
710!
711END MODULE surf_landice_mod
712
713
714
Note: See TracBrowser for help on using the repository browser.