source: trunk/mesoscale/LMDZ.MARS/libf_gcm/phymars/vdifc.F @ 92

Last change on this file since 92 was 57, checked in by aslmd, 15 years ago

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