source: LMDZ5/trunk/libf/phylmd/thermcell_alim.F90 @ 2402

Last change on this file since 2402 was 2392, checked in by fhourdin, 9 years ago

Retour a 1+1=2. Probleme dans thermcell_alim.
Bug fixing

File size: 3.9 KB
Line 
1!
2! $Id: thermcell_plume.F90 2311 2015-06-25 07:45:24Z emillour $
3!
4      SUBROUTINE thermcell_alim(flag,ngrid,klev,ztv,d_temp,zlev,alim_star,lalim)
5IMPLICIT NONE
6
7!--------------------------------------------------------------------------
8! FH : 2015/11/06
9! thermcell_alim: calcule la distribution verticale de l'alimentation
10! laterale a la base des panaches thermiques
11!--------------------------------------------------------------------------
12
13#include "YOMCST.h"
14#include "YOETHF.h"
15#include "FCTTRE.h"
16#include "thermcell.h"
17
18!      fort(10) ptimestep,ztv,zthl,po,zl,rhobarz,zlev,pplev,pphi,zpspsk,f0
19      INTEGER, INTENT(IN) :: ngrid,klev
20      REAL, INTENT(IN) :: ztv(ngrid,klev)
21      REAL, INTENT(IN) :: d_temp(ngrid)
22      REAL, INTENT(IN) :: zlev(ngrid,klev+1)
23      REAL, INTENT(OUT) :: alim_star(ngrid,klev)
24      INTEGER, INTENT(OUT) :: lalim(ngrid)
25      INTEGER, INTENT(IN) :: flag
26
27      REAL :: alim_star_tot(ngrid),zi(ngrid),zh(ngrid)
28      REAL :: zlay(ngrid,klev)
29      REAL ztv_parcel
30
31      INTEGER ig,l
32
33      REAL h,z,falim
34      falim(h,z)=0.2*((z-h)**5+h**5)
35
36
37!===================================================================
38
39   lalim(:)=1
40   alim_star_tot(:)=0.
41
42   IF (ngrid==1) PRINT*,'NEW ALIM flag=',flag
43
44!-------------------------------------------------------------------------
45! Definition de l'alimentation a l'origine dans thermcell_init
46!-------------------------------------------------------------------------
47   IF (flag==0) THEN ! CMIP5 version
48      do l=1,klev-1
49         do ig=1,ngrid
50            if (ztv(ig,l)> ztv(ig,l+1) .and. ztv(ig,1)>=ztv(ig,l) ) then
51               alim_star(ig,l)=MAX((ztv(ig,l)-ztv(ig,l+1)),0.)  &
52     &                       *sqrt(zlev(ig,l+1))
53               lalim(ig)=l+1
54               alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
55            endif
56         enddo
57      enddo
58      do l=1,klev
59         do ig=1,ngrid
60            if (alim_star_tot(ig) > 1.e-10 ) then
61               alim_star(ig,l)=alim_star(ig,l)/alim_star_tot(ig)
62            endif
63         enddo
64      enddo
65      alim_star_tot(:)=1.
66
67!-------------------------------------------------------------------------
68! Nouvelle definition avec possibilite d'introduire un DT en surface
69! On suppose que la forme du profile d'alimentation scale avec la hauteur
70! d'inversion calculée avec une particule partant de la premieere couche
71
72! Fonction  f(z) = z ( h - z ) , avec h = zi/3
73! On utilise l'integralle
74! Int_0^z f(z') dz' = z^2 ( h/2 - z/3 ) = falim(h,z)
75! Pour calculer l'alimentation des couches
76!-------------------------------------------------------------------------
77   ELSE
78! Computing inversion height zi and zh=zi/3.
79      zi(:)=0.
80! Il faut recalculer zlay qui n'est pas dispo dans thermcell_plume
81! A changer eventuellement.
82      do l=1,klev
83         zlay(:,l)=0.5*(zlev(:,l)+zlev(:,l+1))
84      enddo
85
86      do l=1,klev-1
87         do ig=1,ngrid
88            ztv_parcel=ztv(ig,1)+d_temp(ig)
89            if (ztv_parcel<ztv(ig,l+1) .and. lalim(ig)==1 ) THEN
90                lalim(ig)=l
91                zi(ig)=zlay(ig,l)+(zlay(ig,l+1)-zlay(ig,l))/(ztv(ig,l+1)-ztv(ig,l))*(ztv_parcel-ztv(ig,l))
92                IF (zi(ig)<0.) STOP
93            endif
94         enddo
95      enddo
96      zh(:)=zi(:)/2.
97      alim_star_tot(:)=0.
98      alim_star(:,:)=0.
99      lalim(:)=0
100      do l=1,klev-1
101         do ig=1,ngrid
102            if (zlev(ig,l+1)<=zh(ig)) THEN
103               alim_star(ig,l)=(falim(zh(ig),zlev(ig,l+1))-falim(zh(ig),zlev(ig,l)))/falim(zh(ig),zh(ig))
104               lalim(ig)=l
105            ELSE IF (zlev(ig,l)<=zh(ig)) THEN
106               alim_star(ig,l)=(falim(zh(ig),zh(ig))-falim(zh(ig),zlev(ig,l)))/falim(zh(ig),zh(ig))
107               lalim(ig)=l
108            ELSE
109               alim_star(ig,l)=0.
110            ENDIF
111         ENDDO
112         alim_star_tot(:)=alim_star_tot(:)+alim_star(:,l)
113      ENDDO
114      IF (ngrid==1) print*,'NEW ALIM CALCUL DE ZI ',alim_star_tot
115      alim_star_tot(:)=1.
116
117   ENDIF
118
119
120RETURN
121END
Note: See TracBrowser for help on using the repository browser.