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

Last change on this file was 5285, checked in by abarral, 3 weeks ago

As discussed internally, remove generic ONLY: ... for new _mod_h modules

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