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

Last change on this file since 4797 was 4724, checked in by evignon, 13 months ago

debut du travail de replayisation des routines de neige soufflee
+ correction bug pour les temperatures positives

  • 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: 21.2 KB
RevLine 
[781]1!
2MODULE surf_landice_mod
3 
4  IMPLICIT NONE
5
6CONTAINS
7!
8!****************************************************************************************
9!
10  SUBROUTINE surf_landice(itime, dtime, knon, knindex, &
[1865]11       rlon, rlat, debut, lafin, &
12       rmu0, lwdownm, albedo, pphi1, &
[888]13       swnet, lwnet, tsurf, p1lay, &
[4523]14       cdragh, cdragm, precip_rain, precip_snow, precip_bs, temp_air, spechum, &
[1067]15       AcoefH, AcoefQ, BcoefH, BcoefQ, &
16       AcoefU, AcoefV, BcoefU, BcoefV, &
[4529]17       AcoefQBS, BcoefQBS, &
[2240]18       ps, u1, v1, gustiness, rugoro, pctsrf, &
[4523]19       snow, qsurf, qsol, qbs1, agesno, &
20       tsoil, z0m, z0h, SFRWL, alb_dir, alb_dif, evap, fluxsens, fluxlat, fluxbs, &
[1067]21       tsurf_new, dflux_s, dflux_l, &
[3900]22       alt, slope, cloudf, &
[1865]23       snowhgt, qsnow, to_ice, sissnow, &
24       alb3, runoff, &
[1067]25       flux_u1, flux_v1)
[781]26
[1067]27    USE dimphy
[3974]28    USE geometry_mod,     ONLY : longitude,latitude
[3900]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
[1067]31    USE cpl_mod,          ONLY : cpl_send_landice_fields
32    USE calcul_fluxs_mod
[4523]33    USE phys_local_var_mod, ONLY : zxrhoslic, zxustartlic
[4414]34    USE phys_output_var_mod, ONLY : snow_o,zfra_o
[2728]35!FC
36    USE ioipsl_getin_p_mod, ONLY : getin_p
[4724]37    USE lmdz_blowing_snow_ini, ONLY : zeta_bs, pbst_bs, prt_bs, iflag_saltation_bs
[2728]38
[3792]39#ifdef CPP_INLANDSIS
40    USE surf_inlandsis_mod,  ONLY : surf_inlandsis
41#endif
42
[1785]43    USE indice_sol_mod
[1067]44
[1785]45!    INCLUDE "indicesol.h"
[781]46    INCLUDE "dimsoil.h"
[793]47    INCLUDE "YOMCST.h"
48    INCLUDE "clesphys.h"
[781]49
50! Input variables
51!****************************************************************************************
52    INTEGER, INTENT(IN)                           :: itime, knon
53    INTEGER, DIMENSION(klon), INTENT(in)          :: knindex
54    REAL, INTENT(in)                              :: dtime
[888]55    REAL, DIMENSION(klon), INTENT(IN)             :: swnet ! net shortwave radiance
56    REAL, DIMENSION(klon), INTENT(IN)             :: lwnet ! net longwave radiance
[781]57    REAL, DIMENSION(klon), INTENT(IN)             :: tsurf
58    REAL, DIMENSION(klon), INTENT(IN)             :: p1lay
[1067]59    REAL, DIMENSION(klon), INTENT(IN)             :: cdragh, cdragm
[4523]60    REAL, DIMENSION(klon), INTENT(IN)             :: precip_rain, precip_snow, precip_bs
[781]61    REAL, DIMENSION(klon), INTENT(IN)             :: temp_air, spechum
[1067]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
[4529]65    REAL, DIMENSION(klon), INTENT(IN)             :: AcoefQBS, BcoefQBS
[781]66    REAL, DIMENSION(klon), INTENT(IN)             :: ps
[4523]67    REAL, DIMENSION(klon), INTENT(IN)             :: u1, v1, gustiness, qbs1
[781]68    REAL, DIMENSION(klon), INTENT(IN)             :: rugoro
69    REAL, DIMENSION(klon,nbsrf), INTENT(IN)       :: pctsrf
70
[1865]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   
[3900]78    REAL, DIMENSION(klon), INTENT(IN)             :: alt   !mean altitude of the grid box 
[1865]79    REAL, DIMENSION(klon), INTENT(IN)             :: slope   !mean slope in grid box 
80    REAL, DIMENSION(klon), INTENT(IN)             :: cloudf  !total cloud fraction
81
[781]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
[2243]91    REAL, DIMENSION(klon), INTENT(OUT)            :: z0m, z0h
[2227]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
[3792]95    REAL, DIMENSION(6), INTENT(IN)                :: SFRWL
96    REAL, DIMENSION(klon,nsw), INTENT(OUT)        :: alb_dir,alb_dif
[2227]97!albedo SB <<<
[781]98    REAL, DIMENSION(klon), INTENT(OUT)            :: evap, fluxsens, fluxlat
[4523]99    REAL, DIMENSION(klon), INTENT(OUT)            :: fluxbs
[888]100    REAL, DIMENSION(klon), INTENT(OUT)            :: tsurf_new
[781]101    REAL, DIMENSION(klon), INTENT(OUT)            :: dflux_s, dflux_l     
[1067]102    REAL, DIMENSION(klon), INTENT(OUT)            :: flux_u1, flux_v1
[781]103
[1865]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
[781]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
[888]117    REAL, DIMENSION(klon)    :: radsol
[3792]118    REAL, DIMENSION(klon)    :: u0, v0, u1_lay, v1_lay, ustar
119    INTEGER                  :: i,j,nt
[3900]120    REAL, DIMENSION(klon)    :: fqfonte,ffonte
[4283]121    REAL, DIMENSION(klon)    :: run_off_lic_frac
[1865]122    REAL, DIMENSION(klon)    :: emis_new                  !Emissivity
123    REAL, DIMENSION(klon)    :: swdown,lwdown
[3900]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
[1865]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
[3792]131    REAL, DIMENSION(klon,nsoilmx) :: tsoil0               !modif
132    REAL                          :: dtis                ! subtimestep
133    LOGICAL                       :: debut_is, lafin_is  ! debut and lafin for inlandsis
[1865]134
135    CHARACTER (len = 20)                      :: modname = 'surf_landice'
136    CHARACTER (len = 80)                      :: abort_message
137
[2728]138
[3900]139    REAL,DIMENSION(klon) :: alb1,alb2
[4523]140    REAL,DIMENSION(klon) :: precip_totsnow, evap_totsnow
[3900]141    REAL, DIMENSION (klon,6) :: alb6
[4672]142    REAL                   :: rho0, rhoice, ustart0,esalt, rhod
[4529]143    REAL                   :: lambdasalt,fluxsalt, csalt, nunu, aa, bb, cc
[4523]144    REAL                   :: tau_dens, tau_dens0, tau_densmin, rhomax, rhohard
[4672]145    REAL, DIMENSION(klon)  :: ws1, rhos, ustart, qsalt, hsalt, snowini
146    REAL                   :: maxerosion ! in kg/s   
147
[781]148! End definition
149!****************************************************************************************
[2728]150!FC
151!FC
152   REAL,SAVE :: alb_vis_sno_lic
153  !$OMP THREADPRIVATE(alb_vis_sno_lic)
154   REAL,SAVE :: alb_nir_sno_lic
155  !$OMP THREADPRIVATE(alb_nir_sno_lic)
156  LOGICAL, SAVE :: firstcall = .TRUE.
157  !$OMP THREADPRIVATE(firstcall)
158
159
[3792]160!FC firtscall initializations
161!******************************************************************************************
[2728]162  IF (firstcall) THEN
163  alb_vis_sno_lic=0.77
164  CALL getin_p('alb_vis_sno_lic',alb_vis_sno_lic)
165           PRINT*, 'alb_vis_sno_lic',alb_vis_sno_lic
166  alb_nir_sno_lic=0.77
167  CALL getin_p('alb_nir_sno_lic',alb_nir_sno_lic)
168           PRINT*, 'alb_nir_sno_lic',alb_nir_sno_lic
[3792]169 
[2728]170  firstcall=.false.
171  ENDIF
[3792]172!******************************************************************************************
173
[781]174! Initialize output variables
[1865]175    alb3(:) = 999999.
[888]176    alb2(:) = 999999.
177    alb1(:) = 999999.
[4523]178    fluxbs(:)=0. 
[1865]179    runoff(:) = 0.
[888]180!****************************************************************************************
181! Calculate total absorbed radiance at surface
182!
183!****************************************************************************************
184    radsol(:) = 0.0
185    radsol(1:knon) = swnet(1:knon) + lwnet(1:knon)
[781]186
187!****************************************************************************************
[4523]188
189!****************************************************************************************
[3792]190!  landice_opt = 0 : soil_model, calcul_flux, fonte_neige, ... 
[3901]191!  landice_opt = 1  : prepare and call INterace Lmdz SISvat (INLANDSIS)
[1865]192!****************************************************************************************
[3792]193
194
195    IF (landice_opt .EQ. 1) THEN
196
[3901]197!****************************************************************************************   
[3792]198! CALL to INLANDSIS interface
199!****************************************************************************************
200#ifdef CPP_INLANDSIS
201
202        debut_is=debut
203        lafin_is=.false.
204        ! Suppose zero surface speed
205        u0(:)            = 0.0
206        v0(:)            = 0.0
207
208
209        CALL calcul_flux_wind(knon, dtime, &
210         u0, v0, u1, v1, gustiness, cdragm, &
211         AcoefU, AcoefV, BcoefU, BcoefV, &
212         p1lay, temp_air, &
213         flux_u1, flux_v1)
214
215       
216       ! Set constants and compute some input for SISVAT
217       ! = 1000 hPa
218       ! and calculate incoming flux for SW and LW interval: swdown, lwdown
219       swdown(:)        = 0.0
220       lwdown(:)        = 0.0
[3900]221       snow_cont_air(:) = 0.  ! the snow content in air is not a prognostic variable of the model     
[3792]222       alb_soil(:)      = 0.4 ! before albedo(:) but here it is the ice albedo that we have to set
223       ustar(:)         = 0.
224       pref             = 100000.       
225       DO i = 1, knon
226          swdown(i)        = swnet(i)/(1-albedo(i))
227          lwdown(i)        = lwdownm(i)
228          wind_velo(i)     = u1(i)**2 + v1(i)**2
229          wind_velo(i)     = wind_velo(i)**0.5
230          pexner(i)        = (p1lay(i)/pref)**(RD/RCPD)
231          dens_air(i)      = p1lay(i)/RD/temp_air(i)  ! dry air density
232          zsl_height(i)    = pphi1(i)/RG     
233          tsoil0(i,:)      = tsoil(i,:) 
234          ustar(i)= (cdragm(i)*(wind_velo(i)**2))**0.5   
235       END DO
236       
237
238
[3900]239        dtis=dtime
[3792]240
[3900]241          IF (lafin) THEN
[3792]242            lafin_is=.true.
243          END IF
244
[3900]245          CALL surf_inlandsis(knon, rlon, rlat, knindex, itime, dtis, debut_is, lafin_is,&
246            rmu0, swdown, lwdown, albedo, pexner, ps, p1lay, precip_rain, precip_snow,   &
247            zsl_height, wind_velo, ustar, temp_air, dens_air, spechum, tsurf,&
248            rugoro, snow_cont_air, alb_soil, alt, slope, cloudf, &
249            radsol, qsol, tsoil0, snow, zfra, snowhgt, qsnow, to_ice, sissnow,agesno,   &
[3792]250            AcoefH, AcoefQ, BcoefH, BcoefQ, cdragm, cdragh, &
[3900]251            run_off_lic, fqfonte, ffonte, evap, erod, fluxsens, fluxlat,dflux_s, dflux_l, &
252            tsurf_new, alb1, alb2, alb3, alb6, &
253            emis_new, z0m, z0h, qsurf)
[3792]254
[3900]255          debut_is=.false.
[3792]256
257
[3900]258        ! Treatment of snow melting and calving
[3792]259
[3900]260        ! for consistency with standard LMDZ, add calving to run_off_lic
261        run_off_lic(:)=run_off_lic(:) + to_ice(:)
262
263        DO i = 1, knon
264           ffonte_global(knindex(i),is_lic)    = ffonte(i)
265           fqfonte_global(knindex(i),is_lic)   = fqfonte(i)! net melting= melting - refreezing
266           fqcalving_global(knindex(i),is_lic) = to_ice(i) ! flux
267           runofflic_global(knindex(i)) = run_off_lic(i)
268        ENDDO
269        ! Here, we assume that the calving term is equal to the to_ice term
270        ! (no ice accumulation)
271
272
[3792]273#else
[3901]274       abort_message='Pb de coherence: landice_opt = 1  mais CPP_INLANDSIS = .false.'
[3792]275       CALL abort_physic(modname,abort_message,1)
276#endif
277
278
279    ELSE
280
281!****************************************************************************************
[781]282! Soil calculations
283!
284!****************************************************************************************
[3780]285
286    ! EV: use calbeta
287    CALL calbeta(dtime, is_lic, knon, snow, qsol, beta, cal, dif_grnd)
288
289
290    ! use soil model and recalculate properly cal
[781]291    IF (soil_model) THEN
[3974]292       CALL soil(dtime, is_lic, knon, snow, tsurf, qsol, &
293        & longitude(knindex(1:knon)), latitude(knindex(1:knon)), tsoil, soilcap, soilflux)
[781]294       cal(1:knon) = RCPD / soilcap(1:knon)
295       radsol(1:knon)  = radsol(1:knon) + soilflux(1:knon)
296    ELSE
297       cal = RCPD * calice
298       WHERE (snow > 0.0) cal = RCPD * calsno
299    ENDIF
300
301
302!****************************************************************************************
303! Calulate fluxes
304!
305!****************************************************************************************
[3792]306!    beta(:) = 1.0
307!    dif_grnd(:) = 0.0
[781]308
[1067]309! Suppose zero surface speed
310    u0(:)=0.0
311    v0(:)=0.0
312    u1_lay(:) = u1(:) - u0(:)
313    v1_lay(:) = v1(:) - v0(:)
314
[781]315    CALL calcul_fluxs(knon, is_lic, dtime, &
[2254]316         tsurf, p1lay, cal, beta, cdragh, cdragh, ps, &
[781]317         precip_rain, precip_snow, snow, qsurf,  &
[2240]318         radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, gustiness, &
[2254]319         1.,AcoefH, AcoefQ, BcoefH, BcoefQ, &
[781]320         tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
321
[1067]322    CALL calcul_flux_wind(knon, dtime, &
[2240]323         u0, v0, u1, v1, gustiness, cdragm, &
[1067]324         AcoefU, AcoefV, BcoefU, BcoefV, &
325         p1lay, temp_air, &
326         flux_u1, flux_v1)
[781]327
328
329!****************************************************************************************
330! Calculate albedo
331!
332!****************************************************************************************
[4523]333 
[781]334    CALL albsno(klon,knon,dtime,agesno(:),alb_neig(:), precip_snow(:)) 
[3780]335
[4523]336
337! EV: following lines are obsolete since we set alb1 and alb2 to constant values
338! I therefore comment them
339!    alb1(1:knon) = alb_neig(1:knon)*zfra(1:knon) + &
340!         0.6 * (1.0-zfra(1:knon))
[781]341!
342!IM: plusieurs choix/tests sur l'albedo des "glaciers continentaux"
[888]343!       alb1(1 : knon)  = 0.6 !IM cf FH/GK
344!       alb1(1 : knon)  = 0.82
345!       alb1(1 : knon)  = 0.77 !211003 Ksta0.77
346!       alb1(1 : knon)  = 0.8 !KstaTER0.8 & LMD_ARMIP5
[781]347!IM: KstaTER0.77 & LMD_ARMIP6   
348
[3780]349! Attantion: alb1 and alb2 are not the same!
[2728]350    alb1(1:knon)  = alb_vis_sno_lic
351    alb2(1:knon)  = alb_nir_sno_lic
[781]352
353
354!****************************************************************************************
355! Rugosity
356!
357!****************************************************************************************
[4245]358    z0m = z0m_landice
359    z0h = z0h_landice
360    !z0m = SQRT(z0m**2+rugoro**2)
[2243]361
[781]362
[1865]363
[4523]364  ! Simple blowing snow param
365  if (ok_bs) then
366       ustart0 = 0.211
367       rhoice = 920.0
[4672]368       rho0 = 300.0
[4523]369       rhomax=450.0
[4672]370       rhohard=450.0
[4523]371       tau_dens0=86400.0*10.  ! 10 days by default, in s
[4672]372       tau_densmin=21600.0    ! 1/4 days according to in situ obs by C. Amory during blowing snow +
373                              ! Marshall et al. 1999 (snow densification during rain)
[4529]374
375       ! computation of threshold friction velocity
376       ! which depends on surface snow density
[4672]377       do j = 1, knon
[4523]378           ! estimation of snow density
379           ! snow density increases with snow age and
[4672]380           ! increases even faster in case of sedimentation of blowing snow or rain
[4673]381           tau_dens=max(tau_densmin, tau_dens0*exp(-abs(precip_bs(j))/pbst_bs - &
382                    abs(precip_rain(j))/prt_bs)*exp(-max(tsurf(j)-272.0,0.)))
[4672]383           rhos(j)=rho0+(rhohard-rho0)*(1.-exp(-agesno(j)*86400.0/tau_dens))
[4523]384           ! blowing snow flux formula used in MAR
[4672]385           ws1(j)=(u1(j)**2+v1(j)**2)**0.5
386           ustar(j)=(cdragm(j)*(u1(j)**2+v1(j)**2))**0.5
387           ustart(j)=ustart0*exp(max(rhoice/rho0-rhoice/rhos(j),0.))*exp(max(0.,rhos(j)-rhomax))
[4523]388           ! we have multiplied by exp to prevent erosion when rhos>rhomax (usefull till
389           ! rhohard<450)
390       enddo
[4529]391       
392       ! computation of qbs at the top of the saltation layer
393       ! two formulations possible
394       ! pay attention that qbs is a mixing ratio and has to be converted
395       ! to specific content
396       
397       if (iflag_saltation_bs .eq. 1) then
398       ! expression from CRYOWRF (Sharma et al. 2022)
399          aa=2.6
400          bb=2.5
401          cc=2.0
402          lambdasalt=0.45
[4672]403          do j =1, knon
404               rhod=p1lay(j)/RD/temp_air(j)
405               nunu=max(ustar(j)/ustart(j),1.e-3)
406               fluxsalt=rhod/RG*(ustar(j)**3)*(1.-nunu**(-2)) * &
[4529]407                        (aa+bb*nunu**(-2)+cc*nunu**(-1)) 
[4672]408               csalt=fluxsalt/(2.8*ustart(j))
409               hsalt(j)=0.08436*ustar(j)**1.27
410               qsalt(j)=1./rhod*csalt*lambdasalt*RG/(max(ustar(j)**2,1E-6)) &
411                       * exp(-lambdasalt*RG*hsalt(j)/max(ustar(j)**2,1E-6))
412               qsalt(j)=max(qsalt(j),0.)
[4529]413          enddo
[4523]414
[4529]415
416       else
417       ! default formulation from MAR model (Amory et al. 2021, Gallee et al. 2001)       
[4672]418          do j=1, knon
419              esalt=1./(3.25*max(ustar(j),0.001))
420              hsalt(j)=0.08436*ustar(j)**1.27
421              qsalt(j)=(max(ustar(j)**2-ustart(j)**2,0.))/(RG*hsalt(j))*esalt
422              !ep=qsalt*cdragm(j)*sqrt(u1(j)**2+v1(j)**2)
[4529]423          enddo
424       endif
425
[4672]426        ! calculation of erosion (flux positive towards the surface here)
[4529]427        ! consistent with implicit resolution of turbulent mixing equation
[4672]428       do j=1, knon
429              i = knindex(j)
430              rhod=p1lay(j)/RD/temp_air(j)
431              ! Nemoto and Nishimura 2004 show that steady-state saltation is achieved within an interval of about 10s
432              ! we thus prevent snowerosion (snow particle transfer from the saltation layer to the first model level)
433              ! to exceed (in intensity) integrated over 10s to exceed the total mass of snow particle in the saltation layer
434              ! (rho*qsalt*hsalt)
435              maxerosion=hsalt(j)*qsalt(j)*rhod/10.
436              fluxbs(j)=rhod*ws1(j)*cdragm(j)*(AcoefQBS(j)-qsalt(j)) &
437                       / (1.-rhod*ws1(j)*cdragm(j)*BcoefQBS(j)*dtime)
438              !fluxbs(j)= rhod*ws1(j)*cdragm(j)*(qbs1(j)-qsalt(j))
439              fluxbs(j)=max(-maxerosion,fluxbs(j))
[4529]440
[4672]441              ! for outputs
442
443              zxustartlic(i) = ustart(j)
444              zxrhoslic(i) = rhos(j)
[4523]445       enddo
446
447  endif
448
449
450
451!****************************************************************************************
[4672]452! Calculate snow amount
[4523]453!   
454!****************************************************************************************
455    IF (ok_bs) THEN
[4526]456      precip_totsnow(:)=precip_snow(:)+precip_bs(:)
457      evap_totsnow(:)=evap(:)-fluxbs(:) ! flux bs is positive towards the surface (snow erosion)
[4523]458    ELSE
[4526]459      precip_totsnow(:)=precip_snow(:)
460      evap_totsnow(:)=evap(:)
[4523]461    ENDIF
[4672]462   
463 
[4523]464    CALL fonte_neige(knon, is_lic, knindex, dtime, &
465         tsurf, precip_rain, precip_totsnow,  &
466         snow, qsol, tsurf_new, evap_totsnow)
[4672]467   
468   
[4523]469    WHERE (snow(1 : knon) .LT. 0.0001) agesno(1 : knon) = 0.                                         
470    zfra(1:knon) = MAX(0.0,MIN(1.0,snow(1:knon)/(snow(1:knon)+10.0))) 
471
472
[3792]473    END IF ! landice_opt
474
475
[781]476!****************************************************************************************
477! Send run-off on land-ice to coupler if coupled ocean.
[3903]478! run_off_lic has been calculated in fonte_neige or surf_inlandsis
[4283]479! If landice_opt>=2, corresponding call is done from surf_land_orchidee
[781]480!****************************************************************************************
[4283]481    IF (type_ocean=='couple' .AND. landice_opt .LT. 2) THEN
482       ! Compress fraction where run_off_lic is active (here all pctsrf(is_lic))
483       run_off_lic_frac(:)=0.0
484       DO j = 1, knon
485          i = knindex(j)
486          run_off_lic_frac(j) = pctsrf(i,is_lic)
487       ENDDO
488
489       CALL cpl_send_landice_fields(itime, knon, knindex, run_off_lic, run_off_lic_frac)
[781]490    ENDIF
[1865]491
492 ! transfer runoff rate [kg/m2/s](!) to physiq for output
493    runoff(1:knon)=run_off_lic(1:knon)/dtime
494
[1403]495       snow_o=0.
496       zfra_o = 0.
497       DO j = 1, knon
498           i = knindex(j)
499           snow_o(i) = snow(j)
500           zfra_o(i) = zfra(j)
501       ENDDO
502
503
[2227]504!albedo SB >>>
505     select case(NSW)
506     case(2)
507       alb_dir(1:knon,1)=alb1(1:knon)
508       alb_dir(1:knon,2)=alb2(1:knon)
509     case(4)
510       alb_dir(1:knon,1)=alb1(1:knon)
511       alb_dir(1:knon,2)=alb2(1:knon)
512       alb_dir(1:knon,3)=alb2(1:knon)
513       alb_dir(1:knon,4)=alb2(1:knon)
514     case(6)
515       alb_dir(1:knon,1)=alb1(1:knon)
516       alb_dir(1:knon,2)=alb1(1:knon)
517       alb_dir(1:knon,3)=alb1(1:knon)
518       alb_dir(1:knon,4)=alb2(1:knon)
519       alb_dir(1:knon,5)=alb2(1:knon)
520       alb_dir(1:knon,6)=alb2(1:knon)
[3900]521
[3901]522       IF ((landice_opt .EQ. 1) .AND. (iflag_albcalc .EQ. 2)) THEN
[3900]523       alb_dir(1:knon,1)=alb6(1:knon,1)
524       alb_dir(1:knon,2)=alb6(1:knon,2)
525       alb_dir(1:knon,3)=alb6(1:knon,3)
526       alb_dir(1:knon,4)=alb6(1:knon,4)
527       alb_dir(1:knon,5)=alb6(1:knon,5)
528       alb_dir(1:knon,6)=alb6(1:knon,6)
529       ENDIF
530
[2227]531     end select
532alb_dif=alb_dir
533!albedo SB <<<
534
535
[781]536  END SUBROUTINE surf_landice
537!
538!****************************************************************************************
539!
540END MODULE surf_landice_mod
541
542
543
Note: See TracBrowser for help on using the repository browser.