source: LMDZ4/trunk/libf/phytherm/phystokenc.F @ 839

Last change on this file since 839 was 814, checked in by Laurent Fairhead, 17 years ago

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