source: LMDZ6/branches/LMDZ_cdrag_LSCE/libf/phylmd/calcratqs.F90 @ 5428

Last change on this file since 5428 was 4669, checked in by Laurent Fairhead, 16 months ago

Merged with trunk revision 4586 corresponding to june 2023 testing
  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
File size: 9.6 KB
Line 
1SUBROUTINE calcratqs(klon,klev,prt_level,lunout,       &
2           iflag_ratqs,iflag_con,iflag_cld_th,pdtphys, &
3           ratqsbas,ratqshaut,ratqsp0,ratqsdp, &
4           tau_ratqs,fact_cldcon,wake_s, wake_deltaq,   &
5           ptconv,ptconvth,clwcon0th, rnebcon0th,       &
6           paprs,pplay,t_seri,q_seri,                   &
7           qtc_cv, sigt_cv, zqsat,             &
8           tke,tke_dissip,lmix,wprime, &
9           t2m,q2m,fm_therm,cell_area,&
10           ratqs,ratqsc,ratqs_inter)
11
12
13USE indice_sol_mod
14USE phys_state_var_mod, ONLY: pctsrf
15USE calcratqs_multi_mod, ONLY: calcratqs_inter, calcratqs_oro, calcratqs_hetero, calcratqs_tke
16
17implicit none
18
19!========================================================================
20! Computation of ratqs, the width of the subrid scale water distribution
21! (normalized by the mean value)
22! Various options controled by flags iflag_con and iflag_ratqs
23! F Hourdin 2012/12/06
24!========================================================================
25
26! Declarations
27
28! Input
29integer,intent(in) :: klon,klev,prt_level,lunout
30integer,intent(in) :: iflag_con,iflag_cld_th,iflag_ratqs
31real,intent(in) :: pdtphys,ratqsbas,ratqshaut,fact_cldcon,tau_ratqs
32real,intent(in) :: ratqsp0, ratqsdp
33real, dimension(klon,klev+1),intent(in) :: paprs,tke,tke_dissip,lmix,wprime
34real, dimension(klon,klev),intent(in) :: pplay,t_seri,q_seri,zqsat,fm_therm, qtc_cv, sigt_cv
35logical, dimension(klon,klev),intent(in) :: ptconv
36real, dimension(klon,klev),intent(in) :: rnebcon0th,clwcon0th
37real, dimension(klon,klev),intent(in) :: wake_deltaq,wake_s
38real, dimension(klon,nbsrf),intent(in) :: t2m,q2m
39real, dimension(klon), intent(in) :: cell_area
40! Output
41real, dimension(klon,klev),intent(inout) :: ratqs,ratqsc,ratqs_inter
42
43logical, dimension(klon,klev),intent(inout) :: ptconvth
44
45! local
46integer i,k
47real, dimension(klon,klev) :: ratqss
48real facteur,zfratqs1,zfratqs2
49real, dimension(klon,klev) :: ratqs_hetero,ratqs_oro,ratqs_tke
50real resol,resolmax,fact
51
52!-------------------------------------------------------------------------
53!  Caclul des ratqs
54!-------------------------------------------------------------------------
55
56!      print*,'calcul des ratqs'
57!   ratqs convectifs a l'ancienne en fonction de q(z=0)-q / q
58!   ----------------
59!   on ecrase le tableau ratqsc calcule par clouds_gno
60      if (iflag_cld_th.eq.1) then
61         do k=1,klev
62         do i=1,klon
63            if(ptconv(i,k)) then
64              ratqsc(i,k)=ratqsbas &
65              +fact_cldcon*(q_seri(i,1)-q_seri(i,k))/q_seri(i,k)
66            else
67               ratqsc(i,k)=0.
68            endif
69         enddo
70         enddo
71
72!-----------------------------------------------------------------------
73!  par nversion de la fonction log normale
74!-----------------------------------------------------------------------
75      else if (iflag_cld_th.eq.4) then
76         ptconvth(:,:)=.false.
77         ratqsc(:,:)=0.
78         if(prt_level.ge.9) print*,'avant clouds_gno thermique'
79         call clouds_gno &
80         (klon,klev,q_seri,zqsat,clwcon0th,ptconvth,ratqsc,rnebcon0th)
81         if(prt_level.ge.9) print*,' CLOUDS_GNO OK'
82       
83       endif
84
85!   ratqs stables
86!   -------------
87
88      if (iflag_ratqs.eq.0) then
89
90! Le cas iflag_ratqs=0 correspond a la version IPCC 2005 du modele.
91         do k=1,klev
92            do i=1, klon
93               ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)* &
94               min((paprs(i,1)-pplay(i,k))/(paprs(i,1)-30000.),1.)
95            enddo
96         enddo
97
98! Pour iflag_ratqs=1 ou 2, le ratqs est constant au dessus de
99! 300 hPa (ratqshaut), varie lineariement en fonction de la pression
100! entre 600 et 300 hPa et est soit constant (ratqsbas) pour iflag_ratqs=1
101! soit lineaire (entre 0 a la surface et ratqsbas) pour iflag_ratqs=2
102! Il s'agit de differents tests dans la phase de reglage du modele
103! avec thermiques.
104
105      else if (iflag_ratqs.eq.1) then
106
107         do k=1,klev
108            do i=1, klon
109               if (pplay(i,k).ge.60000.) then
110                  ratqss(i,k)=ratqsbas
111               else if ((pplay(i,k).ge.30000.).and.(pplay(i,k).lt.60000.)) then
112                  ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*(60000.-pplay(i,k))/(60000.-30000.)
113               else
114                  ratqss(i,k)=ratqshaut
115               endif
116            enddo
117         enddo
118
119      else if (iflag_ratqs.eq.2) then
120
121         do k=1,klev
122            do i=1, klon
123               if (pplay(i,k).ge.60000.) then
124                  ratqss(i,k)=ratqsbas*(paprs(i,1)-pplay(i,k))/(paprs(i,1)-60000.)
125               else if ((pplay(i,k).ge.30000.).and.(pplay(i,k).lt.60000.)) then
126                    ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*(60000.-pplay(i,k))/(60000.-30000.)
127               else
128                    ratqss(i,k)=ratqshaut
129               endif
130            enddo
131         enddo
132
133      else if (iflag_ratqs==3) then
134         do k=1,klev
135           ratqss(:,k)=ratqsbas+(ratqshaut-ratqsbas) &
136           *min( ((paprs(:,1)-pplay(:,k))/70000.)**2 , 1. )
137         enddo
138
139      else if (iflag_ratqs==4) then
140         do k=1,klev
141           ratqss(:,k)=ratqsbas+0.5*(ratqshaut-ratqsbas) &
142!          *( tanh( (50000.-pplay(:,k))/20000.) + 1.)
143           *( tanh( (ratqsp0-pplay(:,k))/ratqsdp) + 1.)
144         enddo
145
146
147      else if (iflag_ratqs==5) then
148! Dependency of ratqs on model resolution
149! Audran, Meryl, Lea, Gwendal and Etienne
150! April 2023
151        resolmax=sqrt(maxval(cell_area))
152         do k=1,klev
153            do i=1,klon
154              resol=sqrt(cell_area(i))
155              fact=sqrt(resol/resolmax)
156              ratqss(i,k)=ratqsbas*fact+0.5*(ratqshaut-ratqsbas)*fact &
157              *( tanh( (ratqsp0-pplay(i,k))/ratqsdp) + 1.)
158           enddo
159         enddo
160
161
162       else if (iflag_ratqs .GT. 9) then
163 
164       ! interactive ratqs calculations that depend on cold pools, orography, surface heterogeneity and small-scale turbulence
165       ! This should help getting a more realistic ratqs in the low and mid troposphere
166       ! We however need a "background" ratqs to account for subgrid distribution of qt (or qt/qs)
167       ! in the high troposphere
168       
169       ! background ratqs and initialisations
170          do k=1,klev
171             do i=1,klon
172              ratqss(i,k)=ratqsbas+0.5*(ratqshaut-ratqsbas) &
173              *( tanh( (ratqsp0-pplay(i,k))/ratqsdp) + 1.)
174              ratqss(i,k)=max(ratqss(i,k),0.0)
175
176              ratqs_hetero(i,k)=0.
177              ratqs_oro(i,k)=0.
178              ratqs_tke(i,k)=0.
179              ratqs_inter(i,k)=0
180             enddo
181          enddo
182     
183          if (iflag_ratqs .EQ. 10) then
184             ! interactive ratqs in presence of cold pools     
185             call calcratqs_inter(klon,klev,iflag_ratqs,pdtphys,ratqsbas,wake_deltaq,wake_s,q_seri,qtc_cv, sigt_cv,ratqs_inter)
186             do k=1,klev
187                do i=1,klon
188                    ratqs_inter(i,k)=ratqs_inter(i,k)-0.5*ratqs_inter(i,k)*(tanh((ratqsp0-pplay(i,k))/ratqsdp)+1.)
189                enddo
190             enddo
191             ratqss=ratqss+ratqs_inter
192          else if (iflag_ratqs .EQ. 11) then
193            ! interactive ratqs with several sources
194            call calcratqs_inter(klon,klev,iflag_ratqs,pdtphys,ratqsbas,wake_deltaq,wake_s,q_seri,qtc_cv, sigt_cv,ratqs_inter)
195             ratqss=ratqss+ratqs_inter
196          else if (iflag_ratqs .EQ. 12) then
197             ! contribution of surface heterogeneities to ratqs
198             call calcratqs_hetero(klon,klev,t2m,q2m,t_seri,q_seri,pplay,paprs,ratqs_hetero)
199             ratqss=ratqss+ratqs_hetero
200          else if (iflag_ratqs .EQ. 13) then
201             ! contribution of ubgrid orography to ratqs
202             call calcratqs_oro(klon,klev,zqsat,t_seri,pplay,paprs,ratqs_oro)
203             ratqss=ratqss+ratqs_oro
204          else if (iflag_ratqs .EQ. 14) then
205             ! effect of subgrid-scale TKE on ratqs (in development)
206             call calcratqs_tke(klon,klev,pdtphys,t_seri,q_seri,zqsat,pplay,paprs,tke,tke_dissip,lmix,wprime,ratqs_tke)     
207             ratqss=ratqss+ratqs_tke
208          endif
209         
210     
211      endif
212
213
214!  ratqs final
215!  -----------
216
217      if (iflag_cld_th.eq.1 .or.iflag_cld_th.eq.2.or.iflag_cld_th.eq.4) then
218
219! On ajoute une constante au ratqsc*2 pour tenir compte de
220! fluctuations turbulentes de petite echelle
221
222         do k=1,klev
223            do i=1,klon
224               if ((fm_therm(i,k).gt.1.e-10)) then
225                  ratqsc(i,k)=sqrt(ratqsc(i,k)**2+0.05**2)
226               endif
227            enddo
228         enddo
229
230!   les ratqs sont une combinaison de ratqss et ratqsc
231       if(prt_level.ge.9) write(lunout,*)'PHYLMD NOUVEAU TAU_RATQS ',tau_ratqs
232
233         if (tau_ratqs>1.e-10) then
234            facteur=exp(-pdtphys/tau_ratqs)
235         else
236            facteur=0.
237         endif
238         ratqs(:,:)=ratqsc(:,:)*(1.-facteur)+ratqs(:,:)*facteur
239!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
240! FH 22/09/2009
241! La ligne ci-dessous faisait osciller le modele et donnait une solution
242! assymptotique bidon et dépendant fortement du pas de temps.
243!        ratqs(:,:)=sqrt(ratqs(:,:)**2+ratqss(:,:)**2)
244!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
245         ratqs(:,:)=max(ratqs(:,:),ratqss(:,:))
246      else if (iflag_cld_th<=6) then
247!   on ne prend que le ratqs stable pour fisrtilp
248         ratqs(:,:)=ratqss(:,:)
249      else
250          zfratqs1=exp(-pdtphys/10800.)
251          zfratqs2=exp(-pdtphys/10800.)
252          do k=1,klev
253             do i=1,klon
254                if (ratqsc(i,k).gt.1.e-10) then
255                   ratqs(i,k)=ratqs(i,k)*zfratqs2+(iflag_cld_th/100.)*ratqsc(i,k)*(1.-zfratqs2)
256                endif
257                ratqs(i,k)=min(ratqs(i,k)*zfratqs1+ratqss(i,k)*(1.-zfratqs1),0.5)
258             enddo
259          enddo
260      endif
261
262
263return
264end
Note: See TracBrowser for help on using the repository browser.