source: trunk/LMDZ.MARS/libf/phymars/vdifc.F @ 226

Last change on this file since 226 was 224, checked in by emillour, 14 years ago

Mars GCM:

Implemented using 'z0' roughness length map (important: 'z0' reference

field is in datafile surface.nc, which has also been updated).

  • made z0 a z0(ngridmx) array and moved 'z0' from 'planete.h' to 'surfdat.h'; added a 'z0_default' (common in surfdat.h) corresponding to the 'control' array value (contole(19) in startfi.nc).
  • adapted 'tabfi.F' to use 'z0_default'.
  • adapted 'phyetat0.F' to look for a 'z0' field in startfi.nc. If 'z0' is not found in the startfi.nc file, then the uniform default value (z0_default) is used.
  • modified 'physdem1.F' to write 'z0' field to restart.nc
  • adapted use of z0() in 'physiq.F' (diagnostic computation of surface stress), 'vdifc.F' and 'vdif_cd.F'.
  • adapted 'dustdevil.F' to use 'z0_default'.
  • 'testphys1d.F' now uses 'z0_default', and the value to use can be set in run.def (with "z0=TheValueYouWant?").
  • modified 'datareadnc.F' to load reference map of 'z0' from surface.nc, and added a 'z0' option in 'newstart.F' to force a uniform value of z0. Note that the use of the z0 map is automatic when using newstart, but only when it loads a start_archive.nc file.
