source: LMDZ4/trunk/libf/phylmd/thermcell_closure.F90 @ 982

Last change on this file since 982 was 972, checked in by lmdzadmin, 16 years ago

Version thermique FH/CRio
Ajout tests cas physiques non pris en comptes et ajout/enleve prints
Nouvelle routine thermcell_flux2.F90
IM

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 2.2 KB
Line 
1      SUBROUTINE thermcell_closure(ngrid,nlay,r_aspect,ptimestep,rho,  &
2     &   zlev,lalim,alim_star,zmax_sec,wmax_sec,zmax,wmax,f,lev_out)
3
4!-------------------------------------------------------------------------
5!thermcell_closure: fermeture, determination de f
6!-------------------------------------------------------------------------
7      IMPLICIT NONE
8
9#include "iniprint.h"
10      INTEGER ngrid,nlay
11      INTEGER ig,k       
12      REAL r_aspect,ptimestep
13      integer lev_out                           ! niveau pour les print
14
15      INTEGER lalim(ngrid)
16      REAL alim_star(ngrid,nlay)
17      REAL rho(ngrid,nlay)
18      REAL zlev(ngrid,nlay)
19      REAL zmax(ngrid),zmax_sec(ngrid)
20      REAL wmax(ngrid),wmax_sec(ngrid)
21      real zdenom
22
23      REAL alim_star2(ngrid)
24
25      REAL f(ngrid)
26
27      do ig=1,ngrid
28         alim_star2(ig)=0.
29      enddo
30      do ig=1,ngrid
31         if (alim_star(ig,1).LT.1.e-10) then
32            f(ig)=0.
33         else   
34             do k=1,lalim(ig)
35                alim_star2(ig)=alim_star2(ig)+alim_star(ig,k)**2  &
36     &                    /(rho(ig,k)*(zlev(ig,k+1)-zlev(ig,k)))
37             enddo
38             zdenom=max(500.,zmax(ig))*r_aspect*alim_star2(ig)
39             if (zdenom<1.e-14) then
40                print*,'ig=',ig
41                print*,'alim_star2',alim_star2(ig)
42                print*,'zmax',zmax(ig)
43                print*,'r_aspect',r_aspect
44                print*,'zdenom',zdenom
45                print*,'alim_star',alim_star(ig,:)
46                print*,'zmax_sec',zmax_sec(ig)
47                print*,'wmax_sec',wmax_sec(ig)
48                stop
49             endif
50             if ((zmax_sec(ig).gt.1.e-10).and.(1.eq.1)) then
51             f(ig)=wmax_sec(ig)/(max(500.,zmax_sec(ig))*r_aspect  &
52     &             *alim_star2(ig))
53!            f(ig)=f(ig)+(f0(ig)-f(ig))*exp((-ptimestep/  &
54!    &                     zmax_sec(ig))*wmax_sec(ig))
55             else
56             f(ig)=wmax(ig)/zdenom
57!            f(ig)=f(ig)+(f0(ig)-f(ig))*exp((-ptimestep/  &
58!     &                     zmax(ig))*wmax(ig))
59             endif
60         endif
61!         f0(ig)=f(ig)
62      enddo
63      if (prt_level.ge.1) print*,'apres fermeture'
64
65!
66      return
67      end
Note: See TracBrowser for help on using the repository browser.