source: LMDZ6/trunk/libf/phylmd/lmdz_ratqs_main.f90 @ 5816

Last change on this file since 5816 was 5816, checked in by rkazeroni, 2 months ago

For GPU porting of clouds_gno and clouds_bigauss routines:

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