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

Last change on this file since 5006 was 4947, checked in by evignon, 7 months ago

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