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

Last change on this file since 5282 was 5282, checked in by abarral, 6 hours ago

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