source: LMDZ4/trunk/libf/phytherm/thermcell_closure.F90 @ 965

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

Mise a jour de la physique avec thermiques avec la version de FH d'aout 2007
LF

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 2.2 KB
RevLine 
[814]1      SUBROUTINE thermcell_closure(ngrid,nlay,r_aspect,ptimestep,rho,  &
2     &   zlev,lalim,alim_star,zmax_sec,wmax_sec,zmax,wmax,f,f0,lev_out)
3
4!-------------------------------------------------------------------------
5!thermcell_closure: fermeture, determination de f
6!-------------------------------------------------------------------------
7      IMPLICIT NONE
8
9      INTEGER ngrid,nlay
10      INTEGER ig,k       
11      REAL r_aspect,ptimestep
12      integer lev_out                           ! niveau pour les print
13
14      INTEGER lalim(ngrid)
15      REAL alim_star(ngrid,nlay)
16      REAL rho(ngrid,nlay)
17      REAL zlev(ngrid,nlay)
18      REAL zmax(ngrid),zmax_sec(ngrid)
19      REAL wmax(ngrid),wmax_sec(ngrid)
[852]20      real zdenom
[814]21
22      REAL alim_star2(ngrid)
23
24      REAL f(ngrid)
25      REAL f0(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
[852]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
[814]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
[852]56             f(ig)=wmax(ig)/zdenom
[814]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 (lev_out.ge.1) print*,'apres fermeture'
64
65!
66      return
67      end
Note: See TracBrowser for help on using the repository browser.