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

Last change on this file since 263 was 263, checked in by aslmd, 14 years ago

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