source: trunk/LMDZ.MARS/libf/phymars/vdifc_mod.F

Last change on this file was 3727, checked in by emillour, 3 months ago

Mars PCM:
Code tidying: put routines in modules, remove useless "return" statements and
remove obsolete and unused scopyi.F
EM

File size: 60.3 KB
RevLine 
[1969]1      MODULE vdifc_mod
2     
3      IMPLICIT NONE
4     
5      CONTAINS
6     
[3115]7      SUBROUTINE vdifc(ngrid,nlay,nsoil,nq,nqsoil,ppopsk,
[3468]8     $                ptimestep,pcapcal,
[38]9     $                pplay,pplev,pzlay,pzlev,pz0,
[3115]10     $                pu,pv,ph,pq,ptsrf,ptsoil,pemis,pqsurf,qsoil,
[3230]11     $                pore_icefraction,pdufi,pdvfi,pdhfi,
12     $                pdqfi,pfluxsrf,
[38]13     $                pdudif,pdvdif,pdhdif,pdtsrf,pq2,
[3203]14     $                pdqdif,pdqsdif,wstar,
[1996]15     $                hfmax,pcondicea_co2microp,sensibFlux,
[3115]16     $                dustliftday,local_time,watercap, dwatercap_dif)
17
[1036]18      use tracer_mod, only: noms, igcm_dust_mass, igcm_dust_number,
19     &                      igcm_dust_submicron, igcm_h2o_vap,
[2826]20     &                      igcm_h2o_ice, alpha_lift, igcm_co2,
[2312]21     &                      igcm_hdo_vap, igcm_hdo_ice,
[1974]22     &                      igcm_stormdust_mass, igcm_stormdust_number
[3111]23      use surfdat_h, only: watercaptag, frost_albedo_threshold, dryness,
24     &                     old_wsublimation_scheme
[2953]25      USE comcstfi_h, ONLY: cpp, r, rcp, g, pi
[1996]26      use watersat_mod, only: watersat
[1242]27      use turb_mod, only: turb_resolved, ustar, tstar
[2160]28      use compute_dtau_mod, only: ti_injection_sol,tf_injection_sol
[2312]29      use hdo_surfex_mod, only: hdo_surfex
[2515]30c      use geometry_mod, only: longitude_deg,latitude_deg !  Joseph
[2409]31      use dust_param_mod, only: doubleq, submicron, lifting
[2932]32      use write_output_mod, only: write_output
[2953]33      use comslope_mod, ONLY: nslope,def_slope,def_slope_mean,
34     &                      subslope_dist,major_slope,iflat
[3008]35      use microphys_h, only: To
[3619]36      use paleoclimate_mod, only: d_coef,h2o_ice_depth,lag_layer,
37     &                            zdqsdif_ssi_tot
[3262]38      use comsoil_h, only: layer, mlayer,adsorption_soil
[3150]39      use vdif_cd_mod, only: vdif_cd
[3727]40      use vdif_kc_mod, only: vdif_kc
41      use yamada4_mod, only: yamada4
[3167]42      use lmdz_call_atke, only: call_atke
[3292]43      use dust_windstress_lift_mod, only: dust_windstress_lift
[3726]44      use callkeys_mod, only: callcond, callrichsl, callyamada4,
45     &                        callatke, tke_heat_flux, water, hdo,
46     &                        co2clouds, rdstorm, dustinjection,
47     &                        latentheat_surfwater
[38]48      IMPLICIT NONE
49
50c=======================================================================
51c
52c   subject:
53c   --------
54c   Turbulent diffusion (mixing) for potential T, U, V and tracer
55c
56c   Shema implicite
57c   On commence par rajouter au variables x la tendance physique
58c   et on resoult en fait:
59c      x(t+1) =  x(t) + dt * (dx/dt)phys(t)  +  dt * (dx/dt)difv(t+1)
60c
61c   author:
62c   ------
63c      Hourdin/Forget/Fournier
64c=======================================================================
65
66c-----------------------------------------------------------------------
67c   declarations:
68c   -------------
69
70c
71c   arguments:
72c   ----------
73
[3115]74      INTEGER,INTENT(IN) :: ngrid,nlay,nsoil,nqsoil
[1036]75      REAL,INTENT(IN) :: ptimestep
76      REAL,INTENT(IN) :: pplay(ngrid,nlay),pplev(ngrid,nlay+1)
77      REAL,INTENT(IN) :: pzlay(ngrid,nlay),pzlev(ngrid,nlay+1)
78      REAL,INTENT(IN) :: pu(ngrid,nlay),pv(ngrid,nlay)
79      REAL,INTENT(IN) :: ph(ngrid,nlay)
[2953]80      REAL,INTENT(IN) :: ptsrf(ngrid,nslope),pemis(ngrid,nslope)
[1036]81      REAL,INTENT(IN) :: pdufi(ngrid,nlay),pdvfi(ngrid,nlay)
82      REAL,INTENT(IN) :: pdhfi(ngrid,nlay)
[2953]83      REAL,INTENT(IN) :: pfluxsrf(ngrid,nslope)
[1036]84      REAL,INTENT(OUT) :: pdudif(ngrid,nlay),pdvdif(ngrid,nlay)
[2953]85      REAL,INTENT(OUT) :: pdtsrf(ngrid,nslope),pdhdif(ngrid,nlay)
86      REAL,INTENT(IN) :: pcapcal(ngrid,nslope)
[1130]87      REAL,INTENT(INOUT) :: pq2(ngrid,nlay+1)
[3115]88      REAL,INTENT(IN) :: ptsoil(ngrid,nsoil,nslope)
[38]89
90c    Argument added for condensation:
[2826]91      REAL,INTENT(IN) :: ppopsk(ngrid,nlay)
[1996]92      REAL,INTENT(IN) :: pcondicea_co2microp(ngrid,nlay)! tendency due to CO2 condensation (kg/kg.s-1)
93     
[1036]94      REAL,INTENT(IN) :: pz0(ngrid) ! surface roughness length (m)
[224]95
[256]96c    Argument added to account for subgrid gustiness :
97
[1944]98      REAL,INTENT(IN) :: wstar(ngrid), hfmax(ngrid)!, zi(ngrid)
[256]99
[38]100c    Traceurs :
[1036]101      integer,intent(in) :: nq
[2953]102      REAL,INTENT(IN) :: pqsurf(ngrid,nq,nslope)
[2515]103      REAL :: zqsurf(ngrid) ! temporary water tracer
[1036]104      real,intent(in) :: pq(ngrid,nlay,nq), pdqfi(ngrid,nlay,nq)
105      real,intent(out) :: pdqdif(ngrid,nlay,nq)
[2953]106      real,intent(out) :: pdqsdif(ngrid,nq,nslope)
[1974]107      REAL,INTENT(in) :: dustliftday(ngrid)
[3115]108      REAL,INTENT(inout) :: qsoil(ngrid,nsoil,nqsoil,nslope) !subsurface tracers
[1974]109      REAL,INTENT(in) :: local_time(ngrid)
[38]110     
111c   local:
112c   ------
113
[1130]114      REAL :: pt(ngrid,nlay)
[473]115 
[3098]116      INTEGER ilev,ig,ilay,nlev,islope,ik,lice
[38]117
[1047]118      REAL z4st,zdplanck(ngrid)
119      REAL zkv(ngrid,nlay+1),zkh(ngrid,nlay+1)
120      REAL zkq(ngrid,nlay+1)
[3165]121      REAL zcdv(ngrid,nslope),zcdh(ngrid,nslope)
[3203]122      REAL :: zcdv_true(ngrid,nslope)
123      REAL :: zcdh_true(ngrid,nslope)         ! drag coeff are used by the LES to recompute u* and hfx
124      REAL :: zcdv_tmp(ngrid),zcdh_tmp(ngrid) ! drag coeffs for the major sub-grid surface
[3165]125      REAL :: zcdv_true_tmp(ngrid),zcdh_true_tmp(ngrid)    ! drag coeffs (computed with wind gustiness for the major sub-grid surface
[1047]126      REAL zu(ngrid,nlay),zv(ngrid,nlay)
127      REAL zh(ngrid,nlay)
128      REAL ztsrf2(ngrid)
129      REAL z1(ngrid),z2(ngrid)
130      REAL za(ngrid,nlay),zb(ngrid,nlay)
131      REAL zb0(ngrid,nlay)
132      REAL zc(ngrid,nlay),zd(ngrid,nlay)
[38]133      REAL zcst1
[1047]134      REAL zu2(ngrid)
[3253]135 
[38]136
137      EXTERNAL SSUM,SCOPY
138      REAL SSUM
[1036]139      LOGICAL,SAVE :: firstcall=.true.
[38]140
[2616]141!$OMP THREADPRIVATE(firstcall)
[626]142
[38]143c     variable added for CO2 condensation:
144c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[1047]145      REAL hh , zhcond(ngrid,nlay)
146      REAL,PARAMETER :: latcond=5.9e5
147      REAL,PARAMETER :: tcond1mb=136.27
148      REAL,SAVE :: acond,bcond
[2616]149
150!$OMP THREADPRIVATE(acond,bcond)
[669]151     
[2515]152c     Subtimestep & implicit treatment of water vapor
153c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[3115]154      REAL zdqsdif_surf(ngrid) ! subtimestep pdqsdif for water ice
[2515]155      REAL ztsrf(ngrid) ! temporary surface temperature in tsub
[2953]156      REAL zdtsrf(ngrid,nslope) ! surface temperature tendancy in tsub
157      REAL surf_h2o_lh(ngrid,nslope) ! Surface h2o latent heat flux
158      REAL zsurf_h2o_lh(ngrid,nslope) ! Tsub surface h2o latent heat flux
[2515]159
[2179]160c     For latent heat release from ground water ice sublimation   
[2515]161c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
162      REAL tsrf_lh(ngrid) ! temporary surface temperature with lh effect
[2179]163      REAL lh ! latent heat, formulation given in the Technical Document:
164              ! "Modeling water ice sublimation under Phoenix-like conditions", Montmessin et al. 2004 
[38]165
166c    Tracers :
167c    ~~~~~~~
168      INTEGER iq
[1047]169      REAL zq(ngrid,nlay,nq)
170      REAL zq1temp(ngrid)
171      REAL rho(ngrid) ! near surface air density
172      REAL qsat(ngrid)
[3253]173 
[38]174
[2953]175      REAL hdoflux(ngrid,nslope)       ! value of vapour flux of HDO
176      REAL hdoflux_meshavg(ngrid)       ! value of vapour flux of HDO
[2934]177!      REAL h2oflux(ngrid)       ! value of vapour flux of H2O
[2312]178      REAL old_h2o_vap(ngrid)   ! traceur d'eau avant traitement
[2953]179      REAL saved_h2o_vap(ngrid)   ! traceur d'eau avant traitement
[2312]180
[38]181      REAL kmixmin
182
[2260]183c    Argument added for surface water ice budget:
[2953]184      REAL,INTENT(IN) :: watercap(ngrid,nslope)
185      REAL,INTENT(OUT) :: dwatercap_dif(ngrid,nslope)
[2260]186
[2515]187c    Subtimestep to compute h2o latent heat flux:
188      REAL :: dtmax = 0.5 ! subtimestep temp criterion
189      INTEGER tsub ! adaptative subtimestep (seconds)
190      REAL subtimestep !ptimestep/nsubtimestep
191      INTEGER nsubtimestep(ngrid) ! number of subtimestep (int)
192
[473]193c    Mass-variation scheme :
194c    ~~~~~~~
195
196      INTEGER j,l
[1047]197      REAL zcondicea(ngrid,nlay)
198      REAL zt(ngrid,nlay),ztcond(ngrid,nlay+1)
199      REAL betam(ngrid,nlay),dmice(ngrid,nlay)
200      REAL pdtc(ngrid,nlay)
201      REAL zhs(ngrid,nlay)
202      REAL,SAVE :: ccond
[473]203
[2616]204!$OMP THREADPRIVATE(ccond)
205
[473]206c     Theta_m formulation for mass-variation scheme :
207c     ~~~~~~~
208
[1047]209      INTEGER,SAVE :: ico2
210      INTEGER llnt(ngrid)
211      REAL,SAVE :: m_co2, m_noco2, A , B
212      REAL vmr_co2(ngrid,nlay)
[3153]213      REAL qco2,mmean(ngrid,nlay)
[473]214
[2616]215!$OMP THREADPRIVATE(ico2,m_co2,m_noco2,A,B)
216
[1047]217      REAL,INTENT(OUT) :: sensibFlux(ngrid)
[660]218
[2312]219!!MARGAUX
220      REAL DoH_vap(ngrid,nlay)
[2953]221!! Sub-grid scale slopes
222      REAL :: pdqsdif_tmp(ngrid,nq) ! Temporary for dust lifting
223      REAL :: watercap_tmp(ngrid)
224      REAL :: zq_slope_vap(ngrid,nlay,nq,nslope)
225      REAL :: zq_tmp_vap(ngrid,nlay,nq)
226      REAL :: ptsrf_tmp(ngrid)
227      REAL :: pqsurf_tmp(ngrid,nq)
228      REAL :: pdqsdif_tmphdo(ngrid,nq)
229      REAL :: qsat_tmp(ngrid)
230      INTEGER :: indmax
[2312]231
[2953]232      character*2 str2
233
[3115]234!! Subsurface exchanges
235      LOGICAL :: exchange             ! boolean to check if exchange between the subsurface and the atmosphere can occurs
236      REAL :: zdqsdif_regolith(ngrid,nslope) ! Flux from subsurface (positive pointing outwards) (kg/m^2/s)
237      REAL zq1temp_regolith(ngrid)    ! Temporary atmospheric mixing ratio after exchange with subsurface (kg / kg)
238      REAL zdqsdif_tot(ngrid)         ! subtimestep pdqsdif for water ice
239      LOGICAL :: writeoutput          ! boolean to say to soilexchange.F if we are at the last iteration and thus if he  can write in the diagsoil
[3230]240      REAL, INTENT(out) :: pore_icefraction(ngrid,nsoil,nslope) ! ice filling fraction in the pores
[3115]241
[3253]242!! Subsurface ice interactions
243      REAL Tice(ngrid,nslope)                  ! subsurface temperature where ice is located (K)
244      REAL qsat_ssi(ngrid,nslope)              ! qsat(Tice) (kg/kg)
245      REAL resist(ngrid,nslope)                ! subsurface ice flux reduction coef (1)
246      REAL zdqsdif_ssi_atm(ngrid,nslope)       ! SSI - atmosphere flux (kg/m^2/s^-1) at a given subtimestep
247      REAL zdqsdif_ssi_frost(ngrid,nslope)     ! SSI - frost flux (kg/m^2/s^-1) at a given subtimestep
248      REAL zdqsdif_ssi_atm_tot(ngrid,nslope)   ! SSI - atmosphere flux (kg/m^2/s^-1) summed over all subtimestep
249      REAL zdqsdif_ssi_frost_tot(ngrid,nslope) ! SSI - frost flux (kg/m^2/s^-1) summed over all subtimestep
250     
[3285]251      REAL :: tol_frost = 2e-3 ! tolerence for frost thicnkess (kg/m^2) to avoid numerical noise effect
[38]252c    ** un petit test de coherence
253c       --------------------------
254
[1779]255      ! AS: OK firstcall absolute
[38]256      IF (firstcall) THEN
257c        To compute: Tcond= 1./(bcond-acond*log(.0095*p)) (p in pascal)
258         bcond=1./tcond1mb
259         acond=r/latcond
[473]260         ccond=cpp/(g*latcond)
[38]261         PRINT*,'In vdifc: Tcond(P=1mb)=',tcond1mb,' Lcond=',latcond
[473]262         PRINT*,'          acond,bcond,ccond',acond,bcond,ccond
[38]263
[473]264
265         ico2=0
266
267c          Prepare Special treatment if one of the tracer is CO2 gas
[1036]268           do iq=1,nq
[473]269             if (noms(iq).eq."co2") then
270                ico2=iq
271                m_co2 = 44.01E-3  ! CO2 molecular mass (kg/mol)
272                m_noco2 = 33.37E-3  ! Non condensible mol mass (kg/mol)
273c               Compute A and B coefficient use to compute
274c               mean molecular mass Mair defined by
275c               1/Mair = q(ico2)/m_co2 + (1-q(ico2))/m_noco2
276c               1/Mair = A*q(ico2) + B
277                A =(1/m_co2 - 1/m_noco2)
278                B=1/m_noco2
279             endif
280           enddo
281
[38]282        firstcall=.false.
283      ENDIF
284
[2953]285      DO ig = 1,ngrid
286       ptsrf_tmp(ig) = ptsrf(ig,major_slope(ig))
287       pqsurf_tmp(ig,:) = pqsurf(ig,:,major_slope(ig))
288      ENDDO
[38]289
290c-----------------------------------------------------------------------
291c    1. initialisation
292c    -----------------
293
294      nlev=nlay+1
295
[1035]296      ! initialize output tendencies to zero:
297      pdudif(1:ngrid,1:nlay)=0
298      pdvdif(1:ngrid,1:nlay)=0
299      pdhdif(1:ngrid,1:nlay)=0
[2953]300      pdtsrf(1:ngrid,1:nslope)=0
301      zdtsrf(1:ngrid,1:nslope)=0
302      surf_h2o_lh(1:ngrid,1:nslope)=0
303      zsurf_h2o_lh(1:ngrid,1:nslope)=0
[1035]304      pdqdif(1:ngrid,1:nlay,1:nq)=0
[2953]305      pdqsdif(1:ngrid,1:nq,1:nslope)=0
306      pdqsdif_tmp(1:ngrid,1:nq)=0
[3115]307      zdqsdif_surf(1:ngrid)=0
[2953]308      dwatercap_dif(1:ngrid,1:nslope)=0
[3115]309      zdqsdif_regolith(1:ngrid,1:nslope)=0
310      zq1temp_regolith(1:ngrid)=0
311      zdqsdif_tot(1:ngrid)=0
[3230]312      pore_icefraction(:,:,:) = 0.
[3253]313      zdqsdif_ssi_atm(:,:) = 0.
314      zdqsdif_ssi_frost(:,:) = 0.
315      zdqsdif_ssi_tot(:,:) = 0.
316      zdqsdif_ssi_atm_tot(:,:) = 0.
317      zdqsdif_ssi_frost_tot(:,:) = 0.
318     
[3262]319           
[38]320c    ** calcul de rho*dz et dt*rho/dz=dt*rho**2 g/dp
321c       avec rho=p/RT=p/ (R Theta) (p/ps)**kappa
322c       ----------------------------------------
323
324      DO ilay=1,nlay
325         DO ig=1,ngrid
326            za(ig,ilay)=(pplev(ig,ilay)-pplev(ig,ilay+1))/g
[473]327! Mass variation scheme:
328            betam(ig,ilay)=-za(ig,ilay)*latcond/(cpp*ppopsk(ig,ilay))
[38]329         ENDDO
330      ENDDO
331
332      zcst1=4.*g*ptimestep/(r*r)
333      DO ilev=2,nlev-1
334         DO ig=1,ngrid
335            zb0(ig,ilev)=pplev(ig,ilev)*
336     s      (pplev(ig,1)/pplev(ig,ilev))**rcp /
337     s      (ph(ig,ilev-1)+ph(ig,ilev))
338            zb0(ig,ilev)=zcst1*zb0(ig,ilev)*zb0(ig,ilev)/
339     s      (pplay(ig,ilev-1)-pplay(ig,ilev))
340         ENDDO
341      ENDDO
342      DO ig=1,ngrid
[2953]343         zb0(ig,1)=ptimestep*pplev(ig,1)/(r*ptsrf_tmp(ig))
[38]344      ENDDO
345
[473]346c     -----------------------------------
[38]347c     Potential Condensation temperature:
348c     -----------------------------------
349
[473]350c     Compute CO2 Volume mixing ratio
351c     -------------------------------
352      if (callcond.and.(ico2.ne.0)) then
353         DO ilev=1,nlay
354            DO ig=1,ngrid
[529]355              qco2=MAX(1.E-30
356     &             ,pq(ig,ilev,ico2)+pdqfi(ig,ilev,ico2)*ptimestep)
[473]357c             Mean air molecular mass = 1/(q(ico2)/m_co2 + (1-q(ico2))/m_noco2)
[3153]358              mmean(ig,ilev)=1/(A*qco2 +B)
359              vmr_co2(ig,ilev) = qco2*mmean(ig,ilev)/m_co2
[473]360            ENDDO
361         ENDDO
362      else
363         DO ilev=1,nlay
364            DO ig=1,ngrid
365              vmr_co2(ig,ilev)=0.95
366            ENDDO
367         ENDDO
368      end if
[38]369
[473]370c     forecast of atmospheric temperature zt and frost temperature ztcond
371c     --------------------------------------------------------------------
[38]372
[473]373      if (callcond) then
374        DO ilev=1,nlay
375          DO ig=1,ngrid
[884]376              ztcond(ig,ilev)=
377     &      1./(bcond-acond*log(.01*vmr_co2(ig,ilev)*pplay(ig,ilev)))
378            if (pplay(ig,ilev).lt.1e-4) ztcond(ig,ilev)=0.0 !mars Monica
[473]379!            zhcond(ig,ilev) =
380!     &  (1./(bcond-acond*log(.0095*pplay(ig,ilev))))/ppopsk(ig,ilev)
381              zhcond(ig,ilev) = ztcond(ig,ilev)/ppopsk(ig,ilev)
382          END DO
383        END DO
[884]384        ztcond(:,nlay+1)=ztcond(:,nlay)
[473]385      else
[884]386         zhcond(:,:) = 0
387         ztcond(:,:) = 0
[473]388      end if
389
390
[38]391c-----------------------------------------------------------------------
392c   2. ajout des tendances physiques
393c      -----------------------------
394
395      DO ilev=1,nlay
396         DO ig=1,ngrid
397            zu(ig,ilev)=pu(ig,ilev)+pdufi(ig,ilev)*ptimestep
398            zv(ig,ilev)=pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep
399            zh(ig,ilev)=ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep
[473]400!            zh(ig,ilev)=max(zh(ig,ilev),zhcond(ig,ilev))
[38]401         ENDDO
402      ENDDO
[2953]403
[2823]404      zq(1:ngrid,1:nlay,1:nq)=pq(1:ngrid,1:nlay,1:nq)+
405     &                        pdqfi(1:ngrid,1:nlay,1:nq)*ptimestep
[38]406
407c-----------------------------------------------------------------------
408c   3. schema de turbulence
409c      --------------------
410
411c    ** source d'energie cinetique turbulente a la surface
412c       (condition aux limites du schema de diffusion turbulente
413c       dans la couche limite
414c       ---------------------
415
[3165]416      CALL vdif_cd(ngrid,nlay,nslope,pz0,g,pzlay,pplay,pu,pv,wstar,
[3325]417     &          ptsrf,ph,mmean(:,1),zq(:,:,igcm_h2o_vap),
418     &          pqsurf(:,igcm_h2o_ice,:), .false.,
[3153]419     &          zcdv_true,zcdh_true)
[256]420
[3165]421      zu2(:)=pu(:,1)*pu(:,1)+pv(:,1)*pv(:,1)
[256]422
[3165]423      DO islope = 1,nslope
[291]424        IF (callrichsl) THEN
[3165]425          zcdv(:,islope)=zcdv_true(:,islope)*sqrt(zu2(:)+
[545]426     &     (log(1.+0.7*wstar(:) + 2.3*wstar(:)**2))**2)
[3165]427          zcdh(:,islope)=zcdh_true(:,islope)*sqrt(zu2(:)+
[545]428     &     (log(1.+0.7*wstar(:) + 2.3*wstar(:)**2))**2)
[284]429        ELSE
[3165]430           zcdv(:,islope)=zcdv_true(:,islope)*sqrt(zu2(:))     ! 1 / bulk aerodynamic momentum conductance
431           zcdh(:,islope)=zcdh_true(:,islope)*sqrt(zu2(:))     ! 1 / bulk aerodynamic heat conductance
[290]432        ENDIF
[3165]433      ENDDO
434      ustar(:) = 0
435      tstar(:) = 0
436      DO ig = 1,ngrid
437          zcdv_tmp(ig) = zcdv(ig,major_slope(ig))
438          zcdh_tmp(ig) = zcdh(ig,major_slope(ig))
439          zcdv_true_tmp(ig) = zcdv_true(ig,major_slope(ig))
440          zcdh_true_tmp(ig) = zcdh_true(ig,major_slope(ig))
441          IF (callrichsl) THEN
442             ustar(ig)=sqrt(zcdv_true(ig,major_slope(ig)))
443     &                *sqrt(zu2(ig)+(log(1.+0.7*wstar(ig) +
444     &                 2.3*wstar(ig)**2))**2)
445             IF (zcdh_true(ig,major_slope(ig)) .ne. 0.) THEN      ! When Cd=Ch=0, u*=t*=0
446                tstar(ig)=(ph(ig,1)-ptsrf_tmp(ig))
447     &                     *zcdh_tmp(ig)/ustar(ig)
448             ENDIF
449          ELSE
450                ustar(ig)=sqrt(zcdv_true(ig,major_slope(ig)))
451     &                    *sqrt(zu2(ig))
452                tstar(ig)=(ph(ig,1)-ptsrf_tmp(ig))
453     &                     *zcdh_true(ig,major_slope(ig))
454     &                     /sqrt(zcdv_true(ig,major_slope(ig)))
455          ENDIF
456      ENDDO
[38]457
458c    ** schema de diffusion turbulente dans la couche limite
459c       ----------------------------------------------------
[555]460      pt(:,:)=ph(:,:)*ppopsk(:,:)
[3167]461      if (callyamada4) then
462         call yamada4(ngrid,nlay,nq,ptimestep,g,r,pplev,pt
[3165]463     s   ,pzlev,pzlay,pu,pv,ph,pq,zcdv_true_tmp,pq2,zkv,zkh,zkq,ustar
[3167]464     s   ,9)     
465     
466      elseif (callatke) then
467         call call_atke(ptimestep,ngrid,nlay,zcdv_true_tmp,
468     s               zcdh_true_tmp,pu(:,1),pv(:,1),ptsrf_tmp,
469     s               pu,pv,pt,zq(:,1,igcm_h2o_vap),pplay,pplev,
470     s               pzlay,pzlev,pq2,zkv(:,1:nlay),zkh(:,1:nlay))
[544]471
[3167]472         zkv(:,nlay+1) = zkv(:,nlay)
473         zkh(:,nlay+1) = zkh(:,nlay)   
474      else
475        call vdif_kc(ngrid,nlay,nq,ptimestep,g,pzlev,pzlay
476     s              ,pu,pv,ph,zcdv_true_tmp
477     s              ,pq2,zkv,zkh,zq)
478     
479      endif
[38]480      if ((doubleq).and.(ngrid.eq.1)) then
481        kmixmin = 80. !80.! minimum eddy mix coeff in 1D
[3167]482        do ilev=2,nlay
[38]483          do ig=1,ngrid
484           zkh(ig,ilev) = max(kmixmin,zkh(ig,ilev))
485           zkv(ig,ilev) = max(kmixmin,zkv(ig,ilev))
486          end do
487        end do
488      end if
489
490c-----------------------------------------------------------------------
491c   4. inversion pour l'implicite sur u
492c      --------------------------------
493
494c    ** l'equation est
495c       u(t+1) =  u(t) + dt * {(du/dt)phys}(t)  +  dt * {(du/dt)difv}(t+1)
496c       avec
497c       /zu/ = u(t) + dt * {(du/dt)phys}(t)   (voir paragraphe 2.)
498c       et
499c       dt * {(du/dt)difv}(t+1) = dt * {(d/dz)[ Ku (du/dz) ]}(t+1)
500c       donc les entrees sont /zcdv/ pour la condition a la limite sol
501c       et /zkv/ = Ku
[2312]502
[2274]503      zb(1:ngrid,2:nlay)=zkv(1:ngrid,2:nlay)*zb0(1:ngrid,2:nlay)
[3165]504      zb(1:ngrid,1)=zcdv_tmp(1:ngrid)*zb0(1:ngrid,1)
[38]505
506      DO ig=1,ngrid
507         z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
508         zc(ig,nlay)=za(ig,nlay)*zu(ig,nlay)*z1(ig)
509         zd(ig,nlay)=zb(ig,nlay)*z1(ig)
510      ENDDO
511
512      DO ilay=nlay-1,1,-1
513         DO ig=1,ngrid
514            z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
515     $         zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
516            zc(ig,ilay)=(za(ig,ilay)*zu(ig,ilay)+
517     $         zb(ig,ilay+1)*zc(ig,ilay+1))*z1(ig)
518            zd(ig,ilay)=zb(ig,ilay)*z1(ig)
519         ENDDO
520      ENDDO
521
522      DO ig=1,ngrid
523         zu(ig,1)=zc(ig,1)
524      ENDDO
525      DO ilay=2,nlay
526         DO ig=1,ngrid
527            zu(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zu(ig,ilay-1)
528         ENDDO
529      ENDDO
530
531c-----------------------------------------------------------------------
532c   5. inversion pour l'implicite sur v
533c      --------------------------------
534
535c    ** l'equation est
536c       v(t+1) =  v(t) + dt * {(dv/dt)phys}(t)  +  dt * {(dv/dt)difv}(t+1)
537c       avec
538c       /zv/ = v(t) + dt * {(dv/dt)phys}(t)   (voir paragraphe 2.)
539c       et
540c       dt * {(dv/dt)difv}(t+1) = dt * {(d/dz)[ Kv (dv/dz) ]}(t+1)
541c       donc les entrees sont /zcdv/ pour la condition a la limite sol
542c       et /zkv/ = Kv
543
544      DO ig=1,ngrid
545         z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
546         zc(ig,nlay)=za(ig,nlay)*zv(ig,nlay)*z1(ig)
547         zd(ig,nlay)=zb(ig,nlay)*z1(ig)
548      ENDDO
549
550      DO ilay=nlay-1,1,-1
551         DO ig=1,ngrid
552            z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
553     $         zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
554            zc(ig,ilay)=(za(ig,ilay)*zv(ig,ilay)+
555     $         zb(ig,ilay+1)*zc(ig,ilay+1))*z1(ig)
556            zd(ig,ilay)=zb(ig,ilay)*z1(ig)
557         ENDDO
558      ENDDO
559
560      DO ig=1,ngrid
561         zv(ig,1)=zc(ig,1)
562      ENDDO
563      DO ilay=2,nlay
564         DO ig=1,ngrid
565            zv(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zv(ig,ilay-1)
566         ENDDO
567      ENDDO
568
569c-----------------------------------------------------------------------
[3162]570c    Using the wind modified by friction for lifting and  sublimation
571c     ----------------------------------------------------------------
572
573!     This is computed above and takes into account surface-atmosphere flux
574!     enhancement by subgrid gustiness and atmospheric-stability related
575!     variations of transfer coefficients.
576!     Calculate Cd again with wind slowed by friction
577c -------------------------------------------
578
[3165]579      CALL vdif_cd(ngrid,nlay,nslope,pz0,g,pzlay,pplay,zu,zv,wstar,
[3325]580     &          ptsrf,ph,mmean(:,1),zq(:,:,igcm_h2o_vap),
581     &          pqsurf(:,igcm_h2o_ice,:), .true.,
[3163]582     &          zcdv_true,zcdh_true)
[3162]583
[3165]584      zu2(:)=zu(:,1)*zu(:,1)+zv(:,1)*zv(:,1)
585
586      DO islope = 1,nslope
587        IF (callrichsl) THEN
588          zcdv(:,islope)=zcdv_true(:,islope)*sqrt(zu2(:)+
[3162]589     &     (log(1.+0.7*wstar(:) + 2.3*wstar(:)**2))**2)
[3165]590          zcdh(:,islope)=zcdh_true(:,islope)*sqrt(zu2(:)+
[3162]591     &     (log(1.+0.7*wstar(:) + 2.3*wstar(:)**2))**2)
592        ELSE
[3165]593           zcdv(:,islope)=zcdv_true(:,islope)*sqrt(zu2(:))     ! 1 / bulk aerodynamic momentum conductance
594           zcdh(:,islope)=zcdh_true(:,islope)*sqrt(zu2(:))     ! 1 / bulk aerodynamic heat conductance
595        ENDIF
596      ENDDO
597      ustar(:) = 0
598      tstar(:) = 0
599      DO ig = 1,ngrid
600          zcdv_tmp(ig) = zcdv(ig,major_slope(ig))
601          zcdh_tmp(ig) = zcdh(ig,major_slope(ig))
602          zcdv_true_tmp(ig) = zcdv_true(ig,major_slope(ig))
603          zcdh_true_tmp(ig) = zcdh_true(ig,major_slope(ig))
604          IF (callrichsl) THEN
605             ustar(ig)=sqrt(zcdv_true(ig,major_slope(ig)))
606     &                *sqrt(zu2(ig)+(log(1.+0.7*wstar(ig) +
607     &                 2.3*wstar(ig)**2))**2)
608             IF (zcdh_true(ig,major_slope(ig)) .ne. 0.) THEN      ! When Cd=Ch=0, u*=t*=0
609                tstar(ig)=(ph(ig,1)-ptsrf_tmp(ig))
610     &                     *zcdh_tmp(ig)/ustar(ig)
611             ENDIF
612          ELSE
613                ustar(ig)=sqrt(zcdv_true(ig,major_slope(ig)))
614     &                    *sqrt(zu2(ig))
615                tstar(ig)=(ph(ig,1)-ptsrf_tmp(ig))
616     &                     *zcdh_true(ig,major_slope(ig))
617     &                     /sqrt(zcdv_true(ig,major_slope(ig)))
618          ENDIF
619      ENDDO
[3162]620         
621
622c-----------------------------------------------------------------------
[38]623c   6. inversion pour l'implicite sur h sans oublier le couplage
624c      avec le sol (conduction)
625c      ------------------------
626
627c    ** l'equation est
628c       h(t+1) =  h(t) + dt * {(dh/dt)phys}(t)  +  dt * {(dh/dt)difv}(t+1)
629c       avec
630c       /zh/ = h(t) + dt * {(dh/dt)phys}(t)   (voir paragraphe 2.)
631c       et
632c       dt * {(dh/dt)difv}(t+1) = dt * {(d/dz)[ Kh (dh/dz) ]}(t+1)
633c       donc les entrees sont /zcdh/ pour la condition de raccord au sol
634c       et /zkh/ = Kh
635c       -------------
636
[473]637c Mass variation scheme:
[2274]638      zb(1:ngrid,2:nlay)=zkh(1:ngrid,2:nlay)*zb0(1:ngrid,2:nlay)
[3165]639      zb(1:ngrid,1)=zcdh_tmp(1:ngrid)*zb0(1:ngrid,1)
[38]640
[473]641c on initialise dm c
642     
643      pdtc(:,:)=0.
644      zt(:,:)=0.
645      dmice(:,:)=0.
[38]646
647c    ** calcul de (d Planck / dT) a la temperature d'interface
648c       ------------------------------------------------------
649
650      z4st=4.*5.67e-8*ptimestep
[544]651      IF (tke_heat_flux .eq. 0.) THEN
[38]652      DO ig=1,ngrid
[2953]653         indmax = major_slope(ig)
654         zdplanck(ig)=z4st*pemis(ig,indmax)*ptsrf(ig,indmax)*
655     &                ptsrf(ig,indmax)*ptsrf(ig,indmax)
[38]656      ENDDO
[544]657      ELSE
658         zdplanck(:)=0.
659      ENDIF
[38]660
[473]661! calcul de zc et zd pour la couche top en prenant en compte le terme
[2080]662! de variation de masse (on fait une boucle pour que \E7a converge)
[473]663
664! Identification des points de grilles qui ont besoin de la correction
665
666      llnt(:)=1
[1236]667      IF (.not.turb_resolved) THEN
[884]668      IF (callcond) THEN
669       DO ig=1,ngrid
[473]670         DO l=1,nlay
671            if(zh(ig,l) .lt. zhcond(ig,l)) then
672               llnt(ig)=300 
673! 200 and 100 do not go beyond month 9 with normal dissipation
674               goto 5
675            endif
676         ENDDO
[884]6775      continue
678       ENDDO
679      ENDIF
[473]680
[529]681      ENDIF
682
[473]683      DO ig=1,ngrid
[2953]684       indmax = major_slope(ig)
[473]685! Initialization of z1 and zd, which do not depend on dmice
686
687      z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
688      zd(ig,nlay)=zb(ig,nlay)*z1(ig)
689
690      DO ilay=nlay-1,1,-1
691          z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
692     $        zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
693          zd(ig,ilay)=zb(ig,ilay)*z1(ig)
694      ENDDO
695
696! Convergence loop
697
698      DO j=1,llnt(ig)
699
700            z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
701            zc(ig,nlay)=za(ig,nlay)*zh(ig,nlay)
702     &      -betam(ig,nlay)*dmice(ig,nlay)
703            zc(ig,nlay)=zc(ig,nlay)*z1(ig)
704!            zd(ig,nlay)=zb(ig,nlay)*z1(ig)
705
706! calcul de zc et zd pour les couches du haut vers le bas
707
708             DO ilay=nlay-1,1,-1
709               z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
710     $            zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
711               zc(ig,ilay)=(za(ig,ilay)*zh(ig,ilay)+
712     $            zb(ig,ilay+1)*zc(ig,ilay+1)-
713     $            betam(ig,ilay)*dmice(ig,ilay))*z1(ig)
714!               zd(ig,ilay)=zb(ig,ilay)*z1(ig)
715            ENDDO
716
[38]717c    ** calcul de la temperature_d'interface et de sa tendance.
718c       on ecrit que la somme des flux est nulle a l'interface
719c       a t + \delta t,
720c       c'est a dire le flux radiatif a {t + \delta t}
721c       + le flux turbulent a {t + \delta t}
722c            qui s'ecrit K (T1-Tsurf) avec T1 = d1 Tsurf + c1
723c            (notation K dt = /cpp*b/)       
724c       + le flux dans le sol a t
725c       + l'evolution du flux dans le sol lorsque la temperature d'interface
726c       passe de sa valeur a t a sa valeur a {t + \delta t}.
727c       ----------------------------------------------------
728
[2953]729         z1(ig)=pcapcal(ig,indmax)*ptsrf(ig,indmax)
730     s     + cpp*zb(ig,1)*zc(ig,1)
731     s     + zdplanck(ig)*ptsrf(ig,indmax)
732     s     + pfluxsrf(ig,indmax)*ptimestep
733         z2(ig)= pcapcal(ig,indmax)+cpp*zb(ig,1)*(1.-zd(ig,1))
734     s     +zdplanck(ig)
[38]735         ztsrf2(ig)=z1(ig)/z2(ig)
[473]736!         pdtsrf(ig)=(ztsrf2(ig)-ptsrf(ig))/ptimestep  !incremented outside loop
737            zhs(ig,1)=zc(ig,1)+zd(ig,1)*ztsrf2(ig)
[38]738
739c    ** et a partir de la temperature au sol on remonte
740c       -----------------------------------------------
741
[473]742            DO ilay=2,nlay
743               zhs(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zhs(ig,ilay-1)
744            ENDDO
745            DO ilay=1,nlay
746               zt(ig,ilay)=zhs(ig,ilay)*ppopsk(ig,ilay)
747            ENDDO
748
749c      Condensation/sublimation in the atmosphere
750c      ------------------------------------------
751c      (computation of zcondicea and dmice)
752
[1996]753      IF (.NOT. co2clouds) then
754       DO l=nlay , 1, -1
[473]755           IF(zt(ig,l).LT.ztcond(ig,l)) THEN
756               pdtc(ig,l)=(ztcond(ig,l) - zt(ig,l))/ptimestep
757               zcondicea(ig,l)=(pplev(ig,l)-pplev(ig,l+1))
758     &                        *ccond*pdtc(ig,l)
759              dmice(ig,l)= dmice(ig,l) + zcondicea(ig,l)*ptimestep
760            END IF
[1996]761       ENDDO
762      ELSE
763       DO l=nlay , 1, -1
764         zcondicea(ig,l)= 0.!pcondicea_co2microp(ig,l)*
765c     &                        (pplev(ig,l) - pplev(ig,l+1))/g
766         dmice(ig,l)= 0.!dmice(ig,l) + zcondicea(ig,l)*ptimestep
767         pdtc(ig,l)=0.
768       ENDDO   
769      ENDIF
[2953]770
[1996]771       ENDDO!of Do j=1,XXX
[2953]772      pdtsrf(ig,indmax)=(ztsrf2(ig)-ptsrf(ig,indmax))/ptimestep
[1996]773      ENDDO   !of Do ig=1,ngrid
[669]774     
[660]775      DO ig=1,ngrid  ! computing sensible heat flux (atm => surface)
776         sensibFlux(ig)=cpp*zb(ig,1)/ptimestep*(zhs(ig,1)-ztsrf2(ig))
777      ENDDO
[473]778
[2953]779c  Now implicit sheme on each sub-grid subslope:
780      IF (nslope.ne.1) then
781      DO islope=1,nslope
782        DO ig=1,ngrid
783          IF(islope.ne.major_slope(ig)) then
784            IF (tke_heat_flux .eq. 0.) THEN
785              zdplanck(ig)=z4st*pemis(ig,islope)*ptsrf(ig,islope)**3
786            ELSE
787              zdplanck(ig) = 0.
788            ENDIF
[3165]789             zb(ig,1)=zcdh(ig,islope)*zb0(ig,1)
790             z1(ig)=pcapcal(ig,islope)*ptsrf(ig,islope)
[2953]791     s       + cpp*zb(ig,1)*zc(ig,1)
792     s       + zdplanck(ig)*ptsrf(ig,islope)
793     s       + pfluxsrf(ig,islope)*ptimestep
794            z2(ig)= pcapcal(ig,islope)+cpp*zb(ig,1)*(1.-zd(ig,1))
795     s       +zdplanck(ig)
796            ztsrf2(ig)=z1(ig)/z2(ig)
797            pdtsrf(ig,islope)=(ztsrf2(ig)-ptsrf(ig,islope))/ptimestep
798          ENDIF ! islope != indmax
799        ENDDO ! ig
800       ENDDO !islope
801       ENDIF !nslope.ne.1
802
[38]803c-----------------------------------------------------------------------
804c   TRACERS
805c   -------
806c       Calcul du flux vertical au bas de la premiere couche (dust) :
807c       -----------------------------------------------------------
[1047]808        do ig=1,ngrid 
[38]809          rho(ig) = zb0(ig,1) /ptimestep
810c          zb(ig,1) = 0.
811        end do
812c       Dust lifting:
813        if (lifting) then
[310]814#ifndef MESOSCALE
[38]815           if (doubleq.AND.submicron) then
816             do ig=1,ngrid
[2826]817c              if(qsurf(ig,igcm_co2).lt.1) then
[2953]818                 pdqsdif_tmp(ig,igcm_dust_mass) =
[38]819     &             -alpha_lift(igcm_dust_mass) 
[2953]820                 pdqsdif_tmp(ig,igcm_dust_number) =
[38]821     &             -alpha_lift(igcm_dust_number) 
[2953]822                 pdqsdif_tmp(ig,igcm_dust_submicron) =
[38]823     &             -alpha_lift(igcm_dust_submicron)
824c              end if
825             end do
826           else if (doubleq) then
[1974]827             if (dustinjection.eq.0) then !injection scheme 0 (old)
828                                          !or 2 (injection in CL)
829              do ig=1,ngrid
[2953]830               if(pqsurf_tmp(ig,igcm_co2).lt.1) then ! pas de soulevement si glace CO2
831                 pdqsdif_tmp(ig,igcm_dust_mass) =
[38]832     &             -alpha_lift(igcm_dust_mass) 
[2953]833                 pdqsdif_tmp(ig,igcm_dust_number) =
[520]834     &             -alpha_lift(igcm_dust_number)
835               end if
[1974]836              end do
837             elseif(dustinjection.eq.1)then ! dust injection scheme = 1 injection from surface
838              do ig=1,ngrid
[2953]839               if(pqsurf_tmp(ig,igcm_co2).lt.1) then ! pas de soulevement si glace CO2
[2160]840                IF((ti_injection_sol.LE.local_time(ig)).and.
841     &                  (local_time(ig).LE.tf_injection_sol)) THEN
[1974]842                  if (rdstorm) then !Rocket dust storm scheme
[2953]843                        pdqsdif_tmp(ig,igcm_stormdust_mass) =
[1974]844     &                          -alpha_lift(igcm_stormdust_mass)
845     &                          *dustliftday(ig)
[2953]846                        pdqsdif_tmp(ig,igcm_stormdust_number) =
[1974]847     &                          -alpha_lift(igcm_stormdust_number)
848     &                          *dustliftday(ig)
[2953]849                        pdqsdif_tmp(ig,igcm_dust_mass)= 0.
850                        pdqsdif_tmp(ig,igcm_dust_number)= 0.
[1974]851                  else
[2953]852                        pdqsdif_tmp(ig,igcm_dust_mass)=
[1974]853     &                          -dustliftday(ig)*
854     &                          alpha_lift(igcm_dust_mass)               
[2953]855                        pdqsdif_tmp(ig,igcm_dust_number)=
[1974]856     &                          -dustliftday(ig)*
857     &                          alpha_lift(igcm_dust_number)
858                  endif
859                  if (submicron) then
[2953]860                        pdqsdif_tmp(ig,igcm_dust_submicron) = 0.
[1974]861                  endif ! if (submicron)
862                ELSE ! outside dust injection time frame
[2953]863                  pdqsdif_tmp(ig,igcm_dust_mass)= 0.
864                  pdqsdif_tmp(ig,igcm_dust_number)= 0.
[2080]865                  if (rdstorm) then     
[2953]866                        pdqsdif_tmp(ig,igcm_stormdust_mass)= 0.
867                        pdqsdif_tmp(ig,igcm_stormdust_number)= 0.
[2080]868                  end if
[1974]869                ENDIF
870               
[2826]871                end if ! of if(qsurf(ig,igcm_co2).lt.1)
[1974]872              end do
873             endif ! end if dustinjection
[38]874           else if (submicron) then
875             do ig=1,ngrid
[2953]876                 pdqsdif_tmp(ig,igcm_dust_submicron) =
[38]877     &             -alpha_lift(igcm_dust_submicron)
878             end do
879           else
[1236]880#endif
[3292]881            call dust_windstress_lift(ngrid,nlay,nq,rho,
882     &                    zcdh_true_tmp,zcdh_tmp,
[2953]883     &                    pqsurf_tmp(:,igcm_co2),pdqsdif_tmp)
[1236]884#ifndef MESOSCALE
[38]885           endif !doubleq.AND.submicron
[310]886#endif
[38]887        else
[2953]888           pdqsdif_tmp(1:ngrid,1:nq) = 0.
[38]889        end if
890
891c       OU calcul de la valeur de q a la surface (water)  :
892c       ----------------------------------------
893
894c      Inversion pour l'implicite sur q
[2515]895c      Cas des traceurs qui ne sont pas h2o_vap
896c      h2o_vap est traite plus loin avec un sous pas de temps
897c      hdo_vap est traite ensuite car dependant de h2o_vap
[38]898c       --------------------------------
[2515]899
900        do iq=1,nq  !for all tracers except water vapor
901          if ((.not. water).or.(.not. iq.eq.igcm_h2o_vap).or.
902     &       (.not. iq.eq.igcm_hdo_vap)) then
903
904
[2274]905          zb(1:ngrid,2:nlay)=zkh(1:ngrid,2:nlay)*zb0(1:ngrid,2:nlay)
[2515]906          zb(1:ngrid,1)=0
[38]907
[2515]908          DO ig=1,ngrid
909               z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
910               zc(ig,nlay)=za(ig,nlay)*zq(ig,nlay,iq)*z1(ig)
911               zd(ig,nlay)=zb(ig,nlay)*z1(ig)
912          ENDDO
913 
914          DO ilay=nlay-1,2,-1
915               DO ig=1,ngrid
916                z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
917     $           zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
918                zc(ig,ilay)=(za(ig,ilay)*zq(ig,ilay,iq)+
919     $           zb(ig,ilay+1)*zc(ig,ilay+1))*z1(ig)
920                zd(ig,ilay)=zb(ig,ilay)*z1(ig)
921               ENDDO
922          ENDDO
923
924            if ((iq.eq.igcm_h2o_ice)
925     $      .or. (hdo.and.(iq.eq.igcm_hdo_ice) )) then
926
927               DO ig=1,ngrid
928                   z1(ig)=1./(za(ig,1)+zb(ig,1)+
929     $              zb(ig,2)*(1.-zd(ig,2)))
930                   zc(ig,1)=(za(ig,1)*zq(ig,1,iq)+
931     $              zb(ig,2)*zc(ig,2)) *z1(ig)  !special case h2o_ice
932               ENDDO
933            else ! every other tracer
934               DO ig=1,ngrid
935                   z1(ig)=1./(za(ig,1)+zb(ig,1)+
936     $              zb(ig,2)*(1.-zd(ig,2)))
937                   zc(ig,1)=(za(ig,1)*zq(ig,1,iq)+
938     $              zb(ig,2)*zc(ig,2) +
[2953]939     $             (-pdqsdif_tmp(ig,iq)) *ptimestep) *z1(ig)  !tracer flux from surface
[2515]940               ENDDO
941            endif !((iq.eq.igcm_h2o_ice)
942c         Starting upward calculations for simple mixing of tracer (dust)
943          DO ig=1,ngrid
944             zq(ig,1,iq)=zc(ig,1)
945             DO ilay=2,nlay
946               zq(ig,ilay,iq)=zc(ig,ilay)+zd(ig,ilay)*zq(ig,ilay-1,iq)
947             ENDDO
948          ENDDO
[2953]949          DO islope = 1,nslope
950           DO ig = 1,ngrid
951            pdqsdif(ig,iq,islope) = pdqsdif_tmp(ig,iq)
952     &            * cos(pi*def_slope_mean(islope)/180.)
953            ENDDO
954          ENDDO
955
[2515]956          endif! ((.not. water).or.(.not. iq.eq.igcm_h2o_vap)) then
957        enddo ! of do iq=1,nq
958
959c --------- h2o_vap --------------------------------
960
[3253]961c Treatement of the water frost
962c We use a subtimestep to take into account the release of latent heat
[3098]963         
[2515]964        do iq=1,nq 
[2312]965          if ((water).and.(iq.eq.igcm_h2o_vap)) then
[2515]966
[2953]967          DO islope = 1,nslope
[2515]968           DO ig=1,ngrid
[3098]969
[2953]970             zqsurf(ig)=pqsurf(ig,igcm_h2o_ice,islope)/
971     &             cos(pi*def_slope_mean(islope)/180.)
972              watercap_tmp(ig) = watercap(ig,islope)/
973     &                       cos(pi*def_slope_mean(islope)/180.)
[2515]974           ENDDO ! ig=1,ngrid
[3253]975           
976c             Computation of the subtimestep
[3115]977              call make_tsub(ngrid,pdtsrf(:,islope),zqsurf,
[2515]978     &                    ptimestep,dtmax,watercaptag,
[3253]979     &                    nsubtimestep)
980c           Calculation for turbulent exchange (rho Cd,h U (qatm - qsat(Tsurf)) with the surface (for ice)
[2515]981c           initialization of ztsrf, which is surface temperature in
982c           the subtimestep.
[2953]983           saved_h2o_vap(:)= zq(:,1,igcm_h2o_vap)   
[2515]984           DO ig=1,ngrid
[3262]985c           nsubtimestep(ig)=1 !for debug
[3253]986           subtimestep = ptimestep/nsubtimestep(ig)
[2953]987            ztsrf(ig)=ptsrf(ig,islope)  !  +pdtsrf(ig)*subtimestep
988            zq_tmp_vap(ig,:,:) =zq(ig,:,:)
[3253]989c           Beginning of the subtimestep
[2515]990            DO tsub=1,nsubtimestep(ig)
[3253]991             if(tsub.eq.nsubtimestep(ig)) writeoutput = .true.
[2515]992             zb(1:ngrid,2:nlay)=zkh(1:ngrid,2:nlay)*zb0(1:ngrid,2:nlay)
993     &                     /float(nsubtimestep(ig))
[3111]994             if(old_wsublimation_scheme) then
[3165]995               zb(1:ngrid,1)=zcdv(1:ngrid,islope)*zb0(1:ngrid,1)
[2515]996     &                     /float(nsubtimestep(ig))
[3111]997             else
[3165]998               zb(1:ngrid,1)=zcdh(1:ngrid,islope)*zb0(1:ngrid,1)
[3111]999     &                     /float(nsubtimestep(ig))
1000             endif
[2274]1001             zb(1:ngrid,1)=dryness(1:ngrid)*zb(1:ngrid,1)
[2515]1002             
1003             z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
[2953]1004             zc(ig,nlay)=za(ig,nlay)*zq_tmp_vap(ig,nlay,iq)*z1(ig)
[2515]1005             zd(ig,nlay)=zb(ig,nlay)*z1(ig)
1006             DO ilay=nlay-1,2,-1
1007                 z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
1008     $            zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
[2953]1009                 zc(ig,ilay)=(za(ig,ilay)*zq_tmp_vap(ig,ilay,iq)+
[2515]1010     $            zb(ig,ilay+1)*zc(ig,ilay+1))*z1(ig)
1011                 zd(ig,ilay)=zb(ig,ilay)*z1(ig)
1012             ENDDO 
1013             z1(ig)=1./(za(ig,1)+zb(ig,1)+
1014     $          zb(ig,2)*(1.-zd(ig,2)))
[2953]1015             zc(ig,1)=(za(ig,1)*zq_tmp_vap(ig,1,iq)+
[2515]1016     $        zb(ig,2)*zc(ig,2)) * z1(ig)
[38]1017
[2531]1018             call watersat(1,ztsrf(ig),pplev(ig,1),qsat(ig))
[2953]1019             old_h2o_vap(ig)=zq_tmp_vap(ig,1,igcm_h2o_vap)
[2515]1020             zd(ig,1)=zb(ig,1)*z1(ig)
1021             zq1temp(ig)=zc(ig,1)+ zd(ig,1)*qsat(ig)
[3111]1022             if(old_wsublimation_scheme) then
[3165]1023                zdqsdif_surf(ig)=rho(ig)*dryness(ig)*zcdv(ig,islope)
[3111]1024     &                           *(zq1temp(ig)-qsat(ig))
1025             else
[3165]1026                zdqsdif_surf(ig)=rho(ig)*dryness(ig)*zcdh(ig,islope)
[3111]1027     &                           *(zq1temp(ig)-qsat(ig))
1028             endif
[2953]1029
[3165]1030             zdqsdif_tot(ig) = zdqsdif_surf(ig)
[3253]1031! --------------------------------------------------------------------------------------------------------------------------------             
1032! We consider here the possible interactions between the subsurface and the atmosphere in the case of the surface is free of frost
1033! when computing the complete adsorption/desorption model
1034! --------------------------------------------------------------------------------------------------------------------------------     
[2587]1035             if(.not.watercaptag(ig)) then
[3115]1036                if (((-(zdqsdif_surf(ig))*
1037     &             subtimestep).gt.zqsurf(ig))
1038     &             .and.(pqsurf(ig,igcm_co2,islope).eq.0.)) then
1039                      exchange = .true.
1040                 else
1041                      exchange = .false.
1042                 endif
1043             else
1044                  exchange = .false.
1045             endif
1046
[3253]1047! --------------------------------------------------------------------------------------------------------------------------------             
1048! If one consider adsorption, all the fluxes to/from the surface/subsurface/atmosphere are computed here
1049! --------------------------------------------------------------------------------------------------------------------------------   
[3115]1050
1051             if (adsorption_soil) then 
1052                 call soilwater(1,nlay,nq,nsoil, nqsoil,
1053     &                     ztsrf(ig),ptsoil(ig,:,islope),subtimestep,
1054     &                     exchange,qsat(ig),zq_tmp_vap(ig,:,:),
1055     &                     za(ig,:),zb(ig,:),zc(ig,:),zd(ig,:),
1056     &                     zdqsdif_surf(ig), zqsurf(ig),
[3262]1057     &                     qsoil(ig,:,:,islope), pplev(ig,1), rho(ig),
[3115]1058     &                     writeoutput,zdqsdif_regolith(ig,islope),
[3230]1059     &                     zq1temp_regolith(ig),
[3262]1060     &                     pore_icefraction(ig,:,islope))
[3115]1061
[3262]1062
[3115]1063                 if(.not.watercaptag(ig)) then
1064                     if (exchange) then
1065                         zq1temp(ig) = zq1temp_regolith(ig)
1066                         zdqsdif_tot(ig)=
1067     &                        -zqsurf(ig)/subtimestep
1068                     else
1069                         zdqsdif_tot(ig) = zdqsdif_surf(ig) +
1070     &                           zdqsdif_regolith(ig,islope) ! boundary condition = qsat, but pdqsdif is calculated to update qsurf (including loss of surface ice to the subsurface)
1071                     endif  ! of "if exchange = true"
1072                 endif  ! of "if not.watercaptag"
[3248]1073              endif ! adsorption         
[3253]1074! --------------------------------------------------------------------------------------------------------------------------------------             
1075! Here we do the same, but without computing the complete adsorpption/desorption. Note that it work only if one does not use adsorption
1076! If no subsurface ice, then the models computes the surface flux/water vapor flux as usual
1077! --------------------------------------------------------------------------------------------------------------------------------------     
1078                     
1079! ------------------------------------------------------------------------------------------------------------------------------------------             
1080! First, We consider here the possible interactions between the subsurface and the atmosphere in the case of the surface is free of frost
1081! ------------------------------------------------------------------------------------------------------------------------------------------
[3115]1082
[3253]1083              if(.not.watercaptag(ig).and.(.not.adsorption_soil)) then
[3400]1084                 if ((((-zdqsdif_tot(ig)*subtimestep)
1085     &                .ge.(zqsurf(ig))))
1086     &              .or.
1087     &                (((-zdqsdif_tot(ig)*subtimestep)
1088     &                .lt.(zqsurf(ig)))
1089     &                .and. ((zqsurf(ig).lt.tol_frost))
1090     &                .and.lag_layer)) then
1091
[3404]1092                    if ((-zdqsdif_tot(ig)*subtimestep)
1093     &                .ge.(zqsurf(ig))) then
1094                          zdqsdif_tot(ig)=-zqsurf(ig)/subtimestep
1095                    endif
1096
[3253]1097                    if((h2o_ice_depth(ig,islope).gt.0)
[3400]1098     &                  .and. (zqsurf(ig) .lt. tol_frost)
[3253]1099     &                  .and.lag_layer) then
1100                  ! Atm <-> subsurface exchange, we need to update the exchange coefficient zb by a factor 1/1+R; R = zice*Cd,h*/D and add the flux from the subsurface
[3400]1101!                         zqsurf(ig) = 0
[3253]1102                         if(old_wsublimation_scheme) then
1103                            resist(ig,islope)=h2o_ice_depth(ig,islope)
1104     &                          *zcdv(ig,islope)/d_coef(ig,islope)
1105                         else
1106                            resist(ig,islope)=h2o_ice_depth(ig,islope)
1107     &                          *zcdh(ig,islope)/d_coef(ig,islope)                         
1108                         endif
1109                         
1110                         zb(ig,1)=zb(ig,1)*1/(1+resist(ig,islope)) ! change zb to account subsurface ice
1111                  ! Now we add the flux from the subsurface ice : rho Cd,h U*(1/1+R) * (q1-qsat_ssi(Tice))
1112                   
1113                        call compute_Tice(nsoil, ptsoil(ig,:,islope),
1114     &                                    ztsrf(ig),
1115     &                                    h2o_ice_depth(ig,islope),
1116     &                                    Tice(ig,islope))        ! compute ice temperature
[3098]1117   
[3253]1118                        call watersat(1,Tice(ig,islope),pplev(ig,1),
1119     &                                qsat_ssi(ig,islope))
1120                 
1121                         qsat_ssi(ig,islope)=ztsrf(ig)/Tice(ig,islope)
1122     &                                       *qsat_ssi(ig,islope)
[3285]1123     
[3262]1124                  ! And now we solve correctly the system     
1125                        z1(ig)=1./(za(ig,1)+zb(ig,1)+
1126     &                         zb(ig,2)*(1.-zd(ig,2)))
1127                        zc(ig,1)=(za(ig,1)*zq_tmp_vap(ig,1,iq)+
1128     &                           zb(ig,2)*zc(ig,2) +
1129     &                           (-zdqsdif_tot(ig)) *subtimestep)
1130     &                           * z1(ig)
1131                        zd(ig,1)=zb(ig,1)*z1(ig)
1132                        zq1temp(ig)=zc(ig,1)+ zd(ig,1)
1133     &                              *qsat_ssi(ig,islope)
1134           
[3273]1135                  ! Flux from the subsurface     
1136                        if(old_wsublimation_scheme) then
1137                           zdqsdif_ssi_atm(ig,islope)=rho(ig)
1138     &                     *dryness(ig)*zcdv(ig,islope)
1139     &                     *1/(1+resist(ig,islope))
1140     &                     *(zq1temp(ig)-qsat_ssi(ig,islope))
1141                        else
1142                           zdqsdif_ssi_atm(ig,islope)=rho(ig)*
1143     &                     *dryness(ig) *zcdh(ig,islope)
1144     &                     *1/(1+resist(ig,islope))
1145     &                     *(zq1temp(ig)-qsat_ssi(ig,islope))                         
1146                        endif
[3098]1147
[3253]1148                    else
1149                  ! No atm <-> subsurface exchange, we do it the usual way
1150                        zdqsdif_tot(ig)=-zqsurf(ig)/subtimestep
1151                        z1(ig)=1./(za(ig,1)+ zb(ig,2)*(1.-zd(ig,2)))
1152                        zc(ig,1)=(za(ig,1)*
1153     &                      zq_tmp_vap(ig,1,igcm_h2o_vap)+
1154     &                      zb(ig,2)*zc(ig,2) +
1155     &                      (-zdqsdif_tot(ig)) *subtimestep) *z1(ig)
[3262]1156                            zq1temp(ig)=zc(ig,1)
[3253]1157                    endif ! Of subsurface <-> atmosphere exchange
1158                 endif ! sublimating more than available frost & surface - frost exchange
1159             endif !if .not.watercaptag(ig)
[3124]1160
[3253]1161! --------------------------------------------------------------------------------------------------------------------------------             
1162! We check possible frost subsurface ice interaction: since there is no subsurface water ice mass reservoir represented (yet),
1163! we do not include their effect on the mass of surface frost.
1164! --------------------------------------------------------------------------------------------------------------------------------
[3124]1165
[3253]1166             if((h2o_ice_depth(ig,islope).gt.0).and.lag_layer
1167     &                      .and.(.not.adsorption_soil)) then
1168             ! First case: still frost at the surface but no watercaptag
1169                 if(((watercaptag(ig))).or.
1170     &                (((-zdqsdif_tot(ig)*subtimestep)
1171     &                .lt.(zqsurf(ig)))
1172     &                .and. (zqsurf(ig).gt.tol_frost))) then
1173                  ! Still frost at the surface:  we consider the possibility to have subsurface <-> frost exchange
1174                  ! The flux between frost and ssi is D/zice *(qsat(Tsurf)-qsat_ssi(Tice))   
1175                    call compute_Tice(nsoil, ptsoil(ig,:,islope),
1176     &                                ztsrf(ig),
1177     &                                h2o_ice_depth(ig,islope),
1178     &                                Tice(ig,islope))          ! compute ice temperature
1179   
1180                    call watersat(1,Tice(ig,islope),pplev(ig,1),
1181     &                            qsat_ssi(ig,islope))
1182                 
1183                    qsat_ssi(ig,islope)=ztsrf(ig)/Tice(ig,islope)
1184     &                                       *qsat_ssi(ig,islope)
1185     
1186                    zdqsdif_ssi_frost(ig,islope)= d_coef(ig,islope)
1187     &                                 /h2o_ice_depth(ig,islope)
1188     &                                 *rho(ig)*dryness(ig)
1189     &                                 *(qsat(ig)-qsat_ssi(ig,islope))
[3124]1190
[3273]1191! Line to comment for now since we don't have a mass of subsurface frost
1192! in our computation (otherwise, we would not conserve the H2O mass in
1193! the system)     
1194                  if ((-zdqsdif_tot(ig)+zdqsdif_ssi_frost(ig,islope))
1195     &            *subtimestep.ge.(zqsurf(ig))) then
1196                     zdqsdif_ssi_frost(ig,islope) =
1197     &                         zqsurf(ig)/subtimestep
1198     &                         +   zdqsdif_tot(ig)
1199                  endif
1200
[3253]1201                    zdqsdif_tot(ig) = zdqsdif_tot(ig) -
[3262]1202     &                        zdqsdif_ssi_frost(ig,islope)
[3253]1203                 endif ! watercaptag or frost at the surface
1204             endif ! interaction frost <-> subsurface ice
1205             
[3098]1206
[2515]1207c             Starting upward calculations for water :
[2953]1208             zq_tmp_vap(ig,1,igcm_h2o_vap)=zq1temp(ig)
[3121]1209c             Take into account the H2O latent heat impact on the surface temperature
[2515]1210              if (latentheat_surfwater) then
1211                lh=(2834.3-0.28*(ztsrf(ig)-To)-
1212     &              0.004*(ztsrf(ig)-To)*(ztsrf(ig)-To))*1.e+3
[3115]1213                zdtsrf(ig,islope)=  zdqsdif_tot(ig)*lh
1214     &                              /pcapcal(ig,islope)
[2515]1215              endif ! (latentheat_surfwater) then
1216
[2587]1217              DO ilay=2,nlay
[2953]1218                zq_tmp_vap(ig,ilay,iq)=zc(ig,ilay)+zd(ig,ilay)
1219     &                              *zq_tmp_vap(ig,ilay-1,iq)
[2587]1220              ENDDO
[2515]1221c             Subtimestep water budget :
[2953]1222              ztsrf(ig) = ztsrf(ig)+(pdtsrf(ig,islope)
1223     &                 + zdtsrf(ig,islope))*subtimestep
[2515]1224              zqsurf(ig)= zqsurf(ig)+(
[3115]1225     &                       zdqsdif_tot(ig))*subtimestep
[3400]1226
[3253]1227              if (zqsurf(ig)<0 .and.(.not.watercaptag(ig))) then
[3134]1228                      zqsurf(ig)=0
1229              endif
[3253]1230              zdqsdif_ssi_atm_tot(ig,islope) =
1231     &                              zdqsdif_ssi_atm_tot(ig,islope)
[3285]1232     &                     + zdqsdif_ssi_atm(ig,islope)*subtimestep
[3253]1233              zdqsdif_ssi_frost_tot(ig,islope) =
1234     &                                zdqsdif_ssi_frost_tot(ig,islope)
[3285]1235     &                     +zdqsdif_ssi_frost(ig,islope)*subtimestep
[2840]1236c             Monitoring instantaneous latent heat flux in W.m-2 :
[2953]1237              zsurf_h2o_lh(ig,islope) = zsurf_h2o_lh(ig,islope)+
1238     &                    (zdtsrf(ig,islope)*pcapcal(ig,islope))
1239     &                    *subtimestep
[2515]1240
1241c             We ensure that surface temperature can't rise above the solid domain if there
1242c             is still ice on the surface (oldschool)
1243               if(zqsurf(ig)
[3115]1244     &           +zdqsdif_tot(ig)*subtimestep
[2953]1245     &           .gt.frost_albedo_threshold) then ! if there is still ice, T cannot exceed To
1246                 zdtsrf(ig,islope) = min(zdtsrf(ig,islope),
1247     &                          (To-ztsrf(ig))/subtimestep) ! ice melt case
1248               endif
[2515]1249
[3253]1250c             End of the subtimestep
[3124]1251            ENDDO ! tsub=1,nsubtimestep
1252
[2587]1253c             Integration of subtimestep temp and water budget :
1254c             (btw could also compute the post timestep temp and ice
1255c             by simply adding the subtimestep trend instead of this)
[2953]1256            surf_h2o_lh(ig,islope)= zsurf_h2o_lh(ig,islope)/ptimestep
1257            pdtsrf(ig,islope)= (ztsrf(ig) -
1258     &                     ptsrf(ig,islope))/ptimestep
1259            pdqsdif(ig,igcm_h2o_ice,islope)=
1260     &                  (zqsurf(ig)- pqsurf(ig,igcm_h2o_ice,islope)/
1261     &                       cos(pi*def_slope_mean(islope)/180.))
1262     &                       /ptimestep
[3262]1263     
[3285]1264            zdqsdif_ssi_atm_tot(ig,islope)=
1265     &           zdqsdif_ssi_atm_tot(ig,islope)/ptimestep
1266            zdqsdif_ssi_frost_tot(ig,islope)=
1267     &           zdqsdif_ssi_frost_tot(ig,islope)/ptimestep
1268                   zdqsdif_ssi_tot(ig,islope) =
[3253]1269     &                        zdqsdif_ssi_atm_tot(ig,islope)
1270     &                      + zdqsdif_ssi_frost_tot(ig,islope)
[2587]1271c if subliming more than qsurf(ice) and on watercaptag, water
1272c sublimates from watercap reservoir (dwatercap_dif is <0)
1273            if(watercaptag(ig)) then
[2953]1274              if ((-pdqsdif(ig,igcm_h2o_ice,islope)*ptimestep)
1275     &           .gt.(pqsurf(ig,igcm_h2o_ice,islope)
1276     &                 /cos(pi*def_slope_mean(islope)/180.))) then
1277                 dwatercap_dif(ig,islope)=
1278     &                     pdqsdif(ig,igcm_h2o_ice,islope)+
1279     &                     (pqsurf(ig,igcm_h2o_ice,islope) /
1280     &          cos(pi*def_slope_mean(islope)/180.))/ptimestep
1281                 pdqsdif(ig,igcm_h2o_ice,islope)=
1282     &                   - (pqsurf(ig,igcm_h2o_ice,islope)/
1283     &          cos(pi*def_slope_mean(islope)/180.))/ptimestep
[2587]1284              endif! ((-pdqsdif(ig)*ptimestep)
1285            endif !(watercaptag(ig)) then
[2953]1286           zq_slope_vap(ig,:,:,islope) = zq_tmp_vap(ig,:,:)
1287           ENDDO ! of DO ig=1,ngrid
1288          ENDDO ! islope
1289c Some grid box averages: interface with the atmosphere
1290       DO ig = 1,ngrid
1291         DO ilay = 1,nlay
1292           zq(ig,ilay,iq) = 0.
1293           DO islope = 1,nslope
1294             zq(ig,ilay,iq) = zq(ig,ilay,iq) +
1295     $           zq_slope_vap(ig,ilay,iq,islope) *
1296     $           subslope_dist(ig,islope)
1297           ENDDO
1298         ENDDO
1299       ENDDO
1300! Recompute values in kg/m^2 slopped
1301        DO ig = 1,ngrid
1302          DO islope = 1,nslope 
1303             pdqsdif(ig,igcm_h2o_ice,islope) =
1304     &      pdqsdif(ig,igcm_h2o_ice,islope)
1305     &      * cos(pi*def_slope_mean(islope)/180.)
1306
1307             dwatercap_dif(ig,islope) =
1308     &      dwatercap_dif(ig,islope)
1309     &      * cos(pi*def_slope_mean(islope)/180.)
1310          ENDDO
1311         ENDDO
1312
[2515]1313          END IF ! of IF ((water).and.(iq.eq.igcm_h2o_vap))
1314
1315c --------- end of h2o_vap ----------------------------
1316
1317c --------- hdo_vap -----------------------------------
1318
1319c    hdo_ice has already been with along h2o_ice
1320c    amongst "normal" tracers (ie not h2o_vap)
1321
1322         if (hdo.and.(iq.eq.igcm_hdo_vap)) then
1323          zb(1:ngrid,2:nlay)=zkh(1:ngrid,2:nlay)*zb0(1:ngrid,2:nlay)
1324          zb(1:ngrid,1)=0
1325
[38]1326          DO ig=1,ngrid
1327               z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
1328               zc(ig,nlay)=za(ig,nlay)*zq(ig,nlay,iq)*z1(ig)
1329               zd(ig,nlay)=zb(ig,nlay)*z1(ig)
1330          ENDDO
[2515]1331
[38]1332          DO ilay=nlay-1,2,-1
1333               DO ig=1,ngrid
1334                z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
1335     $           zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
1336                zc(ig,ilay)=(za(ig,ilay)*zq(ig,ilay,iq)+
1337     $           zb(ig,ilay+1)*zc(ig,ilay+1))*z1(ig)
1338                zd(ig,ilay)=zb(ig,ilay)*z1(ig)
1339               ENDDO
1340          ENDDO
[2953]1341          hdoflux_meshavg(:) = 0.
1342          DO islope = 1,nslope
[38]1343
[2953]1344            pdqsdif_tmphdo(:,:) = pdqsdif(:,:,islope)
1345     &      /cos(pi*def_slope_mean(islope)/180.)
1346
1347            call watersat(ngrid,pdtsrf(:,islope)*ptimestep +
1348     &                     ptsrf(:,islope),pplev(:,1),qsat_tmp)
1349
[2312]1350            CALL hdo_surfex(ngrid,nlay,nq,ptimestep,
[2953]1351     &                      zt,pplay,zq,pqsurf(:,:,islope),
1352     &                      saved_h2o_vap,qsat_tmp,
1353     &                      pdqsdif_tmphdo,
1354     & dwatercap_dif(:,islope)/cos(pi*def_slope_mean(islope)/180.),
1355     &                      hdoflux(:,islope))
1356
1357            pdqsdif(:,:,islope) = pdqsdif_tmphdo(:,:) *
1358     &               cos(pi*def_slope_mean(islope)/180.)
1359            DO ig = 1,ngrid
1360              hdoflux_meshavg(ig) = hdoflux_meshavg(ig) +
1361     &                     hdoflux(ig,islope)*subslope_dist(ig,islope)
1362
1363            ENDDO !ig = 1,ngrid
1364          ENDDO !islope = 1,nslope
1365
[2312]1366            DO ig=1,ngrid
1367                z1(ig)=1./(za(ig,1)+zb(ig,1)+
1368     $           zb(ig,2)*(1.-zd(ig,2)))
1369                zc(ig,1)=(za(ig,1)*zq(ig,1,iq)+
1370     $         zb(ig,2)*zc(ig,2) +
[2953]1371     $        (-hdoflux_meshavg(ig)) *ptimestep) *z1(ig)  !tracer flux from surface
[2312]1372            ENDDO
1373
[38]1374            DO ig=1,ngrid
[2515]1375              zq(ig,1,iq)=zc(ig,1)
1376              DO ilay=2,nlay
1377                zq(ig,ilay,iq)=zc(ig,ilay)+zd(ig,ilay)*zq(ig,ilay-1,iq)
1378              ENDDO
[38]1379            ENDDO
[2515]1380         endif ! (hdo.and.(iq.eq.igcm_hdo_vap))
[2312]1381
[2515]1382c --------- end of hdo ----------------------------
[38]1383
[2515]1384        enddo ! of do iq=1,nq
[38]1385
[2515]1386c --------- end of tracers  ----------------------------
[2312]1387
[2932]1388         call write_output("surf_h2o_lh",
[3107]1389     &                     "Ground ice latent heat flux",
1390     &                     "W.m-2",surf_h2o_lh(:,iflat))
[3124]1391
[3253]1392         call write_output('zdqsdif_ssi_frost_tot',
1393     &     'Flux between frost and subsurface ice','kg.m-2.s-1',
1394     &                zdqsdif_ssi_frost_tot(:,iflat))
1395
1396         call write_output('zdqsdif_ssi_atm_tot',
1397     &     'Flux between atmosphere and subsurface ice','kg.m-2.s-1',
1398     &                zdqsdif_ssi_atm_tot(:,iflat))
1399
1400         call write_output('zdqsdif_ssi_tot',
1401     &     'Total flux echange with subsurface ice','kg.m-2.s-1',
1402     &                zdqsdif_ssi_tot(:,iflat))
1403
1404
[2312]1405C       Diagnostic output for HDO
[2934]1406!        if (hdo) then
1407!            CALL write_output('hdoflux',
1408!     &                       'hdoflux',
[2953]1409!     &                       ' ',hdoflux_meshavg(:)) 
[2934]1410!            CALL write_output('h2oflux',
1411!     &                       'h2oflux',
1412!     &                       ' ',h2oflux(:))
1413!        endif
[2312]1414
[38]1415c-----------------------------------------------------------------------
1416c   8. calcul final des tendances de la diffusion verticale
1417c      ----------------------------------------------------
1418
1419      DO ilev = 1, nlay
1420         DO ig=1,ngrid
1421            pdudif(ig,ilev)=(    zu(ig,ilev)-
1422     $      (pu(ig,ilev)+pdufi(ig,ilev)*ptimestep)    )/ptimestep
1423            pdvdif(ig,ilev)=(    zv(ig,ilev)-
1424     $      (pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep)    )/ptimestep
[473]1425            hh = ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep
1426     $  + (latcond*dmice(ig,ilev)/cpp)/ppopsk(ig,ilev)
1427            pdhdif(ig,ilev)=( zhs(ig,ilev)- hh )/ptimestep
[38]1428         ENDDO
1429      ENDDO
1430
[2823]1431      pdqdif(1:ngrid,1:nlay,1:nq)=(zq(1:ngrid,1:nlay,1:nq)-
1432     &                             (pq(1:ngrid,1:nlay,1:nq)
1433     &                              +pdqfi(1:ngrid,1:nlay,1:nq)
1434     &                                    *ptimestep))/ptimestep
[38]1435
[1036]1436      END SUBROUTINE vdifc
[1969]1437
[2312]1438c====================================
1439
[2515]1440      SUBROUTINE make_tsub(naersize,dtsurf,qsurf,ptimestep,
1441     $                     dtmax,watercaptag,ntsub)
[2312]1442
[2515]1443c Pas de temps adaptatif en estimant le taux de sublimation
1444c et en adaptant avec un critere "dtmax" du chauffage a accomoder
1445c dtmax est regle empiriquement (pour l'instant) a 0.5 K
[3348]1446      use surfdat_h, only: frost_albedo_threshold
[2515]1447
1448      integer,intent(in) :: naersize
1449      real,intent(in) :: dtsurf(naersize)
1450      real,intent(in) :: qsurf(naersize)
1451      logical,intent(in) :: watercaptag(naersize)
1452      real,intent(in) :: ptimestep
1453      real,intent(in)  :: dtmax
1454      real  :: ztsub(naersize)
1455      integer  :: i
1456      integer,intent(out) :: ntsub(naersize)
1457
1458      do i=1,naersize
[3348]1459c     If not on permanent ice or thin frost layer :
1460c     do not use a subtimestep (ntsub=1)
1461        if ((qsurf(i).lt.frost_albedo_threshold).and.
[2515]1462     &      (.not.watercaptag(i))) then
1463          ntsub(i) = 1
1464        else
1465          ztsub(i) = ptimestep * dtsurf(i) / dtmax
1466          ntsub(i) = ceiling(abs(ztsub(i)))
1467        endif ! (qsurf(i).eq.0) then
1468c     
1469c       write(78,*), dtsurf*ptimestep, dtsurf, ntsub
1470      enddo! 1=1,ngrid
1471
1472
1473
1474      END SUBROUTINE make_tsub
[3253]1475     
1476     
1477c====================================
1478
1479      SUBROUTINE compute_Tice(nsoil, ptsoil, ptsurf, ice_depth, Tice)
1480
1481c Compute subsurface ice temperature by interpolating the temperatures between the two adjacent cells.
1482      use comsoil_h, only: layer, mlayer
1483
1484      implicit none
1485      integer,intent(in) :: nsoil            ! Number of soil layers
1486      real,intent(in) :: ptsoil(nsoil)       ! Soil temperature (K)
1487      real,intent(in) :: ptsurf              ! Soil temperature (K)
1488      real,intent(in) :: ice_depth           ! Ice depth (m)
1489      real,intent(out) :: Tice               ! Ice temperatures (K)
1490
1491c Local
1492      integer :: ik                          ! Loop variables
1493      integer :: indexice                    ! Index of the ice
1494
1495c Code:
1496      indexice = -1
1497      if(ice_depth.lt.mlayer(0)) then
1498          indexice = 0.
1499      else
1500          do ik = 0,nsoil-2 ! go through all the layers to find the ice locations
1501             if((mlayer(ik).le.ice_depth).and.
1502     &          (mlayer(ik+1).gt.ice_depth)) then
1503                    indexice = ik+1
1504                    exit
1505             endif
1506          enddo
1507      endif
1508     
1509      if(indexice.lt.0) then
1510          call abort_physic("vdifc - compute Tice",
1511     &         "subsurface ice is below the last soil layer",1)
1512      else
1513          if(indexice .ge. 1) then ! Linear inteprolation between soil temperature
1514              Tice = (ptsoil(indexice)-ptsoil(indexice+1))
1515     &               /(mlayer(indexice-1)-mlayer(indexice))
1516     &               *(ice_depth-mlayer(indexice)) + ptsoil(indexice+1)
1517          else ! Linear inteprolation between the 1st soil temperature and the surface temperature
1518              Tice = (ptsoil(1) - ptsurf)/mlayer(0)
1519     &               *(ice_depth-mlayer(0)) + ptsoil(1)
1520          endif ! index ice >=0
1521      endif !indexice <0
1522     
1523
1524      END SUBROUTINE compute_Tice
[1969]1525      END MODULE vdifc_mod
Note: See TracBrowser for help on using the repository browser.