source: LMDZ4/trunk/libf/phylmd/phystokenc.F @ 1291

Last change on this file since 1291 was 1146, checked in by Laurent Fairhead, 16 years ago

Réintegration dans le tronc des modifications issues de la branche LMDZ-dev
comprises entre la révision 1074 et 1145
Validation: une simulation de 1 jour en séquentiel sur PC donne les mêmes
résultats entre la trunk et la dev
LF

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 13.3 KB
RevLine 
[524]1!
2c
3c
4      SUBROUTINE phystokenc (
5     I                   nlon,nlev,pdtphys,rlon,rlat,
6     I                   pt,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
[541]7     I                   pfm_therm,pentr_therm,
[1067]8     I                   cdragh, pcoefh,yu1,yv1,ftsol,pctsrf,
[541]9     I                   frac_impa,frac_nucl,
[524]10     I                   pphis,paire,dtime,itap)
11      USE ioipsl
[766]12      USE dimphy
[1146]13      USE infotrac, ONLY : nqtot
[766]14      USE iophy
[524]15      IMPLICIT none
16
17c======================================================================
18c Auteur(s) FH
19c Objet: Moniteur general des tendances traceurs
20c
21
22c======================================================================
23#include "dimensions.h"
24#include "tracstoke.h"
25#include "indicesol.h"
26#include "control.h"
27c======================================================================
28
29c Arguments:
30c
31c   EN ENTREE:
32c   ==========
33c
34c   divers:
35c   -------
36c
37      integer nlon ! nombre de points horizontaux
38      integer nlev ! nombre de couches verticales
39      real pdtphys ! pas d'integration pour la physique (seconde)
40c
[541]41      integer physid, itap
42      save physid
[766]43c$OMP THREADPRIVATE(physid)
[541]44      integer ndex2d(iim*(jjm+1)),ndex3d(iim*(jjm+1)*klev)
[524]45
46c   convection:
47c   -----------
48c
49      REAL pmfu(klon,klev)  ! flux de masse dans le panache montant
50      REAL pmfd(klon,klev)  ! flux de masse dans le panache descendant
51      REAL pen_u(klon,klev) ! flux entraine dans le panache montant
52      REAL pde_u(klon,klev) ! flux detraine dans le panache montant
53      REAL pen_d(klon,klev) ! flux entraine dans le panache descendant
54      REAL pde_d(klon,klev) ! flux detraine dans le panache descendant
[766]55      real pt(klon,klev)
56      REAL,allocatable,save :: t(:,:)
57c$OMP THREADPRIVATE(t)
[524]58c
59      REAL rlon(klon), rlat(klon), dtime
60      REAL zx_tmp_3d(iim,jjm+1,klev),zx_tmp_2d(iim,jjm+1)
61
62c   Couche limite:
63c   --------------
64c
[1067]65      REAL cdragh(klon)          ! cdrag
66      REAL pcoefh(klon,klev)     ! coeff melange CL
67      REAL pcoefh_buf(klon,klev) ! coeff melange CL + cdrag
[524]68      REAL yv1(klon)
69      REAL yu1(klon),pphis(klon),paire(klon)
[541]70
71c   Les Thermiques : (Abderr 25 11 02)
72c   ---------------
73      REAL pfm_therm(klon,klev+1)
[766]74      real fm_therm1(klon,klev)
[541]75      REAL pentr_therm(klon,klev)
[766]76   
77      REAL,allocatable,save :: entr_therm(:,:)
78      REAL,allocatable,save :: fm_therm(:,:)
79c$OMP THREADPRIVATE(entr_therm)
80c$OMP THREADPRIVATE(fm_therm)
[524]81c
82c   Lessivage:
83c   ----------
84c
[541]85      REAL frac_impa(klon,klev)
86      REAL frac_nucl(klon,klev)
[524]87c
88c Arguments necessaires pour les sources et puits de traceur
89C
90      real ftsol(klon,nbsrf)  ! Temperature du sol (surf)(Kelvin)
91      real pctsrf(klon,nbsrf) ! Pourcentage de sol f(nature du sol)
92c======================================================================
93c
94      INTEGER i, k
95c
[766]96      REAL,allocatable,save :: mfu(:,:)  ! flux de masse dans le panache montant
97      REAL,allocatable,save :: mfd(:,:)  ! flux de masse dans le panache descendant
98      REAL,allocatable,save :: en_u(:,:) ! flux entraine dans le panache montant
99      REAL,allocatable,save :: de_u(:,:) ! flux detraine dans le panache montant
100      REAL,allocatable,save :: en_d(:,:) ! flux entraine dans le panache descendant
101      REAL,allocatable,save :: de_d(:,:) ! flux detraine dans le panache descendant
102      REAL,allocatable,save :: coefh(:,:) ! flux detraine dans le panache descendant
[524]103
[766]104      REAL,allocatable,save :: pyu1(:)
105      REAL,allocatable,save :: pyv1(:)
106      REAL,allocatable,save :: pftsol(:,:)
107      REAL,allocatable,save :: ppsrf(:,:)
108c$OMP THREADPRIVATE(mfu,mfd,en_u,de_u,en_d,de_d,coefh)
109c$OMP THREADPRIVATE(pyu1,pyv1,pftsol,ppsrf)
[524]110      real pftsol1(klon),pftsol2(klon),pftsol3(klon),pftsol4(klon)
111      real ppsrf1(klon),ppsrf2(klon),ppsrf3(klon),ppsrf4(klon)
112
113      REAL dtcum
114
115      integer iadvtr,irec
116      real zmin,zmax
[541]117      logical ok_sync
118 
[766]119      save dtcum
[524]120      save iadvtr,irec
[766]121c$OMP THREADPRIVATE(dtcum,iadvtr,irec)
[524]122      data iadvtr,irec/0,1/
[766]123      logical,save :: first=.true.
124c$OMP THREADPRIVATE(first)
[524]125c
126c   Couche limite:
127c======================================================================
128
[1067]129c Dans le meme vecteur on recombine le drag et les coeff d'echange
130      pcoefh_buf(:,1)      = cdragh(:)
131      pcoefh_buf(:,2:klev) = pcoefh(:,2:klev)
132
[541]133      ok_sync = .true.
134        print*,'Dans phystokenc.F'
[524]135      print*,'iadvtr= ',iadvtr
136      print*,'istphy= ',istphy
137      print*,'istdyn= ',istdyn
138
[766]139      if (first) then
140     
141        allocate( t(klon,klev))
142        allocate( mfu(klon,klev)) 
143        allocate( mfd(klon,klev)) 
144        allocate( en_u(klon,klev))
145        allocate( de_u(klon,klev))
146        allocate( en_d(klon,klev))
147        allocate( de_d(klon,klev))
148        allocate( coefh(klon,klev))
149        allocate( entr_therm(klon,klev))
150        allocate( fm_therm(klon,klev))
151        allocate( pyu1(klon))
152        allocate( pyv1(klon))
153        allocate( pftsol(klon,nbsrf))
154        allocate( ppsrf(klon,nbsrf))
155 
156        first=.false.
157      endif
158     
[524]159      IF (iadvtr.eq.0) THEN
160       
161        CALL initphysto('phystoke',
[1146]162     . rlon,rlat,dtime, dtime*istphy,dtime*istphy,nqtot,physid)
[524]163       
164        write(*,*) 'apres initphysto ds phystokenc'
165
166       
167      ENDIF
168c
[541]169      ndex2d = 0
170      ndex3d = 0
171      i=itap
[766]172cym      CALL gr_fi_ecrit(1,klon,iim,jjm+1,pphis,zx_tmp_2d)
173      CALL histwrite_phy(physid,"phis",i,pphis)
[541]174c
175      i=itap
[766]176cym      CALL gr_fi_ecrit(1,klon,iim,jjm+1,paire,zx_tmp_2d)
177      CALL histwrite_phy(physid,"aire",i,paire)
[541]178
[524]179      iadvtr=iadvtr+1
180c
181      if (mod(iadvtr,istphy).eq.1.or.istphy.eq.1) then
182        print*,'reinitialisation des champs cumules
183     s          a iadvtr=',iadvtr
184         do k=1,klev
185            do i=1,klon
186               mfu(i,k)=0.
187               mfd(i,k)=0.
188               en_u(i,k)=0.
189               de_u(i,k)=0.
190               en_d(i,k)=0.
191               de_d(i,k)=0.
192               coefh(i,k)=0.
[541]193                t(i,k)=0.
194                fm_therm(i,k)=0.
195               entr_therm(i,k)=0.
[524]196            enddo
197         enddo
198         do i=1,klon
199            pyv1(i)=0.
200            pyu1(i)=0.
201         end do
202         do k=1,nbsrf
203             do i=1,klon
204               pftsol(i,k)=0.
205               ppsrf(i,k)=0.
206            enddo
207         enddo
208
209         dtcum=0.
210      endif
211
212      do k=1,klev
213         do i=1,klon
214            mfu(i,k)=mfu(i,k)+pmfu(i,k)*pdtphys
215            mfd(i,k)=mfd(i,k)+pmfd(i,k)*pdtphys
216            en_u(i,k)=en_u(i,k)+pen_u(i,k)*pdtphys
217            de_u(i,k)=de_u(i,k)+pde_u(i,k)*pdtphys
218            en_d(i,k)=en_d(i,k)+pen_d(i,k)*pdtphys
219            de_d(i,k)=de_d(i,k)+pde_d(i,k)*pdtphys
[1067]220            coefh(i,k)=coefh(i,k)+pcoefh_buf(i,k)*pdtphys
[541]221                t(i,k)=t(i,k)+pt(i,k)*pdtphys
222       fm_therm(i,k)=fm_therm(i,k)+pfm_therm(i,k)*pdtphys
223       entr_therm(i,k)=entr_therm(i,k)+pentr_therm(i,k)*pdtphys
[524]224         enddo
225      enddo
226         do i=1,klon
227            pyv1(i)=pyv1(i)+yv1(i)*pdtphys
228            pyu1(i)=pyu1(i)+yu1(i)*pdtphys
229         end do
230         do k=1,nbsrf
231             do i=1,klon
232               pftsol(i,k)=pftsol(i,k)+ftsol(i,k)*pdtphys
233               ppsrf(i,k)=ppsrf(i,k)+pctsrf(i,k)*pdtphys
234            enddo
235         enddo
236
237      dtcum=dtcum+pdtphys
[541]238
239      IF(mod(iadvtr,istphy).eq.0) THEN
[524]240c
241c   normalisation par le temps cumule
242         do k=1,klev
243            do i=1,klon
244               mfu(i,k)=mfu(i,k)/dtcum
245               mfd(i,k)=mfd(i,k)/dtcum
246               en_u(i,k)=en_u(i,k)/dtcum
247               de_u(i,k)=de_u(i,k)/dtcum
248               en_d(i,k)=en_d(i,k)/dtcum
249               de_d(i,k)=de_d(i,k)/dtcum
250               coefh(i,k)=coefh(i,k)/dtcum
[541]251c Unitel a enlever
252              t(i,k)=t(i,k)/dtcum       
253               fm_therm(i,k)=fm_therm(i,k)/dtcum
254               entr_therm(i,k)=entr_therm(i,k)/dtcum
[524]255            enddo
256         enddo
257         do i=1,klon
258            pyv1(i)=pyv1(i)/dtcum
259            pyu1(i)=pyu1(i)/dtcum
260         end do
[541]261         do k=1,nbsrf
[524]262             do i=1,klon
263               pftsol(i,k)=pftsol(i,k)/dtcum
264               pftsol1(i) = pftsol(i,1)
265               pftsol2(i) = pftsol(i,2)
266               pftsol3(i) = pftsol(i,3)
267               pftsol4(i) = pftsol(i,4)
268
[541]269               ppsrf(i,k)=ppsrf(i,k)/dtcum
[524]270               ppsrf1(i) = ppsrf(i,1)
271               ppsrf2(i) = ppsrf(i,2)
272               ppsrf3(i) = ppsrf(i,3)
273               ppsrf4(i) = ppsrf(i,4)
274
275            enddo
[541]276         enddo
[524]277c
278c   ecriture des champs
279c
280         irec=irec+1
281
282ccccc
[766]283cym         CALL gr_fi_ecrit(klev,klon,iim,jjm+1, t, zx_tmp_3d)
284         CALL histwrite_phy(physid,"t",itap,t)
[524]285
[766]286cym         CALL gr_fi_ecrit(klev,klon,iim,jjm+1, mfu, zx_tmp_3d)
287      CALL histwrite_phy(physid,"mfu",itap,mfu)
288cym     CALL gr_fi_ecrit(klev,klon,iim,jjm+1, mfd, zx_tmp_3d)
289      CALL histwrite_phy(physid,"mfd",itap,mfd)
290cym        CALL gr_fi_ecrit(klev,klon,iim,jjm+1, en_u, zx_tmp_3d)
291      CALL histwrite_phy(physid,"en_u",itap,en_u)
292cym        CALL gr_fi_ecrit(klev,klon,iim,jjm+1, de_u, zx_tmp_3d)
293      CALL histwrite_phy(physid,"de_u",itap,de_u)
294cym        CALL gr_fi_ecrit(klev,klon,iim,jjm+1, en_d, zx_tmp_3d)
295      CALL histwrite_phy(physid,"en_d",itap,en_d)
296cym        CALL gr_fi_ecrit(klev,klon,iim,jjm+1, de_d, zx_tmp_3d)       
297      CALL histwrite_phy(physid,"de_d",itap,de_d)
298cym        CALL gr_fi_ecrit(klev,klon,iim,jjm+1, coefh, zx_tmp_3d)         
299      CALL histwrite_phy(physid,"coefh",itap,coefh)     
[541]300
301c ajou...
302        do k=1,klev
303           do i=1,klon
304         fm_therm1(i,k)=fm_therm(i,k)   
305           enddo
306        enddo
307
[766]308cym      CALL gr_fi_ecrit(klev,klon,iim,jjm+1, fm_therm1, zx_tmp_3d)
309      CALL histwrite_phy(physid,"fm_th",itap,fm_therm1)
[541]310c
[766]311cym      CALL gr_fi_ecrit(klev,klon,iim,jjm+1, entr_therm, zx_tmp_3d)
312      CALL histwrite_phy(physid,"en_th",itap,entr_therm)
[524]313cccc
[766]314cym       CALL gr_fi_ecrit(klev,klon,iim,jjm+1,frac_impa,zx_tmp_3d)
315        CALL histwrite_phy(physid,"frac_impa",itap,frac_impa)
[524]316
[766]317cym        CALL gr_fi_ecrit(klev,klon,iim,jjm+1,frac_nucl,zx_tmp_3d)
318        CALL histwrite_phy(physid,"frac_nucl",itap,frac_nucl)
[541]319 
[766]320cym        CALL gr_fi_ecrit(1, klon,iim,jjm+1, pyu1,zx_tmp_2d)
321      CALL histwrite_phy(physid,"pyu1",itap,pyu1)
[541]322       
[766]323cym     CALL gr_fi_ecrit(1, klon,iim,jjm+1, pyv1,zx_tmp_2d)
324      CALL histwrite_phy(physid,"pyv1",itap,pyv1)
[541]325       
[766]326cym     CALL gr_fi_ecrit(1,klon,iim,jjm+1, pftsol1, zx_tmp_2d)
327      CALL histwrite_phy(physid,"ftsol1",itap,pftsol1)
328cym         CALL gr_fi_ecrit(1,klon,iim,jjm+1, pftsol2, zx_tmp_2d)
329      CALL histwrite_phy(physid,"ftsol2",itap,pftsol2)
330cym          CALL gr_fi_ecrit(1,klon,iim,jjm+1, pftsol3, zx_tmp_2d)
331      CALL histwrite_phy(physid,"ftsol3",itap,pftsol3)
332cym         CALL gr_fi_ecrit(1,klon,iim,jjm+1, pftsol4, zx_tmp_2d)
333      CALL histwrite_phy(physid,"ftsol4",itap,pftsol4)
[524]334
[766]335cym        CALL gr_fi_ecrit(1,klon,iim,jjm+1, ppsrf1, zx_tmp_2d)
336      CALL histwrite_phy(physid,"psrf1",itap,ppsrf1)
337cym        CALL gr_fi_ecrit(1,klon,iim,jjm+1, ppsrf2, zx_tmp_2d)
338      CALL histwrite_phy(physid,"psrf2",itap,ppsrf2)
339cym        CALL gr_fi_ecrit(1,klon,iim,jjm+1, ppsrf3, zx_tmp_2d)
340      CALL histwrite_phy(physid,"psrf3",itap,ppsrf3)
341cym        CALL gr_fi_ecrit(1,klon,iim,jjm+1, ppsrf4, zx_tmp_2d)
342      CALL histwrite_phy(physid,"psrf4",itap,ppsrf4)
[524]343
[766]344c$OMP MASTER
[541]345      if (ok_sync) call histsync(physid)
[766]346c$OMP END MASTER
[541]347c     if (ok_sync) call histsync
348       
[524]349c
[541]350cAA Test sur la valeur des coefficients de lessivage
[524]351c
352         zmin=1e33
353         zmax=-1e33
354         do k=1,klev
355            do i=1,klon
356                  zmax=max(zmax,frac_nucl(i,k))
357                  zmin=min(zmin,frac_nucl(i,k))
358            enddo
359         enddo
360         Print*,'------ coefs de lessivage (min et max) --------'
361         Print*,'facteur de nucleation ',zmin,zmax
362         zmin=1e33
363         zmax=-1e33
364         do k=1,klev
365            do i=1,klon
366                  zmax=max(zmax,frac_impa(i,k))
367                  zmin=min(zmin,frac_impa(i,k))
368            enddo
369         enddo
370         Print*,'facteur d impaction ',zmin,zmax
371
[541]372      ENDIF
[524]373
[541]374c   reinitialisation des champs cumules
375        go to 768
376      if (mod(iadvtr,istphy).eq.1) then
377         do k=1,klev
378            do i=1,klon
379               mfu(i,k)=0.
380               mfd(i,k)=0.
381               en_u(i,k)=0.
382               de_u(i,k)=0.
383               en_d(i,k)=0.
384               de_d(i,k)=0.
385               coefh(i,k)=0.
386               t(i,k)=0.
387               fm_therm(i,k)=0.
388               entr_therm(i,k)=0.
389            enddo
390         enddo
391         do i=1,klon
392            pyv1(i)=0.
393            pyu1(i)=0.
394         end do
395         do k=1,nbsrf
396             do i=1,klon
397               pftsol(i,k)=0.
398               ppsrf(i,k)=0.
399            enddo
400         enddo
[524]401
[541]402         dtcum=0.
403      endif
404
405      do k=1,klev
406         do i=1,klon
407            mfu(i,k)=mfu(i,k)+pmfu(i,k)*pdtphys
408            mfd(i,k)=mfd(i,k)+pmfd(i,k)*pdtphys
409            en_u(i,k)=en_u(i,k)+pen_u(i,k)*pdtphys
410            de_u(i,k)=de_u(i,k)+pde_u(i,k)*pdtphys
411            en_d(i,k)=en_d(i,k)+pen_d(i,k)*pdtphys
412            de_d(i,k)=de_d(i,k)+pde_d(i,k)*pdtphys
[1067]413            coefh(i,k)=coefh(i,k)+pcoefh_buf(i,k)*pdtphys
[541]414                t(i,k)=t(i,k)+pt(i,k)*pdtphys
415       fm_therm(i,k)=fm_therm(i,k)+pfm_therm(i,k)*pdtphys
416       entr_therm(i,k)=entr_therm(i,k)+pentr_therm(i,k)*pdtphys
417         enddo
418      enddo
419         do i=1,klon
420            pyv1(i)=pyv1(i)+yv1(i)*pdtphys
421            pyu1(i)=pyu1(i)+yu1(i)*pdtphys
422         end do
423         do k=1,nbsrf
424             do i=1,klon
425               pftsol(i,k)=pftsol(i,k)+ftsol(i,k)*pdtphys
426               ppsrf(i,k)=ppsrf(i,k)+pctsrf(i,k)*pdtphys
427            enddo
428         enddo
429
430      dtcum=dtcum+pdtphys
431768   continue
432
[524]433      RETURN
434      END
Note: See TracBrowser for help on using the repository browser.