File size: 20.4 KB
Line 
1      SUBROUTINE vdifc(ngrid,nlay,nq,co2ice,ppopsk,
2     $                ptimestep,pcapcal,lecrit,
3     $                pplay,pplev,pzlay,pzlev,pz0,
4     $                pu,pv,ph,pq,ptsrf,pemis,pqsurf,
5     $                pdufi,pdvfi,pdhfi,pdqfi,pfluxsrf,
6     $                pdudif,pdvdif,pdhdif,pdtsrf,pq2,
7     $                pdqdif,pdqsdif)
8      IMPLICIT NONE
9
10c=======================================================================
11c
12c   subject:
13c   --------
14c   Turbulent diffusion (mixing) for potential T, U, V and tracer
15c
16c   Shema implicite
17c   On commence par rajouter au variables x la tendance physique
18c   et on resoult en fait:
19c      x(t+1) =  x(t) + dt * (dx/dt)phys(t)  +  dt * (dx/dt)difv(t+1)
20c
21c   author:
22c   ------
23c      Hourdin/Forget/Fournier
24c=======================================================================
25
26c-----------------------------------------------------------------------
27c   declarations:
28c   -------------
29
30#include "dimensions.h"
31#include "dimphys.h"
32#include "comcstfi.h"
33#include "callkeys.h"
34#include "surfdat.h"
35#include "comgeomfi.h"
36#include "tracer.h"
37
38#include "watercap.h"
39c
40c   arguments:
41c   ----------
42
43      INTEGER ngrid,nlay
44      REAL ptimestep
45      REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1)
46      REAL pzlay(ngrid,nlay),pzlev(ngrid,nlay+1)
47      REAL pu(ngrid,nlay),pv(ngrid,nlay),ph(ngrid,nlay)
48      REAL ptsrf(ngrid),pemis(ngrid)
49      REAL pdufi(ngrid,nlay),pdvfi(ngrid,nlay),pdhfi(ngrid,nlay)
50      REAL pfluxsrf(ngrid)
51      REAL pdudif(ngrid,nlay),pdvdif(ngrid,nlay),pdhdif(ngrid,nlay)
52      REAL pdtsrf(ngrid),pcapcal(ngrid)
53      REAL pq2(ngrid,nlay+1)
54
55c    Argument added for condensation:
56      REAL co2ice (ngrid), ppopsk(ngrid,nlay)
57      logical lecrit
58
59      REAL pz0(ngridmx) ! surface roughness length (m)
60
61c    Traceurs :
62      integer nq
63      REAL pqsurf(ngrid,nq)
64      real pq(ngrid,nlay,nq), pdqfi(ngrid,nlay,nq)
65      real pdqdif(ngrid,nlay,nq)
66      real pdqsdif(ngrid,nq)
67     
68c   local:
69c   ------
70
71      INTEGER ilev,ig,ilay,nlev
72
73      REAL z4st,zdplanck(ngridmx)
74      REAL zkv(ngridmx,nlayermx+1),zkh(ngridmx,nlayermx+1)
75      REAL zcdv(ngridmx),zcdh(ngridmx)
76      REAL zcdv_true(ngridmx),zcdh_true(ngridmx)
77      REAL zu(ngridmx,nlayermx),zv(ngridmx,nlayermx)
78      REAL zh(ngridmx,nlayermx)
79      REAL ztsrf2(ngridmx)
80      REAL z1(ngridmx),z2(ngridmx)
81      REAL za(ngridmx,nlayermx),zb(ngridmx,nlayermx)
82      REAL zb0(ngridmx,nlayermx)
83      REAL zc(ngridmx,nlayermx),zd(ngridmx,nlayermx)
84      REAL zcst1
85      REAL zu2
86
87      EXTERNAL SSUM,SCOPY
88      REAL SSUM
89      LOGICAL firstcall
90      SAVE firstcall
91
92c     variable added for CO2 condensation:
93c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
94      REAL hh , zhcond(ngridmx,nlayermx)
95      REAL latcond,tcond1mb
96      REAL acond,bcond
97      SAVE acond,bcond
98      DATA latcond,tcond1mb/5.9e5,136.27/
99
100c    Tracers :
101c    ~~~~~~~
102      INTEGER iq
103      REAL zq(ngridmx,nlayermx,nqmx)
104      REAL zq1temp(ngridmx)
105      REAL rho(ngridmx) ! near surface air density
106      REAL qsat(ngridmx)
107      DATA firstcall/.true./
108
109      REAL kmixmin
110
111c    ** un petit test de coherence
112c       --------------------------
113
114      IF (firstcall) THEN
115         IF(ngrid.NE.ngridmx) THEN
116            PRINT*,'STOP dans vdifc'
117            PRINT*,'probleme de dimensions :'
118            PRINT*,'ngrid  =',ngrid
119            PRINT*,'ngridmx  =',ngridmx
120            STOP
121         ENDIF
122c        To compute: Tcond= 1./(bcond-acond*log(.0095*p)) (p in pascal)
123         bcond=1./tcond1mb
124         acond=r/latcond
125         PRINT*,'In vdifc: Tcond(P=1mb)=',tcond1mb,' Lcond=',latcond
126         PRINT*,'          acond,bcond',acond,bcond
127
128        firstcall=.false.
129      ENDIF
130
131
132
133
134
135c-----------------------------------------------------------------------
136c    1. initialisation
137c    -----------------
138
139      nlev=nlay+1
140
141c    ** calcul de rho*dz et dt*rho/dz=dt*rho**2 g/dp
142c       avec rho=p/RT=p/ (R Theta) (p/ps)**kappa
143c       ----------------------------------------
144
145      DO ilay=1,nlay
146         DO ig=1,ngrid
147            za(ig,ilay)=(pplev(ig,ilay)-pplev(ig,ilay+1))/g
148         ENDDO
149      ENDDO
150
151      zcst1=4.*g*ptimestep/(r*r)
152      DO ilev=2,nlev-1
153         DO ig=1,ngrid
154            zb0(ig,ilev)=pplev(ig,ilev)*
155     s      (pplev(ig,1)/pplev(ig,ilev))**rcp /
156     s      (ph(ig,ilev-1)+ph(ig,ilev))
157            zb0(ig,ilev)=zcst1*zb0(ig,ilev)*zb0(ig,ilev)/
158     s      (pplay(ig,ilev-1)-pplay(ig,ilev))
159         ENDDO
160      ENDDO
161      DO ig=1,ngrid
162                 zb0(ig,1)=ptimestep*pplev(ig,1)/(r*ptsrf(ig))
163      ENDDO
164
165c    ** diagnostique pour l'initialisation
166c       ----------------------------------
167
168      IF(lecrit) THEN
169         ig=ngrid/2+1
170         PRINT*,'Pression (mbar) ,altitude (km),u,v,theta, rho dz'
171         DO ilay=1,nlay
172            WRITE(*,'(6f11.5)')
173     s      .01*pplay(ig,ilay),.001*pzlay(ig,ilay),
174     s      pu(ig,ilay),pv(ig,ilay),ph(ig,ilay),za(ig,ilay)
175         ENDDO
176         PRINT*,'Pression (mbar) ,altitude (km),zb'
177         DO ilev=1,nlay
178            WRITE(*,'(3f15.7)')
179     s      .01*pplev(ig,ilev),.001*pzlev(ig,ilev),
180     s      zb0(ig,ilev)
181         ENDDO
182      ENDIF
183
184c     Potential Condensation temperature:
185c     -----------------------------------
186
187c     if (callcond) then
188c       DO ilev=1,nlay
189c         DO ig=1,ngrid
190c           zhcond(ig,ilev) =
191c    &  (1./(bcond-acond*log(.0095*pplay(ig,ilev))))/ppopsk(ig,ilev)
192c         END DO
193c       END DO
194c     else
195        call zerophys(ngrid*nlay,zhcond)
196c     end if
197
198
199c-----------------------------------------------------------------------
200c   2. ajout des tendances physiques
201c      -----------------------------
202
203      DO ilev=1,nlay
204         DO ig=1,ngrid
205            zu(ig,ilev)=pu(ig,ilev)+pdufi(ig,ilev)*ptimestep
206            zv(ig,ilev)=pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep
207            zh(ig,ilev)=ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep
208            zh(ig,ilev)=max(zh(ig,ilev),zhcond(ig,ilev))
209         ENDDO
210      ENDDO
211      if(tracer) then
212        DO iq =1, nq
213         DO ilev=1,nlay
214           DO ig=1,ngrid
215              zq(ig,ilev,iq)=pq(ig,ilev,iq)+pdqfi(ig,ilev,iq)*ptimestep
216           ENDDO
217         ENDDO
218        ENDDO
219      end if
220
221c-----------------------------------------------------------------------
222c   3. schema de turbulence
223c      --------------------
224
225c    ** source d'energie cinetique turbulente a la surface
226c       (condition aux limites du schema de diffusion turbulente
227c       dans la couche limite
228c       ---------------------
229
230      CALL vdif_cd( ngrid,nlay,pz0,g,pzlay,pu,pv,ptsrf,ph
231     &             ,zcdv_true,zcdh_true)
232      DO ig=1,ngrid
233        zu2=pu(ig,1)*pu(ig,1)+pv(ig,1)*pv(ig,1)
234        zcdv(ig)=zcdv_true(ig)*sqrt(zu2)
235        zcdh(ig)=zcdh_true(ig)*sqrt(zu2)
236      ENDDO
237
238c    ** schema de diffusion turbulente dans la couche limite
239c       ----------------------------------------------------
240
241        CALL vdif_kc(ptimestep,g,pzlev,pzlay
242     &              ,pu,pv,ph,zcdv_true
243     &              ,pq2,zkv,zkh)
244
245      if ((doubleq).and.(ngrid.eq.1)) then
246        kmixmin = 80. !80.! minimum eddy mix coeff in 1D
247        do ilev=1,nlay
248          do ig=1,ngrid
249           zkh(ig,ilev) = max(kmixmin,zkh(ig,ilev))
250           zkv(ig,ilev) = max(kmixmin,zkv(ig,ilev))
251          end do
252        end do
253      end if
254
255c    ** diagnostique pour le schema de turbulence
256c       -----------------------------------------
257
258      IF(lecrit) THEN
259         PRINT*
260         PRINT*,'Diagnostic for the vertical turbulent mixing'
261         PRINT*,'Cd for momentum and potential temperature'
262
263         PRINT*,zcdv(ngrid/2+1),zcdh(ngrid/2+1)
264         PRINT*,'Mixing coefficient for momentum and pot.temp.'
265         DO ilev=1,nlay
266            PRINT*,zkv(ngrid/2+1,ilev),zkh(ngrid/2+1,ilev)
267         ENDDO
268      ENDIF
269
270
271
272
273c-----------------------------------------------------------------------
274c   4. inversion pour l'implicite sur u
275c      --------------------------------
276
277c    ** l'equation est
278c       u(t+1) =  u(t) + dt * {(du/dt)phys}(t)  +  dt * {(du/dt)difv}(t+1)
279c       avec
280c       /zu/ = u(t) + dt * {(du/dt)phys}(t)   (voir paragraphe 2.)
281c       et
282c       dt * {(du/dt)difv}(t+1) = dt * {(d/dz)[ Ku (du/dz) ]}(t+1)
283c       donc les entrees sont /zcdv/ pour la condition a la limite sol
284c       et /zkv/ = Ku
285 
286      CALL multipl((nlay-1)*ngrid,zkv(1,2),zb0(1,2),zb(1,2))
287      CALL multipl(ngrid,zcdv,zb0,zb)
288
289      DO ig=1,ngrid
290         z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
291         zc(ig,nlay)=za(ig,nlay)*zu(ig,nlay)*z1(ig)
292         zd(ig,nlay)=zb(ig,nlay)*z1(ig)
293      ENDDO
294
295      DO ilay=nlay-1,1,-1
296         DO ig=1,ngrid
297            z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
298     $         zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
299            zc(ig,ilay)=(za(ig,ilay)*zu(ig,ilay)+
300     $         zb(ig,ilay+1)*zc(ig,ilay+1))*z1(ig)
301            zd(ig,ilay)=zb(ig,ilay)*z1(ig)
302         ENDDO
303      ENDDO
304
305      DO ig=1,ngrid
306         zu(ig,1)=zc(ig,1)
307      ENDDO
308      DO ilay=2,nlay
309         DO ig=1,ngrid
310            zu(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zu(ig,ilay-1)
311         ENDDO
312      ENDDO
313
314
315
316
317
318c-----------------------------------------------------------------------
319c   5. inversion pour l'implicite sur v
320c      --------------------------------
321
322c    ** l'equation est
323c       v(t+1) =  v(t) + dt * {(dv/dt)phys}(t)  +  dt * {(dv/dt)difv}(t+1)
324c       avec
325c       /zv/ = v(t) + dt * {(dv/dt)phys}(t)   (voir paragraphe 2.)
326c       et
327c       dt * {(dv/dt)difv}(t+1) = dt * {(d/dz)[ Kv (dv/dz) ]}(t+1)
328c       donc les entrees sont /zcdv/ pour la condition a la limite sol
329c       et /zkv/ = Kv
330
331      DO ig=1,ngrid
332         z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
333         zc(ig,nlay)=za(ig,nlay)*zv(ig,nlay)*z1(ig)
334         zd(ig,nlay)=zb(ig,nlay)*z1(ig)
335      ENDDO
336
337      DO ilay=nlay-1,1,-1
338         DO ig=1,ngrid
339            z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
340     $         zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
341            zc(ig,ilay)=(za(ig,ilay)*zv(ig,ilay)+
342     $         zb(ig,ilay+1)*zc(ig,ilay+1))*z1(ig)
343            zd(ig,ilay)=zb(ig,ilay)*z1(ig)
344         ENDDO
345      ENDDO
346
347      DO ig=1,ngrid
348         zv(ig,1)=zc(ig,1)
349      ENDDO
350      DO ilay=2,nlay
351         DO ig=1,ngrid
352            zv(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zv(ig,ilay-1)
353         ENDDO
354      ENDDO
355
356
357
358
359
360c-----------------------------------------------------------------------
361c   6. inversion pour l'implicite sur h sans oublier le couplage
362c      avec le sol (conduction)
363c      ------------------------
364
365c    ** l'equation est
366c       h(t+1) =  h(t) + dt * {(dh/dt)phys}(t)  +  dt * {(dh/dt)difv}(t+1)
367c       avec
368c       /zh/ = h(t) + dt * {(dh/dt)phys}(t)   (voir paragraphe 2.)
369c       et
370c       dt * {(dh/dt)difv}(t+1) = dt * {(d/dz)[ Kh (dh/dz) ]}(t+1)
371c       donc les entrees sont /zcdh/ pour la condition de raccord au sol
372c       et /zkh/ = Kh
373c       -------------
374
375      CALL multipl((nlay-1)*ngrid,zkh(1,2),zb0(1,2),zb(1,2))
376      CALL multipl(ngrid,zcdh,zb0,zb)
377
378      DO ig=1,ngrid
379         z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
380         zc(ig,nlay)=za(ig,nlay)*zh(ig,nlay)*z1(ig)
381         zd(ig,nlay)=zb(ig,nlay)*z1(ig)
382      ENDDO
383
384      DO ilay=nlay-1,1,-1
385         DO ig=1,ngrid
386            z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
387     $         zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
388            zc(ig,ilay)=(za(ig,ilay)*zh(ig,ilay)+
389     $         zb(ig,ilay+1)*zc(ig,ilay+1))*z1(ig)
390            zd(ig,ilay)=zb(ig,ilay)*z1(ig)
391         ENDDO
392      ENDDO
393
394c    ** calcul de (d Planck / dT) a la temperature d'interface
395c       ------------------------------------------------------
396
397      z4st=4.*5.67e-8*ptimestep
398      DO ig=1,ngrid
399         zdplanck(ig)=z4st*pemis(ig)*ptsrf(ig)*ptsrf(ig)*ptsrf(ig)
400      ENDDO
401
402c    ** calcul de la temperature_d'interface et de sa tendance.
403c       on ecrit que la somme des flux est nulle a l'interface
404c       a t + \delta t,
405c       c'est a dire le flux radiatif a {t + \delta t}
406c       + le flux turbulent a {t + \delta t}
407c            qui s'ecrit K (T1-Tsurf) avec T1 = d1 Tsurf + c1
408c            (notation K dt = /cpp*b/)       
409c       + le flux dans le sol a t
410c       + l'evolution du flux dans le sol lorsque la temperature d'interface
411c       passe de sa valeur a t a sa valeur a {t + \delta t}.
412c       ----------------------------------------------------
413
414      DO ig=1,ngrid
415         z1(ig)=pcapcal(ig)*ptsrf(ig)+cpp*zb(ig,1)*zc(ig,1)
416     s     +zdplanck(ig)*ptsrf(ig)+ pfluxsrf(ig)*ptimestep
417         z2(ig)= pcapcal(ig)+cpp*zb(ig,1)*(1.-zd(ig,1))+zdplanck(ig)
418         ztsrf2(ig)=z1(ig)/z2(ig)
419         pdtsrf(ig)=(ztsrf2(ig)-ptsrf(ig))/ptimestep
420
421c        Modif speciale CO2 condensation:
422c        tconds = 1./(bcond-acond*log(.0095*pplev(ig,1)))
423c        if ((callcond).and.
424c    &      ((co2ice(ig).ne.0).or.(ztsrf2(ig).lt.tconds)))then
425c           zh(ig,1)=zc(ig,1)+zd(ig,1)*tconds
426c        else
427            zh(ig,1)=zc(ig,1)+zd(ig,1)*ztsrf2(ig)
428c        end if
429      ENDDO
430
431c    ** et a partir de la temperature au sol on remonte
432c       -----------------------------------------------
433
434      DO ilay=2,nlay
435         DO ig=1,ngrid
436            hh = max( zh(ig,ilay-1) , zhcond(ig,ilay-1) ) ! modif co2cond
437            zh(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*hh
438         ENDDO
439      ENDDO
440
441
442c-----------------------------------------------------------------------
443c   TRACERS
444c   -------
445
446      if(tracer) then
447
448c     Using the wind modified by friction for lifting and  sublimation
449c     ----------------------------------------------------------------
450        DO ig=1,ngrid
451          zu2=zu(ig,1)*zu(ig,1)+zv(ig,1)*zv(ig,1)
452          zcdv(ig)=zcdv_true(ig)*sqrt(zu2)
453          zcdh(ig)=zcdh_true(ig)*sqrt(zu2)
454        ENDDO
455
456c       Calcul du flux vertical au bas de la premiere couche (dust) :
457c       -----------------------------------------------------------
458        do ig=1,ngridmx 
459          rho(ig) = zb0(ig,1) /ptimestep
460c          zb(ig,1) = 0.
461        end do
462c       Dust lifting:
463        if (lifting) then
464           if (doubleq.AND.submicron) then
465             do ig=1,ngrid
466c              if(co2ice(ig).lt.1) then
467                 pdqsdif(ig,igcm_dust_mass) =
468     &             -alpha_lift(igcm_dust_mass) 
469                 pdqsdif(ig,igcm_dust_number) =
470     &             -alpha_lift(igcm_dust_number) 
471                 pdqsdif(ig,igcm_dust_submicron) =
472     &             -alpha_lift(igcm_dust_submicron)
473c              end if
474             end do
475           else if (doubleq) then
476             do ig=1,ngrid
477           !!! soulevement constant
478                 pdqsdif(ig,igcm_dust_mass) =
479     &             -alpha_lift(igcm_dust_mass) 
480                 pdqsdif(ig,igcm_dust_number) =
481     &             -alpha_lift(igcm_dust_number) 
482             end do
483           else if (submicron) then
484             do ig=1,ngrid
485                 pdqsdif(ig,igcm_dust_submicron) =
486     &             -alpha_lift(igcm_dust_submicron)
487             end do
488           else
489            call dustlift(ngrid,nlay,nq,rho,zcdh_true,zcdh,co2ice,
490     &                    pdqsdif)
491           endif !doubleq.AND.submicron
492        else
493           pdqsdif(1:ngrid,1:nq) = 0.
494        end if
495
496c       OU calcul de la valeur de q a la surface (water)  :
497c       ----------------------------------------
498        if (water) then
499            call watersat(ngridmx,ptsrf,pplev(1,1),qsat)
500        end if
501
502c      Inversion pour l'implicite sur q
503c       --------------------------------
504        do iq=1,nq
505          CALL multipl((nlay-1)*ngrid,zkh(1,2),zb0(1,2),zb(1,2))
506
507          if ((water).and.(iq.eq.igcm_h2o_vap)) then
508c            This line is required to account for turbulent transport
509c            from surface (e.g. ice) to mid-layer of atmosphere:
510             CALL multipl(ngrid,zcdv,zb0,zb(1,1))
511             CALL multipl(ngrid,dryness,zb(1,1),zb(1,1))
512          else ! (re)-initialize zb(:,1)
513             zb(1:ngrid,1)=0
514          end if
515
516          DO ig=1,ngrid
517               z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
518               zc(ig,nlay)=za(ig,nlay)*zq(ig,nlay,iq)*z1(ig)
519               zd(ig,nlay)=zb(ig,nlay)*z1(ig)
520          ENDDO
521 
522          DO ilay=nlay-1,2,-1
523               DO ig=1,ngrid
524                z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
525     $           zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
526                zc(ig,ilay)=(za(ig,ilay)*zq(ig,ilay,iq)+
527     $           zb(ig,ilay+1)*zc(ig,ilay+1))*z1(ig)
528                zd(ig,ilay)=zb(ig,ilay)*z1(ig)
529               ENDDO
530          ENDDO
531
532          if (water.and.(iq.eq.igcm_h2o_ice)) then
533            ! special case for water ice tracer: do not include
534            ! h2o ice tracer from surface (which is set when handling
535            ! h2o vapour case (see further down).
536            DO ig=1,ngrid
537                z1(ig)=1./(za(ig,1)+zb(ig,1)+
538     $           zb(ig,2)*(1.-zd(ig,2)))
539                zc(ig,1)=(za(ig,1)*zq(ig,1,iq)+
540     $         zb(ig,2)*zc(ig,2)) *z1(ig)
541            ENDDO
542          else ! general case
543            DO ig=1,ngrid
544                z1(ig)=1./(za(ig,1)+zb(ig,1)+
545     $           zb(ig,2)*(1.-zd(ig,2)))
546                zc(ig,1)=(za(ig,1)*zq(ig,1,iq)+
547     $         zb(ig,2)*zc(ig,2) +
548     $        (-pdqsdif(ig,iq)) *ptimestep) *z1(ig)  !tracer flux from surface
549            ENDDO
550          endif ! of if (water.and.(iq.eq.igcm_h2o_ice))
551 
552          IF ((water).and.(iq.eq.igcm_h2o_vap)) then
553c           Calculation for turbulent exchange with the surface (for ice)
554            DO ig=1,ngrid
555              zd(ig,1)=zb(ig,1)*z1(ig)
556              zq1temp(ig)=zc(ig,1)+ zd(ig,1)*qsat(ig)
557
558              pdqsdif(ig,igcm_h2o_ice)=rho(ig)*dryness(ig)*zcdv(ig)
559     &                       *(zq1temp(ig)-qsat(ig))
560c             write(*,*)'flux vers le sol=',pdqsdif(ig,nq)
561            END DO
562
563            DO ig=1,ngrid
564              if(.not.watercaptag(ig)) then
565                if ((-pdqsdif(ig,igcm_h2o_ice)*ptimestep)
566     &             .gt.pqsurf(ig,igcm_h2o_ice)) then
567c                 write(*,*)'on sublime plus que qsurf!'
568                  pdqsdif(ig,igcm_h2o_ice)=
569     &                         -pqsurf(ig,igcm_h2o_ice)/ptimestep
570c                 write(*,*)'flux vers le sol=',pdqsdif(ig,nq)
571                  z1(ig)=1./(za(ig,1)+ zb(ig,2)*(1.-zd(ig,2)))
572                  zc(ig,1)=(za(ig,1)*zq(ig,1,igcm_h2o_vap)+
573     $            zb(ig,2)*zc(ig,2) +
574     $            (-pdqsdif(ig,igcm_h2o_ice)) *ptimestep) *z1(ig)
575                  zq1temp(ig)=zc(ig,1)
576                endif   
577              endif ! if (.not.watercaptag(ig))
578c             Starting upward calculations for water :
579               zq(ig,1,igcm_h2o_vap)=zq1temp(ig)
580            ENDDO ! of DO ig=1,ngrid
581          ELSE
582c           Starting upward calculations for simple mixing of tracer (dust)
583            DO ig=1,ngrid
584               zq(ig,1,iq)=zc(ig,1)
585            ENDDO
586          END IF ! of IF ((water).and.(iq.eq.igcm_h2o_vap))
587
588          DO ilay=2,nlay
589             DO ig=1,ngrid
590                zq(ig,ilay,iq)=zc(ig,ilay)+zd(ig,ilay)*zq(ig,ilay-1,iq)
591             ENDDO
592          ENDDO
593        enddo ! of do iq=1,nq
594      end if ! of if(tracer)
595
596c-----------------------------------------------------------------------
597c   8. calcul final des tendances de la diffusion verticale
598c      ----------------------------------------------------
599
600      DO ilev = 1, nlay
601         DO ig=1,ngrid
602            pdudif(ig,ilev)=(    zu(ig,ilev)-
603     $      (pu(ig,ilev)+pdufi(ig,ilev)*ptimestep)    )/ptimestep
604            pdvdif(ig,ilev)=(    zv(ig,ilev)-
605     $      (pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep)    )/ptimestep
606            hh = max(ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep ,
607     $           zhcond(ig,ilev))        ! modif co2cond
608            pdhdif(ig,ilev)=( zh(ig,ilev)- hh )/ptimestep
609         ENDDO
610      ENDDO
611
612   
613      if (tracer) then
614        DO iq = 1, nq
615          DO ilev = 1, nlay
616            DO ig=1,ngrid
617              pdqdif(ig,ilev,iq)=(zq(ig,ilev,iq)-
618     $      (pq(ig,ilev,iq) + pdqfi(ig,ilev,iq)*ptimestep))/ptimestep
619            ENDDO
620          ENDDO
621        ENDDO
622      end if
623
624c    ** diagnostique final
625c       ------------------
626
627      IF(lecrit) THEN
628         PRINT*,'In vdif'
629         PRINT*,'Ts (t) and Ts (t+st)'
630         WRITE(*,'(a10,3a15)')
631     s   'theta(t)','theta(t+dt)','u(t)','u(t+dt)'
632         PRINT*,ptsrf(ngrid/2+1),ztsrf2(ngrid/2+1)
633         DO ilev=1,nlay
634            WRITE(*,'(4f15.7)')
635     s      ph(ngrid/2+1,ilev),zh(ngrid/2+1,ilev),
636     s      pu(ngrid/2+1,ilev),zu(ngrid/2+1,ilev)
637
638         ENDDO
639      ENDIF
640
641      if (calltherm .and. outptherm) then
642      if (ngrid .eq. 1) then
643        call WRITEDIAGFI(ngrid,'zh','zh inside vdifc',
644     &                       'SI',1,ph(:,:)+pdhfi(:,:)*ptimestep)
645        call WRITEDIAGFI(ngrid,'zkh','zkh',
646     &                       'SI',1,zkh)
647      endif
648      endif
649
650      RETURN
651      END
Note: See TracBrowser for help on using the repository browser.