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

Last change on this file since 2613 was 2593, checked in by mvals, 3 years ago

Mars GCM:
hdo_surfex_mod.F, vdifc_mod.F: addings to account for the effect of kinetics in the fractionation by condensation of HDO at the surface
MV

File size: 38.4 KB
RevLine 
[1969]1      MODULE vdifc_mod
2     
3      IMPLICIT NONE
4     
5      CONTAINS
6     
[38]7      SUBROUTINE vdifc(ngrid,nlay,nq,co2ice,ppopsk,
8     $                ptimestep,pcapcal,lecrit,
9     $                pplay,pplev,pzlay,pzlev,pz0,
10     $                pu,pv,ph,pq,ptsrf,pemis,pqsurf,
11     $                pdufi,pdvfi,pdhfi,pdqfi,pfluxsrf,
12     $                pdudif,pdvdif,pdhdif,pdtsrf,pq2,
[660]13     $                pdqdif,pdqsdif,wstar,zcdv_true,zcdh_true,
[1996]14     $                hfmax,pcondicea_co2microp,sensibFlux,
[2260]15     $                dustliftday,local_time,watercap, dwatercap_dif)
[1236]16
[1036]17      use tracer_mod, only: noms, igcm_dust_mass, igcm_dust_number,
18     &                      igcm_dust_submicron, igcm_h2o_vap,
[1974]19     &                      igcm_h2o_ice, alpha_lift,
[2312]20     &                      igcm_hdo_vap, igcm_hdo_ice,
[1974]21     &                      igcm_stormdust_mass, igcm_stormdust_number
[1224]22      use surfdat_h, only: watercaptag, frost_albedo_threshold, dryness
[1969]23      USE comcstfi_h, ONLY: cpp, r, rcp, g
[1996]24      use watersat_mod, only: watersat
[1242]25      use turb_mod, only: turb_resolved, ustar, tstar
[2160]26      use compute_dtau_mod, only: ti_injection_sol,tf_injection_sol
[2312]27      use hdo_surfex_mod, only: hdo_surfex
[2515]28c      use geometry_mod, only: longitude_deg,latitude_deg !  Joseph
[2409]29      use dust_param_mod, only: doubleq, submicron, lifting
[2312]30
[38]31      IMPLICIT NONE
32
33c=======================================================================
34c
35c   subject:
36c   --------
37c   Turbulent diffusion (mixing) for potential T, U, V and tracer
38c
39c   Shema implicite
40c   On commence par rajouter au variables x la tendance physique
41c   et on resoult en fait:
42c      x(t+1) =  x(t) + dt * (dx/dt)phys(t)  +  dt * (dx/dt)difv(t+1)
43c
44c   author:
45c   ------
46c      Hourdin/Forget/Fournier
47c=======================================================================
48
49c-----------------------------------------------------------------------
50c   declarations:
51c   -------------
52
[1944]53      include "callkeys.h"
54      include "microphys.h"
[38]55
56c
57c   arguments:
58c   ----------
59
[1036]60      INTEGER,INTENT(IN) :: ngrid,nlay
61      REAL,INTENT(IN) :: ptimestep
62      REAL,INTENT(IN) :: pplay(ngrid,nlay),pplev(ngrid,nlay+1)
63      REAL,INTENT(IN) :: pzlay(ngrid,nlay),pzlev(ngrid,nlay+1)
64      REAL,INTENT(IN) :: pu(ngrid,nlay),pv(ngrid,nlay)
65      REAL,INTENT(IN) :: ph(ngrid,nlay)
66      REAL,INTENT(IN) :: ptsrf(ngrid),pemis(ngrid)
67      REAL,INTENT(IN) :: pdufi(ngrid,nlay),pdvfi(ngrid,nlay)
68      REAL,INTENT(IN) :: pdhfi(ngrid,nlay)
69      REAL,INTENT(IN) :: pfluxsrf(ngrid)
70      REAL,INTENT(OUT) :: pdudif(ngrid,nlay),pdvdif(ngrid,nlay)
71      REAL,INTENT(OUT) :: pdtsrf(ngrid),pdhdif(ngrid,nlay)
[1130]72      REAL,INTENT(IN) :: pcapcal(ngrid)
73      REAL,INTENT(INOUT) :: pq2(ngrid,nlay+1)
[38]74
75c    Argument added for condensation:
[1036]76      REAL,INTENT(IN) :: co2ice (ngrid), ppopsk(ngrid,nlay)
77      logical,INTENT(IN) :: lecrit
[1996]78      REAL,INTENT(IN) :: pcondicea_co2microp(ngrid,nlay)! tendency due to CO2 condensation (kg/kg.s-1)
79     
[1036]80      REAL,INTENT(IN) :: pz0(ngrid) ! surface roughness length (m)
[224]81
[256]82c    Argument added to account for subgrid gustiness :
83
[1944]84      REAL,INTENT(IN) :: wstar(ngrid), hfmax(ngrid)!, zi(ngrid)
[256]85
[38]86c    Traceurs :
[1036]87      integer,intent(in) :: nq
88      REAL,INTENT(IN) :: pqsurf(ngrid,nq)
[2515]89      REAL :: zqsurf(ngrid) ! temporary water tracer
[1036]90      real,intent(in) :: pq(ngrid,nlay,nq), pdqfi(ngrid,nlay,nq)
91      real,intent(out) :: pdqdif(ngrid,nlay,nq)
[1974]92      real,intent(out) :: pdqsdif(ngrid,nq)
93      REAL,INTENT(in) :: dustliftday(ngrid)
94      REAL,INTENT(in) :: local_time(ngrid)
[38]95     
96c   local:
97c   ------
98
[1130]99      REAL :: pt(ngrid,nlay)
[473]100 
[38]101      INTEGER ilev,ig,ilay,nlev
102
[1047]103      REAL z4st,zdplanck(ngrid)
104      REAL zkv(ngrid,nlay+1),zkh(ngrid,nlay+1)
105      REAL zkq(ngrid,nlay+1)
106      REAL zcdv(ngrid),zcdh(ngrid)
107      REAL zcdv_true(ngrid),zcdh_true(ngrid)    ! drag coeff are used by the LES to recompute u* and hfx
108      REAL zu(ngrid,nlay),zv(ngrid,nlay)
109      REAL zh(ngrid,nlay)
110      REAL ztsrf2(ngrid)
111      REAL z1(ngrid),z2(ngrid)
112      REAL za(ngrid,nlay),zb(ngrid,nlay)
113      REAL zb0(ngrid,nlay)
114      REAL zc(ngrid,nlay),zd(ngrid,nlay)
[38]115      REAL zcst1
[1047]116      REAL zu2(ngrid)
[38]117
118      EXTERNAL SSUM,SCOPY
119      REAL SSUM
[1036]120      LOGICAL,SAVE :: firstcall=.true.
[38]121
[626]122
[38]123c     variable added for CO2 condensation:
124c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[1047]125      REAL hh , zhcond(ngrid,nlay)
126      REAL,PARAMETER :: latcond=5.9e5
127      REAL,PARAMETER :: tcond1mb=136.27
128      REAL,SAVE :: acond,bcond
[669]129     
[2515]130c     Subtimestep & implicit treatment of water vapor
131c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
132      REAL zdqsdif(ngrid) ! subtimestep pdqsdif for water ice
133      REAL ztsrf(ngrid) ! temporary surface temperature in tsub
134      REAL zdtsrf(ngrid) ! surface temperature tendancy in tsub
135
[2179]136c     For latent heat release from ground water ice sublimation   
[2515]137c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
138      REAL tsrf_lh(ngrid) ! temporary surface temperature with lh effect
[2179]139      REAL lh ! latent heat, formulation given in the Technical Document:
140              ! "Modeling water ice sublimation under Phoenix-like conditions", Montmessin et al. 2004 
[38]141
142c    Tracers :
143c    ~~~~~~~
144      INTEGER iq
[1047]145      REAL zq(ngrid,nlay,nq)
146      REAL zq1temp(ngrid)
147      REAL rho(ngrid) ! near surface air density
148      REAL qsat(ngrid)
[38]149
[2312]150      REAL hdoflux(ngrid)       ! value of vapour flux of HDO
151      REAL h2oflux(ngrid)       ! value of vapour flux of H2O
152      REAL old_h2o_vap(ngrid)   ! traceur d'eau avant traitement
153
[38]154      REAL kmixmin
155
[2260]156c    Argument added for surface water ice budget:
157      REAL,INTENT(IN) :: watercap(ngrid)
158      REAL,INTENT(OUT) :: dwatercap_dif(ngrid)
159
[2515]160c    Subtimestep to compute h2o latent heat flux:
161      REAL :: dtmax = 0.5 ! subtimestep temp criterion
162      INTEGER tsub ! adaptative subtimestep (seconds)
163      REAL subtimestep !ptimestep/nsubtimestep
164      INTEGER nsubtimestep(ngrid) ! number of subtimestep (int)
165
[473]166c    Mass-variation scheme :
167c    ~~~~~~~
168
169      INTEGER j,l
[1047]170      REAL zcondicea(ngrid,nlay)
171      REAL zt(ngrid,nlay),ztcond(ngrid,nlay+1)
172      REAL betam(ngrid,nlay),dmice(ngrid,nlay)
173      REAL pdtc(ngrid,nlay)
174      REAL zhs(ngrid,nlay)
175      REAL,SAVE :: ccond
[473]176
177c     Theta_m formulation for mass-variation scheme :
178c     ~~~~~~~
179
[1047]180      INTEGER,SAVE :: ico2
181      INTEGER llnt(ngrid)
182      REAL,SAVE :: m_co2, m_noco2, A , B
183      REAL vmr_co2(ngrid,nlay)
[473]184      REAL qco2,mmean
185
[1047]186      REAL,INTENT(OUT) :: sensibFlux(ngrid)
[660]187
[2312]188!!MARGAUX
189      REAL DoH_vap(ngrid,nlay)
190
[38]191c    ** un petit test de coherence
192c       --------------------------
193
[1779]194      ! AS: OK firstcall absolute
[38]195      IF (firstcall) THEN
196c        To compute: Tcond= 1./(bcond-acond*log(.0095*p)) (p in pascal)
197         bcond=1./tcond1mb
198         acond=r/latcond
[473]199         ccond=cpp/(g*latcond)
[38]200         PRINT*,'In vdifc: Tcond(P=1mb)=',tcond1mb,' Lcond=',latcond
[473]201         PRINT*,'          acond,bcond,ccond',acond,bcond,ccond
[38]202
[473]203
204         ico2=0
205
206         if (tracer) then
207c          Prepare Special treatment if one of the tracer is CO2 gas
[1036]208           do iq=1,nq
[473]209             if (noms(iq).eq."co2") then
210                ico2=iq
211                m_co2 = 44.01E-3  ! CO2 molecular mass (kg/mol)
212                m_noco2 = 33.37E-3  ! Non condensible mol mass (kg/mol)
213c               Compute A and B coefficient use to compute
214c               mean molecular mass Mair defined by
215c               1/Mair = q(ico2)/m_co2 + (1-q(ico2))/m_noco2
216c               1/Mair = A*q(ico2) + B
217                A =(1/m_co2 - 1/m_noco2)
218                B=1/m_noco2
219             endif
220           enddo
221         end if
222
[38]223        firstcall=.false.
224      ENDIF
225
226
227c-----------------------------------------------------------------------
228c    1. initialisation
229c    -----------------
230
231      nlev=nlay+1
232
[1035]233      ! initialize output tendencies to zero:
234      pdudif(1:ngrid,1:nlay)=0
235      pdvdif(1:ngrid,1:nlay)=0
236      pdhdif(1:ngrid,1:nlay)=0
237      pdtsrf(1:ngrid)=0
[2515]238      zdtsrf(1:ngrid)=0
[1035]239      pdqdif(1:ngrid,1:nlay,1:nq)=0
240      pdqsdif(1:ngrid,1:nq)=0
[2515]241      zdqsdif(1:ngrid)=0
[2260]242      dwatercap_dif(1:ngrid)=0
[1035]243
[38]244c    ** calcul de rho*dz et dt*rho/dz=dt*rho**2 g/dp
245c       avec rho=p/RT=p/ (R Theta) (p/ps)**kappa
246c       ----------------------------------------
247
248      DO ilay=1,nlay
249         DO ig=1,ngrid
250            za(ig,ilay)=(pplev(ig,ilay)-pplev(ig,ilay+1))/g
[473]251! Mass variation scheme:
252            betam(ig,ilay)=-za(ig,ilay)*latcond/(cpp*ppopsk(ig,ilay))
[38]253         ENDDO
254      ENDDO
255
256      zcst1=4.*g*ptimestep/(r*r)
257      DO ilev=2,nlev-1
258         DO ig=1,ngrid
259            zb0(ig,ilev)=pplev(ig,ilev)*
260     s      (pplev(ig,1)/pplev(ig,ilev))**rcp /
261     s      (ph(ig,ilev-1)+ph(ig,ilev))
262            zb0(ig,ilev)=zcst1*zb0(ig,ilev)*zb0(ig,ilev)/
263     s      (pplay(ig,ilev-1)-pplay(ig,ilev))
264         ENDDO
265      ENDDO
266      DO ig=1,ngrid
[473]267         zb0(ig,1)=ptimestep*pplev(ig,1)/(r*ptsrf(ig))
[38]268      ENDDO
269
270c    ** diagnostique pour l'initialisation
271c       ----------------------------------
272
273      IF(lecrit) THEN
274         ig=ngrid/2+1
275         PRINT*,'Pression (mbar) ,altitude (km),u,v,theta, rho dz'
276         DO ilay=1,nlay
277            WRITE(*,'(6f11.5)')
278     s      .01*pplay(ig,ilay),.001*pzlay(ig,ilay),
279     s      pu(ig,ilay),pv(ig,ilay),ph(ig,ilay),za(ig,ilay)
280         ENDDO
281         PRINT*,'Pression (mbar) ,altitude (km),zb'
282         DO ilev=1,nlay
283            WRITE(*,'(3f15.7)')
284     s      .01*pplev(ig,ilev),.001*pzlev(ig,ilev),
285     s      zb0(ig,ilev)
286         ENDDO
287      ENDIF
288
[473]289c     -----------------------------------
[38]290c     Potential Condensation temperature:
291c     -----------------------------------
292
[473]293c     Compute CO2 Volume mixing ratio
294c     -------------------------------
295      if (callcond.and.(ico2.ne.0)) then
296         DO ilev=1,nlay
297            DO ig=1,ngrid
[529]298              qco2=MAX(1.E-30
299     &             ,pq(ig,ilev,ico2)+pdqfi(ig,ilev,ico2)*ptimestep)
[473]300c             Mean air molecular mass = 1/(q(ico2)/m_co2 + (1-q(ico2))/m_noco2)
301              mmean=1/(A*qco2 +B)
302              vmr_co2(ig,ilev) = qco2*mmean/m_co2
303            ENDDO
304         ENDDO
305      else
306         DO ilev=1,nlay
307            DO ig=1,ngrid
308              vmr_co2(ig,ilev)=0.95
309            ENDDO
310         ENDDO
311      end if
[38]312
[473]313c     forecast of atmospheric temperature zt and frost temperature ztcond
314c     --------------------------------------------------------------------
[38]315
[473]316      if (callcond) then
317        DO ilev=1,nlay
318          DO ig=1,ngrid
[884]319              ztcond(ig,ilev)=
320     &      1./(bcond-acond*log(.01*vmr_co2(ig,ilev)*pplay(ig,ilev)))
321            if (pplay(ig,ilev).lt.1e-4) ztcond(ig,ilev)=0.0 !mars Monica
[473]322!            zhcond(ig,ilev) =
323!     &  (1./(bcond-acond*log(.0095*pplay(ig,ilev))))/ppopsk(ig,ilev)
324              zhcond(ig,ilev) = ztcond(ig,ilev)/ppopsk(ig,ilev)
325          END DO
326        END DO
[884]327        ztcond(:,nlay+1)=ztcond(:,nlay)
[473]328      else
[884]329         zhcond(:,:) = 0
330         ztcond(:,:) = 0
[473]331      end if
332
333
[38]334c-----------------------------------------------------------------------
335c   2. ajout des tendances physiques
336c      -----------------------------
337
338      DO ilev=1,nlay
339         DO ig=1,ngrid
340            zu(ig,ilev)=pu(ig,ilev)+pdufi(ig,ilev)*ptimestep
341            zv(ig,ilev)=pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep
342            zh(ig,ilev)=ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep
[473]343!            zh(ig,ilev)=max(zh(ig,ilev),zhcond(ig,ilev))
[38]344         ENDDO
345      ENDDO
346      if(tracer) then
347        DO iq =1, nq
348         DO ilev=1,nlay
349           DO ig=1,ngrid
350              zq(ig,ilev,iq)=pq(ig,ilev,iq)+pdqfi(ig,ilev,iq)*ptimestep
351           ENDDO
352         ENDDO
353        ENDDO
354      end if
355
356c-----------------------------------------------------------------------
357c   3. schema de turbulence
358c      --------------------
359
360c    ** source d'energie cinetique turbulente a la surface
361c       (condition aux limites du schema de diffusion turbulente
362c       dans la couche limite
363c       ---------------------
364
[499]365      CALL vdif_cd(ngrid,nlay,pz0,g,pzlay,pu,pv,wstar,ptsrf,ph
[1236]366     &             ,zcdv_true,zcdh_true)
[256]367
[290]368        zu2(:)=pu(:,1)*pu(:,1)+pv(:,1)*pv(:,1)
[256]369
[291]370        IF (callrichsl) THEN
[545]371          zcdv(:)=zcdv_true(:)*sqrt(zu2(:)+
372     &     (log(1.+0.7*wstar(:) + 2.3*wstar(:)**2))**2)
373          zcdh(:)=zcdh_true(:)*sqrt(zu2(:)+
374     &     (log(1.+0.7*wstar(:) + 2.3*wstar(:)**2))**2)
[496]375
[1242]376           ustar(:)=sqrt(zcdv_true(:))*sqrt(zu2(:)+
[545]377     &     (log(1.+0.7*wstar(:) + 2.3*wstar(:)**2))**2)
378
[1242]379           tstar(:)=0.
[545]380           DO ig=1,ngrid
381              IF (zcdh_true(ig) .ne. 0.) THEN      ! When Cd=Ch=0, u*=t*=0
[1242]382                 tstar(ig)=(ph(ig,1)-ptsrf(ig))*zcdh(ig)/ustar(ig)
[545]383              ENDIF
384           ENDDO
385
[284]386        ELSE
[545]387           zcdv(:)=zcdv_true(:)*sqrt(zu2(:))     ! 1 / bulk aerodynamic momentum conductance
388           zcdh(:)=zcdh_true(:)*sqrt(zu2(:))     ! 1 / bulk aerodynamic heat conductance
[1242]389           ustar(:)=sqrt(zcdv_true(:))*sqrt(zu2(:))
390           tstar(:)=(ph(:,1)-ptsrf(:))*zcdh_true(:)/sqrt(zcdv_true(:))
[290]391        ENDIF
[38]392
[529]393! Some usefull diagnostics for the new surface layer parametrization :
[290]394
[1130]395!        call WRITEDIAGFI(ngrid,'vdifc_zcdv_true',
[256]396!     &              'momentum drag','no units',
397!     &                         2,zcdv_true)
[1130]398!        call WRITEDIAGFI(ngrid,'vdifc_zcdh_true',
[256]399!     &              'heat drag','no units',
400!     &                         2,zcdh_true)
[1130]401!        call WRITEDIAGFI(ngrid,'vdifc_ust',
[473]402!     &              'friction velocity','m/s',2,ust)
[1130]403!       call WRITEDIAGFI(ngrid,'vdifc_tst',
[473]404!     &              'friction temperature','K',2,tst)
[1130]405!        call WRITEDIAGFI(ngrid,'vdifc_zcdv',
[268]406!     &              'aerodyn momentum conductance','m/s',
[256]407!     &                         2,zcdv)
[1130]408!        call WRITEDIAGFI(ngrid,'vdifc_zcdh',
[268]409!     &              'aerodyn heat conductance','m/s',
[256]410!     &                         2,zcdh)
411
[38]412c    ** schema de diffusion turbulente dans la couche limite
413c       ----------------------------------------------------
[1240]414       IF (.not. callyamada4) THEN
[529]415
[1130]416       CALL vdif_kc(ngrid,nlay,nq,ptimestep,g,pzlev,pzlay
[38]417     &              ,pu,pv,ph,zcdv_true
[325]418     &              ,pq2,zkv,zkh,zq)
[529]419
[544]420      ELSE
[38]421
[555]422      pt(:,:)=ph(:,:)*ppopsk(:,:)
[1036]423      CALL yamada4(ngrid,nlay,nq,ptimestep,g,r,pplev,pt
[1242]424     s   ,pzlev,pzlay,pu,pv,ph,pq,zcdv_true,pq2,zkv,zkh,zkq,ustar
[652]425     s   ,9)
[544]426      ENDIF
427
[38]428      if ((doubleq).and.(ngrid.eq.1)) then
429        kmixmin = 80. !80.! minimum eddy mix coeff in 1D
430        do ilev=1,nlay
431          do ig=1,ngrid
432           zkh(ig,ilev) = max(kmixmin,zkh(ig,ilev))
433           zkv(ig,ilev) = max(kmixmin,zkv(ig,ilev))
434          end do
435        end do
436      end if
437
438c    ** diagnostique pour le schema de turbulence
439c       -----------------------------------------
440
441      IF(lecrit) THEN
442         PRINT*
443         PRINT*,'Diagnostic for the vertical turbulent mixing'
444         PRINT*,'Cd for momentum and potential temperature'
445
446         PRINT*,zcdv(ngrid/2+1),zcdh(ngrid/2+1)
447         PRINT*,'Mixing coefficient for momentum and pot.temp.'
448         DO ilev=1,nlay
449            PRINT*,zkv(ngrid/2+1,ilev),zkh(ngrid/2+1,ilev)
450         ENDDO
451      ENDIF
452
453
454
455
456c-----------------------------------------------------------------------
457c   4. inversion pour l'implicite sur u
458c      --------------------------------
459
460c    ** l'equation est
461c       u(t+1) =  u(t) + dt * {(du/dt)phys}(t)  +  dt * {(du/dt)difv}(t+1)
462c       avec
463c       /zu/ = u(t) + dt * {(du/dt)phys}(t)   (voir paragraphe 2.)
464c       et
465c       dt * {(du/dt)difv}(t+1) = dt * {(d/dz)[ Ku (du/dz) ]}(t+1)
466c       donc les entrees sont /zcdv/ pour la condition a la limite sol
467c       et /zkv/ = Ku
[2312]468
[2274]469      zb(1:ngrid,2:nlay)=zkv(1:ngrid,2:nlay)*zb0(1:ngrid,2:nlay)
470      zb(1:ngrid,1)=zcdv(1:ngrid)*zb0(1:ngrid,1)
[38]471
472      DO ig=1,ngrid
473         z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
474         zc(ig,nlay)=za(ig,nlay)*zu(ig,nlay)*z1(ig)
475         zd(ig,nlay)=zb(ig,nlay)*z1(ig)
476      ENDDO
477
478      DO ilay=nlay-1,1,-1
479         DO ig=1,ngrid
480            z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
481     $         zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
482            zc(ig,ilay)=(za(ig,ilay)*zu(ig,ilay)+
483     $         zb(ig,ilay+1)*zc(ig,ilay+1))*z1(ig)
484            zd(ig,ilay)=zb(ig,ilay)*z1(ig)
485         ENDDO
486      ENDDO
487
488      DO ig=1,ngrid
489         zu(ig,1)=zc(ig,1)
490      ENDDO
491      DO ilay=2,nlay
492         DO ig=1,ngrid
493            zu(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zu(ig,ilay-1)
494         ENDDO
495      ENDDO
496
497
498
499
500
501c-----------------------------------------------------------------------
502c   5. inversion pour l'implicite sur v
503c      --------------------------------
504
505c    ** l'equation est
506c       v(t+1) =  v(t) + dt * {(dv/dt)phys}(t)  +  dt * {(dv/dt)difv}(t+1)
507c       avec
508c       /zv/ = v(t) + dt * {(dv/dt)phys}(t)   (voir paragraphe 2.)
509c       et
510c       dt * {(dv/dt)difv}(t+1) = dt * {(d/dz)[ Kv (dv/dz) ]}(t+1)
511c       donc les entrees sont /zcdv/ pour la condition a la limite sol
512c       et /zkv/ = Kv
513
514      DO ig=1,ngrid
515         z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
516         zc(ig,nlay)=za(ig,nlay)*zv(ig,nlay)*z1(ig)
517         zd(ig,nlay)=zb(ig,nlay)*z1(ig)
518      ENDDO
519
520      DO ilay=nlay-1,1,-1
521         DO ig=1,ngrid
522            z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
523     $         zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
524            zc(ig,ilay)=(za(ig,ilay)*zv(ig,ilay)+
525     $         zb(ig,ilay+1)*zc(ig,ilay+1))*z1(ig)
526            zd(ig,ilay)=zb(ig,ilay)*z1(ig)
527         ENDDO
528      ENDDO
529
530      DO ig=1,ngrid
531         zv(ig,1)=zc(ig,1)
532      ENDDO
533      DO ilay=2,nlay
534         DO ig=1,ngrid
535            zv(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zv(ig,ilay-1)
536         ENDDO
537      ENDDO
538
539
540
541
542
543c-----------------------------------------------------------------------
544c   6. inversion pour l'implicite sur h sans oublier le couplage
545c      avec le sol (conduction)
546c      ------------------------
547
548c    ** l'equation est
549c       h(t+1) =  h(t) + dt * {(dh/dt)phys}(t)  +  dt * {(dh/dt)difv}(t+1)
550c       avec
551c       /zh/ = h(t) + dt * {(dh/dt)phys}(t)   (voir paragraphe 2.)
552c       et
553c       dt * {(dh/dt)difv}(t+1) = dt * {(d/dz)[ Kh (dh/dz) ]}(t+1)
554c       donc les entrees sont /zcdh/ pour la condition de raccord au sol
555c       et /zkh/ = Kh
556c       -------------
557
[473]558c Mass variation scheme:
[2274]559      zb(1:ngrid,2:nlay)=zkh(1:ngrid,2:nlay)*zb0(1:ngrid,2:nlay)
560      zb(1:ngrid,1)=zcdh(1:ngrid)*zb0(1:ngrid,1)
[38]561
[473]562c on initialise dm c
563     
564      pdtc(:,:)=0.
565      zt(:,:)=0.
566      dmice(:,:)=0.
[38]567
568c    ** calcul de (d Planck / dT) a la temperature d'interface
569c       ------------------------------------------------------
570
571      z4st=4.*5.67e-8*ptimestep
[544]572      IF (tke_heat_flux .eq. 0.) THEN
[38]573      DO ig=1,ngrid
574         zdplanck(ig)=z4st*pemis(ig)*ptsrf(ig)*ptsrf(ig)*ptsrf(ig)
575      ENDDO
[544]576      ELSE
577         zdplanck(:)=0.
578      ENDIF
[38]579
[473]580! calcul de zc et zd pour la couche top en prenant en compte le terme
[2080]581! de variation de masse (on fait une boucle pour que \E7a converge)
[473]582
583! Identification des points de grilles qui ont besoin de la correction
584
585      llnt(:)=1
[1236]586      IF (.not.turb_resolved) THEN
[884]587      IF (callcond) THEN
588       DO ig=1,ngrid
[473]589         DO l=1,nlay
590            if(zh(ig,l) .lt. zhcond(ig,l)) then
591               llnt(ig)=300 
592! 200 and 100 do not go beyond month 9 with normal dissipation
593               goto 5
594            endif
595         ENDDO
[884]5965      continue
597       ENDDO
598      ENDIF
[473]599
[529]600      ENDIF
601
[473]602      DO ig=1,ngrid
603
604! Initialization of z1 and zd, which do not depend on dmice
605
606      z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
607      zd(ig,nlay)=zb(ig,nlay)*z1(ig)
608
609      DO ilay=nlay-1,1,-1
610          z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
611     $        zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
612          zd(ig,ilay)=zb(ig,ilay)*z1(ig)
613      ENDDO
614
615! Convergence loop
616
617      DO j=1,llnt(ig)
618
619            z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
620            zc(ig,nlay)=za(ig,nlay)*zh(ig,nlay)
621     &      -betam(ig,nlay)*dmice(ig,nlay)
622            zc(ig,nlay)=zc(ig,nlay)*z1(ig)
623!            zd(ig,nlay)=zb(ig,nlay)*z1(ig)
624
625! calcul de zc et zd pour les couches du haut vers le bas
626
627             DO ilay=nlay-1,1,-1
628               z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
629     $            zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
630               zc(ig,ilay)=(za(ig,ilay)*zh(ig,ilay)+
631     $            zb(ig,ilay+1)*zc(ig,ilay+1)-
632     $            betam(ig,ilay)*dmice(ig,ilay))*z1(ig)
633!               zd(ig,ilay)=zb(ig,ilay)*z1(ig)
634            ENDDO
635
[38]636c    ** calcul de la temperature_d'interface et de sa tendance.
637c       on ecrit que la somme des flux est nulle a l'interface
638c       a t + \delta t,
639c       c'est a dire le flux radiatif a {t + \delta t}
640c       + le flux turbulent a {t + \delta t}
641c            qui s'ecrit K (T1-Tsurf) avec T1 = d1 Tsurf + c1
642c            (notation K dt = /cpp*b/)       
643c       + le flux dans le sol a t
644c       + l'evolution du flux dans le sol lorsque la temperature d'interface
645c       passe de sa valeur a t a sa valeur a {t + \delta t}.
646c       ----------------------------------------------------
647
648         z1(ig)=pcapcal(ig)*ptsrf(ig)+cpp*zb(ig,1)*zc(ig,1)
649     s     +zdplanck(ig)*ptsrf(ig)+ pfluxsrf(ig)*ptimestep
650         z2(ig)= pcapcal(ig)+cpp*zb(ig,1)*(1.-zd(ig,1))+zdplanck(ig)
651         ztsrf2(ig)=z1(ig)/z2(ig)
[473]652!         pdtsrf(ig)=(ztsrf2(ig)-ptsrf(ig))/ptimestep  !incremented outside loop
653            zhs(ig,1)=zc(ig,1)+zd(ig,1)*ztsrf2(ig)
[38]654
655c    ** et a partir de la temperature au sol on remonte
656c       -----------------------------------------------
657
[473]658            DO ilay=2,nlay
659               zhs(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zhs(ig,ilay-1)
660            ENDDO
661            DO ilay=1,nlay
662               zt(ig,ilay)=zhs(ig,ilay)*ppopsk(ig,ilay)
663            ENDDO
664
665c      Condensation/sublimation in the atmosphere
666c      ------------------------------------------
667c      (computation of zcondicea and dmice)
668
[1996]669      IF (.NOT. co2clouds) then
670       DO l=nlay , 1, -1
[473]671           IF(zt(ig,l).LT.ztcond(ig,l)) THEN
672               pdtc(ig,l)=(ztcond(ig,l) - zt(ig,l))/ptimestep
673               zcondicea(ig,l)=(pplev(ig,l)-pplev(ig,l+1))
674     &                        *ccond*pdtc(ig,l)
675              dmice(ig,l)= dmice(ig,l) + zcondicea(ig,l)*ptimestep
676            END IF
[1996]677       ENDDO
678      ELSE
679       DO l=nlay , 1, -1
680         zcondicea(ig,l)= 0.!pcondicea_co2microp(ig,l)*
681c     &                        (pplev(ig,l) - pplev(ig,l+1))/g
682         dmice(ig,l)= 0.!dmice(ig,l) + zcondicea(ig,l)*ptimestep
683         pdtc(ig,l)=0.
684       ENDDO   
685      ENDIF
686     
687       ENDDO!of Do j=1,XXX
[38]688
[1996]689      ENDDO   !of Do ig=1,ngrid
[38]690
[473]691      pdtsrf(:)=(ztsrf2(:)-ptsrf(:))/ptimestep
[669]692     
[660]693      DO ig=1,ngrid  ! computing sensible heat flux (atm => surface)
694         sensibFlux(ig)=cpp*zb(ig,1)/ptimestep*(zhs(ig,1)-ztsrf2(ig))
695      ENDDO
[473]696
[38]697c-----------------------------------------------------------------------
698c   TRACERS
699c   -------
700
701      if(tracer) then
702
703c     Using the wind modified by friction for lifting and  sublimation
704c     ----------------------------------------------------------------
705
[529]706!     This is computed above and takes into account surface-atmosphere flux
707!     enhancement by subgrid gustiness and atmospheric-stability related
708!     variations of transfer coefficients.
709
710!        DO ig=1,ngrid
711!          zu2(ig)=zu(ig,1)*zu(ig,1)+zv(ig,1)*zv(ig,1)
712!          zcdv(ig)=zcdv_true(ig)*sqrt(zu2(ig))
713!          zcdh(ig)=zcdh_true(ig)*sqrt(zu2(ig))
714!        ENDDO
715
[38]716c       Calcul du flux vertical au bas de la premiere couche (dust) :
717c       -----------------------------------------------------------
[1047]718        do ig=1,ngrid 
[38]719          rho(ig) = zb0(ig,1) /ptimestep
720c          zb(ig,1) = 0.
721        end do
722c       Dust lifting:
723        if (lifting) then
[310]724#ifndef MESOSCALE
[38]725           if (doubleq.AND.submicron) then
726             do ig=1,ngrid
727c              if(co2ice(ig).lt.1) then
728                 pdqsdif(ig,igcm_dust_mass) =
729     &             -alpha_lift(igcm_dust_mass) 
730                 pdqsdif(ig,igcm_dust_number) =
731     &             -alpha_lift(igcm_dust_number) 
732                 pdqsdif(ig,igcm_dust_submicron) =
733     &             -alpha_lift(igcm_dust_submicron)
734c              end if
735             end do
736           else if (doubleq) then
[1974]737             if (dustinjection.eq.0) then !injection scheme 0 (old)
738                                          !or 2 (injection in CL)
739              do ig=1,ngrid
[1455]740               if(co2ice(ig).lt.1) then ! pas de soulevement si glace CO2
[38]741                 pdqsdif(ig,igcm_dust_mass) =
742     &             -alpha_lift(igcm_dust_mass) 
743                 pdqsdif(ig,igcm_dust_number) =
[520]744     &             -alpha_lift(igcm_dust_number)
745               end if
[1974]746              end do
747             elseif(dustinjection.eq.1)then ! dust injection scheme = 1 injection from surface
748              do ig=1,ngrid
749               if(co2ice(ig).lt.1) then ! pas de soulevement si glace CO2
[2160]750                IF((ti_injection_sol.LE.local_time(ig)).and.
751     &                  (local_time(ig).LE.tf_injection_sol)) THEN
[1974]752                  if (rdstorm) then !Rocket dust storm scheme
753                        pdqsdif(ig,igcm_stormdust_mass) =
754     &                          -alpha_lift(igcm_stormdust_mass)
755     &                          *dustliftday(ig)
756                        pdqsdif(ig,igcm_stormdust_number) =
757     &                          -alpha_lift(igcm_stormdust_number)
758     &                          *dustliftday(ig)
759                        pdqsdif(ig,igcm_dust_mass)= 0.
760                        pdqsdif(ig,igcm_dust_number)= 0.
761                  else
762                        pdqsdif(ig,igcm_dust_mass)=
763     &                          -dustliftday(ig)*
764     &                          alpha_lift(igcm_dust_mass)               
765                        pdqsdif(ig,igcm_dust_number)=
766     &                          -dustliftday(ig)*
767     &                          alpha_lift(igcm_dust_number)
768                  endif
769                  if (submicron) then
770                        pdqsdif(ig,igcm_dust_submicron) = 0.
771                  endif ! if (submicron)
772                ELSE ! outside dust injection time frame
[2080]773                  pdqsdif(ig,igcm_dust_mass)= 0.
774                  pdqsdif(ig,igcm_dust_number)= 0.
775                  if (rdstorm) then     
[1974]776                        pdqsdif(ig,igcm_stormdust_mass)= 0.
777                        pdqsdif(ig,igcm_stormdust_number)= 0.
[2080]778                  end if
[1974]779                ENDIF
780               
781                end if ! of if(co2ice(ig).lt.1)
782              end do
783             endif ! end if dustinjection
[38]784           else if (submicron) then
785             do ig=1,ngrid
786                 pdqsdif(ig,igcm_dust_submicron) =
787     &             -alpha_lift(igcm_dust_submicron)
788             end do
789           else
[1236]790#endif
[38]791            call dustlift(ngrid,nlay,nq,rho,zcdh_true,zcdh,co2ice,
792     &                    pdqsdif)
[1236]793#ifndef MESOSCALE
[38]794           endif !doubleq.AND.submicron
[310]795#endif
[38]796        else
797           pdqsdif(1:ngrid,1:nq) = 0.
798        end if
799
800c       OU calcul de la valeur de q a la surface (water)  :
801c       ----------------------------------------
802
803c      Inversion pour l'implicite sur q
[2515]804c      Cas des traceurs qui ne sont pas h2o_vap
805c      h2o_vap est traite plus loin avec un sous pas de temps
806c      hdo_vap est traite ensuite car dependant de h2o_vap
[38]807c       --------------------------------
[2515]808
809        do iq=1,nq  !for all tracers except water vapor
810          if ((.not. water).or.(.not. iq.eq.igcm_h2o_vap).or.
811     &       (.not. iq.eq.igcm_hdo_vap)) then
812
813
[2274]814          zb(1:ngrid,2:nlay)=zkh(1:ngrid,2:nlay)*zb0(1:ngrid,2:nlay)
[2515]815          zb(1:ngrid,1)=0
[38]816
[2515]817          DO ig=1,ngrid
818               z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
819               zc(ig,nlay)=za(ig,nlay)*zq(ig,nlay,iq)*z1(ig)
820               zd(ig,nlay)=zb(ig,nlay)*z1(ig)
821          ENDDO
822 
823          DO ilay=nlay-1,2,-1
824               DO ig=1,ngrid
825                z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
826     $           zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
827                zc(ig,ilay)=(za(ig,ilay)*zq(ig,ilay,iq)+
828     $           zb(ig,ilay+1)*zc(ig,ilay+1))*z1(ig)
829                zd(ig,ilay)=zb(ig,ilay)*z1(ig)
830               ENDDO
831          ENDDO
832
833            if ((iq.eq.igcm_h2o_ice)
834     $      .or. (hdo.and.(iq.eq.igcm_hdo_ice) )) then
835
836               DO ig=1,ngrid
837                   z1(ig)=1./(za(ig,1)+zb(ig,1)+
838     $              zb(ig,2)*(1.-zd(ig,2)))
839                   zc(ig,1)=(za(ig,1)*zq(ig,1,iq)+
840     $              zb(ig,2)*zc(ig,2)) *z1(ig)  !special case h2o_ice
841               ENDDO
842            else ! every other tracer
843               DO ig=1,ngrid
844                   z1(ig)=1./(za(ig,1)+zb(ig,1)+
845     $              zb(ig,2)*(1.-zd(ig,2)))
846                   zc(ig,1)=(za(ig,1)*zq(ig,1,iq)+
847     $              zb(ig,2)*zc(ig,2) +
848     $             (-pdqsdif(ig,iq)) *ptimestep) *z1(ig)  !tracer flux from surface
849               ENDDO
850            endif !((iq.eq.igcm_h2o_ice)
851c         Starting upward calculations for simple mixing of tracer (dust)
852          DO ig=1,ngrid
853             zq(ig,1,iq)=zc(ig,1)
854             DO ilay=2,nlay
855               zq(ig,ilay,iq)=zc(ig,ilay)+zd(ig,ilay)*zq(ig,ilay-1,iq)
856             ENDDO
857          ENDDO
858          endif! ((.not. water).or.(.not. iq.eq.igcm_h2o_vap)) then
859        enddo ! of do iq=1,nq
860
861c --------- h2o_vap --------------------------------
862
863
864c      Traitement de la vapeur d'eau h2o_vap
865c      Utilisation d'un sous pas de temps afin
866c      de decrire le flux de chaleur latente
867
868
869        do iq=1,nq 
[2312]870          if ((water).and.(iq.eq.igcm_h2o_vap)) then
[2515]871
872
873           DO ig=1,ngrid
874             zqsurf(ig)=pqsurf(ig,igcm_h2o_ice)
875           ENDDO ! ig=1,ngrid
876
877c          make_tsub : sous pas de temps adaptatif
878c          la subroutine est a la fin du fichier
879
880           call make_tsub(ngrid,pdtsrf,zqsurf,
881     &                    ptimestep,dtmax,watercaptag,
882     &                    nsubtimestep)
883
884c           Calculation for turbulent exchange with the surface (for ice)
885c           initialization of ztsrf, which is surface temperature in
886c           the subtimestep.
887           DO ig=1,ngrid
888            subtimestep = ptimestep/nsubtimestep(ig)
889            ztsrf(ig)=ptsrf(ig)  !  +pdtsrf(ig)*subtimestep
890
891c           Debut du sous pas de temps
892
893            DO tsub=1,nsubtimestep(ig)
894
895c           C'est parti !
896
897             zb(1:ngrid,2:nlay)=zkh(1:ngrid,2:nlay)*zb0(1:ngrid,2:nlay)
898     &                     /float(nsubtimestep(ig))
[2274]899             zb(1:ngrid,1)=zcdv(1:ngrid)*zb0(1:ngrid,1)
[2515]900     &                     /float(nsubtimestep(ig))
[2274]901             zb(1:ngrid,1)=dryness(1:ngrid)*zb(1:ngrid,1)
[2515]902             
903             z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
904             zc(ig,nlay)=za(ig,nlay)*zq(ig,nlay,iq)*z1(ig)
905             zd(ig,nlay)=zb(ig,nlay)*z1(ig)
906             DO ilay=nlay-1,2,-1
907                 z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
908     $            zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
909                 zc(ig,ilay)=(za(ig,ilay)*zq(ig,ilay,iq)+
910     $            zb(ig,ilay+1)*zc(ig,ilay+1))*z1(ig)
911                 zd(ig,ilay)=zb(ig,ilay)*z1(ig)
912             ENDDO 
913             z1(ig)=1./(za(ig,1)+zb(ig,1)+
914     $          zb(ig,2)*(1.-zd(ig,2)))
915             zc(ig,1)=(za(ig,1)*zq(ig,1,iq)+
916     $        zb(ig,2)*zc(ig,2)) * z1(ig)
[38]917
[2531]918             call watersat(1,ztsrf(ig),pplev(ig,1),qsat(ig))
[2515]919             old_h2o_vap(ig)=zq(ig,1,igcm_h2o_vap)
920             zd(ig,1)=zb(ig,1)*z1(ig)
921             zq1temp(ig)=zc(ig,1)+ zd(ig,1)*qsat(ig)
922
923             zdqsdif(ig)=rho(ig)*dryness(ig)*zcdv(ig)
924     &                      *(zq1temp(ig)-qsat(ig))
[2587]925c             write(*,*)'subliming more than available frost:  qsurf!'
926             if(.not.watercaptag(ig)) then
927               if ((-zdqsdif(ig)*subtimestep)
[2515]928     &            .gt.(zqsurf(ig))) then
[2587]929c             pdqsdif > 0 : ice condensing
930c             pdqsdif < 0 : ice subliming
931c             write(*,*) "subliming more than available frost:  qsurf!"
[2515]932                  zdqsdif(ig)=
933     &                        -zqsurf(ig)/subtimestep
934c                write(*,*)'flux vers le sol=',pdqsdif(ig,nq)
935                 z1(ig)=1./(za(ig,1)+ zb(ig,2)*(1.-zd(ig,2)))
936                 zc(ig,1)=(za(ig,1)*zq(ig,1,igcm_h2o_vap)+
[2587]937     $           zb(ig,2)*zc(ig,2) +
938     $           (-zdqsdif(ig)) *subtimestep) *z1(ig)
[2515]939                 zq1temp(ig)=zc(ig,1)
940               endif  !if .not.watercaptag(ig)
941             endif ! if sublim more than surface
942
943c             Starting upward calculations for water :
944c             Actualisation de h2o_vap dans le premier niveau
945             zq(ig,1,igcm_h2o_vap)=zq1temp(ig)
[2587]946
947c    Take into account the H2O latent heat impact on the surface temperature
[2515]948              if (latentheat_surfwater) then
949                lh=(2834.3-0.28*(ztsrf(ig)-To)-
950     &              0.004*(ztsrf(ig)-To)*(ztsrf(ig)-To))*1.e+3
[2587]951                zdtsrf(ig)=  zdqsdif(ig)*lh /pcapcal(ig)
[2515]952              endif ! (latentheat_surfwater) then
953
[2587]954              DO ilay=2,nlay
955                zq(ig,ilay,iq)=zc(ig,ilay)+zd(ig,ilay)*zq(ig,ilay-1,iq)
956              ENDDO
[2515]957
958c             Subtimestep water budget :
959
960              ztsrf(ig) = ztsrf(ig)+(pdtsrf(ig) + zdtsrf(ig))
961     &                          *subtimestep
962              zqsurf(ig)= zqsurf(ig)+(
963     &                       zdqsdif(ig))*subtimestep
964
965
966c             We ensure that surface temperature can't rise above the solid domain if there
967c             is still ice on the surface (oldschool)
968               if(zqsurf(ig)
969     &           +zdqsdif(ig)*subtimestep
970     &           .gt.frost_albedo_threshold) ! if there is still ice, T cannot exceed To
971     &           zdtsrf(ig) = min(zdtsrf(ig),(To-ztsrf(ig))/subtimestep) ! ice melt case
972
973
[2587]974
[2515]975c             Fin du sous pas de temps
976            ENDDO ! tsub=1,nsubtimestep
977
[2587]978c             Integration of subtimestep temp and water budget :
979c             (btw could also compute the post timestep temp and ice
980c             by simply adding the subtimestep trend instead of this)
[2515]981            pdtsrf(ig)= (ztsrf(ig) -
982     &                     ptsrf(ig))/ptimestep
983            pdqsdif(ig,igcm_h2o_ice)=
984     &                  (zqsurf(ig)- pqsurf(ig,igcm_h2o_ice))/ptimestep
985
[2587]986c if subliming more than qsurf(ice) and on watercaptag, water
987c sublimates from watercap reservoir (dwatercap_dif is <0)
988            if(watercaptag(ig)) then
989              if ((-pdqsdif(ig,igcm_h2o_ice)*ptimestep)
990     &           .gt.(pqsurf(ig,igcm_h2o_ice))) then
991                 dwatercap_dif(ig)=pdqsdif(ig,igcm_h2o_ice)+
992     &                     pqsurf(ig,igcm_h2o_ice)/ptimestep
993                 pdqsdif(ig,igcm_h2o_ice)=
994     &                   - pqsurf(ig,igcm_h2o_ice)/ptimestep
995              endif! ((-pdqsdif(ig)*ptimestep)
996            endif !(watercaptag(ig)) then
997
[2515]998           ENDDO ! of DO ig=1,ngrid
999          END IF ! of IF ((water).and.(iq.eq.igcm_h2o_vap))
1000
1001c --------- end of h2o_vap ----------------------------
1002
1003c --------- hdo_vap -----------------------------------
1004
1005c    hdo_ice has already been with along h2o_ice
1006c    amongst "normal" tracers (ie not h2o_vap)
1007
1008         if (hdo.and.(iq.eq.igcm_hdo_vap)) then
1009          zb(1:ngrid,2:nlay)=zkh(1:ngrid,2:nlay)*zb0(1:ngrid,2:nlay)
1010          zb(1:ngrid,1)=0
1011
[38]1012          DO ig=1,ngrid
1013               z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
1014               zc(ig,nlay)=za(ig,nlay)*zq(ig,nlay,iq)*z1(ig)
1015               zd(ig,nlay)=zb(ig,nlay)*z1(ig)
1016          ENDDO
[2515]1017
[38]1018          DO ilay=nlay-1,2,-1
1019               DO ig=1,ngrid
1020                z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
1021     $           zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
1022                zc(ig,ilay)=(za(ig,ilay)*zq(ig,ilay,iq)+
1023     $           zb(ig,ilay+1)*zc(ig,ilay+1))*z1(ig)
1024                zd(ig,ilay)=zb(ig,ilay)*z1(ig)
1025               ENDDO
1026          ENDDO
1027
[2312]1028            CALL hdo_surfex(ngrid,nlay,nq,ptimestep,
[2593]1029     &                      zt,pplay,zq,pqsurf,
1030     &                      old_h2o_vap,qsat,pdqsdif,dwatercap_dif,
[2312]1031     &                      hdoflux)
1032            DO ig=1,ngrid
1033                z1(ig)=1./(za(ig,1)+zb(ig,1)+
1034     $           zb(ig,2)*(1.-zd(ig,2)))
1035                zc(ig,1)=(za(ig,1)*zq(ig,1,iq)+
1036     $         zb(ig,2)*zc(ig,2) +
1037     $        (-hdoflux(ig)) *ptimestep) *z1(ig)  !tracer flux from surface
1038            ENDDO
1039
[38]1040            DO ig=1,ngrid
[2515]1041              zq(ig,1,iq)=zc(ig,1)
1042              DO ilay=2,nlay
1043                zq(ig,ilay,iq)=zc(ig,ilay)+zd(ig,ilay)*zq(ig,ilay-1,iq)
1044              ENDDO
[38]1045            ENDDO
[2515]1046         endif ! (hdo.and.(iq.eq.igcm_hdo_vap))
[2312]1047
[2515]1048c --------- end of hdo ----------------------------
[38]1049
[2515]1050        enddo ! of do iq=1,nq
1051      end if ! of if(tracer)
[38]1052
[2515]1053c --------- end of tracers  ----------------------------
[2312]1054
[2260]1055
[2312]1056C       Diagnostic output for HDO
1057        if (hdo) then
1058            CALL WRITEDIAGFI(ngrid,'hdoflux',
1059     &                       'hdoflux',
1060     &                       ' ',2,hdoflux) 
1061            CALL WRITEDIAGFI(ngrid,'h2oflux',
1062     &                       'h2oflux',
1063     &                       ' ',2,h2oflux)
1064        endif
1065
[38]1066c-----------------------------------------------------------------------
1067c   8. calcul final des tendances de la diffusion verticale
1068c      ----------------------------------------------------
1069
1070      DO ilev = 1, nlay
1071         DO ig=1,ngrid
1072            pdudif(ig,ilev)=(    zu(ig,ilev)-
1073     $      (pu(ig,ilev)+pdufi(ig,ilev)*ptimestep)    )/ptimestep
1074            pdvdif(ig,ilev)=(    zv(ig,ilev)-
1075     $      (pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep)    )/ptimestep
[473]1076            hh = ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep
1077     $  + (latcond*dmice(ig,ilev)/cpp)/ppopsk(ig,ilev)
1078            pdhdif(ig,ilev)=( zhs(ig,ilev)- hh )/ptimestep
[38]1079         ENDDO
1080      ENDDO
1081
1082      if (tracer) then
1083        DO iq = 1, nq
1084          DO ilev = 1, nlay
1085            DO ig=1,ngrid
1086              pdqdif(ig,ilev,iq)=(zq(ig,ilev,iq)-
1087     $      (pq(ig,ilev,iq) + pdqfi(ig,ilev,iq)*ptimestep))/ptimestep
1088            ENDDO
1089          ENDDO
1090        ENDDO
1091      end if
1092
1093c    ** diagnostique final
1094c       ------------------
1095
1096      IF(lecrit) THEN
1097         PRINT*,'In vdif'
1098         PRINT*,'Ts (t) and Ts (t+st)'
1099         WRITE(*,'(a10,3a15)')
1100     s   'theta(t)','theta(t+dt)','u(t)','u(t+dt)'
1101         PRINT*,ptsrf(ngrid/2+1),ztsrf2(ngrid/2+1)
1102         DO ilev=1,nlay
1103            WRITE(*,'(4f15.7)')
[473]1104     s      ph(ngrid/2+1,ilev),zhs(ngrid/2+1,ilev),
[38]1105     s      pu(ngrid/2+1,ilev),zu(ngrid/2+1,ilev)
1106
1107         ENDDO
1108      ENDIF
1109
[1036]1110      END SUBROUTINE vdifc
[1969]1111
[2312]1112c====================================
1113
[2515]1114      SUBROUTINE make_tsub(naersize,dtsurf,qsurf,ptimestep,
1115     $                     dtmax,watercaptag,ntsub)
[2312]1116
[2515]1117c Pas de temps adaptatif en estimant le taux de sublimation
1118c et en adaptant avec un critere "dtmax" du chauffage a accomoder
1119c dtmax est regle empiriquement (pour l'instant) a 0.5 K
1120
1121      integer,intent(in) :: naersize
1122      real,intent(in) :: dtsurf(naersize)
1123      real,intent(in) :: qsurf(naersize)
1124      logical,intent(in) :: watercaptag(naersize)
1125      real,intent(in) :: ptimestep
1126      real,intent(in)  :: dtmax
1127      real  :: ztsub(naersize)
1128      integer  :: i
1129      integer,intent(out) :: ntsub(naersize)
1130
1131      do i=1,naersize
1132        if ((qsurf(i).eq.0).and.
1133     &      (.not.watercaptag(i))) then
1134          ntsub(i) = 1
1135        else
1136          ztsub(i) = ptimestep * dtsurf(i) / dtmax
1137          ntsub(i) = ceiling(abs(ztsub(i)))
1138        endif ! (qsurf(i).eq.0) then
1139c     
1140c       write(78,*), dtsurf*ptimestep, dtsurf, ntsub
1141      enddo! 1=1,ngrid
1142
1143
1144
1145      END SUBROUTINE make_tsub
[1969]1146      END MODULE vdifc_mod
Note: See TracBrowser for help on using the repository browser.