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

Last change on this file since 1619 was 1455, checked in by aslmd, 9 years ago

MESOSCALE. can now use same callphys than GCM. necessary changes are in run.def. modified initracer to be compliant with lifting differences. and also made GCM 29 levels by default.
File size: 28.2 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,wstar,zcdv_true,zcdh_true,
8     $                hfmax,sensibFlux)
9
10      use tracer_mod, only: noms, igcm_dust_mass, igcm_dust_number,
11     &                      igcm_dust_submicron, igcm_h2o_vap,
12     &                      igcm_h2o_ice, alpha_lift
13      use surfdat_h, only: watercaptag, frost_albedo_threshold, dryness
14      USE comcstfi_h
15      use turb_mod, only: turb_resolved, ustar, tstar
16      IMPLICIT NONE
17
18c=======================================================================
19c
20c   subject:
21c   --------
22c   Turbulent diffusion (mixing) for potential T, U, V and tracer
23c
24c   Shema implicite
25c   On commence par rajouter au variables x la tendance physique
26c   et on resoult en fait:
27c      x(t+1) =  x(t) + dt * (dx/dt)phys(t)  +  dt * (dx/dt)difv(t+1)
28c
29c   author:
30c   ------
31c      Hourdin/Forget/Fournier
32c=======================================================================
33
34c-----------------------------------------------------------------------
35c   declarations:
36c   -------------
37
38#include "callkeys.h"
39#include "microphys.h"
40
41c
42c   arguments:
43c   ----------
44
45      INTEGER,INTENT(IN) :: ngrid,nlay
46      REAL,INTENT(IN) :: ptimestep
47      REAL,INTENT(IN) :: pplay(ngrid,nlay),pplev(ngrid,nlay+1)
48      REAL,INTENT(IN) :: pzlay(ngrid,nlay),pzlev(ngrid,nlay+1)
49      REAL,INTENT(IN) :: pu(ngrid,nlay),pv(ngrid,nlay)
50      REAL,INTENT(IN) :: ph(ngrid,nlay)
51      REAL,INTENT(IN) :: ptsrf(ngrid),pemis(ngrid)
52      REAL,INTENT(IN) :: pdufi(ngrid,nlay),pdvfi(ngrid,nlay)
53      REAL,INTENT(IN) :: pdhfi(ngrid,nlay)
54      REAL,INTENT(IN) :: pfluxsrf(ngrid)
55      REAL,INTENT(OUT) :: pdudif(ngrid,nlay),pdvdif(ngrid,nlay)
56      REAL,INTENT(OUT) :: pdtsrf(ngrid),pdhdif(ngrid,nlay)
57      REAL,INTENT(IN) :: pcapcal(ngrid)
58      REAL,INTENT(INOUT) :: pq2(ngrid,nlay+1)
59
60c    Argument added for condensation:
61      REAL,INTENT(IN) :: co2ice (ngrid), ppopsk(ngrid,nlay)
62      logical,INTENT(IN) :: lecrit
63
64      REAL,INTENT(IN) :: pz0(ngrid) ! surface roughness length (m)
65
66c    Argument added to account for subgrid gustiness :
67
68      REAL wstar(ngrid), hfmax(ngrid)!, zi(ngrid)
69
70c    Traceurs :
71      integer,intent(in) :: nq
72      REAL,INTENT(IN) :: pqsurf(ngrid,nq)
73      real,intent(in) :: pq(ngrid,nlay,nq), pdqfi(ngrid,nlay,nq)
74      real,intent(out) :: pdqdif(ngrid,nlay,nq)
75      real,intent(out) :: pdqsdif(ngrid,nq)
76     
77c   local:
78c   ------
79
80      REAL :: pt(ngrid,nlay)
81 
82      INTEGER ilev,ig,ilay,nlev
83
84      REAL z4st,zdplanck(ngrid)
85      REAL zkv(ngrid,nlay+1),zkh(ngrid,nlay+1)
86      REAL zkq(ngrid,nlay+1)
87      REAL zcdv(ngrid),zcdh(ngrid)
88      REAL zcdv_true(ngrid),zcdh_true(ngrid)    ! drag coeff are used by the LES to recompute u* and hfx
89      REAL zu(ngrid,nlay),zv(ngrid,nlay)
90      REAL zh(ngrid,nlay)
91      REAL ztsrf2(ngrid)
92      REAL z1(ngrid),z2(ngrid)
93      REAL za(ngrid,nlay),zb(ngrid,nlay)
94      REAL zb0(ngrid,nlay)
95      REAL zc(ngrid,nlay),zd(ngrid,nlay)
96      REAL zcst1
97      REAL zu2(ngrid)
98
99      EXTERNAL SSUM,SCOPY
100      REAL SSUM
101      LOGICAL,SAVE :: firstcall=.true.
102
103
104c     variable added for CO2 condensation:
105c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
106      REAL hh , zhcond(ngrid,nlay)
107      REAL,PARAMETER :: latcond=5.9e5
108      REAL,PARAMETER :: tcond1mb=136.27
109      REAL,SAVE :: acond,bcond
110     
111c     For latent heat release from ground ice sublimation   
112!      REAL tsrf_lw(ngrid)
113!      REAL alpha
114      REAL T1,T2
115      SAVE T1,T2
116      DATA T1,T2/-604.3,1080.7/ ! zeros of latent heat equation for ice
117
118c    Tracers :
119c    ~~~~~~~
120      INTEGER iq
121      REAL zq(ngrid,nlay,nq)
122      REAL zq1temp(ngrid)
123      REAL rho(ngrid) ! near surface air density
124      REAL qsat(ngrid)
125
126      REAL kmixmin
127
128c    Mass-variation scheme :
129c    ~~~~~~~
130
131      INTEGER j,l
132      REAL zcondicea(ngrid,nlay)
133      REAL zt(ngrid,nlay),ztcond(ngrid,nlay+1)
134      REAL betam(ngrid,nlay),dmice(ngrid,nlay)
135      REAL pdtc(ngrid,nlay)
136      REAL zhs(ngrid,nlay)
137      REAL,SAVE :: ccond
138
139c     Theta_m formulation for mass-variation scheme :
140c     ~~~~~~~
141
142      INTEGER,SAVE :: ico2
143      INTEGER llnt(ngrid)
144      REAL,SAVE :: m_co2, m_noco2, A , B
145      REAL vmr_co2(ngrid,nlay)
146      REAL qco2,mmean
147
148      REAL,INTENT(OUT) :: sensibFlux(ngrid)
149
150c    ** un petit test de coherence
151c       --------------------------
152
153      IF (firstcall) THEN
154c        To compute: Tcond= 1./(bcond-acond*log(.0095*p)) (p in pascal)
155         bcond=1./tcond1mb
156         acond=r/latcond
157         ccond=cpp/(g*latcond)
158         PRINT*,'In vdifc: Tcond(P=1mb)=',tcond1mb,' Lcond=',latcond
159         PRINT*,'          acond,bcond,ccond',acond,bcond,ccond
160
161
162         ico2=0
163
164         if (tracer) then
165c          Prepare Special treatment if one of the tracer is CO2 gas
166           do iq=1,nq
167             if (noms(iq).eq."co2") then
168                ico2=iq
169                m_co2 = 44.01E-3  ! CO2 molecular mass (kg/mol)
170                m_noco2 = 33.37E-3  ! Non condensible mol mass (kg/mol)
171c               Compute A and B coefficient use to compute
172c               mean molecular mass Mair defined by
173c               1/Mair = q(ico2)/m_co2 + (1-q(ico2))/m_noco2
174c               1/Mair = A*q(ico2) + B
175                A =(1/m_co2 - 1/m_noco2)
176                B=1/m_noco2
177             endif
178           enddo
179         end if
180
181        firstcall=.false.
182      ENDIF
183
184
185c-----------------------------------------------------------------------
186c    1. initialisation
187c    -----------------
188
189      nlev=nlay+1
190
191      ! initialize output tendencies to zero:
192      pdudif(1:ngrid,1:nlay)=0
193      pdvdif(1:ngrid,1:nlay)=0
194      pdhdif(1:ngrid,1:nlay)=0
195      pdtsrf(1:ngrid)=0
196      pdqdif(1:ngrid,1:nlay,1:nq)=0
197      pdqsdif(1:ngrid,1:nq)=0
198
199c    ** calcul de rho*dz et dt*rho/dz=dt*rho**2 g/dp
200c       avec rho=p/RT=p/ (R Theta) (p/ps)**kappa
201c       ----------------------------------------
202
203      DO ilay=1,nlay
204         DO ig=1,ngrid
205            za(ig,ilay)=(pplev(ig,ilay)-pplev(ig,ilay+1))/g
206! Mass variation scheme:
207            betam(ig,ilay)=-za(ig,ilay)*latcond/(cpp*ppopsk(ig,ilay))
208         ENDDO
209      ENDDO
210
211      zcst1=4.*g*ptimestep/(r*r)
212      DO ilev=2,nlev-1
213         DO ig=1,ngrid
214            zb0(ig,ilev)=pplev(ig,ilev)*
215     s      (pplev(ig,1)/pplev(ig,ilev))**rcp /
216     s      (ph(ig,ilev-1)+ph(ig,ilev))
217            zb0(ig,ilev)=zcst1*zb0(ig,ilev)*zb0(ig,ilev)/
218     s      (pplay(ig,ilev-1)-pplay(ig,ilev))
219         ENDDO
220      ENDDO
221      DO ig=1,ngrid
222         zb0(ig,1)=ptimestep*pplev(ig,1)/(r*ptsrf(ig))
223      ENDDO
224
225c    ** diagnostique pour l'initialisation
226c       ----------------------------------
227
228      IF(lecrit) THEN
229         ig=ngrid/2+1
230         PRINT*,'Pression (mbar) ,altitude (km),u,v,theta, rho dz'
231         DO ilay=1,nlay
232            WRITE(*,'(6f11.5)')
233     s      .01*pplay(ig,ilay),.001*pzlay(ig,ilay),
234     s      pu(ig,ilay),pv(ig,ilay),ph(ig,ilay),za(ig,ilay)
235         ENDDO
236         PRINT*,'Pression (mbar) ,altitude (km),zb'
237         DO ilev=1,nlay
238            WRITE(*,'(3f15.7)')
239     s      .01*pplev(ig,ilev),.001*pzlev(ig,ilev),
240     s      zb0(ig,ilev)
241         ENDDO
242      ENDIF
243
244c     -----------------------------------
245c     Potential Condensation temperature:
246c     -----------------------------------
247
248c     Compute CO2 Volume mixing ratio
249c     -------------------------------
250      if (callcond.and.(ico2.ne.0)) then
251         DO ilev=1,nlay
252            DO ig=1,ngrid
253              qco2=MAX(1.E-30
254     &             ,pq(ig,ilev,ico2)+pdqfi(ig,ilev,ico2)*ptimestep)
255c             Mean air molecular mass = 1/(q(ico2)/m_co2 + (1-q(ico2))/m_noco2)
256              mmean=1/(A*qco2 +B)
257              vmr_co2(ig,ilev) = qco2*mmean/m_co2
258            ENDDO
259         ENDDO
260      else
261         DO ilev=1,nlay
262            DO ig=1,ngrid
263              vmr_co2(ig,ilev)=0.95
264            ENDDO
265         ENDDO
266      end if
267
268c     forecast of atmospheric temperature zt and frost temperature ztcond
269c     --------------------------------------------------------------------
270
271      if (callcond) then
272        DO ilev=1,nlay
273          DO ig=1,ngrid
274              ztcond(ig,ilev)=
275     &      1./(bcond-acond*log(.01*vmr_co2(ig,ilev)*pplay(ig,ilev)))
276            if (pplay(ig,ilev).lt.1e-4) ztcond(ig,ilev)=0.0 !mars Monica
277!            zhcond(ig,ilev) =
278!     &  (1./(bcond-acond*log(.0095*pplay(ig,ilev))))/ppopsk(ig,ilev)
279              zhcond(ig,ilev) = ztcond(ig,ilev)/ppopsk(ig,ilev)
280          END DO
281        END DO
282        ztcond(:,nlay+1)=ztcond(:,nlay)
283      else
284         zhcond(:,:) = 0
285         ztcond(:,:) = 0
286      end if
287
288
289c-----------------------------------------------------------------------
290c   2. ajout des tendances physiques
291c      -----------------------------
292
293      DO ilev=1,nlay
294         DO ig=1,ngrid
295            zu(ig,ilev)=pu(ig,ilev)+pdufi(ig,ilev)*ptimestep
296            zv(ig,ilev)=pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep
297            zh(ig,ilev)=ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep
298!            zh(ig,ilev)=max(zh(ig,ilev),zhcond(ig,ilev))
299         ENDDO
300      ENDDO
301      if(tracer) then
302        DO iq =1, nq
303         DO ilev=1,nlay
304           DO ig=1,ngrid
305              zq(ig,ilev,iq)=pq(ig,ilev,iq)+pdqfi(ig,ilev,iq)*ptimestep
306           ENDDO
307         ENDDO
308        ENDDO
309      end if
310
311c-----------------------------------------------------------------------
312c   3. schema de turbulence
313c      --------------------
314
315c    ** source d'energie cinetique turbulente a la surface
316c       (condition aux limites du schema de diffusion turbulente
317c       dans la couche limite
318c       ---------------------
319
320      CALL vdif_cd(ngrid,nlay,pz0,g,pzlay,pu,pv,wstar,ptsrf,ph
321     &             ,zcdv_true,zcdh_true)
322
323        zu2(:)=pu(:,1)*pu(:,1)+pv(:,1)*pv(:,1)
324
325        IF (callrichsl) THEN
326          zcdv(:)=zcdv_true(:)*sqrt(zu2(:)+
327     &     (log(1.+0.7*wstar(:) + 2.3*wstar(:)**2))**2)
328          zcdh(:)=zcdh_true(:)*sqrt(zu2(:)+
329     &     (log(1.+0.7*wstar(:) + 2.3*wstar(:)**2))**2)
330
331           ustar(:)=sqrt(zcdv_true(:))*sqrt(zu2(:)+
332     &     (log(1.+0.7*wstar(:) + 2.3*wstar(:)**2))**2)
333
334           tstar(:)=0.
335           DO ig=1,ngrid
336              IF (zcdh_true(ig) .ne. 0.) THEN      ! When Cd=Ch=0, u*=t*=0
337                 tstar(ig)=(ph(ig,1)-ptsrf(ig))*zcdh(ig)/ustar(ig)
338              ENDIF
339           ENDDO
340
341        ELSE
342           zcdv(:)=zcdv_true(:)*sqrt(zu2(:))     ! 1 / bulk aerodynamic momentum conductance
343           zcdh(:)=zcdh_true(:)*sqrt(zu2(:))     ! 1 / bulk aerodynamic heat conductance
344           ustar(:)=sqrt(zcdv_true(:))*sqrt(zu2(:))
345           tstar(:)=(ph(:,1)-ptsrf(:))*zcdh_true(:)/sqrt(zcdv_true(:))
346        ENDIF
347
348! Some usefull diagnostics for the new surface layer parametrization :
349
350!        call WRITEDIAGFI(ngrid,'vdifc_zcdv_true',
351!     &              'momentum drag','no units',
352!     &                         2,zcdv_true)
353!        call WRITEDIAGFI(ngrid,'vdifc_zcdh_true',
354!     &              'heat drag','no units',
355!     &                         2,zcdh_true)
356!        call WRITEDIAGFI(ngrid,'vdifc_ust',
357!     &              'friction velocity','m/s',2,ust)
358!       call WRITEDIAGFI(ngrid,'vdifc_tst',
359!     &              'friction temperature','K',2,tst)
360!        call WRITEDIAGFI(ngrid,'vdifc_zcdv',
361!     &              'aerodyn momentum conductance','m/s',
362!     &                         2,zcdv)
363!        call WRITEDIAGFI(ngrid,'vdifc_zcdh',
364!     &              'aerodyn heat conductance','m/s',
365!     &                         2,zcdh)
366
367c    ** schema de diffusion turbulente dans la couche limite
368c       ----------------------------------------------------
369       IF (.not. callyamada4) THEN
370
371       CALL vdif_kc(ngrid,nlay,nq,ptimestep,g,pzlev,pzlay
372     &              ,pu,pv,ph,zcdv_true
373     &              ,pq2,zkv,zkh,zq)
374
375      ELSE
376
377      pt(:,:)=ph(:,:)*ppopsk(:,:)
378      CALL yamada4(ngrid,nlay,nq,ptimestep,g,r,pplev,pt
379     s   ,pzlev,pzlay,pu,pv,ph,pq,zcdv_true,pq2,zkv,zkh,zkq,ustar
380     s   ,9)
381      ENDIF
382
383      if ((doubleq).and.(ngrid.eq.1)) then
384        kmixmin = 80. !80.! minimum eddy mix coeff in 1D
385        do ilev=1,nlay
386          do ig=1,ngrid
387           zkh(ig,ilev) = max(kmixmin,zkh(ig,ilev))
388           zkv(ig,ilev) = max(kmixmin,zkv(ig,ilev))
389          end do
390        end do
391      end if
392
393c    ** diagnostique pour le schema de turbulence
394c       -----------------------------------------
395
396      IF(lecrit) THEN
397         PRINT*
398         PRINT*,'Diagnostic for the vertical turbulent mixing'
399         PRINT*,'Cd for momentum and potential temperature'
400
401         PRINT*,zcdv(ngrid/2+1),zcdh(ngrid/2+1)
402         PRINT*,'Mixing coefficient for momentum and pot.temp.'
403         DO ilev=1,nlay
404            PRINT*,zkv(ngrid/2+1,ilev),zkh(ngrid/2+1,ilev)
405         ENDDO
406      ENDIF
407
408
409
410
411c-----------------------------------------------------------------------
412c   4. inversion pour l'implicite sur u
413c      --------------------------------
414
415c    ** l'equation est
416c       u(t+1) =  u(t) + dt * {(du/dt)phys}(t)  +  dt * {(du/dt)difv}(t+1)
417c       avec
418c       /zu/ = u(t) + dt * {(du/dt)phys}(t)   (voir paragraphe 2.)
419c       et
420c       dt * {(du/dt)difv}(t+1) = dt * {(d/dz)[ Ku (du/dz) ]}(t+1)
421c       donc les entrees sont /zcdv/ pour la condition a la limite sol
422c       et /zkv/ = Ku
423 
424      CALL multipl((nlay-1)*ngrid,zkv(1,2),zb0(1,2),zb(1,2))
425      CALL multipl(ngrid,zcdv,zb0,zb)
426
427      DO ig=1,ngrid
428         z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
429         zc(ig,nlay)=za(ig,nlay)*zu(ig,nlay)*z1(ig)
430         zd(ig,nlay)=zb(ig,nlay)*z1(ig)
431      ENDDO
432
433      DO ilay=nlay-1,1,-1
434         DO ig=1,ngrid
435            z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
436     $         zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
437            zc(ig,ilay)=(za(ig,ilay)*zu(ig,ilay)+
438     $         zb(ig,ilay+1)*zc(ig,ilay+1))*z1(ig)
439            zd(ig,ilay)=zb(ig,ilay)*z1(ig)
440         ENDDO
441      ENDDO
442
443      DO ig=1,ngrid
444         zu(ig,1)=zc(ig,1)
445      ENDDO
446      DO ilay=2,nlay
447         DO ig=1,ngrid
448            zu(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zu(ig,ilay-1)
449         ENDDO
450      ENDDO
451
452
453
454
455
456c-----------------------------------------------------------------------
457c   5. inversion pour l'implicite sur v
458c      --------------------------------
459
460c    ** l'equation est
461c       v(t+1) =  v(t) + dt * {(dv/dt)phys}(t)  +  dt * {(dv/dt)difv}(t+1)
462c       avec
463c       /zv/ = v(t) + dt * {(dv/dt)phys}(t)   (voir paragraphe 2.)
464c       et
465c       dt * {(dv/dt)difv}(t+1) = dt * {(d/dz)[ Kv (dv/dz) ]}(t+1)
466c       donc les entrees sont /zcdv/ pour la condition a la limite sol
467c       et /zkv/ = Kv
468
469      DO ig=1,ngrid
470         z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
471         zc(ig,nlay)=za(ig,nlay)*zv(ig,nlay)*z1(ig)
472         zd(ig,nlay)=zb(ig,nlay)*z1(ig)
473      ENDDO
474
475      DO ilay=nlay-1,1,-1
476         DO ig=1,ngrid
477            z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
478     $         zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
479            zc(ig,ilay)=(za(ig,ilay)*zv(ig,ilay)+
480     $         zb(ig,ilay+1)*zc(ig,ilay+1))*z1(ig)
481            zd(ig,ilay)=zb(ig,ilay)*z1(ig)
482         ENDDO
483      ENDDO
484
485      DO ig=1,ngrid
486         zv(ig,1)=zc(ig,1)
487      ENDDO
488      DO ilay=2,nlay
489         DO ig=1,ngrid
490            zv(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zv(ig,ilay-1)
491         ENDDO
492      ENDDO
493
494
495
496
497
498c-----------------------------------------------------------------------
499c   6. inversion pour l'implicite sur h sans oublier le couplage
500c      avec le sol (conduction)
501c      ------------------------
502
503c    ** l'equation est
504c       h(t+1) =  h(t) + dt * {(dh/dt)phys}(t)  +  dt * {(dh/dt)difv}(t+1)
505c       avec
506c       /zh/ = h(t) + dt * {(dh/dt)phys}(t)   (voir paragraphe 2.)
507c       et
508c       dt * {(dh/dt)difv}(t+1) = dt * {(d/dz)[ Kh (dh/dz) ]}(t+1)
509c       donc les entrees sont /zcdh/ pour la condition de raccord au sol
510c       et /zkh/ = Kh
511c       -------------
512
513c Mass variation scheme:
514      CALL multipl((nlay-1)*ngrid,zkh(1,2),zb0(1,2),zb(1,2))
515      CALL multipl(ngrid,zcdh,zb0,zb)
516
517c on initialise dm c
518     
519      pdtc(:,:)=0.
520      zt(:,:)=0.
521      dmice(:,:)=0.
522
523c    ** calcul de (d Planck / dT) a la temperature d'interface
524c       ------------------------------------------------------
525
526      z4st=4.*5.67e-8*ptimestep
527      IF (tke_heat_flux .eq. 0.) THEN
528      DO ig=1,ngrid
529         zdplanck(ig)=z4st*pemis(ig)*ptsrf(ig)*ptsrf(ig)*ptsrf(ig)
530      ENDDO
531      ELSE
532         zdplanck(:)=0.
533      ENDIF
534
535! calcul de zc et zd pour la couche top en prenant en compte le terme
536! de variation de masse (on fait une boucle pour que ça converge)
537
538! Identification des points de grilles qui ont besoin de la correction
539
540      llnt(:)=1
541      IF (.not.turb_resolved) THEN
542      IF (callcond) THEN
543       DO ig=1,ngrid
544         DO l=1,nlay
545            if(zh(ig,l) .lt. zhcond(ig,l)) then
546               llnt(ig)=300 
547! 200 and 100 do not go beyond month 9 with normal dissipation
548               goto 5
549            endif
550         ENDDO
5515      continue
552       ENDDO
553      ENDIF
554
555      ENDIF
556
557      DO ig=1,ngrid
558
559! Initialization of z1 and zd, which do not depend on dmice
560
561      z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
562      zd(ig,nlay)=zb(ig,nlay)*z1(ig)
563
564      DO ilay=nlay-1,1,-1
565          z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
566     $        zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
567          zd(ig,ilay)=zb(ig,ilay)*z1(ig)
568      ENDDO
569
570! Convergence loop
571
572      DO j=1,llnt(ig)
573
574            z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
575            zc(ig,nlay)=za(ig,nlay)*zh(ig,nlay)
576     &      -betam(ig,nlay)*dmice(ig,nlay)
577            zc(ig,nlay)=zc(ig,nlay)*z1(ig)
578!            zd(ig,nlay)=zb(ig,nlay)*z1(ig)
579
580! calcul de zc et zd pour les couches du haut vers le bas
581
582             DO ilay=nlay-1,1,-1
583               z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
584     $            zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
585               zc(ig,ilay)=(za(ig,ilay)*zh(ig,ilay)+
586     $            zb(ig,ilay+1)*zc(ig,ilay+1)-
587     $            betam(ig,ilay)*dmice(ig,ilay))*z1(ig)
588!               zd(ig,ilay)=zb(ig,ilay)*z1(ig)
589            ENDDO
590
591c    ** calcul de la temperature_d'interface et de sa tendance.
592c       on ecrit que la somme des flux est nulle a l'interface
593c       a t + \delta t,
594c       c'est a dire le flux radiatif a {t + \delta t}
595c       + le flux turbulent a {t + \delta t}
596c            qui s'ecrit K (T1-Tsurf) avec T1 = d1 Tsurf + c1
597c            (notation K dt = /cpp*b/)       
598c       + le flux dans le sol a t
599c       + l'evolution du flux dans le sol lorsque la temperature d'interface
600c       passe de sa valeur a t a sa valeur a {t + \delta t}.
601c       ----------------------------------------------------
602
603         z1(ig)=pcapcal(ig)*ptsrf(ig)+cpp*zb(ig,1)*zc(ig,1)
604     s     +zdplanck(ig)*ptsrf(ig)+ pfluxsrf(ig)*ptimestep
605         z2(ig)= pcapcal(ig)+cpp*zb(ig,1)*(1.-zd(ig,1))+zdplanck(ig)
606         ztsrf2(ig)=z1(ig)/z2(ig)
607!         pdtsrf(ig)=(ztsrf2(ig)-ptsrf(ig))/ptimestep  !incremented outside loop
608            zhs(ig,1)=zc(ig,1)+zd(ig,1)*ztsrf2(ig)
609
610c    ** et a partir de la temperature au sol on remonte
611c       -----------------------------------------------
612
613            DO ilay=2,nlay
614               zhs(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zhs(ig,ilay-1)
615            ENDDO
616            DO ilay=1,nlay
617               zt(ig,ilay)=zhs(ig,ilay)*ppopsk(ig,ilay)
618            ENDDO
619
620c      Condensation/sublimation in the atmosphere
621c      ------------------------------------------
622c      (computation of zcondicea and dmice)
623
624      zcondicea(ig,:)=0.
625      pdtc(ig,:)=0.
626
627      DO l=nlay , 1, -1
628           IF(zt(ig,l).LT.ztcond(ig,l)) THEN
629               pdtc(ig,l)=(ztcond(ig,l) - zt(ig,l))/ptimestep
630               zcondicea(ig,l)=(pplev(ig,l)-pplev(ig,l+1))
631     &                        *ccond*pdtc(ig,l)
632              dmice(ig,l)= dmice(ig,l) + zcondicea(ig,l)*ptimestep
633            END IF
634      ENDDO
635
636         ENDDO    !of Do j=1,XXX
637
638      ENDDO   !of Do ig=1,nlay
639
640      pdtsrf(:)=(ztsrf2(:)-ptsrf(:))/ptimestep
641     
642      DO ig=1,ngrid  ! computing sensible heat flux (atm => surface)
643         sensibFlux(ig)=cpp*zb(ig,1)/ptimestep*(zhs(ig,1)-ztsrf2(ig))
644      ENDDO
645
646c-----------------------------------------------------------------------
647c   TRACERS
648c   -------
649
650      if(tracer) then
651
652c     Using the wind modified by friction for lifting and  sublimation
653c     ----------------------------------------------------------------
654
655!     This is computed above and takes into account surface-atmosphere flux
656!     enhancement by subgrid gustiness and atmospheric-stability related
657!     variations of transfer coefficients.
658
659!        DO ig=1,ngrid
660!          zu2(ig)=zu(ig,1)*zu(ig,1)+zv(ig,1)*zv(ig,1)
661!          zcdv(ig)=zcdv_true(ig)*sqrt(zu2(ig))
662!          zcdh(ig)=zcdh_true(ig)*sqrt(zu2(ig))
663!        ENDDO
664
665c       Calcul du flux vertical au bas de la premiere couche (dust) :
666c       -----------------------------------------------------------
667        do ig=1,ngrid 
668          rho(ig) = zb0(ig,1) /ptimestep
669c          zb(ig,1) = 0.
670        end do
671c       Dust lifting:
672        if (lifting) then
673#ifndef MESOSCALE
674           if (doubleq.AND.submicron) then
675             do ig=1,ngrid
676c              if(co2ice(ig).lt.1) then
677                 pdqsdif(ig,igcm_dust_mass) =
678     &             -alpha_lift(igcm_dust_mass) 
679                 pdqsdif(ig,igcm_dust_number) =
680     &             -alpha_lift(igcm_dust_number) 
681                 pdqsdif(ig,igcm_dust_submicron) =
682     &             -alpha_lift(igcm_dust_submicron)
683c              end if
684             end do
685           else if (doubleq) then
686             do ig=1,ngrid
687               if(co2ice(ig).lt.1) then ! pas de soulevement si glace CO2
688                 pdqsdif(ig,igcm_dust_mass) =
689     &             -alpha_lift(igcm_dust_mass) 
690                 pdqsdif(ig,igcm_dust_number) =
691     &             -alpha_lift(igcm_dust_number)
692               end if
693             end do
694           else if (submicron) then
695             do ig=1,ngrid
696                 pdqsdif(ig,igcm_dust_submicron) =
697     &             -alpha_lift(igcm_dust_submicron)
698             end do
699           else
700#endif
701            call dustlift(ngrid,nlay,nq,rho,zcdh_true,zcdh,co2ice,
702     &                    pdqsdif)
703#ifndef MESOSCALE
704           endif !doubleq.AND.submicron
705#endif
706        else
707           pdqsdif(1:ngrid,1:nq) = 0.
708        end if
709
710c       OU calcul de la valeur de q a la surface (water)  :
711c       ----------------------------------------
712        if (water) then
713            call watersat(ngrid,ptsrf,pplev(1,1),qsat)
714        end if
715
716c      Inversion pour l'implicite sur q
717c       --------------------------------
718        do iq=1,nq
719          CALL multipl((nlay-1)*ngrid,zkh(1,2),zb0(1,2),zb(1,2))
720
721          if ((water).and.(iq.eq.igcm_h2o_vap)) then
722c            This line is required to account for turbulent transport
723c            from surface (e.g. ice) to mid-layer of atmosphere:
724             CALL multipl(ngrid,zcdv,zb0,zb(1,1))
725             CALL multipl(ngrid,dryness,zb(1,1),zb(1,1))
726          else ! (re)-initialize zb(:,1)
727             zb(1:ngrid,1)=0
728          end if
729
730          DO ig=1,ngrid
731               z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
732               zc(ig,nlay)=za(ig,nlay)*zq(ig,nlay,iq)*z1(ig)
733               zd(ig,nlay)=zb(ig,nlay)*z1(ig)
734          ENDDO
735 
736          DO ilay=nlay-1,2,-1
737               DO ig=1,ngrid
738                z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
739     $           zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
740                zc(ig,ilay)=(za(ig,ilay)*zq(ig,ilay,iq)+
741     $           zb(ig,ilay+1)*zc(ig,ilay+1))*z1(ig)
742                zd(ig,ilay)=zb(ig,ilay)*z1(ig)
743               ENDDO
744          ENDDO
745
746          if (water.and.(iq.eq.igcm_h2o_ice)) then
747            ! special case for water ice tracer: do not include
748            ! h2o ice tracer from surface (which is set when handling
749           ! h2o vapour case (see further down).
750            DO ig=1,ngrid
751                z1(ig)=1./(za(ig,1)+zb(ig,1)+
752     $           zb(ig,2)*(1.-zd(ig,2)))
753                zc(ig,1)=(za(ig,1)*zq(ig,1,iq)+
754     $         zb(ig,2)*zc(ig,2)) *z1(ig)
755            ENDDO
756          else ! general case
757            DO ig=1,ngrid
758                z1(ig)=1./(za(ig,1)+zb(ig,1)+
759     $           zb(ig,2)*(1.-zd(ig,2)))
760                zc(ig,1)=(za(ig,1)*zq(ig,1,iq)+
761     $         zb(ig,2)*zc(ig,2) +
762     $        (-pdqsdif(ig,iq)) *ptimestep) *z1(ig)  !tracer flux from surface
763            ENDDO
764          endif ! of if (water.and.(iq.eq.igcm_h2o_ice))
765 
766          IF ((water).and.(iq.eq.igcm_h2o_vap)) then
767c           Calculation for turbulent exchange with the surface (for ice)
768            DO ig=1,ngrid
769              zd(ig,1)=zb(ig,1)*z1(ig)
770              zq1temp(ig)=zc(ig,1)+ zd(ig,1)*qsat(ig)
771
772              pdqsdif(ig,igcm_h2o_ice)=rho(ig)*dryness(ig)*zcdv(ig)
773     &                       *(zq1temp(ig)-qsat(ig))
774c             write(*,*)'flux vers le sol=',pdqsdif(ig,nq)
775            END DO
776
777            DO ig=1,ngrid
778              if(.not.watercaptag(ig)) then
779                if ((-pdqsdif(ig,igcm_h2o_ice)*ptimestep)
780     &             .gt.pqsurf(ig,igcm_h2o_ice)) then
781c                 write(*,*)'on sublime plus que qsurf!'
782                  pdqsdif(ig,igcm_h2o_ice)=
783     &                         -pqsurf(ig,igcm_h2o_ice)/ptimestep
784c                 write(*,*)'flux vers le sol=',pdqsdif(ig,nq)
785                  z1(ig)=1./(za(ig,1)+ zb(ig,2)*(1.-zd(ig,2)))
786                  zc(ig,1)=(za(ig,1)*zq(ig,1,igcm_h2o_vap)+
787     $            zb(ig,2)*zc(ig,2) +
788     $            (-pdqsdif(ig,igcm_h2o_ice)) *ptimestep) *z1(ig)
789                  zq1temp(ig)=zc(ig,1)
790                endif   
791              endif ! if (.not.watercaptag(ig))
792c             Starting upward calculations for water :
793              zq(ig,1,igcm_h2o_vap)=zq1temp(ig)
794             
795!c             Take into account H2O latent heat in surface energy budget
796!c             We solve dT/dt = (2834.3-0.28*(T-To)-0.004*(T-To)^2)*1e3*iceflux/cpp
797!              tsrf_lw(ig) = ptsrf(ig) + pdtsrf(ig) *ptimestep
798!             
799!              alpha = exp(-4*abs(T1-T2)*pdqsdif(ig,igcm_h2o_ice)
800!     &            *ptimestep/pcapcal(ig))
801!
802!              tsrf_lw(ig) = (tsrf_lw(ig)*(T2-alpha*T1)+T1*T2*(alpha-1))
803!     &                      /(tsrf_lw(ig)*(1-alpha)+alpha*T2-T1)  ! surface temperature at t+1
804!
805!              pdtsrf(ig) = (tsrf_lw(ig)-ptsrf(ig))/ptimestep
806
807               if(pqsurf(ig,igcm_h2o_ice)
808     &           +pdqsdif(ig,igcm_h2o_ice)*ptimestep
809     &           .gt.frost_albedo_threshold) ! if there is still ice, T cannot exceed To
810     &           pdtsrf(ig) = min(pdtsrf(ig),(To-ptsrf(ig))/ptimestep) ! ice melt case
811     
812            ENDDO ! of DO ig=1,ngrid
813          ELSE
814c           Starting upward calculations for simple mixing of tracer (dust)
815            DO ig=1,ngrid
816               zq(ig,1,iq)=zc(ig,1)
817            ENDDO
818          END IF ! of IF ((water).and.(iq.eq.igcm_h2o_vap))
819
820          DO ilay=2,nlay
821             DO ig=1,ngrid
822               zq(ig,ilay,iq)=zc(ig,ilay)+zd(ig,ilay)*zq(ig,ilay-1,iq)
823             ENDDO
824          ENDDO
825        enddo ! of do iq=1,nq
826      end if ! of if(tracer)
827     
828
829c-----------------------------------------------------------------------
830c   8. calcul final des tendances de la diffusion verticale
831c      ----------------------------------------------------
832
833      DO ilev = 1, nlay
834         DO ig=1,ngrid
835            pdudif(ig,ilev)=(    zu(ig,ilev)-
836     $      (pu(ig,ilev)+pdufi(ig,ilev)*ptimestep)    )/ptimestep
837            pdvdif(ig,ilev)=(    zv(ig,ilev)-
838     $      (pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep)    )/ptimestep
839            hh = ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep
840     $  + (latcond*dmice(ig,ilev)/cpp)/ppopsk(ig,ilev)
841            pdhdif(ig,ilev)=( zhs(ig,ilev)- hh )/ptimestep
842         ENDDO
843      ENDDO
844
845      if (tracer) then
846        DO iq = 1, nq
847          DO ilev = 1, nlay
848            DO ig=1,ngrid
849              pdqdif(ig,ilev,iq)=(zq(ig,ilev,iq)-
850     $      (pq(ig,ilev,iq) + pdqfi(ig,ilev,iq)*ptimestep))/ptimestep
851            ENDDO
852          ENDDO
853        ENDDO
854      end if
855
856c    ** diagnostique final
857c       ------------------
858
859      IF(lecrit) THEN
860         PRINT*,'In vdif'
861         PRINT*,'Ts (t) and Ts (t+st)'
862         WRITE(*,'(a10,3a15)')
863     s   'theta(t)','theta(t+dt)','u(t)','u(t+dt)'
864         PRINT*,ptsrf(ngrid/2+1),ztsrf2(ngrid/2+1)
865         DO ilev=1,nlay
866            WRITE(*,'(4f15.7)')
867     s      ph(ngrid/2+1,ilev),zhs(ngrid/2+1,ilev),
868     s      pu(ngrid/2+1,ilev),zu(ngrid/2+1,ilev)
869
870         ENDDO
871      ENDIF
872
873      RETURN
874      END SUBROUTINE vdifc
Note: See TracBrowser for help on using the repository browser.