source: trunk/mars/libf/phymars/vdifc.F @ 73

Last change on this file since 73 was 38, checked in by emillour, 15 years ago

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