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

Last change on this file since 1357 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
Line 
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,
7     I                   pfm_therm,pentr_therm,
8     I                   cdragh, pcoefh,yu1,yv1,ftsol,pctsrf,
9     I                   frac_impa,frac_nucl,
10     I                   pphis,paire,dtime,itap)
11      USE ioipsl
12      USE dimphy
13      USE infotrac, ONLY : nqtot
14      USE iophy
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
41      integer physid, itap
42      save physid
43c$OMP THREADPRIVATE(physid)
44      integer ndex2d(iim*(jjm+1)),ndex3d(iim*(jjm+1)*klev)
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
55      real pt(klon,klev)
56      REAL,allocatable,save :: t(:,:)
57c$OMP THREADPRIVATE(t)
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
65      REAL cdragh(klon)          ! cdrag
66      REAL pcoefh(klon,klev)     ! coeff melange CL
67      REAL pcoefh_buf(klon,klev) ! coeff melange CL + cdrag
68      REAL yv1(klon)
69      REAL yu1(klon),pphis(klon),paire(klon)
70
71c   Les Thermiques : (Abderr 25 11 02)
72c   ---------------
73      REAL pfm_therm(klon,klev+1)
74      real fm_therm1(klon,klev)
75      REAL pentr_therm(klon,klev)
76   
77      REAL,allocatable,save :: entr_therm(:,:)
78      REAL,allocatable,save :: fm_therm(:,:)
79c$OMP THREADPRIVATE(entr_therm)
80c$OMP THREADPRIVATE(fm_therm)
81c
82c   Lessivage:
83c   ----------
84c
85      REAL frac_impa(klon,klev)
86      REAL frac_nucl(klon,klev)
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
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
103
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)
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
117      logical ok_sync
118 
119      save dtcum
120      save iadvtr,irec
121c$OMP THREADPRIVATE(dtcum,iadvtr,irec)
122      data iadvtr,irec/0,1/
123      logical,save :: first=.true.
124c$OMP THREADPRIVATE(first)
125c
126c   Couche limite:
127c======================================================================
128
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
133      ok_sync = .true.
134        print*,'Dans phystokenc.F'
135      print*,'iadvtr= ',iadvtr
136      print*,'istphy= ',istphy
137      print*,'istdyn= ',istdyn
138
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     
159      IF (iadvtr.eq.0) THEN
160       
161        CALL initphysto('phystoke',
162     . rlon,rlat,dtime, dtime*istphy,dtime*istphy,nqtot,physid)
163       
164        write(*,*) 'apres initphysto ds phystokenc'
165
166       
167      ENDIF
168c
169      ndex2d = 0
170      ndex3d = 0
171      i=itap
172cym      CALL gr_fi_ecrit(1,klon,iim,jjm+1,pphis,zx_tmp_2d)
173      CALL histwrite_phy(physid,"phis",i,pphis)
174c
175      i=itap
176cym      CALL gr_fi_ecrit(1,klon,iim,jjm+1,paire,zx_tmp_2d)
177      CALL histwrite_phy(physid,"aire",i,paire)
178
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.
193                t(i,k)=0.
194                fm_therm(i,k)=0.
195               entr_therm(i,k)=0.
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
220            coefh(i,k)=coefh(i,k)+pcoefh_buf(i,k)*pdtphys
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
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
238
239      IF(mod(iadvtr,istphy).eq.0) THEN
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
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
255            enddo
256         enddo
257         do i=1,klon
258            pyv1(i)=pyv1(i)/dtcum
259            pyu1(i)=pyu1(i)/dtcum
260         end do
261         do k=1,nbsrf
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
269               ppsrf(i,k)=ppsrf(i,k)/dtcum
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
276         enddo
277c
278c   ecriture des champs
279c
280         irec=irec+1
281
282ccccc
283cym         CALL gr_fi_ecrit(klev,klon,iim,jjm+1, t, zx_tmp_3d)
284         CALL histwrite_phy(physid,"t",itap,t)
285
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)     
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
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)
310c
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)
313cccc
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)
316
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)
319 
320cym        CALL gr_fi_ecrit(1, klon,iim,jjm+1, pyu1,zx_tmp_2d)
321      CALL histwrite_phy(physid,"pyu1",itap,pyu1)
322       
323cym     CALL gr_fi_ecrit(1, klon,iim,jjm+1, pyv1,zx_tmp_2d)
324      CALL histwrite_phy(physid,"pyv1",itap,pyv1)
325       
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)
334
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)
343
344c$OMP MASTER
345      if (ok_sync) call histsync(physid)
346c$OMP END MASTER
347c     if (ok_sync) call histsync
348       
349c
350cAA Test sur la valeur des coefficients de lessivage
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
372      ENDIF
373
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
401
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
413            coefh(i,k)=coefh(i,k)+pcoefh_buf(i,k)*pdtphys
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
433      RETURN
434      END
Note: See TracBrowser for help on using the repository browser.