source: LMDZ6/trunk/libf/phylmdiso/surf_landice_mod.F90 @ 4722

Last change on this file since 4722 was 4673, checked in by evignon, 16 months ago

correction commit precedent

File size: 24.8 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 surface_data,     ONLY : type_ocean, calice, calsno, landice_opt, iflag_albcalc
34    USE fonte_neige_mod,  ONLY : fonte_neige,run_off_lic,fqcalving_global,ffonte_global,fqfonte_global,runofflic_global
35    USE cpl_mod,          ONLY : cpl_send_landice_fields
36    USE calcul_fluxs_mod
37    USE phys_local_var_mod, ONLY : zxrhoslic, zxustartlic
38    USE phys_output_var_mod, ONLY : snow_o,zfra_o
39#ifdef ISO   
40    USE fonte_neige_mod,  ONLY : xtrun_off_lic
41    USE infotrac_phy, ONLY : ntiso,niso
42    USE isotopes_routines_mod, ONLY: calcul_iso_surf_lic_vectall
43#ifdef ISOVERIF
44    use isotopes_mod, ONLY: iso_eau,ridicule
45    use isotopes_verif_mod
46#endif
47#endif
48    USE geometry_mod,     ONLY : longitude,latitude
49
50!FC
51    USE ioipsl_getin_p_mod, ONLY : getin_p
52    USE blowing_snow_ini_mod, ONLY : zeta_bs, pbst_bs, prt_bs, iflag_saltation_bs
53
54#ifdef CPP_INLANDSIS
55    USE surf_inlandsis_mod,  ONLY : surf_inlandsis
56#endif
57
58    USE indice_sol_mod
59
60
61!    INCLUDE "indicesol.h"
62    INCLUDE "dimsoil.h"
63    INCLUDE "YOMCST.h"
64    INCLUDE "clesphys.h"
65
66! Input variables
67!****************************************************************************************
68    INTEGER, INTENT(IN)                           :: itime, knon
69    INTEGER, DIMENSION(klon), INTENT(in)          :: knindex
70    REAL, INTENT(in)                              :: dtime
71    REAL, DIMENSION(klon), INTENT(IN)             :: swnet ! net shortwave radiance
72    REAL, DIMENSION(klon), INTENT(IN)             :: lwnet ! net longwave radiance
73    REAL, DIMENSION(klon), INTENT(IN)             :: tsurf
74    REAL, DIMENSION(klon), INTENT(IN)             :: p1lay
75    REAL, DIMENSION(klon), INTENT(IN)             :: cdragh, cdragm
76    REAL, DIMENSION(klon), INTENT(IN)             :: precip_rain, precip_snow, precip_bs
77    REAL, DIMENSION(klon), INTENT(IN)             :: temp_air, spechum
78    REAL, DIMENSION(klon), INTENT(IN)             :: AcoefH, AcoefQ
79    REAL, DIMENSION(klon), INTENT(IN)             :: BcoefH, BcoefQ
80    REAL, DIMENSION(klon), INTENT(IN)             :: AcoefU, AcoefV, BcoefU, BcoefV
81    REAL, DIMENSION(klon), INTENT(IN)             :: AcoefQBS, BcoefQBS
82    REAL, DIMENSION(klon), INTENT(IN)             :: ps
83    REAL, DIMENSION(klon), INTENT(IN)             :: u1, v1, gustiness, qbs1
84    REAL, DIMENSION(klon), INTENT(IN)             :: rugoro
85    REAL, DIMENSION(klon,nbsrf), INTENT(IN)       :: pctsrf
86#ifdef ISO
87    REAL, DIMENSION(ntiso,klon), INTENT(IN)       :: xtprecip_rain, xtprecip_snow
88    REAL, DIMENSION(ntiso,klon), INTENT(IN)       :: xtspechum
89#endif
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! Output variables
113!****************************************************************************************
114    REAL, DIMENSION(klon), INTENT(OUT)            :: qsurf
115    REAL, DIMENSION(klon), INTENT(OUT)            :: z0m, z0h
116!albedo SB >>>
117!    REAL, DIMENSION(klon), INTENT(OUT)            :: alb1  ! new albedo in visible SW interval
118!    REAL, DIMENSION(klon), INTENT(OUT)            :: alb2  ! new albedo in near IR interval
119    REAL, DIMENSION(6), INTENT(IN)                :: SFRWL
120    REAL, DIMENSION(klon,nsw), INTENT(OUT)        :: alb_dir,alb_dif
121!albedo SB <<<
122    REAL, DIMENSION(klon), INTENT(OUT)            :: evap, fluxsens, fluxlat
123    REAL, DIMENSION(klon), INTENT(OUT)            :: fluxbs
124    REAL, DIMENSION(klon), INTENT(OUT)            :: tsurf_new
125    REAL, DIMENSION(klon), INTENT(OUT)            :: dflux_s, dflux_l     
126    REAL, DIMENSION(klon), INTENT(OUT)            :: flux_u1, flux_v1
127
128    REAL, DIMENSION(klon), INTENT(OUT)           :: alb3
129    REAL, DIMENSION(klon), INTENT(OUT)           :: qsnow   !column water in snow [kg/m2]
130    REAL, DIMENSION(klon), INTENT(OUT)           :: snowhgt !Snow height (m)
131    REAL, DIMENSION(klon), INTENT(OUT)           :: to_ice
132    REAL, DIMENSION(klon), INTENT(OUT)           :: sissnow
133    REAL, DIMENSION(klon), INTENT(OUT)           :: runoff  !Land ice runoff
134#ifdef ISO
135    REAL, DIMENSION(ntiso,klon), INTENT(OUT)     :: xtevap     
136!    real, DIMENSION(niso,klon) :: xtrun_off_lic_0_diag ! est une variable globale de
137!    fonte_neige
138#endif
139 
140
141! Local variables
142!****************************************************************************************
143    REAL, DIMENSION(klon)    :: soilcap, soilflux
144    REAL, DIMENSION(klon)    :: cal, beta, dif_grnd
145    REAL, DIMENSION(klon)    :: zfra, alb_neig
146    REAL, DIMENSION(klon)    :: radsol
147    REAL, DIMENSION(klon)    :: u0, v0, u1_lay, v1_lay, ustar
148    INTEGER                  :: i,j,nt
149    REAL, DIMENSION(klon)    :: fqfonte,ffonte
150    REAL, DIMENSION(klon)    :: run_off_lic_frac
151#ifdef ISO       
152      real, parameter :: t_coup = 273.15
153      real, dimension(klon) :: fqfonte_diag
154      real, dimension(klon) :: fq_fonte_diag
155      real, dimension(klon) ::  snow_evap_diag
156      real, dimension(klon) ::  fqcalving_diag
157      real max_eau_sol_diag 
158      real, dimension(klon) ::  runoff_diag
159      real, dimension(klon) ::    run_off_lic_diag
160      real ::  coeff_rel_diag
161      integer ixt
162      REAL, DIMENSION(niso,klon) :: xtsnow_prec,xtsol_prec
163      REAL, DIMENSION(klon) :: snow_prec,qsol_prec
164!      real, DIMENSION(klon) :: run_off_lic_0_diag
165#endif
166
167    REAL, DIMENSION(klon)    :: emis_new                  !Emissivity
168    REAL, DIMENSION(klon)    :: swdown,lwdown
169    REAL, DIMENSION(klon)    :: precip_snow_adv, snow_adv !Snow Drift precip./advection (not used in inlandsis)
170    REAL, DIMENSION(klon)    :: erod                      !erosion of surface snow (flux, kg/m2/s like evap)
171    REAL, DIMENSION(klon)    :: zsl_height, wind_velo     !surface layer height, wind spd
172    REAL, DIMENSION(klon)    :: dens_air,  snow_cont_air  !air density; snow content air
173    REAL, DIMENSION(klon)    :: alb_soil                  !albedo of underlying ice
174    REAL, DIMENSION(klon)    :: pexner                    !Exner potential
175    REAL                     :: pref
176    REAL, DIMENSION(klon,nsoilmx) :: tsoil0               !modif
177    REAL                          :: dtis                ! subtimestep
178    LOGICAL                       :: debut_is, lafin_is  ! debut and lafin for inlandsis
179
180    CHARACTER (len = 20)                      :: modname = 'surf_landice'
181    CHARACTER (len = 80)                      :: abort_message
182
183
184    REAL,DIMENSION(klon) :: alb1,alb2
185    REAL,DIMENSION(klon) :: precip_totsnow, evap_totsnow
186    REAL, DIMENSION (klon,6) :: alb6
187    REAL                   :: rho0, rhoice, ustart0,esalt, rhod
188    REAL                   :: lambdasalt,fluxsalt, csalt, nunu, aa, bb, cc
189    REAL                   :: tau_dens, tau_dens0, tau_densmin, rhomax, rhohard
190    REAL, DIMENSION(klon)  :: ws1, rhos, ustart, qsalt, hsalt, snowini
191    REAL                   :: maxerosion ! in kg/s
192
193! End definition
194!****************************************************************************************
195!FC
196!FC
197   REAL,SAVE :: alb_vis_sno_lic
198  !$OMP THREADPRIVATE(alb_vis_sno_lic)
199   REAL,SAVE :: alb_nir_sno_lic
200  !$OMP THREADPRIVATE(alb_nir_sno_lic)
201  LOGICAL, SAVE :: firstcall = .TRUE.
202  !$OMP THREADPRIVATE(firstcall)
203
204
205!FC firtscall initializations
206!******************************************************************************************
207#ifdef ISO
208#ifdef ISOVERIF
209!     write(*,*) 'surf_land_ice 1499'   
210     do i=1,knon
211        if (iso_eau.gt.0) then
212             call iso_verif_egalite_choix(xtsnow(iso_eau,i),snow(i), &
213    &            'surf_land_ice 126',errmax,errmaxrel)
214        endif !if (iso_eau.gt.0) then     
215      enddo !do i=1,knon 
216#endif
217#endif
218
219  IF (firstcall) THEN
220  alb_vis_sno_lic=0.77
221  CALL getin_p('alb_vis_sno_lic',alb_vis_sno_lic)
222           PRINT*, 'alb_vis_sno_lic',alb_vis_sno_lic
223  alb_nir_sno_lic=0.77
224  CALL getin_p('alb_nir_sno_lic',alb_nir_sno_lic)
225           PRINT*, 'alb_nir_sno_lic',alb_nir_sno_lic
226 
227  firstcall=.false.
228  ENDIF
229!******************************************************************************************
230
231! Initialize output variables
232    alb3(:) = 999999.
233    alb2(:) = 999999.
234    alb1(:) = 999999.
235    fluxbs(:)=0. 
236    runoff(:) = 0.
237!****************************************************************************************
238! Calculate total absorbed radiance at surface
239!
240!****************************************************************************************
241    radsol(:) = 0.0
242    radsol(1:knon) = swnet(1:knon) + lwnet(1:knon)
243
244!****************************************************************************************
245
246!****************************************************************************************
247!  landice_opt = 0 : soil_model, calcul_flux, fonte_neige, ... 
248!  landice_opt = 1  : prepare and call INterace Lmdz SISvat (INLANDSIS)
249!****************************************************************************************
250
251
252    IF (landice_opt .EQ. 1) THEN
253
254!****************************************************************************************   
255! CALL to INLANDSIS interface
256!****************************************************************************************
257#ifdef CPP_INLANDSIS
258
259#ifdef ISO
260        CALL abort_gcm('surf_landice 235','isotopes pas dans INLANDSIS',1)
261#endif
262
263        debut_is=debut
264        lafin_is=.false.
265        ! Suppose zero surface speed
266        u0(:)            = 0.0
267        v0(:)            = 0.0
268
269
270        CALL calcul_flux_wind(knon, dtime, &
271         u0, v0, u1, v1, gustiness, cdragm, &
272         AcoefU, AcoefV, BcoefU, BcoefV, &
273         p1lay, temp_air, &
274         flux_u1, flux_v1)
275
276       
277       ! Set constants and compute some input for SISVAT
278       ! = 1000 hPa
279       ! and calculate incoming flux for SW and LW interval: swdown, lwdown
280       swdown(:)        = 0.0
281       lwdown(:)        = 0.0
282       snow_cont_air(:) = 0.  ! the snow content in air is not a prognostic variable of the model     
283       alb_soil(:)      = 0.4 ! before albedo(:) but here it is the ice albedo that we have to set
284       ustar(:)         = 0.
285       pref             = 100000.       
286       DO i = 1, knon
287          swdown(i)        = swnet(i)/(1-albedo(i))
288          lwdown(i)        = lwdownm(i)
289          wind_velo(i)     = u1(i)**2 + v1(i)**2
290          wind_velo(i)     = wind_velo(i)**0.5
291          pexner(i)        = (p1lay(i)/pref)**(RD/RCPD)
292          dens_air(i)      = p1lay(i)/RD/temp_air(i)  ! dry air density
293          zsl_height(i)    = pphi1(i)/RG     
294          tsoil0(i,:)      = tsoil(i,:) 
295          ustar(i)= (cdragm(i)*(wind_velo(i)**2))**0.5   
296       END DO
297       
298
299
300        dtis=dtime
301
302          IF (lafin) THEN
303            lafin_is=.true.
304          END IF
305
306          CALL surf_inlandsis(knon, rlon, rlat, knindex, itime, dtis, debut_is, lafin_is,&
307            rmu0, swdown, lwdown, albedo, pexner, ps, p1lay, precip_rain, precip_snow,   &
308            zsl_height, wind_velo, ustar, temp_air, dens_air, spechum, tsurf,&
309            rugoro, snow_cont_air, alb_soil, alt, slope, cloudf, &
310            radsol, qsol, tsoil0, snow, zfra, snowhgt, qsnow, to_ice, sissnow,agesno,   &
311            AcoefH, AcoefQ, BcoefH, BcoefQ, cdragm, cdragh, &
312            run_off_lic, fqfonte, ffonte, evap, erod, fluxsens, fluxlat,dflux_s, dflux_l, &
313            tsurf_new, alb1, alb2, alb3, alb6, &
314            emis_new, z0m, z0h, qsurf)
315
316          debut_is=.false.
317
318
319        ! Treatment of snow melting and calving
320
321        ! for consistency with standard LMDZ, add calving to run_off_lic
322        run_off_lic(:)=run_off_lic(:) + to_ice(:)
323
324        DO i = 1, knon
325           ffonte_global(knindex(i),is_lic)    = ffonte(i)
326           fqfonte_global(knindex(i),is_lic)   = fqfonte(i)! net melting= melting - refreezing
327           fqcalving_global(knindex(i),is_lic) = to_ice(i) ! flux
328           runofflic_global(knindex(i)) = run_off_lic(i)
329        ENDDO
330        ! Here, we assume that the calving term is equal to the to_ice term
331        ! (no ice accumulation)
332
333
334#else
335       abort_message='Pb de coherence: landice_opt = 1  mais CPP_INLANDSIS = .false.'
336       CALL abort_physic(modname,abort_message,1)
337#endif
338
339
340    ELSE
341
342!****************************************************************************************
343! Soil calculations
344!
345!****************************************************************************************
346
347    ! EV: use calbeta
348    CALL calbeta(dtime, is_lic, knon, snow, qsol, beta, cal, dif_grnd)
349
350
351    ! use soil model and recalculate properly cal
352    IF (soil_model) THEN
353       CALL soil(dtime, is_lic, knon, snow, tsurf, qsol, &
354        & longitude(knindex(1:knon)), latitude(knindex(1:knon)), tsoil, soilcap, soilflux)
355       cal(1:knon) = RCPD / soilcap(1:knon)
356       radsol(1:knon)  = radsol(1:knon) + soilflux(1:knon)
357    ELSE
358       cal = RCPD * calice
359       WHERE (snow > 0.0) cal = RCPD * calsno
360    ENDIF
361
362
363!****************************************************************************************
364! Calulate fluxes
365!
366!****************************************************************************************
367!    beta(:) = 1.0
368!    dif_grnd(:) = 0.0
369
370! Suppose zero surface speed
371    u0(:)=0.0
372    v0(:)=0.0
373    u1_lay(:) = u1(:) - u0(:)
374    v1_lay(:) = v1(:) - v0(:)
375
376    CALL calcul_fluxs(knon, is_lic, dtime, &
377         tsurf, p1lay, cal, beta, cdragh, cdragh, ps, &
378         precip_rain, precip_snow, snow, qsurf,  &
379         radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, gustiness, &
380         1.,AcoefH, AcoefQ, BcoefH, BcoefQ, &
381         tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
382
383
384#ifdef ISO
385   ! verif
386#ifdef ISOVERIF
387     !write(*,*) 'surf_land_ice 1499'   
388     do i=1,knon
389       if (iso_eau.gt.0) then
390           if (snow(i).gt.ridicule) then
391             call iso_verif_egalite_choix(xtsnow(iso_eau,i),snow(i), &
392    &            'surf_land_ice 1151',errmax,errmaxrel)
393            endif !if ((snow(i).gt.ridicule)) then
394        endif !if (iso_eau.gt.0) then     
395      enddo !do i=1,knon 
396#endif
397!#ifdef ISOVERIF
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!#ifdef ISO
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    CALL albsno(klon,knon,dtime,agesno(:),alb_neig(:), precip_snow(:)) 
424
425
426! EV: following lines are obsolete since we set alb1 and alb2 to constant values
427! I therefore comment them
428!    alb1(1:knon) = alb_neig(1:knon)*zfra(1:knon) + &
429!         0.6 * (1.0-zfra(1:knon))
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    z0m = z0m_landice
448    z0h = z0h_landice
449    !z0m = SQRT(z0m**2+rugoro**2)
450
451
452
453  ! Simple blowing snow param
454  if (ok_bs) then
455       ustart0 = 0.211
456       rhoice = 920.0
457       rho0 = 300.0
458       rhomax=450.0
459       rhohard=450.0
460       tau_dens0=86400.0*10.  ! 10 days by default, in s
461       tau_densmin=21600.0    ! 1/4 days according to in situ obs by C. Amory during blowing snow +
462                              ! Marshall et al. 1999 (snow densification during rain)
463
464       ! computation of threshold friction velocity
465       ! which depends on surface snow density
466       do j = 1, knon
467           ! estimation of snow density
468           ! snow density increases with snow age and
469           ! increases even faster in case of sedimentation of blowing snow or rain
470           tau_dens=max(tau_densmin, tau_dens0*exp(-abs(precip_bs(j))/pbst_bs &
471                   -abs(precip_rain(j))/prt_bs) *exp(-max(tsurf(j)-272.0,0.)))
472           rhos(j)=rho0+(rhohard-rho0)*(1.-exp(-agesno(j)*86400.0/tau_dens))
473           ! blowing snow flux formula used in MAR
474           ws1(j)=(u1(j)**2+v1(j)**2)**0.5
475           ustar(j)=(cdragm(j)*(u1(j)**2+v1(j)**2))**0.5
476           ustart(j)=ustart0*exp(max(rhoice/rho0-rhoice/rhos(j),0.))*exp(max(0.,rhos(j)-rhomax))
477           ! we have multiplied by exp to prevent erosion when rhos>rhomax (usefull till
478           ! rhohard<450)
479       enddo
480       
481       ! computation of qbs at the top of the saltation layer
482       ! two formulations possible
483       ! pay attention that qbs is a mixing ratio and has to be converted
484       ! to specific content
485       
486       if (iflag_saltation_bs .eq. 1) then
487       ! expression from CRYOWRF (Sharma et al. 2022)
488          aa=2.6
489          bb=2.5
490          cc=2.0
491          lambdasalt=0.45
492          do j =1, knon
493               rhod=p1lay(j)/RD/temp_air(j)
494               nunu=max(ustar(j)/ustart(j),1.e-3)
495               fluxsalt=rhod/RG*(ustar(j)**3)*(1.-nunu**(-2)) * &
496                        (aa+bb*nunu**(-2)+cc*nunu**(-1)) 
497               csalt=fluxsalt/(2.8*ustart(j))
498               hsalt(j)=0.08436*ustar(j)**1.27
499               qsalt(j)=1./rhod*csalt*lambdasalt*RG/(max(ustar(j)**2,1E-6)) &
500                       * exp(-lambdasalt*RG*hsalt(j)/max(ustar(j)**2,1E-6))
501               qsalt(j)=max(qsalt(j),0.)
502          enddo
503
504
505       else
506       ! default formulation from MAR model (Amory et al. 2021, Gallee et al. 2001)       
507          do j=1, knon
508              esalt=1./(3.25*max(ustar(j),0.001))
509              hsalt(j)=0.08436*ustar(j)**1.27
510              qsalt(j)=(max(ustar(j)**2-ustart(j)**2,0.))/(RG*hsalt(j))*esalt
511              !ep=qsalt*cdragm(j)*sqrt(u1(j)**2+v1(j)**2)
512          enddo
513       endif
514
515        ! calculation of erosion (flux positive towards the surface here)
516        ! consistent with implicit resolution of turbulent mixing equation
517       do j=1, knon
518              i = knindex(j)
519              rhod=p1lay(j)/RD/temp_air(j)
520              ! Nemoto and Nishimura 2004 show that steady-state saltation is achieved within an interval of about 10s
521              ! we thus prevent snowerosion (snow particle transfer from the saltation layer to the first model level)
522              ! to exceed (in intensity) integrated over 10s to exceed the total mass of snow particle in the saltation layer
523              ! (rho*qsalt*hsalt)
524              maxerosion=hsalt(j)*qsalt(j)*rhod/10.
525              fluxbs(j)=rhod*ws1(j)*cdragm(j)*(AcoefQBS(j)-qsalt(j)) &
526                       / (1.-rhod*ws1(j)*cdragm(j)*BcoefQBS(j)*dtime)
527              !fluxbs(j)= rhod*ws1(j)*cdragm(j)*(qbs1(j)-qsalt(j))
528              fluxbs(j)=max(-maxerosion,fluxbs(j))
529
530              ! for outputs
531
532              zxustartlic(i) = ustart(j)
533              zxrhoslic(i) = rhos(j)
534       enddo
535
536  endif
537
538
539
540!****************************************************************************************
541! Calculate surface snow amount
542!   
543!****************************************************************************************
544    IF (ok_bs) THEN
545      precip_totsnow(:)=precip_snow(:)+precip_bs(:)
546      evap_totsnow(:)=evap(:)-fluxbs(:) ! flux bs is positive towards the surface (snow erosion)
547    ELSE
548      precip_totsnow(:)=precip_snow(:)
549      evap_totsnow(:)=evap(:)
550    ENDIF
551
552    CALL fonte_neige(knon, is_lic, knindex, dtime, &
553         tsurf, precip_rain, precip_totsnow, &
554         snow, qsol, tsurf_new, evap_totsnow &
555#ifdef ISO   
556     & ,fq_fonte_diag,fqfonte_diag,snow_evap_diag,fqcalving_diag   &
557     & ,max_eau_sol_diag,runoff_diag,run_off_lic_diag,coeff_rel_diag   &
558#endif
559     &   )
560
561
562#ifdef ISO
563#ifdef ISOVERIF
564      do i=1,knon 
565       if (iso_eau.gt.0) then 
566          call iso_verif_egalite_choix(Rland_ice(iso_eau,i),1.0, &
567     &         'surf_landice_mod 217',errmax,errmaxrel)
568       endif !if (iso_eau.gt.0) then   
569      enddo !do i=1,knon 
570#endif
571!#ifdef ISOVERIF
572
573       call calcul_iso_surf_lic_vectall(klon,knon, &
574     &           evap,snow_evap_diag,Tsurf_new,snow, &
575     &           fq_fonte_diag,fqfonte_diag,dtime,t_coup, &
576     &           precip_snow,xtprecip_snow,precip_rain,xtprecip_rain, snow_prec,xtsnow_prec, &
577     &           xtspechum,spechum,ps,Rland_ice, &
578     &           xtevap,xtsnow,fqcalving_diag, &
579     &           knindex,is_lic,run_off_lic_diag,coeff_rel_diag &
580     &  )
581
582!        call fonte_neige_export_xtrun_off_lic_0(knon,xtrun_off_lic_0_diag)
583
584#endif
585!#ifdef ISO
586
587
588    WHERE (snow(1 : knon) .LT. 0.0001) agesno(1 : knon) = 0.                                         
589    zfra(1:knon) = MAX(0.0,MIN(1.0,snow(1:knon)/(snow(1:knon)+10.0))) 
590
591
592    END IF ! landice_opt
593
594
595!****************************************************************************************
596! Send run-off on land-ice to coupler if coupled ocean.
597! run_off_lic has been calculated in fonte_neige or surf_inlandsis
598! If landice_opt>=2, corresponding call is done from surf_land_orchidee
599!****************************************************************************************
600    IF (type_ocean=='couple' .AND. landice_opt .LT. 2) THEN
601       ! Compress fraction where run_off_lic is active (here all pctsrf(is_lic))
602       run_off_lic_frac(:)=0.0
603       DO j = 1, knon
604          i = knindex(j)
605          run_off_lic_frac(j) = pctsrf(i,is_lic)
606       ENDDO
607
608       CALL cpl_send_landice_fields(itime, knon, knindex, run_off_lic, run_off_lic_frac)
609    ENDIF
610
611 ! transfer runoff rate [kg/m2/s](!) to physiq for output
612    runoff(1:knon)=run_off_lic(1:knon)/dtime
613
614       snow_o=0.
615       zfra_o = 0.
616       DO j = 1, knon
617           i = knindex(j)
618           snow_o(i) = snow(j)
619           zfra_o(i) = zfra(j)
620       ENDDO
621
622
623!albedo SB >>>
624     select case(NSW)
625     case(2)
626       alb_dir(1:knon,1)=alb1(1:knon)
627       alb_dir(1:knon,2)=alb2(1:knon)
628     case(4)
629       alb_dir(1:knon,1)=alb1(1:knon)
630       alb_dir(1:knon,2)=alb2(1:knon)
631       alb_dir(1:knon,3)=alb2(1:knon)
632       alb_dir(1:knon,4)=alb2(1:knon)
633     case(6)
634       alb_dir(1:knon,1)=alb1(1:knon)
635       alb_dir(1:knon,2)=alb1(1:knon)
636       alb_dir(1:knon,3)=alb1(1:knon)
637       alb_dir(1:knon,4)=alb2(1:knon)
638       alb_dir(1:knon,5)=alb2(1:knon)
639       alb_dir(1:knon,6)=alb2(1:knon)
640
641       IF ((landice_opt .EQ. 1) .AND. (iflag_albcalc .EQ. 2)) THEN
642       alb_dir(1:knon,1)=alb6(1:knon,1)
643       alb_dir(1:knon,2)=alb6(1:knon,2)
644       alb_dir(1:knon,3)=alb6(1:knon,3)
645       alb_dir(1:knon,4)=alb6(1:knon,4)
646       alb_dir(1:knon,5)=alb6(1:knon,5)
647       alb_dir(1:knon,6)=alb6(1:knon,6)
648       ENDIF
649
650     end select
651alb_dif=alb_dir
652!albedo SB <<<
653
654
655  END SUBROUTINE surf_landice
656!
657!****************************************************************************************
658!
659END MODULE surf_landice_mod
660
661
662
Note: See TracBrowser for help on using the repository browser.