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

Last change on this file since 2947 was 2934, checked in by emillour, 21 months ago

Mars PCM:
A first step in cleaning up outputs and adding field definitions for XIOS.
EM

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