source: LMDZ6/trunk/libf/phylmd/lmdz_ratqs_main.F90 @ 4706

Last change on this file since 4706 was 4613, checked in by fhourdin, 17 months ago

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