Context Navigation

physiq.F @ 1378

Last change on this file since 1378 was 1378, checked in by musat, 14 years ago
Oups forgot some changes concerning the dynamical declaration of tracers IM
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 123.5 KB

Line
1	! $Id: physiq.F 1378 2010-05-06 10:57:26Z musat $
2	!
3	c#define IO_DEBUG
4
5	SUBROUTINE physiq (nlon,nlev,
6	. debut,lafin,jD_cur, jH_cur,pdtphys,
7	. paprs,pplay,pphi,pphis,presnivs,clesphy0,
8	. u,v,t,qx,
9	. flxmass_w,
10	. d_u, d_v, d_t, d_qx, d_ps
11	. , dudyn
12	. , PVteta)
13
14	USE ioipsl, only: histbeg, histvert, histdef, histend, histsync,
15	$ histwrite, ju2ymds, ymds2ju, ioget_year_len
16	USE comgeomphy
17	USE phys_cal_mod
18	USE write_field_phy
19	USE dimphy
20	USE infotrac
21	USE mod_grid_phy_lmdz
22	USE mod_phys_lmdz_para
23	USE iophy
24	USE misc_mod, mydebug=>debug
25	USE vampir
26	USE pbl_surface_mod, ONLY : pbl_surface
27	USE change_srf_frac_mod
28	USE surface_data, ONLY : type_ocean, ok_veget
29	USE phys_local_var_mod ! Variables internes non sauvegardees de la physique
30	USE phys_state_var_mod ! Variables sauvegardees de la physique
31	USE phys_output_var_mod ! Variables pour les ecritures des sorties
32	USE fonte_neige_mod, ONLY : fonte_neige_get_vars
33	USE phys_output_mod
34	use open_climoz_m, only: open_climoz ! ozone climatology from a file
35	use regr_pr_av_m, only: regr_pr_av
36	use netcdf95, only: nf95_close
37	cIM for NMC files
38	use netcdf, only: nf90_fill_real
39	use mod_phys_lmdz_mpi_data, only: is_mpi_root
40	USE aero_mod
41	use ozonecm_m, only: ozonecm ! ozone of J.-F. Royer
42	use conf_phys_m, only: conf_phys
43	use radlwsw_m, only: radlwsw
44	USE control_mod
45
46
47	IMPLICIT none
48	c======================================================================
49	c
50	c Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
51	c
52	c Objet: Moniteur general de la physique du modele
53	cAA Modifications quant aux traceurs :
54	cAA - uniformisation des parametrisations ds phytrac
55	cAA - stockage des moyennes des champs necessaires
56	cAA en mode traceur off-line
57	c======================================================================
58	c CLEFS CPP POUR LES IO
59	c =====================
60	#define histNMC
61	c#define histISCCP
62	c======================================================================
63	c modif ( P. Le Van , 12/10/98 )
64	c
65	c Arguments:
66	c
67	c nlon----input-I-nombre de points horizontaux
68	c nlev----input-I-nombre de couches verticales, doit etre egale a klev
69	c debut---input-L-variable logique indiquant le premier passage
70	c lafin---input-L-variable logique indiquant le dernier passage
71	c jD_cur -R-jour courant a l'appel de la physique (jour julien)
72	c jH_cur -R-heure courante a l'appel de la physique (jour julien)
73	c pdtphys-input-R-pas d'integration pour la physique (seconde)
74	c paprs---input-R-pression pour chaque inter-couche (en Pa)
75	c pplay---input-R-pression pour le mileu de chaque couche (en Pa)
76	c pphi----input-R-geopotentiel de chaque couche (g z) (reference sol)
77	c pphis---input-R-geopotentiel du sol
78	c presnivs-input_R_pressions approximat. des milieux couches ( en PA)
79	c u-------input-R-vitesse dans la direction X (de O a E) en m/s
80	c v-------input-R-vitesse Y (de S a N) en m/s
81	c t-------input-R-temperature (K)
82	c qx------input-R-humidite specifique (kg/kg) et d'autres traceurs
83	c d_t_dyn-input-R-tendance dynamique pour "t" (K/s)
84	c d_q_dyn-input-R-tendance dynamique pour "q" (kg/kg/s)
85	c flxmass_w -input-R- flux de masse verticale
86	c d_u-----output-R-tendance physique de "u" (m/s/s)
87	c d_v-----output-R-tendance physique de "v" (m/s/s)
88	c d_t-----output-R-tendance physique de "t" (K/s)
89	c d_qx----output-R-tendance physique de "qx" (kg/kg/s)
90	c d_ps----output-R-tendance physique de la pression au sol
91	cIM
92	c PVteta--output-R-vorticite potentielle a des thetas constantes
93	c======================================================================
94	#include "dimensions.h"
95	integer jjmp1
96	parameter (jjmp1=jjm+1-1/jjm)
97	integer iip1
98	parameter (iip1=iim+1)
99
100	#include "regdim.h"
101	#include "indicesol.h"
102	#include "dimsoil.h"
103	#include "clesphys.h"
104	#include "temps.h"
105	#include "iniprint.h"
106	#include "thermcell.h"
107	c======================================================================
108	LOGICAL ok_cvl ! pour activer le nouveau driver pour convection KE
109	PARAMETER (ok_cvl=.TRUE.)
110	LOGICAL ok_gust ! pour activer l'effet des gust sur flux surface
111	PARAMETER (ok_gust=.FALSE.)
112	integer iflag_radia ! active ou non le rayonnement (MPL)
113	save iflag_radia
114	c$OMP THREADPRIVATE(iflag_radia)
115	c======================================================================
116	LOGICAL check ! Verifier la conservation du modele en eau
117	PARAMETER (check=.FALSE.)
118	LOGICAL ok_stratus ! Ajouter artificiellement les stratus
119	PARAMETER (ok_stratus=.FALSE.)
120	c======================================================================
121	REAL amn, amx
122	INTEGER igout
123	c======================================================================
124	c Clef controlant l'activation du cycle diurne:
125	ccc LOGICAL cycle_diurne
126	ccc PARAMETER (cycle_diurne=.FALSE.)
127	c======================================================================
128	c Modele thermique du sol, a activer pour le cycle diurne:
129	ccc LOGICAL soil_model
130	ccc PARAMETER (soil_model=.FALSE.)
131	c======================================================================
132	c Dans les versions precedentes, l'eau liquide nuageuse utilisee dans
133	c le calcul du rayonnement est celle apres la precipitation des nuages.
134	c Si cette cle new_oliq est activee, ce sera une valeur moyenne entre
135	c la condensation et la precipitation. Cette cle augmente les impacts
136	c radiatifs des nuages.
137	ccc LOGICAL new_oliq
138	ccc PARAMETER (new_oliq=.FALSE.)
139	c======================================================================
140	c Clefs controlant deux parametrisations de l'orographie:
141	cc LOGICAL ok_orodr
142	ccc PARAMETER (ok_orodr=.FALSE.)
143	ccc LOGICAL ok_orolf
144	ccc PARAMETER (ok_orolf=.FALSE.)
145	c======================================================================
146	LOGICAL ok_journe ! sortir le fichier journalier
147	save ok_journe
148	c$OMP THREADPRIVATE(ok_journe)
149	c
150	LOGICAL ok_mensuel ! sortir le fichier mensuel
151	save ok_mensuel
152	c$OMP THREADPRIVATE(ok_mensuel)
153	c
154	LOGICAL ok_instan ! sortir le fichier instantane
155	save ok_instan
156	c$OMP THREADPRIVATE(ok_instan)
157	c
158	LOGICAL ok_LES ! sortir le fichier LES
159	save ok_LES
160	c$OMP THREADPRIVATE(ok_LES)
161	c
162	LOGICAL ok_region ! sortir le fichier regional
163	PARAMETER (ok_region=.FALSE.)
164	c======================================================================
165	real weak_inversion(klon),dthmin(klon)
166	real seuil_inversion
167	save seuil_inversion
168	c$OMP THREADPRIVATE(seuil_inversion)
169	integer iflag_ratqs
170	save iflag_ratqs
171	c$OMP THREADPRIVATE(iflag_ratqs)
172	REAL lambda_th(klon,klev),zz,znum,zden
173	REAL wmax_th(klon)
174	REAL zmax_th(klon)
175	REAL tau_overturning_th(klon)
176
177	integer lmax_th(klon)
178	integer limbas(klon)
179	real ratqscth(klon,klev)
180	real ratqsdiff(klon,klev)
181	real zqsatth(klon,klev)
182
183	c======================================================================
184	c
185	INTEGER ivap ! indice de traceurs pour vapeur d'eau
186	PARAMETER (ivap=1)
187	INTEGER iliq ! indice de traceurs pour eau liquide
188	PARAMETER (iliq=2)
189
190	c
191	c
192	c Variables argument:
193	c
194	INTEGER nlon
195	INTEGER nlev
196	REAL, intent(in):: jD_cur, jH_cur
197
198	REAL pdtphys
199	LOGICAL debut, lafin
200	REAL paprs(klon,klev+1)
201	REAL pplay(klon,klev)
202	REAL pphi(klon,klev)
203	REAL pphis(klon)
204	REAL presnivs(klev)
205	REAL znivsig(klev)
206	real pir
207
208	REAL u(klon,klev)
209	REAL v(klon,klev)
210	REAL t(klon,klev),theta(klon,klev)
211	REAL qx(klon,klev,nqtot)
212	REAL flxmass_w(klon,klev)
213	REAL omega(klon,klev) ! vitesse verticale en Pa/s
214	REAL d_u(klon,klev)
215	REAL d_v(klon,klev)
216	REAL d_t(klon,klev)
217	REAL d_qx(klon,klev,nqtot)
218	REAL d_ps(klon)
219	real da(klon,klev),phi(klon,klev,klev),mp(klon,klev)
220	!IM definition dynamique o_trac dans phys_output_open
221	type(ctrl_out) :: o_trac(nqtot)
222	c
223	cIM Amip2 PV a theta constante
224	c
225	INTEGER nbteta
226	PARAMETER(nbteta=3)
227	CHARACTER*3 ctetaSTD(nbteta)
228	DATA ctetaSTD/'350','380','405'/
229	SAVE ctetaSTD
230	c$OMP THREADPRIVATE(ctetaSTD)
231	REAL rtetaSTD(nbteta)
232	DATA rtetaSTD/350., 380., 405./
233	SAVE rtetaSTD
234	c$OMP THREADPRIVATE(rtetaSTD)
235	c
236	REAL PVteta(klon,nbteta)
237	REAL zx_tmp_3dte(iim,jjmp1,nbteta)
238	c
239	cMI Amip2 PV a theta constante
240
241	cym INTEGER klevp1, klevm1
242	cym PARAMETER(klevp1=klev+1,klevm1=klev-1)
243	cym#include "raddim.h"
244	c
245	c
246	cIM Amip2
247	c variables a une pression donnee
248	c
249	real rlevSTD(nlevSTD)
250	DATA rlevSTD/100000., 92500., 85000., 70000.,
251	.60000., 50000., 40000., 30000., 25000., 20000.,
252	.15000., 10000., 7000., 5000., 3000., 2000., 1000./
253	SAVE rlevstd
254	c$OMP THREADPRIVATE(rlevstd)
255	CHARACTER*4 clevSTD(nlevSTD)
256	DATA clevSTD/'1000','925 ','850 ','700 ','600 ',
257	.'500 ','400 ','300 ','250 ','200 ','150 ','100 ',
258	.'70 ','50 ','30 ','20 ','10 '/
259	SAVE clevSTD
260	c$OMP THREADPRIVATE(clevSTD)
261	c
262	CHARACTER*4 bb2
263	CHARACTER*2 bb3
264	c
265
266	real twriteSTD(klon,nlevSTD,nfiles)
267	real qwriteSTD(klon,nlevSTD,nfiles)
268	real rhwriteSTD(klon,nlevSTD,nfiles)
269	real phiwriteSTD(klon,nlevSTD,nfiles)
270	real uwriteSTD(klon,nlevSTD,nfiles)
271	real vwriteSTD(klon,nlevSTD,nfiles)
272	real wwriteSTD(klon,nlevSTD,nfiles)
273	cIM for NMC files
274	REAL geo500(klon)
275	real :: rlevSTD3(nlevSTD3)
276	DATA rlevSTD3/85000., 50000., 25000./
277	SAVE rlevSTD3
278	c$OMP THREADPRIVATE(rlevSTD3)
279	real :: rlevSTD8(nlevSTD8)
280	DATA rlevSTD8/100000., 85000., 70000., 50000., 25000., 10000.,
281	$ 5000., 1000./
282	SAVE rlevSTD8
283	c$OMP THREADPRIVATE(rlevSTD8)
284	real twriteSTD3(klon,nlevSTD3)
285	real qwriteSTD3(klon,nlevSTD3)
286	real rhwriteSTD3(klon,nlevSTD3)
287	real phiwriteSTD3(klon,nlevSTD3)
288	real uwriteSTD3(klon,nlevSTD3)
289	real vwriteSTD3(klon,nlevSTD3)
290	real wwriteSTD3(klon,nlevSTD3)
291	c
292	real tnondefSTD8(klon,nlevSTD8)
293	real twriteSTD8(klon,nlevSTD8)
294	real qwriteSTD8(klon,nlevSTD8)
295	real rhwriteSTD8(klon,nlevSTD8)
296	real phiwriteSTD8(klon,nlevSTD8)
297	real uwriteSTD8(klon,nlevSTD8)
298	real vwriteSTD8(klon,nlevSTD8)
299	real wwriteSTD8(klon,nlevSTD8)
300	c
301	c plevSTD3 END
302	c
303	c nout : niveau de output des variables a une pression donnee
304	logical oknondef(klon,nlevSTD,nout)
305	c
306	c les produits uvSTD, vqSTD, .., T2STD sont calcules
307	c a partir des valeurs instantannees toutes les 6 h
308	c qui sont moyennees sur le mois
309	c
310	#include "radopt.h"
311	c
312	c
313	c prw: precipitable water
314	real prw(klon)
315
316	REAL convliq(klon,klev) ! eau liquide nuageuse convective
317	REAL convfra(klon,klev) ! fraction nuageuse convective
318
319	REAL cldl_c(klon),cldm_c(klon),cldh_c(klon) !nuages bas, moyen et haut
320	REAL cldt_c(klon),cldq_c(klon) !nuage total, eau liquide integree
321	REAL cldl_s(klon),cldm_s(klon),cldh_s(klon) !nuages bas, moyen et haut
322	REAL cldt_s(klon),cldq_s(klon) !nuage total, eau liquide integree
323
324	INTEGER linv, kp1
325	c flwp, fiwp = Liquid Water Path & Ice Water Path (kg/m2)
326	c flwc, fiwc = Liquid Water Content & Ice Water Content (kg/kg)
327	REAL flwp(klon), fiwp(klon)
328	REAL flwc(klon,klev), fiwc(klon,klev)
329	REAL flwp_c(klon), fiwp_c(klon)
330	REAL flwc_c(klon,klev), fiwc_c(klon,klev)
331	REAL flwp_s(klon), fiwp_s(klon)
332	REAL flwc_s(klon,klev), fiwc_s(klon,klev)
333
334	cIM ISCCP simulator v3.4
335	c dans clesphys.h top_height, overlap
336	cv3.4
337	INTEGER debug, debugcol
338	cym INTEGER npoints
339	cym PARAMETER(npoints=klon)
340	c
341	INTEGER sunlit(klon) !sunlit=1 if day; sunlit=0 if night
342	INTEGER nregISCtot
343	PARAMETER(nregISCtot=1)
344	c
345	c imin_debut, nbpti, jmin_debut, nbptj : parametres pour sorties sur 1 region rectangulaire
346	c y compris pour 1 point
347	c imin_debut : indice minimum de i; nbpti : nombre de points en direction i (longitude)
348	c jmin_debut : indice minimum de j; nbptj : nombre de points en direction j (latitude)
349	INTEGER imin_debut, nbpti
350	INTEGER jmin_debut, nbptj
351	cIM parametres ISCCP BEG
352	INTEGER nbapp_isccp
353	! INTEGER nbapp_isccp,isccppas
354	! PARAMETER(isccppas=6) !appel du simulateurs tous les 6pas de temps de la physique
355	! !i.e. toutes les 3 heures
356	INTEGER n
357	INTEGER ifreq_isccp(napisccp), freqin_pdt(napisccp)
358	DATA ifreq_isccp/3/
359	SAVE ifreq_isccp
360	c$OMP THREADPRIVATE(ifreq_isccp)
361	CHARACTER*5 typinout(napisccp)
362	DATA typinout/'i3od'/
363	SAVE typinout
364	c$OMP THREADPRIVATE(typinout)
365	cIM verif boxptop BEG
366	CHARACTER*1 verticaxe(napisccp)
367	DATA verticaxe/'1'/
368	SAVE verticaxe
369	c$OMP THREADPRIVATE(verticaxe)
370	cIM verif boxptop END
371	INTEGER nvlev(napisccp)
372	c INTEGER nvlev
373	REAL t1, aa
374	REAL seed_re(klon,napisccp)
375	cym !!!! A voir plus tard
376	cym INTEGER iphy(iim,jjmp1)
377	cIM parametres ISCCP END
378	c
379	c ncol = nb. de sous-colonnes pour chaque maille du GCM
380	c ncolmx = No. max. de sous-colonnes pour chaque maille du GCM
381	c INTEGER ncol(napisccp), ncolmx, seed(klon,napisccp)
382	INTEGER,SAVE :: ncol(napisccp)
383	c$OMP THREADPRIVATE(ncol)
384	INTEGER ncolmx, seed(klon,napisccp)
385	REAL nbsunlit(nregISCtot,klon,napisccp) !nbsunlit : moyenne de sunlit
386	c PARAMETER(ncolmx=1500)
387	PARAMETER(ncolmx=300)
388	c
389	cIM verif boxptop BEG
390	REAL vertlev(ncolmx,napisccp)
391	cIM verif boxptop END
392	c
393	REAL,SAVE :: tautab_omp(0:255),tautab(0:255)
394	INTEGER,SAVE :: invtau_omp(-20:45000),invtau(-20:45000)
395	c$OMP THREADPRIVATE(tautab,invtau)
396	REAL emsfc_lw
397	PARAMETER(emsfc_lw=0.99)
398	c REAL ran0 ! type for random number fuction
399	c
400	REAL cldtot(klon,klev)
401	c variables de haut en bas pour le simulateur ISCCP
402	REAL dtau_s(klon,klev) !tau nuages startiformes
403	REAL dtau_c(klon,klev) !tau nuages convectifs
404	REAL dem_s(klon,klev) !emissivite nuages startiformes
405	REAL dem_c(klon,klev) !emissivite nuages convectifs
406	c
407	c variables de haut en bas pour le simulateur ISCCP
408	REAL pfull(klon,klev)
409	REAL phalf(klon,klev+1)
410	REAL qv(klon,klev)
411	REAL cc(klon,klev)
412	REAL conv(klon,klev)
413	REAL dtau_sH2B(klon,klev)
414	REAL dtau_cH2B(klon,klev)
415	REAL at(klon,klev)
416	REAL dem_sH2B(klon,klev)
417	REAL dem_cH2B(klon,klev)
418	c
419	INTEGER kmax, lmax, lmax3
420	PARAMETER(kmax=8, lmax=8, lmax3=3)
421	INTEGER kmaxm1, lmaxm1
422	PARAMETER(kmaxm1=kmax-1, lmaxm1=lmax-1)
423	INTEGER iimx7, jjmx7, jjmp1x7
424	PARAMETER(iimx7=iimkmaxm1, jjmx7=jjmlmaxm1,
425	.jjmp1x7=jjmp1*lmaxm1)
426	c
427	c output from ISCCP simulator
428	REAL fq_isccp(klon,kmaxm1,lmaxm1,napisccp)
429	REAL fq_is_true(klon,kmaxm1,lmaxm1,napisccp)
430	REAL totalcldarea(klon,napisccp)
431	REAL meanptop(klon,napisccp)
432	REAL meantaucld(klon,napisccp)
433	REAL boxtau(klon,ncolmx,napisccp)
434	REAL boxptop(klon,ncolmx,napisccp)
435	REAL zx_tmp_fi3d_bx(klon,ncolmx)
436	REAL zx_tmp_3d_bx(iim,jjmp1,ncolmx)
437	c
438	REAL cld_fi3d(klon,lmax3)
439	REAL cld_3d(iim,jjmp1,lmax3)
440	c
441	INTEGER iw, iwmax
442	REAL wmin, pas_w
443	c PARAMETER(wmin=-100.,pas_w=10.,iwmax=30)
444	cIM 051005 PARAMETER(wmin=-200.,pas_w=10.,iwmax=40)
445	PARAMETER(wmin=-100.,pas_w=10.,iwmax=20)
446	REAL o500(klon)
447	c
448
449	c sorties ISCCP
450
451	integer nid_isccp
452	save nid_isccp
453	c$OMP THREADPRIVATE(nid_isccp)
454
455	REAL zx_tau(kmaxm1), zx_pc(lmaxm1), zx_o500(iwmax)
456	DATA zx_tau/0.0, 0.3, 1.3, 3.6, 9.4, 23., 60./
457	SAVE zx_tau
458	DATA zx_pc/180., 310., 440., 560., 680., 800., 1000./
459	SAVE zx_pc
460	c$OMP THREADPRIVATE(zx_tau,zx_pc)
461	c cldtopres pression au sommet des nuages
462	REAL cldtopres(lmaxm1), cldtopres3(lmax3)
463	DATA cldtopres/180., 310., 440., 560., 680., 800., 1000./
464	DATA cldtopres3/440., 680., 1000./
465	SAVE cldtopres,cldtopres3
466	c$OMP THREADPRIVATE(cldtopres,cldtopres3)
467	cIM 051005 BEG
468	INTEGER komega, nhoriRD
469
470	c taulev: numero du niveau de tau dans les sorties ISCCP
471	CHARACTER *4 taulev(kmaxm1)
472	c DATA taulev/'tau1','tau2','tau3','tau4','tau5','tau6','tau7'/
473	DATA taulev/'tau0','tau1','tau2','tau3','tau4','tau5','tau6'/
474	CHARACTER *3 pclev(lmaxm1)
475	DATA pclev/'pc1','pc2','pc3','pc4','pc5','pc6','pc7'/
476	SAVE taulev,pclev
477	c$OMP THREADPRIVATE(taulev,pclev)
478	c
479	c cnameisccp
480	CHARACTER *29 cnameisccp(lmaxm1,kmaxm1)
481	cIM bad 151205 DATA cnameisccp/'pc< 50hPa, tau< 0.3',
482	DATA cnameisccp/'pc= 50-180hPa, tau< 0.3',
483	. 'pc= 180-310hPa, tau< 0.3',
484	. 'pc= 310-440hPa, tau< 0.3',
485	. 'pc= 440-560hPa, tau< 0.3',
486	. 'pc= 560-680hPa, tau< 0.3',
487	. 'pc= 680-800hPa, tau< 0.3',
488	. 'pc= 800-1000hPa, tau< 0.3',
489	. 'pc= 50-180hPa, tau= 0.3-1.3',
490	. 'pc= 180-310hPa, tau= 0.3-1.3',
491	. 'pc= 310-440hPa, tau= 0.3-1.3',
492	. 'pc= 440-560hPa, tau= 0.3-1.3',
493	. 'pc= 560-680hPa, tau= 0.3-1.3',
494	. 'pc= 680-800hPa, tau= 0.3-1.3',
495	. 'pc= 800-1000hPa, tau= 0.3-1.3',
496	. 'pc= 50-180hPa, tau= 1.3-3.6',
497	. 'pc= 180-310hPa, tau= 1.3-3.6',
498	. 'pc= 310-440hPa, tau= 1.3-3.6',
499	. 'pc= 440-560hPa, tau= 1.3-3.6',
500	. 'pc= 560-680hPa, tau= 1.3-3.6',
501	. 'pc= 680-800hPa, tau= 1.3-3.6',
502	. 'pc= 800-1000hPa, tau= 1.3-3.6',
503	. 'pc= 50-180hPa, tau= 3.6-9.4',
504	. 'pc= 180-310hPa, tau= 3.6-9.4',
505	. 'pc= 310-440hPa, tau= 3.6-9.4',
506	. 'pc= 440-560hPa, tau= 3.6-9.4',
507	. 'pc= 560-680hPa, tau= 3.6-9.4',
508	. 'pc= 680-800hPa, tau= 3.6-9.4',
509	. 'pc= 800-1000hPa, tau= 3.6-9.4',
510	. 'pc= 50-180hPa, tau= 9.4-23',
511	. 'pc= 180-310hPa, tau= 9.4-23',
512	. 'pc= 310-440hPa, tau= 9.4-23',
513	. 'pc= 440-560hPa, tau= 9.4-23',
514	. 'pc= 560-680hPa, tau= 9.4-23',
515	. 'pc= 680-800hPa, tau= 9.4-23',
516	. 'pc= 800-1000hPa, tau= 9.4-23',
517	. 'pc= 50-180hPa, tau= 23-60',
518	. 'pc= 180-310hPa, tau= 23-60',
519	. 'pc= 310-440hPa, tau= 23-60',
520	. 'pc= 440-560hPa, tau= 23-60',
521	. 'pc= 560-680hPa, tau= 23-60',
522	. 'pc= 680-800hPa, tau= 23-60',
523	. 'pc= 800-1000hPa, tau= 23-60',
524	. 'pc= 50-180hPa, tau> 60.',
525	. 'pc= 180-310hPa, tau> 60.',
526	. 'pc= 310-440hPa, tau> 60.',
527	. 'pc= 440-560hPa, tau> 60.',
528	. 'pc= 560-680hPa, tau> 60.',
529	. 'pc= 680-800hPa, tau> 60.',
530	. 'pc= 800-1000hPa, tau> 60.'/
531	SAVE cnameisccp
532	c$OMP THREADPRIVATE(cnameisccp)
533	c
534	c REAL zx_lonx7(iimx7), zx_latx7(jjmp1x7)
535	c INTEGER nhorix7
536	cIM: region='3d' <==> sorties en global
537	CHARACTER*3 region
538	PARAMETER(region='3d')
539	c
540	cIM ISCCP simulator v3.4
541	c
542	logical ok_hf
543	c
544	integer nid_hf, nid_hf3d
545	save ok_hf, nid_hf, nid_hf3d
546	c$OMP THREADPRIVATE(ok_hf, nid_hf, nid_hf3d)
547	c QUESTION : noms de variables ?
548
549	INTEGER longcles
550	PARAMETER ( longcles = 20 )
551	REAL clesphy0( longcles )
552	c
553	c Variables propres a la physique
554	INTEGER itap
555	SAVE itap ! compteur pour la physique
556	c$OMP THREADPRIVATE(itap)
557	c
558	real slp(klon) ! sea level pressure
559	c
560	REAL fevap(klon,nbsrf)
561	REAL fluxlat(klon,nbsrf)
562	c
563	REAL qsol(klon)
564	REAL,save :: solarlong0
565	c$OMP THREADPRIVATE(solarlong0)
566
567	c
568	c Parametres de l'Orographie a l'Echelle Sous-Maille (OESM):
569	c
570	cIM 141004 REAL zulow(klon),zvlow(klon),zustr(klon), zvstr(klon)
571	REAL zulow(klon),zvlow(klon)
572	c
573	INTEGER igwd,idx(klon),itest(klon)
574	c
575	REAL agesno(klon,nbsrf)
576	c
577	c REAL,allocatable,save :: run_off_lic_0(:)
578	cc$OMP THREADPRIVATE(run_off_lic_0)
579	cym SAVE run_off_lic_0
580	cKE43
581	c Variables liees a la convection de K. Emanuel (sb):
582	c
583	REAL bas, top ! cloud base and top levels
584	SAVE bas
585	SAVE top
586	c$OMP THREADPRIVATE(bas, top)
587
588	REAL wdn(klon), tdn(klon), qdn(klon)
589	c
590	c=================================================================================================
591	cCR04.12.07: on ajoute les nouvelles variables du nouveau schema de convection avec poches froides
592	c Variables liées à la poche froide (jyg)
593
594	REAL mip(klon,klev) ! mass flux shed by the adiab ascent at each level
595	REAL Vprecip(klon,klev+1) ! precipitation vertical profile
596	c
597	REAL wape_prescr, fip_prescr
598	INTEGER it_wape_prescr
599	SAVE wape_prescr, fip_prescr, it_wape_prescr
600	c$OMP THREADPRIVATE(wape_prescr, fip_prescr, it_wape_prescr)
601	c
602	c variables supplementaires de concvl
603	REAL Tconv(klon,klev)
604	REAL ment(klon,klev,klev),sij(klon,klev,klev)
605	REAL dd_t(klon,klev),dd_q(klon,klev)
606
607	real, save :: alp_bl_prescr=0.
608	real, save :: ale_bl_prescr=0.
609
610	real, save :: ale_max=100.
611	real, save :: alp_max=2.
612
613	c$OMP THREADPRIVATE(alp_bl_prescr,ale_bl_prescr)
614	c$OMP THREADPRIVATE(ale_max,alp_max)
615
616	real ale_wake(klon)
617	real alp_wake(klon)
618	cRC
619	c Variables liées à la poche froide (jyg et rr)
620	c Version diagnostique pour l'instant : pas de rétroaction sur la convection
621
622	REAL t_wake(klon,klev),q_wake(klon,klev) ! wake pour la convection
623
624	REAL wake_dth(klon,klev) ! wake : temp pot difference
625
626	REAL wake_d_deltat_gw(klon,klev)! wake : delta T tendency due to Gravity Wave (/s)
627	REAL wake_omgbdth(klon,klev) ! Wake : flux of Delta_Theta transported by LS omega
628	REAL wake_dp_omgb(klon,klev) ! Wake : vertical gradient of large scale omega
629	REAL wake_dtKE(klon,klev) ! Wake : differential heating (wake - unpertubed) CONV
630	REAL wake_dqKE(klon,klev) ! Wake : differential moistening (wake - unpertubed) CONV
631	REAL wake_dtPBL(klon,klev) ! Wake : differential heating (wake - unpertubed) PBL
632	REAL wake_dqPBL(klon,klev) ! Wake : differential moistening (wake - unpertubed) PBL
633	REAL wake_omg(klon,klev) ! Wake : velocity difference (wake - unpertubed)
634	REAL wake_ddeltat(klon,klev),wake_ddeltaq(klon,klev)
635	REAL wake_dp_deltomg(klon,klev) ! Wake : gradient vertical de wake_omg
636	REAL wake_spread(klon,klev) ! spreading term in wake_delt
637	c
638	cpourquoi y'a pas de save??
639	REAL wake_h(klon) ! Wake : hauteur de la poche froide
640	c
641	INTEGER wake_k(klon) ! Wake sommet
642	c
643	REAL t_undi(klon,klev) ! temperature moyenne dans la zone non perturbee
644	REAL q_undi(klon,klev) ! humidite moyenne dans la zone non perturbee
645	c
646	cjyg
647	ccc REAL wake_pe(klon) ! Wake potential energy - WAPE
648
649	REAL wake_gfl(klon) ! Gust Front Length
650	REAL wake_dens(klon)
651	c
652	c
653	REAL dt_dwn(klon,klev)
654	REAL dq_dwn(klon,klev)
655	REAL wdt_PBL(klon,klev)
656	REAL udt_PBL(klon,klev)
657	REAL wdq_PBL(klon,klev)
658	REAL udq_PBL(klon,klev)
659	REAL M_dwn(klon,klev)
660	REAL M_up(klon,klev)
661	REAL dt_a(klon,klev)
662	REAL dq_a(klon,klev)
663	REAL, SAVE :: alp_offset
664	c$OMP THREADPRIVATE(alp_offset)
665
666	c
667	cRR:fin declarations poches froides
668	c=======================================================================================================
669
670	REAL zw2(klon,klev+1)
671	REAL fraca(klon,klev+1)
672
673	c Variables locales pour la couche limite (al1):
674	c
675	cAl1 REAL pblh(klon) ! Hauteur de couche limite
676	cAl1 SAVE pblh
677	c34EK
678	c
679	c Variables locales:
680	c
681	REAL cdragh(klon) ! drag coefficient pour T and Q
682	REAL cdragm(klon) ! drag coefficient pour vent
683	cAA
684	cAA Pour phytrac
685	cAA
686	REAL coefh(klon,klev) ! coef d'echange pour phytrac, valable pour 2<=k<=klev
687	REAL u1(klon) ! vents dans la premiere couche U
688	REAL v1(klon) ! vents dans la premiere couche V
689
690	REAL zxffonte(klon), zxfqcalving(klon),zxfqfonte(klon)
691
692	c@$$ LOGICAL offline ! Controle du stockage ds "physique"
693	c@$$ PARAMETER (offline=.false.)
694	c@$$ INTEGER physid
695	REAL frac_impa(klon,klev) ! fractions d'aerosols lessivees (impaction)
696	REAL frac_nucl(klon,klev) ! idem (nucleation)
697	INTEGER :: iii
698	REAL :: calday
699
700	cIM cf FH pour Tiedtke 080604
701	REAL rain_tiedtke(klon),snow_tiedtke(klon)
702	c
703	cIM 050204 END
704	REAL evap(klon), devap(klon) ! evaporation et sa derivee
705	REAL sens(klon), dsens(klon) ! chaleur sensible et sa derivee
706
707	REAL bils(klon) ! bilan de chaleur au sol
708	REAL wfbilo(klon,nbsrf) ! bilan d'eau, pour chaque
709	C ! type de sous-surface et pondere par la fraction
710	REAL wfbils(klon,nbsrf) ! bilan de chaleur au sol, pour chaque
711	C ! type de sous-surface et pondere par la fraction
712	REAL slab_wfbils(klon) ! bilan de chaleur au sol pour le cas de slab, sur les points d'ocean
713
714	REAL fder(klon)
715	REAL ve(klon) ! integr. verticale du transport meri. de l'energie
716	REAL vq(klon) ! integr. verticale du transport meri. de l'eau
717	REAL ue(klon) ! integr. verticale du transport zonal de l'energie
718	REAL uq(klon) ! integr. verticale du transport zonal de l'eau
719	c
720	REAL frugs(klon,nbsrf)
721	REAL zxrugs(klon) ! longueur de rugosite
722	c
723	c Conditions aux limites
724	c
725	!
726	REAL :: day_since_equinox
727	! Date de l'equinoxe de printemps
728	INTEGER, parameter :: mth_eq=3, day_eq=21
729	REAL :: jD_eq
730
731	LOGICAL, parameter :: new_orbit = .true.
732
733	c
734	INTEGER lmt_pas
735	SAVE lmt_pas ! frequence de mise a jour
736	c$OMP THREADPRIVATE(lmt_pas)
737	real zmasse(klon, llm)
738	C (column-density of mass of air in a cell, in kg m-2)
739	real, parameter:: dobson_u = 2.1415e-05 ! Dobson unit, in kg m-2
740
741	cIM sorties
742	REAL un_jour
743	PARAMETER(un_jour=86400.)
744	c======================================================================
745	c
746	c Declaration des procedures appelees
747	c
748	EXTERNAL angle ! calculer angle zenithal du soleil
749	EXTERNAL alboc ! calculer l'albedo sur ocean
750	EXTERNAL ajsec ! ajustement sec
751	EXTERNAL conlmd ! convection (schema LMD)
752	cKE43
753	EXTERNAL conema3 ! convect4.3
754	EXTERNAL fisrtilp ! schema de condensation a grande echelle (pluie)
755	cAA
756	EXTERNAL fisrtilp_tr ! schema de condensation a grande echelle (pluie)
757	c ! stockage des coefficients necessaires au
758	c ! lessivage OFF-LINE et ON-LINE
759	EXTERNAL hgardfou ! verifier les temperatures
760	EXTERNAL nuage ! calculer les proprietes radiatives
761	CC EXTERNAL o3cm ! initialiser l'ozone
762	EXTERNAL orbite ! calculer l'orbite terrestre
763	EXTERNAL phyetat0 ! lire l'etat initial de la physique
764	EXTERNAL phyredem ! ecrire l'etat de redemarrage de la physique
765	EXTERNAL suphel ! initialiser certaines constantes
766	EXTERNAL transp ! transport total de l'eau et de l'energie
767	EXTERNAL ecribina ! ecrire le fichier binaire global
768	EXTERNAL ecribins ! ecrire le fichier binaire global
769	EXTERNAL ecrirega ! ecrire le fichier binaire regional
770	EXTERNAL ecriregs ! ecrire le fichier binaire regional
771	cIM
772	EXTERNAL haut2bas !variables de haut en bas
773	INTEGER lnblnk1
774	EXTERNAL lnblnk1 !enleve les blancs a la fin d'une variable de type
775	!caracter
776	EXTERNAL ini_undefSTD !initialise a 0 une variable a 1 niveau de pression
777	EXTERNAL undefSTD !somme les valeurs definies d'1 var a 1 niveau de pression
778	c EXTERNAL moy_undefSTD !moyenne d'1 var a 1 niveau de pression
779	c EXTERNAL moyglo_aire !moyenne globale d'1 var ponderee par l'aire de la maille (moyglo_pondaire)
780	c !par la masse/airetot (moyglo_pondaima) et la vraie masse (moyglo_pondmass)
781	c
782	c Variables locales
783	c
784	REAL rhcl(klon,klev) ! humiditi relative ciel clair
785	REAL dialiq(klon,klev) ! eau liquide nuageuse
786	REAL diafra(klon,klev) ! fraction nuageuse
787	REAL cldliq(klon,klev) ! eau liquide nuageuse
788	REAL cldfra(klon,klev) ! fraction nuageuse
789	REAL cldtau(klon,klev) ! epaisseur optique
790	REAL cldemi(klon,klev) ! emissivite infrarouge
791	c
792	CXXX PB
793	REAL fluxq(klon,klev, nbsrf) ! flux turbulent d'humidite
794	REAL fluxt(klon,klev, nbsrf) ! flux turbulent de chaleur
795	REAL fluxu(klon,klev, nbsrf) ! flux turbulent de vitesse u
796	REAL fluxv(klon,klev, nbsrf) ! flux turbulent de vitesse v
797	c
798	REAL zxfluxt(klon, klev)
799	REAL zxfluxq(klon, klev)
800	REAL zxfluxu(klon, klev)
801	REAL zxfluxv(klon, klev)
802	CXXX
803	c
804	REAL fsollw(klon, nbsrf) ! bilan flux IR pour chaque sous surface
805	REAL fsolsw(klon, nbsrf) ! flux solaire absorb. pour chaque sous surface
806	c Le rayonnement n'est pas calcule tous les pas, il faut donc
807	c sauvegarder les sorties du rayonnement
808	cym SAVE heat,cool,albpla,topsw,toplw,solsw,sollw,sollwdown
809	cym SAVE sollwdownclr, toplwdown, toplwdownclr
810	cym SAVE topsw0,toplw0,solsw0,sollw0, heat0, cool0
811	c
812	INTEGER itaprad
813	SAVE itaprad
814	c$OMP THREADPRIVATE(itaprad)
815	c
816	REAL conv_q(klon,klev) ! convergence de l'humidite (kg/kg/s)
817	REAL conv_t(klon,klev) ! convergence de la temperature(K/s)
818	c
819	REAL cldl(klon),cldm(klon),cldh(klon) !nuages bas, moyen et haut
820	REAL cldt(klon),cldq(klon) !nuage total, eau liquide integree
821	c
822	REAL zxtsol(klon), zxqsurf(klon), zxsnow(klon), zxfluxlat(klon)
823	REAL zxsnow_dummy(klon)
824	c
825	REAL dist, rmu0(klon), fract(klon)
826	REAL zdtime, zlongi
827	c
828	CHARACTER*2 str2
829	CHARACTER*2 iqn
830	c
831	REAL qcheck
832	REAL z_avant(klon), z_apres(klon), z_factor(klon)
833	LOGICAL zx_ajustq
834	c
835	REAL za, zb
836	REAL zx_t, zx_qs, zdelta, zcor, zfra, zlvdcp, zlsdcp
837	real zqsat(klon,klev)
838	INTEGER i, k, iq, ig, j, nsrf, ll, l, iiq, iff
839	REAL t_coup
840	PARAMETER (t_coup=234.0)
841	c
842	REAL zphi(klon,klev)
843	cym A voir plus tard !!
844	cym REAL zx_relief(iim,jjmp1)
845	cym REAL zx_aire(iim,jjmp1)
846	c
847	c Grandeurs de sorties
848	REAL s_pblh(klon), s_lcl(klon), s_capCL(klon)
849	REAL s_oliqCL(klon), s_cteiCL(klon), s_pblt(klon)
850	REAL s_therm(klon), s_trmb1(klon), s_trmb2(klon)
851	REAL s_trmb3(klon)
852	cKE43
853	c Variables locales pour la convection de K. Emanuel (sb):
854	c
855	REAL upwd(klon,klev) ! saturated updraft mass flux
856	REAL dnwd(klon,klev) ! saturated downdraft mass flux
857	REAL dnwd0(klon,klev) ! unsaturated downdraft mass flux
858	REAL tvp(klon,klev) ! virtual temp of lifted parcel
859	CHARACTER*40 capemaxcels !max(CAPE)
860
861	REAL rflag(klon) ! flag fonctionnement de convect
862	INTEGER iflagctrl(klon) ! flag fonctionnement de convect
863	c -- convect43:
864	INTEGER ntra ! nb traceurs pour convect4.3
865	REAL pori_con(klon) ! pressure at the origin level of lifted parcel
866	REAL plcl_con(klon),dtma_con(klon),dtlcl_con(klon)
867	REAL dtvpdt1(klon,klev), dtvpdq1(klon,klev)
868	REAL dplcldt(klon), dplcldr(klon)
869	c? . condm_con(klon,klev),conda_con(klon,klev),
870	c? . mr_con(klon,klev),ep_con(klon,klev)
871	c? . ,sadiab(klon,klev),wadiab(klon,klev)
872	c --
873	c34EK
874	c
875	c Variables du changement
876	c
877	c con: convection
878	c lsc: condensation a grande echelle (Large-Scale-Condensation)
879	c ajs: ajustement sec
880	c eva: evaporation de l'eau liquide nuageuse
881	c vdf: couche limite (Vertical DiFfusion)
882	REAL rneb(klon,klev)
883
884	! tendance nulles
885	REAL du0(klon,klev),dv0(klon,klev),dq0(klon,klev),dql0(klon,klev)
886
887	c
888	*********************************************************
889	* declarations
890
891	*********************************************************
892	cIM 081204 END
893	c
894	REAL pmfu(klon,klev), pmfd(klon,klev)
895	REAL pen_u(klon,klev), pen_d(klon,klev)
896	REAL pde_u(klon,klev), pde_d(klon,klev)
897	INTEGER kcbot(klon), kctop(klon), kdtop(klon)
898	REAL pmflxr(klon,klev+1), pmflxs(klon,klev+1)
899	REAL prfl(klon,klev+1), psfl(klon,klev+1)
900	c
901	REAL rain_lsc(klon)
902	REAL snow_lsc(klon)
903	c
904	REAL ratqss(klon,klev),ratqsc(klon,klev)
905	real ratqsbas,ratqshaut,tau_ratqs
906	save ratqsbas,ratqshaut,tau_ratqs
907	c$OMP THREADPRIVATE(ratqsbas,ratqshaut,tau_ratqs)
908	real zpt_conv(klon,klev)
909
910	c Parametres lies au nouveau schema de nuages (SB, PDF)
911	real fact_cldcon
912	real facttemps
913	logical ok_newmicro
914	save ok_newmicro
915	real ref_liq(klon,klev), ref_ice(klon,klev)
916	c$OMP THREADPRIVATE(ok_newmicro)
917	save fact_cldcon,facttemps
918	c$OMP THREADPRIVATE(fact_cldcon,facttemps)
919	real facteur
920
921	integer iflag_cldcon
922	save iflag_cldcon
923	c$OMP THREADPRIVATE(iflag_cldcon)
924	logical ptconv(klon,klev)
925	cIM cf. AM 081204 BEG
926	logical ptconvth(klon,klev)
927	cIM cf. AM 081204 END
928	c
929	c Variables liees a l'ecriture de la bande histoire physique
930	c
931	c======================================================================
932	c
933	cIM cf. AM 081204 BEG
934	c declarations pour sortir sur une sous-region
935	integer imin_ins,imax_ins,jmin_ins,jmax_ins
936	save imin_ins,imax_ins,jmin_ins,jmax_ins
937	c$OMP THREADPRIVATE(imin_ins,imax_ins,jmin_ins,jmax_ins)
938	c real lonmin_ins,lonmax_ins,latmin_ins
939	c s ,latmax_ins
940	c data lonmin_ins,lonmax_ins,latmin_ins
941	c s ,latmax_ins/
942	c valeurs initiales s -5.,20.,41.,55./
943	c s 100.,130.,-20.,20./
944	c s -180.,180.,-90.,90./
945	c======================================================================
946	cIM cf. AM 081204 END
947
948	c
949	integer itau_w ! pas de temps ecriture = itap + itau_phy
950	c
951	c
952	c Variables locales pour effectuer les appels en serie
953	c
954	REAL zx_rh(klon,klev)
955	cIM RH a 2m (la surface)
956	REAL rh2m(klon), qsat2m(klon)
957	REAL tpot(klon), tpote(klon)
958	REAL Lheat
959
960	INTEGER length
961	PARAMETER ( length = 100 )
962	REAL tabcntr0( length )
963	c
964	INTEGER ndex2d(iimjjmp1),ndex3d(iimjjmp1*klev)
965	cIM
966	INTEGER ndex2d1(iwmax)
967	c
968	cIM AMIP2 BEG
969	REAL moyglo, mountor
970	cIM 141004 BEG
971	REAL zustrdr(klon), zvstrdr(klon)
972	REAL zustrli(klon), zvstrli(klon)
973	REAL zustrph(klon), zvstrph(klon)
974	REAL zustrhi(klon), zvstrhi(klon)
975	REAL aam, torsfc
976	cIM 141004 END
977	cIM 190504 BEG
978	INTEGER ij, imp1jmp1
979	PARAMETER(imp1jmp1=(iim+1)*jjmp1)
980	cym A voir plus tard
981	REAL zx_tmp(imp1jmp1), airedyn(iim+1,jjmp1)
982	REAL padyn(iim+1,jjmp1,klev+1)
983	REAL dudyn(iim+1,jjmp1,klev)
984	REAL rlatdyn(iim+1,jjmp1)
985	cIM 190504 END
986	LOGICAL ok_msk
987	REAL msk(klon)
988	cIM
989	REAL airetot, pi
990	cym A voir plus tard
991	cym REAL zm_wo(jjmp1, klev)
992	cIM AMIP2 END
993	c
994	REAL zx_tmp_fi2d(klon) ! variable temporaire grille physique
995	REAL zx_tmp_fi3d(klon,klev) ! variable temporaire pour champs 3D
996	c#ifdef histNMC
997	cym A voir plus tard !!!!
998	cym REAL zx_tmp_NC(iim,jjmp1,nlevSTD)
999	REAL zx_tmp_fiNC(klon,nlevSTD)
1000	c#endif
1001	REAL(KIND=8) zx_tmp2_fi3d(klon,klev) ! variable temporaire pour champs 3D
1002	REAL zx_tmp_2d(iim,jjmp1), zx_tmp_3d(iim,jjmp1,klev)
1003	REAL zx_lon(iim,jjmp1), zx_lat(iim,jjmp1)
1004	cIM for NMC files
1005	REAL missing_val
1006	REAL, SAVE :: freq_moyNMC(nout)
1007	c$OMP THREADPRIVATE(freq_moyNMC)
1008	c
1009	INTEGER nid_day, nid_mth, nid_ins, nid_mthnmc, nid_daynmc
1010	INTEGER nid_hfnmc, nid_day_seri, nid_ctesGCM
1011	SAVE nid_day, nid_mth, nid_ins, nid_mthnmc, nid_daynmc
1012	SAVE nid_hfnmc, nid_day_seri, nid_ctesGCM
1013	c$OMP THREADPRIVATE(nid_day, nid_mth, nid_ins)
1014	c$OMP THREADPRIVATE(nid_mthnmc, nid_daynmc, nid_hfnmc)
1015	c$OMP THREADPRIVATE(nid_day_seri,nid_ctesGCM)
1016	c
1017	cIM 280405 BEG
1018	INTEGER nid_bilKPins, nid_bilKPave
1019	SAVE nid_bilKPins, nid_bilKPave
1020	c$OMP THREADPRIVATE(nid_bilKPins, nid_bilKPave)
1021	c
1022	REAL ve_lay(klon,klev) ! transport meri. de l'energie a chaque niveau vert.
1023	REAL vq_lay(klon,klev) ! transport meri. de l'eau a chaque niveau vert.
1024	REAL ue_lay(klon,klev) ! transport zonal de l'energie a chaque niveau vert.
1025	REAL uq_lay(klon,klev) ! transport zonal de l'eau a chaque niveau vert.
1026	c
1027	cIM 280405 END
1028	c
1029	INTEGER nhori, nvert, nvert1, nvert3
1030	REAL zsto, zsto1, zsto2
1031	REAL zstophy, zstorad, zstohf, zstoday, zstomth, zout
1032	REAL zcals(napisccp), zcalh(napisccp), zoutj(napisccp)
1033	REAL zout_isccp(napisccp)
1034	SAVE zcals, zcalh, zoutj, zout_isccp
1035	c$OMP THREADPRIVATE(zcals, zcalh, zoutj, zout_isccp)
1036
1037	real zjulian
1038	save zjulian
1039	c$OMP THREADPRIVATE(zjulian)
1040
1041	character*20 modname
1042	character*80 abort_message
1043	logical ok_sync
1044	real date0
1045	integer idayref
1046
1047	C essai writephys
1048	integer fid_day, fid_mth, fid_ins
1049	parameter (fid_ins = 1, fid_day = 2, fid_mth = 3)
1050	integer prof2d_on, prof3d_on, prof2d_av, prof3d_av
1051	parameter (prof2d_on = 1, prof3d_on = 2,
1052	. prof2d_av = 3, prof3d_av = 4)
1053	character*30 nom_fichier
1054	character*10 varname
1055	character*40 vartitle
1056	character*20 varunits
1057	C Variables liees au bilan d'energie et d'enthalpi
1058	REAL ztsol(klon)
1059	REAL h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
1060	$ , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
1061	SAVE h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
1062	$ , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
1063	c$OMP THREADPRIVATE(h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot,
1064	c$OMP+ h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot)
1065	REAL d_h_vcol, d_h_dair, d_qt, d_qw, d_ql, d_qs, d_ec
1066	REAL d_h_vcol_phy
1067	REAL fs_bound, fq_bound
1068	SAVE d_h_vcol_phy
1069	c$OMP THREADPRIVATE(d_h_vcol_phy)
1070	REAL zero_v(klon)
1071	CHARACTER*15 ztit
1072	INTEGER ip_ebil ! PRINT level for energy conserv. diag.
1073	SAVE ip_ebil
1074	DATA ip_ebil/0/
1075	c$OMP THREADPRIVATE(ip_ebil)
1076	INTEGER if_ebil ! level for energy conserv. dignostics
1077	SAVE if_ebil
1078	c$OMP THREADPRIVATE(if_ebil)
1079	c+jld ec_conser
1080	REAL ZRCPD
1081	c-jld ec_conser
1082	REAL t2m(klon,nbsrf) ! temperature a 2m
1083	REAL q2m(klon,nbsrf) ! humidite a 2m
1084
1085	cIM: t2m, q2m, u10m, v10m et t2mincels, t2maxcels
1086	REAL zt2m(klon), zq2m(klon) !temp., hum. 2m moyenne s/ 1 maille
1087	REAL zu10m(klon), zv10m(klon) !vents a 10m moyennes s/1 maille
1088	CHARACTER*40 t2mincels, t2maxcels !t2m min., t2m max
1089	CHARACTER*40 tinst, tave, typeval
1090	REAL cldtaupi(klon,klev) ! Cloud optical thickness for pre-industrial (pi) aerosols
1091
1092	REAL re(klon, klev) ! Cloud droplet effective radius
1093	REAL fl(klon, klev) ! denominator of re
1094
1095	REAL re_top(klon), fl_top(klon) ! CDR at top of liquid water clouds
1096
1097	! Aerosol optical properties
1098	CHARACTER*4, DIMENSION(naero_grp) :: rfname
1099	REAL, DIMENSION(klon) :: aerindex ! POLDER aerosol index
1100	REAL, DIMENSION(klon,klev) :: mass_solu_aero ! total mass concentration for all soluble aerosols[ug/m3]
1101	REAL, DIMENSION(klon,klev) :: mass_solu_aero_pi ! - " - (pre-industrial value)
1102	INTEGER :: naero ! aerosol species
1103
1104	! Parameters
1105	LOGICAL ok_ade, ok_aie ! Apply aerosol (in)direct effects or not
1106	REAL bl95_b0, bl95_b1 ! Parameter in Boucher and Lohmann (1995)
1107	SAVE ok_ade, ok_aie, bl95_b0, bl95_b1
1108	c$OMP THREADPRIVATE(ok_ade, ok_aie, bl95_b0, bl95_b1)
1109	LOGICAL, SAVE :: aerosol_couple ! true : calcul des aerosols dans INCA
1110	! false : lecture des aerosol dans un fichier
1111	c$OMP THREADPRIVATE(aerosol_couple)
1112	INTEGER, SAVE :: flag_aerosol
1113	c$OMP THREADPRIVATE(flag_aerosol)
1114	LOGICAL, SAVE :: new_aod
1115	c$OMP THREADPRIVATE(new_aod)
1116
1117	c
1118	c Declaration des constantes et des fonctions thermodynamiques
1119	c
1120	LOGICAL,SAVE :: first=.true.
1121	c$OMP THREADPRIVATE(first)
1122
1123	integer iunit
1124
1125	integer, save:: read_climoz ! read ozone climatology
1126	C Allowed values are 0, 1 and 2
1127	C 0: do not read an ozone climatology
1128	C 1: read a single ozone climatology that will be used day and night
1129	C 2: read two ozone climatologies, the average day and night
1130	C climatology and the daylight climatology
1131
1132	integer, save:: ncid_climoz ! NetCDF file containing ozone climatologies
1133
1134	real, pointer, save:: press_climoz(:)
1135	! edges of pressure intervals for ozone climatologies, in Pa, in strictly
1136	! ascending order
1137	cIM for NMC files en parallele
1138	c$OMP THREADPRIVATE(read_climoz, ncid_climoz, press_climoz)
1139
1140	integer, save:: co3i = 0
1141	! time index in NetCDF file of current ozone fields
1142	c$OMP THREADPRIVATE(co3i)
1143
1144	integer ro3i
1145	! required time index in NetCDF file for the ozone fields, between 1
1146	! and 360
1147
1148	#include "YOMCST.h"
1149	#include "YOETHF.h"
1150	#include "FCTTRE.h"
1151	cIM 100106 BEG : pouvoir sortir les ctes de la physique
1152	#include "conema3.h"
1153	#include "fisrtilp.h"
1154	#include "nuage.h"
1155	#include "compbl.h"
1156	cIM 100106 END : pouvoir sortir les ctes de la physique
1157	c
1158	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1159	c Declarations pour Simulateur COSP
1160	c============================================================
1161	real :: mr_ozone(klon,klev)
1162	cIM for NMC files
1163	missing_val=nf90_fill_real
1164	c======================================================================
1165	! Gestion calendrier : mise a jour du module phys_cal_mod
1166	!
1167	CALL phys_cal_update(jD_cur,jH_cur)
1168
1169	c======================================================================
1170	! Ecriture eventuelle d'un profil verticale en entree de la physique.
1171	! Utilise notamment en 1D mais peut etre active egalement en 3D
1172	! en imposant la valeur de igout.
1173	c======================================================================
1174
1175	if (prt_level.ge.1) then
1176	igout=klon/2+1/klon
1177	write(lunout,*) 'DEBUT DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
1178	write(lunout,*)
1179	s 'nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur,pdtphys'
1180	write(lunout,*)
1181	s nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur,pdtphys
1182
1183	write(lunout,*) 'paprs, play, phi, u, v, t'
1184	do k=1,klev
1185	write(lunout,*) paprs(igout,k),pplay(igout,k),pphi(igout,k),
1186	s u(igout,k),v(igout,k),t(igout,k)
1187	enddo
1188	write(lunout,*) 'ovap (g/kg), oliq (g/kg)'
1189	do k=1,klev
1190	write(lunout,) qx(igout,k,1)1000,qx(igout,k,2)*1000.
1191	enddo
1192	endif
1193
1194	c======================================================================
1195
1196	cym => necessaire pour iflag_con != 2
1197	pmfd(:,:) = 0.
1198	pen_u(:,:) = 0.
1199	pen_d(:,:) = 0.
1200	pde_d(:,:) = 0.
1201	pde_u(:,:) = 0.
1202	aam=0.
1203
1204	torsfc=0.
1205	forall (k=1: llm) zmasse(:, k) = (paprs(:, k)-paprs(:, k+1)) / rg
1206
1207	if (first) then
1208
1209	cCR:nvelles variables convection/poches froides
1210
1211	print*, '================================================='
1212	print*, 'Allocation des variables locales et sauvegardees'
1213	call phys_local_var_init
1214	c
1215	pasphys=pdtphys
1216	c appel a la lecture du run.def physique
1217	call conf_phys(ok_journe, ok_mensuel,
1218	. ok_instan, ok_hf,
1219	. ok_LES,
1220	. solarlong0,seuil_inversion,
1221	. fact_cldcon, facttemps,ok_newmicro,iflag_radia,
1222	. iflag_cldcon,iflag_ratqs,ratqsbas,ratqshaut,tau_ratqs,
1223	. ok_ade, ok_aie, aerosol_couple,
1224	. flag_aerosol, new_aod,
1225	. bl95_b0, bl95_b1,
1226	. iflag_thermals,nsplit_thermals,tau_thermals,
1227	. iflag_thermals_ed,iflag_thermals_optflux,
1228	c nv flags pour la convection et les poches froides
1229	. iflag_coupl,iflag_clos,iflag_wake, read_climoz,
1230	& alp_offset)
1231	call phys_state_var_init(read_climoz)
1232	call phys_output_var_init
1233
1234	print*, '================================================='
1235	cIM for NMC files
1236	cIM freq_moyNMC = frequences auxquelles on moyenne les champs accumules
1237	cIM sur les niveaux de pression standard du NMC
1238	DO n=1, nout
1239	freq_moyNMC(n)=freq_outNMC(n)/freq_calNMC(n)
1240	ENDDO
1241	c
1242	cIM beg
1243	dnwd0=0.0
1244	ftd=0.0
1245	fqd=0.0
1246	cin=0.
1247	cym Attention pbase pas initialise dans concvl !!!!
1248	pbase=0
1249	paire_ter(:)=0.
1250	cIM 180608
1251	c pmflxr=0.
1252	c pmflxs=0.
1253	itau_con=0
1254	first=.false.
1255
1256	endif ! first
1257
1258	modname = 'physiq'
1259	cIM
1260	IF (ip_ebil_phy.ge.1) THEN
1261	DO i=1,klon
1262	zero_v(i)=0.
1263	END DO
1264	END IF
1265	ok_sync=.TRUE.
1266
1267	IF (debut) THEN
1268	CALL suphel ! initialiser constantes et parametres phys.
1269	ENDIF
1270
1271	if(prt_level.ge.1) print*,'CONVERGENCE PHYSIQUE THERM 1 '
1272
1273
1274	c======================================================================
1275	! Gestion calendrier : mise a jour du module phys_cal_mod
1276	!
1277	c CALL phys_cal_update(jD_cur,jH_cur)
1278
1279	c
1280	c Si c'est le debut, il faut initialiser plusieurs choses
1281	c ********
1282	c
1283	IF (debut) THEN
1284	!rv
1285	cCRinitialisation de wght_th et lalim_conv pour la definition de la couche alimentation
1286	cde la convection a partir des caracteristiques du thermique
1287	wght_th(:,:)=1.
1288	lalim_conv(:)=1
1289	cRC
1290	u10m(:,:)=0.
1291	v10m(:,:)=0.
1292	rain_con(:)=0.
1293	snow_con(:)=0.
1294	topswai(:)=0.
1295	topswad(:)=0.
1296	solswai(:)=0.
1297	solswad(:)=0.
1298
1299	lambda_th(:,:)=0.
1300	wmax_th(:)=0.
1301	tau_overturning_th(:)=0.
1302
1303	IF (config_inca /= 'none') THEN
1304	! jg : initialisation jusqu'au ces variables sont dans restart
1305	ccm(:,:,:) = 0.
1306	tau_aero(:,:,:,:) = 0.
1307	piz_aero(:,:,:,:) = 0.
1308	cg_aero(:,:,:,:) = 0.
1309	END IF
1310
1311	rnebcon0(:,:) = 0.0
1312	clwcon0(:,:) = 0.0
1313	rnebcon(:,:) = 0.0
1314	clwcon(:,:) = 0.0
1315
1316	cIM
1317	IF (ip_ebil_phy.ge.1) d_h_vcol_phy=0.
1318	c
1319	print*,'iflag_coupl,iflag_clos,iflag_wake',
1320	. iflag_coupl,iflag_clos,iflag_wake
1321	print*,'CYCLE_DIURNE', cycle_diurne
1322	c
1323	IF (iflag_con.EQ.2.AND.iflag_cldcon.GT.-1) THEN
1324	abort_message = 'Tiedtke needs iflag_cldcon=-2 or -1'
1325	CALL abort_gcm (modname,abort_message,1)
1326	ENDIF
1327	c
1328	IF(ok_isccp.AND.iflag_con.LE.2) THEN
1329	abort_message = 'ISCCP-like outputs may be available for KE
1330	.(iflag_con >= 3); for Tiedtke (iflag_con=-2) put ok_isccp=n'
1331	CALL abort_gcm (modname,abort_message,1)
1332	ENDIF
1333	c
1334	c Initialiser les compteurs:
1335	c
1336	itap = 0
1337	itaprad = 0
1338
1339	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1340	!! Un petit travail à faire ici.
1341	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1342
1343	if (iflag_pbl>1) then
1344	PRINT*, "Using method MELLOR&YAMADA"
1345	endif
1346
1347	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1348	! FH 2008/05/02 changement lie a la lecture de nbapp_rad dans phylmd plutot que
1349	! dyn3d
1350	! Attention : la version precedente n'etait pas tres propre.
1351	! Il se peut qu'il faille prendre une valeur differente de nbapp_rad
1352	! pour obtenir le meme resultat.
1353	dtime=pdtphys
1354	radpas = NINT( 86400./dtime/nbapp_rad)
1355	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1356
1357	CALL phyetat0 ("startphy.nc",clesphy0,tabcntr0)
1358	cIM begin
1359	print*,'physiq: clwcon rnebcon ratqs',clwcon(1,1),rnebcon(1,1)
1360	$,ratqs(1,1)
1361	cIM end
1362
1363
1364
1365	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1366	c
1367	C on remet le calendrier a zero
1368	c
1369	IF (raz_date .eq. 1) THEN
1370	itau_phy = 0
1371	ENDIF
1372
1373	cIM cf. AM 081204 BEG
1374	PRINT*,'cycle_diurne3 =',cycle_diurne
1375	cIM cf. AM 081204 END
1376	c
1377	CALL printflag( tabcntr0,radpas,ok_journe,
1378	, ok_instan, ok_region )
1379	c
1380	IF (ABS(dtime-pdtphys).GT.0.001) THEN
1381	WRITE(lunout,*) 'Pas physique n est pas correct',dtime,
1382	. pdtphys
1383	abort_message='Pas physique n est pas correct '
1384	! call abort_gcm(modname,abort_message,1)
1385	dtime=pdtphys
1386	ENDIF
1387	IF (nlon .NE. klon) THEN
1388	WRITE(lunout,*)'nlon et klon ne sont pas coherents', nlon,
1389	. klon
1390	abort_message='nlon et klon ne sont pas coherents'
1391	call abort_gcm(modname,abort_message,1)
1392	ENDIF
1393	IF (nlev .NE. klev) THEN
1394	WRITE(lunout,*)'nlev et klev ne sont pas coherents', nlev,
1395	. klev
1396	abort_message='nlev et klev ne sont pas coherents'
1397	call abort_gcm(modname,abort_message,1)
1398	ENDIF
1399	c
1400	IF (dtime*REAL(radpas).GT.21600..AND.cycle_diurne) THEN
1401	WRITE(lunout,*)'Nbre d appels au rayonnement insuffisant'
1402	WRITE(lunout,*)"Au minimum 4 appels par jour si cycle diurne"
1403	abort_message='Nbre d appels au rayonnement insuffisant'
1404	call abort_gcm(modname,abort_message,1)
1405	ENDIF
1406	WRITE(lunout,*)"Clef pour la convection, iflag_con=", iflag_con
1407	WRITE(lunout,*)"Clef pour le driver de la convection, ok_cvl=",
1408	. ok_cvl
1409	c
1410	cKE43
1411	c Initialisation pour la convection de K.E. (sb):
1412	IF (iflag_con.GE.3) THEN
1413
1414	WRITE(lunout,)"** Convection de Kerry Emanuel 4.3 "
1415	WRITE(lunout,*)
1416	. "On va utiliser le melange convectif des traceurs qui"
1417	WRITE(lunout,*)"est calcule dans convect4.3"
1418	WRITE(lunout,*)" !!! penser aux logical flags de phytrac"
1419
1420	DO i = 1, klon
1421	ema_cbmf(i) = 0.
1422	ema_pcb(i) = 0.
1423	ema_pct(i) = 0.
1424	ema_workcbmf(i) = 0.
1425	ENDDO
1426	cIM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>BEG
1427	DO i = 1, klon
1428	ibas_con(i) = 1
1429	itop_con(i) = 1
1430	ENDDO
1431	cIM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>END
1432	c===============================================================================
1433	cCR:04.12.07: initialisations poches froides
1434	c Controle de ALE et ALP pour la fermeture convective (jyg)
1435	if (iflag_wake.eq.1) then
1436	CALL ini_wake(0.,0.,it_wape_prescr,wape_prescr,fip_prescr
1437	s ,alp_bl_prescr, ale_bl_prescr)
1438	c 11/09/06 rajout initialisation ALE et ALP du wake et PBL(YU)
1439	c print*,'apres ini_wake iflag_cldcon=', iflag_cldcon
1440	endif
1441
1442	do i = 1,klon
1443	Ale_bl(i)=0.
1444	Alp_bl(i)=0.
1445	enddo
1446
1447	c================================================================================
1448
1449	ENDIF
1450
1451	DO i=1,klon
1452	rugoro(i) = f_rugoro * MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
1453	ENDDO
1454
1455	c34EK
1456	IF (ok_orodr) THEN
1457
1458	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1459	! FH sans doute a enlever de finitivement ou, si on le garde, l'activer
1460	! justement quand ok_orodr = false.
1461	! ce rugoro est utilise par la couche limite et fait double emploi
1462	! avec les paramétrisations spécifiques de Francois Lott.
1463	! DO i=1,klon
1464	! rugoro(i) = MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
1465	! ENDDO
1466	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1467	IF (ok_strato) THEN
1468	CALL SUGWD_strato(klon,klev,paprs,pplay)
1469	ELSE
1470	CALL SUGWD(klon,klev,paprs,pplay)
1471	ENDIF
1472
1473	DO i=1,klon
1474	zuthe(i)=0.
1475	zvthe(i)=0.
1476	if(zstd(i).gt.10.)then
1477	zuthe(i)=(1.-zgam(i))*cos(zthe(i))
1478	zvthe(i)=(1.-zgam(i))*sin(zthe(i))
1479	endif
1480	ENDDO
1481	ENDIF
1482	c
1483	c
1484	lmt_pas = NINT(86400./dtime * 1.0) ! tous les jours
1485	WRITE(lunout,*)'La frequence de lecture surface est de ',
1486	. lmt_pas
1487	c
1488	cIM 030306 END
1489
1490	capemaxcels = 't_max(X)'
1491	t2mincels = 't_min(X)'
1492	t2maxcels = 't_max(X)'
1493	tinst = 'inst(X)'
1494	tave = 'ave(X)'
1495	cIM cf. AM 081204 BEG
1496	write(lunout,*)'AVANT HIST IFLAG_CON=',iflag_con
1497	cIM cf. AM 081204 END
1498	c
1499	c=============================================================
1500	c Initialisation des sorties
1501	c=============================================================
1502
1503	#ifdef CPP_IOIPSL
1504
1505	c$OMP MASTER
1506	call phys_output_open(jjmp1,nlevSTD,clevSTD,nbteta,
1507	& ctetaSTD,dtime,ok_veget,
1508	& type_ocean,iflag_pbl,ok_mensuel,ok_journe,
1509	& ok_hf,ok_instan,ok_LES,ok_ade,ok_aie,
1510	& read_climoz, new_aod, aerosol_couple,
1511	& o_trac)
1512	c$OMP END MASTER
1513	c$OMP BARRIER
1514
1515	#ifdef histISCCP
1516	#include "ini_histISCCP.h"
1517	#endif
1518
1519	#ifdef histNMC
1520	#include "ini_histhfNMC.h"
1521	#include "ini_histdayNMC.h"
1522	#include "ini_histmthNMC.h"
1523	#endif
1524
1525	#include "ini_histday_seri.h"
1526
1527	#include "ini_paramLMDZ_phy.h"
1528
1529	#endif
1530
1531	cIM 250308bad guide ecrit_hf2mth = 30*1/ecrit_hf
1532	ecrit_hf2mth = ecrit_mth/ecrit_hf
1533
1534	ecrit_hf = ecrit_hf * un_jour
1535	cIM
1536	IF(ecrit_day.LE.1.) THEN
1537	ecrit_day = ecrit_day * un_jour !en secondes
1538	ENDIF
1539	cIM
1540	ecrit_mth = ecrit_mth * un_jour
1541	ecrit_ins = ecrit_ins * un_jour
1542	ecrit_reg = ecrit_reg * un_jour
1543	ecrit_tra = ecrit_tra * un_jour
1544	ecrit_ISCCP = ecrit_ISCCP * un_jour
1545	ecrit_LES = ecrit_LES * un_jour
1546	c
1547	PRINT*,'physiq ecrit_ hf day mth reg tra ISCCP hf2mth',
1548	. ecrit_hf,ecrit_day,ecrit_mth,ecrit_reg,ecrit_tra,ecrit_ISCCP,
1549	. ecrit_hf2mth
1550	cIM 030306 END
1551
1552
1553	cXXXPB Positionner date0 pour initialisation de ORCHIDEE
1554	date0 = jD_ref
1555	WRITE(,) 'physiq date0 : ',date0
1556	c
1557	c
1558	c
1559	c Prescrire l'ozone dans l'atmosphere
1560	c
1561	c
1562	cc DO i = 1, klon
1563	cc DO k = 1, klev
1564	cc CALL o3cm (paprs(i,k)/100.,paprs(i,k+1)/100., wo(i,k),20)
1565	cc ENDDO
1566	cc ENDDO
1567	c
1568	IF (config_inca /= 'none') THEN
1569	#ifdef INCA
1570	CALL VTe(VTphysiq)
1571	CALL VTb(VTinca)
1572	! iii = MOD(NINT(xjour),360)
1573	! calday = REAL(iii) + jH_cur
1574	calday = REAL(days_elapsed) + jH_cur
1575	WRITE(lunout,*) 'initial time chemini', days_elapsed, calday
1576
1577	CALL chemini(
1578	$ rg,
1579	$ ra,
1580	$ airephy,
1581	$ rlat,
1582	$ rlon,
1583	$ presnivs,
1584	$ calday,
1585	$ klon,
1586	$ nqtot,
1587	$ pdtphys,
1588	$ annee_ref,
1589	$ day_ref,
1590	$ itau_phy)
1591
1592	CALL VTe(VTinca)
1593	CALL VTb(VTphysiq)
1594	#endif
1595	END IF
1596	c
1597	c
1598	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1599	! Nouvelle initialisation pour le rayonnement RRTM
1600	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1601
1602	call iniradia(klon,klev,paprs(1,1:klev+1))
1603
1604	C$omp single
1605	if (read_climoz >= 1) then
1606	call open_climoz(ncid_climoz, press_climoz)
1607	END IF
1608	C$omp end single
1609	ENDIF
1610	!
1611	! ************** Fin de IF ( debut ) *************
1612	!
1613	!
1614	! Incrementer le compteur de la physique
1615	!
1616	itap = itap + 1
1617	!
1618	! Update fraction of the sub-surfaces (pctsrf) and
1619	! initialize, where a new fraction has appeared, all variables depending
1620	! on the surface fraction.
1621	!
1622	CALL change_srf_frac(itap, dtime, days_elapsed+1,
1623	* pctsrf, falb1, falb2, ftsol, u10m, v10m, pbl_tke)
1624
1625	! Tendances bidons pour les processus qui n'affectent pas certaines
1626	! variables.
1627	du0(:,:)=0.
1628	dv0(:,:)=0.
1629	dq0(:,:)=0.
1630	dql0(:,:)=0.
1631	c
1632	c Mettre a zero des variables de sortie (pour securite)
1633	c
1634	DO i = 1, klon
1635	d_ps(i) = 0.0
1636	ENDDO
1637	DO k = 1, klev
1638	DO i = 1, klon
1639	d_t(i,k) = 0.0
1640	d_u(i,k) = 0.0
1641	d_v(i,k) = 0.0
1642	ENDDO
1643	ENDDO
1644	DO iq = 1, nqtot
1645	DO k = 1, klev
1646	DO i = 1, klon
1647	d_qx(i,k,iq) = 0.0
1648	ENDDO
1649	ENDDO
1650	ENDDO
1651	da(:,:)=0.
1652	mp(:,:)=0.
1653	phi(:,:,:)=0.
1654	c
1655	c Ne pas affecter les valeurs entrees de u, v, h, et q
1656	c
1657	DO k = 1, klev
1658	DO i = 1, klon
1659	t_seri(i,k) = t(i,k)
1660	u_seri(i,k) = u(i,k)
1661	v_seri(i,k) = v(i,k)
1662	q_seri(i,k) = qx(i,k,ivap)
1663	ql_seri(i,k) = qx(i,k,iliq)
1664	qs_seri(i,k) = 0.
1665	ENDDO
1666	ENDDO
1667	IF (nqtot.GE.3) THEN
1668	DO iq = 3, nqtot
1669	DO k = 1, klev
1670	DO i = 1, klon
1671	tr_seri(i,k,iq-2) = qx(i,k,iq)
1672	ENDDO
1673	ENDDO
1674	ENDDO
1675	ELSE
1676	DO k = 1, klev
1677	DO i = 1, klon
1678	tr_seri(i,k,1) = 0.0
1679	ENDDO
1680	ENDDO
1681	ENDIF
1682	C
1683	DO i = 1, klon
1684	ztsol(i) = 0.
1685	ENDDO
1686	DO nsrf = 1, nbsrf
1687	DO i = 1, klon
1688	ztsol(i) = ztsol(i) + ftsol(i,nsrf)*pctsrf(i,nsrf)
1689	ENDDO
1690	ENDDO
1691	cIM
1692	IF (ip_ebil_phy.ge.1) THEN
1693	ztit='after dynamic'
1694	CALL diagetpq(airephy,ztit,ip_ebil_phy,1,1,dtime
1695	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1696	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1697	C Comme les tendances de la physique sont ajoute dans la dynamique,
1698	C on devrait avoir que la variation d'entalpie par la dynamique
1699	C est egale a la variation de la physique au pas de temps precedent.
1700	C Donc la somme de ces 2 variations devrait etre nulle.
1701	call diagphy(airephy,ztit,ip_ebil_phy
1702	e , zero_v, zero_v, zero_v, zero_v, zero_v
1703	e , zero_v, zero_v, zero_v, ztsol
1704	e , d_h_vcol+d_h_vcol_phy, d_qt, 0.
1705	s , fs_bound, fq_bound )
1706	END IF
1707
1708	c Diagnostiquer la tendance dynamique
1709	c
1710	IF (ancien_ok) THEN
1711	DO k = 1, klev
1712	DO i = 1, klon
1713	d_u_dyn(i,k) = (u_seri(i,k)-u_ancien(i,k))/dtime
1714	d_v_dyn(i,k) = (v_seri(i,k)-v_ancien(i,k))/dtime
1715	d_t_dyn(i,k) = (t_seri(i,k)-t_ancien(i,k))/dtime
1716	d_q_dyn(i,k) = (q_seri(i,k)-q_ancien(i,k))/dtime
1717	ENDDO
1718	ENDDO
1719	ELSE
1720	DO k = 1, klev
1721	DO i = 1, klon
1722	d_u_dyn(i,k) = 0.0
1723	d_v_dyn(i,k) = 0.0
1724	d_t_dyn(i,k) = 0.0
1725	d_q_dyn(i,k) = 0.0
1726	ENDDO
1727	ENDDO
1728	ancien_ok = .TRUE.
1729	ENDIF
1730	c
1731	c Ajouter le geopotentiel du sol:
1732	c
1733	DO k = 1, klev
1734	DO i = 1, klon
1735	zphi(i,k) = pphi(i,k) + pphis(i)
1736	ENDDO
1737	ENDDO
1738	c
1739	c Verifier les temperatures
1740	c
1741	cIM BEG
1742	IF (check) THEN
1743	amn=MIN(ftsol(1,is_ter),1000.)
1744	amx=MAX(ftsol(1,is_ter),-1000.)
1745	DO i=2, klon
1746	amn=MIN(ftsol(i,is_ter),amn)
1747	amx=MAX(ftsol(i,is_ter),amx)
1748	ENDDO
1749	c
1750	PRINT*,' debut avant hgardfou min max ftsol',itap,amn,amx
1751	ENDIF !(check) THEN
1752	cIM END
1753	c
1754	CALL hgardfou(t_seri,ftsol,'debutphy')
1755	c
1756	cIM BEG
1757	IF (check) THEN
1758	amn=MIN(ftsol(1,is_ter),1000.)
1759	amx=MAX(ftsol(1,is_ter),-1000.)
1760	DO i=2, klon
1761	amn=MIN(ftsol(i,is_ter),amn)
1762	amx=MAX(ftsol(i,is_ter),amx)
1763	ENDDO
1764	c
1765	PRINT*,' debut apres hgardfou min max ftsol',itap,amn,amx
1766	ENDIF !(check) THEN
1767	cIM END
1768	c
1769	c Mettre en action les conditions aux limites (albedo, sst, etc.).
1770	c Prescrire l'ozone et calculer l'albedo sur l'ocean.
1771	c
1772	if (read_climoz >= 1) then
1773	C Ozone from a file
1774	! Update required ozone index:
1775	ro3i = int((days_elapsed + jh_cur - jh_1jan)
1776	$ / ioget_year_len(year_cur) * 360.) + 1
1777	if (ro3i == 361) ro3i = 360
1778	C (This should never occur, except perhaps because of roundup
1779	C error. See documentation.)
1780	if (ro3i /= co3i) then
1781	C Update ozone field:
1782	if (read_climoz == 1) then
1783	call regr_pr_av(ncid_climoz, (/"tro3"/), julien=ro3i,
1784	$ press_in_edg=press_climoz, paprs=paprs, v3=wo)
1785	else
1786	C read_climoz == 2
1787	call regr_pr_av(ncid_climoz,
1788	$ (/"tro3 ", "tro3_daylight"/),
1789	$ julien=ro3i, press_in_edg=press_climoz, paprs=paprs,
1790	$ v3=wo)
1791	end if
1792	! Convert from mole fraction of ozone to column density of ozone in a
1793	! cell, in kDU:
1794	forall (l = 1: read_climoz) wo(:, :, l) = wo(:, :, l)
1795	$ * rmo3 / rmd * zmasse / dobson_u / 1e3
1796	C (By regridding ozone values for LMDZ only once every 360th of
1797	C year, we have already neglected the variation of pressure in one
1798	C 360th of year. So do not recompute "wo" at each time step even if
1799	C "zmasse" changes a little.)
1800	co3i = ro3i
1801	end if
1802	elseif (MOD(itap-1,lmt_pas) == 0) THEN
1803	C Once per day, update ozone from Royer:
1804	wo(:, :, 1) = ozonecm(rlat, paprs, rjour=real(days_elapsed+1))
1805	ENDIF
1806	c
1807	c Re-evaporer l'eau liquide nuageuse
1808	c
1809	DO k = 1, klev ! re-evaporation de l'eau liquide nuageuse
1810	DO i = 1, klon
1811	zlvdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1812	c zlsdcp=RLSTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1813	zlsdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1814	zdelta = MAX(0.,SIGN(1.,RTT-t_seri(i,k)))
1815	zb = MAX(0.0,ql_seri(i,k))
1816	za = - MAX(0.0,ql_seri(i,k))
1817	. * (zlvdcp(1.-zdelta)+zlsdcpzdelta)
1818	t_seri(i,k) = t_seri(i,k) + za
1819	q_seri(i,k) = q_seri(i,k) + zb
1820	ql_seri(i,k) = 0.0
1821	d_t_eva(i,k) = za
1822	d_q_eva(i,k) = zb
1823	ENDDO
1824	ENDDO
1825	cIM
1826	IF (ip_ebil_phy.ge.2) THEN
1827	ztit='after reevap'
1828	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,1,dtime
1829	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1830	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1831	call diagphy(airephy,ztit,ip_ebil_phy
1832	e , zero_v, zero_v, zero_v, zero_v, zero_v
1833	e , zero_v, zero_v, zero_v, ztsol
1834	e , d_h_vcol, d_qt, d_ec
1835	s , fs_bound, fq_bound )
1836	C
1837	END IF
1838
1839	c
1840	c=========================================================================
1841	! Calculs de l'orbite.
1842	! Necessaires pour le rayonnement et la surface (calcul de l'albedo).
1843	! doit donc etre placé avant radlwsw et pbl_surface
1844
1845	! calcul selon la routine utilisee pour les planetes
1846	if (new_orbit) then
1847	call ymds2ju(year_cur, mth_eq, day_eq,0., jD_eq)
1848	day_since_equinox = (jD_cur + jH_cur) - jD_eq
1849	! day_since_equinox = (jD_cur) - jD_eq
1850	call solarlong(day_since_equinox, zlongi, dist)
1851	else
1852	! calcul selon la routine utilisee pour l'AR4
1853	! choix entre calcul de la longitude solaire vraie ou valeur fixee a
1854	! solarlong0
1855	if (solarlong0<-999.) then
1856	CALL orbite(REAL(days_elapsed+1),zlongi,dist)
1857	else
1858	zlongi=solarlong0 ! longitude solaire vraie
1859	dist=1. ! distance au soleil / moyenne
1860	endif
1861	endif
1862	if(prt_level.ge.1) &
1863	& write(lunout,*)'Longitude solaire ',zlongi,solarlong0,dist
1864
1865	! Avec ou sans cycle diurne
1866	IF (cycle_diurne) THEN
1867	zdtime=dtime*REAL(radpas) ! pas de temps du rayonnement (s)
1868	CALL zenang(zlongi,jH_cur,zdtime,rlat,rlon,rmu0,fract)
1869	ELSE
1870	CALL angle(zlongi, rlat, fract, rmu0)
1871	ENDIF
1872
1873	if (mydebug) then
1874	call writefield_phy('u_seri',u_seri,llm)
1875	call writefield_phy('v_seri',v_seri,llm)
1876	call writefield_phy('t_seri',t_seri,llm)
1877	call writefield_phy('q_seri',q_seri,llm)
1878	endif
1879
1880	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1881	c Appel au pbl_surface : Planetary Boudary Layer et Surface
1882	c Cela implique tous les interactions des sous-surfaces et la partie diffusion
1883	c turbulent du couche limit.
1884	c
1885	c Certains varibales de sorties de pbl_surface sont utiliser que pour
1886	c ecriture des fihiers hist_XXXX.nc, ces sont :
1887	c qsol, zq2m, s_pblh, s_lcl,
1888	c s_capCL, s_oliqCL, s_cteiCL,s_pblT,
1889	c s_therm, s_trmb1, s_trmb2, s_trmb3,
1890	c zxrugs, zu10m, zv10m, fder,
1891	c zxqsurf, rh2m, zxfluxu, zxfluxv,
1892	c frugs, agesno, fsollw, fsolsw,
1893	c d_ts, fevap, fluxlat, t2m,
1894	c wfbils, wfbilo, fluxt, fluxu, fluxv,
1895	c
1896	c Certains ne sont pas utiliser du tout :
1897	c dsens, devap, zxsnow, zxfluxt, zxfluxq, q2m, fluxq
1898	c
1899
1900	CALL pbl_surface(
1901	e dtime, date0, itap, days_elapsed+1,
1902	e debut, lafin,
1903	e rlon, rlat, rugoro, rmu0,
1904	e rain_fall, snow_fall, solsw, sollw,
1905	e t_seri, q_seri, u_seri, v_seri,
1906	e pplay, paprs, pctsrf,
1907	+ ftsol, falb1, falb2, u10m, v10m,
1908	s sollwdown, cdragh, cdragm, u1, v1,
1909	s albsol1, albsol2, sens, evap,
1910	s zxtsol, zxfluxlat, zt2m, qsat2m,
1911	s d_t_vdf, d_q_vdf, d_u_vdf, d_v_vdf,
1912	s coefh, slab_wfbils,
1913	d qsol, zq2m, s_pblh, s_lcl,
1914	d s_capCL, s_oliqCL, s_cteiCL,s_pblT,
1915	d s_therm, s_trmb1, s_trmb2, s_trmb3,
1916	d zxrugs, zu10m, zv10m, fder,
1917	d zxqsurf, rh2m, zxfluxu, zxfluxv,
1918	d frugs, agesno, fsollw, fsolsw,
1919	d d_ts, fevap, fluxlat, t2m,
1920	d wfbils, wfbilo, fluxt, fluxu, fluxv,
1921	- dsens, devap, zxsnow,
1922	- zxfluxt, zxfluxq, q2m, fluxq, pbl_tke )
1923
1924
1925	!-----------------------------------------------------------------------------------------
1926	! ajout des tendances de la diffusion turbulente
1927	CALL add_phys_tend(d_u_vdf,d_v_vdf,d_t_vdf,d_q_vdf,dql0,'vdf')
1928	!-----------------------------------------------------------------------------------------
1929
1930	if (mydebug) then
1931	call writefield_phy('u_seri',u_seri,llm)
1932	call writefield_phy('v_seri',v_seri,llm)
1933	call writefield_phy('t_seri',t_seri,llm)
1934	call writefield_phy('q_seri',q_seri,llm)
1935	endif
1936
1937
1938	IF (ip_ebil_phy.ge.2) THEN
1939	ztit='after surface_main'
1940	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
1941	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1942	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1943	call diagphy(airephy,ztit,ip_ebil_phy
1944	e , zero_v, zero_v, zero_v, zero_v, sens
1945	e , evap , zero_v, zero_v, ztsol
1946	e , d_h_vcol, d_qt, d_ec
1947	s , fs_bound, fq_bound )
1948	END IF
1949
1950	c =================================================================== c
1951	c Calcul de Qsat
1952
1953	DO k = 1, klev
1954	DO i = 1, klon
1955	zx_t = t_seri(i,k)
1956	IF (thermcep) THEN
1957	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
1958	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
1959	zx_qs = MIN(0.5,zx_qs)
1960	zcor = 1./(1.-retv*zx_qs)
1961	zx_qs = zx_qs*zcor
1962	ELSE
1963	IF (zx_t.LT.t_coup) THEN
1964	zx_qs = qsats(zx_t)/pplay(i,k)
1965	ELSE
1966	zx_qs = qsatl(zx_t)/pplay(i,k)
1967	ENDIF
1968	ENDIF
1969	zqsat(i,k)=zx_qs
1970	ENDDO
1971	ENDDO
1972
1973	if (prt_level.ge.1) then
1974	write(lunout,*) 'L qsat (g/kg) avant clouds_gno'
1975	write(lunout,'(i4,f15.4)') (k,1000.*zqsat(igout,k),k=1,klev)
1976	endif
1977	c
1978	c Appeler la convection (au choix)
1979	c
1980	DO k = 1, klev
1981	DO i = 1, klon
1982	conv_q(i,k) = d_q_dyn(i,k)
1983	. + d_q_vdf(i,k)/dtime
1984	conv_t(i,k) = d_t_dyn(i,k)
1985	. + d_t_vdf(i,k)/dtime
1986	ENDDO
1987	ENDDO
1988	IF (check) THEN
1989	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
1990	WRITE(lunout,*) "avantcon=", za
1991	ENDIF
1992	zx_ajustq = .FALSE.
1993	IF (iflag_con.EQ.2) zx_ajustq=.TRUE.
1994	IF (zx_ajustq) THEN
1995	DO i = 1, klon
1996	z_avant(i) = 0.0
1997	ENDDO
1998	DO k = 1, klev
1999	DO i = 1, klon
2000	z_avant(i) = z_avant(i) + (q_seri(i,k)+ql_seri(i,k))
2001	. *(paprs(i,k)-paprs(i,k+1))/RG
2002	ENDDO
2003	ENDDO
2004	ENDIF
2005
2006	c Calcule de vitesse verticale a partir de flux de masse verticale
2007	DO k = 1, klev
2008	DO i = 1, klon
2009	omega(i,k) = RG*flxmass_w(i,k) / airephy(i)
2010	END DO
2011	END DO
2012	if (prt_level.ge.1) write(lunout,*) 'omega(igout, :) = ',
2013	$ omega(igout, :)
2014
2015	IF (iflag_con.EQ.1) THEN
2016	abort_message ='reactiver le call conlmd dans physiq.F'
2017	CALL abort_gcm (modname,abort_message,1)
2018	c CALL conlmd (dtime, paprs, pplay, t_seri, q_seri, conv_q,
2019	c . d_t_con, d_q_con,
2020	c . rain_con, snow_con, ibas_con, itop_con)
2021	ELSE IF (iflag_con.EQ.2) THEN
2022	CALL conflx(dtime, paprs, pplay, t_seri, q_seri,
2023	e conv_t, conv_q, -evap, omega,
2024	s d_t_con, d_q_con, rain_con, snow_con,
2025	s pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
2026	s kcbot, kctop, kdtop, pmflxr, pmflxs)
2027	d_u_con = 0.
2028	d_v_con = 0.
2029
2030	WHERE (rain_con < 0.) rain_con = 0.
2031	WHERE (snow_con < 0.) snow_con = 0.
2032	DO i = 1, klon
2033	ibas_con(i) = klev+1 - kcbot(i)
2034	itop_con(i) = klev+1 - kctop(i)
2035	ENDDO
2036	ELSE IF (iflag_con.GE.3) THEN
2037	c nb of tracers for the KE convection:
2038	c MAF la partie traceurs est faite dans phytrac
2039	c on met ntra=1 pour limiter les appels mais on peut
2040	c supprimer les calculs / ftra.
2041	ntra = 1
2042
2043	c=====================================================================================
2044	cajout pour la parametrisation des poches froides:
2045	ccalcul de t_wake et t_undi: si pas de poches froides, t_wake=t_undi=t_seri
2046	do k=1,klev
2047	do i=1,klon
2048	if (iflag_wake.eq.1) then
2049	t_wake(i,k) = t_seri(i,k)
2050	. +(1-wake_s(i))*wake_deltat(i,k)
2051	q_wake(i,k) = q_seri(i,k)
2052	. +(1-wake_s(i))*wake_deltaq(i,k)
2053	t_undi(i,k) = t_seri(i,k)
2054	. -wake_s(i)*wake_deltat(i,k)
2055	q_undi(i,k) = q_seri(i,k)
2056	. -wake_s(i)*wake_deltaq(i,k)
2057	else
2058	t_wake(i,k) = t_seri(i,k)
2059	q_wake(i,k) = q_seri(i,k)
2060	t_undi(i,k) = t_seri(i,k)
2061	q_undi(i,k) = q_seri(i,k)
2062	endif
2063	enddo
2064	enddo
2065
2066	cc-- Calcul de l'energie disponible ALE (J/kg) et de la puissance disponible ALP (W/m2)
2067	cc-- pour le soulevement des particules dans le modele convectif
2068	c
2069	do i = 1,klon
2070	ALE(i) = 0.
2071	ALP(i) = 0.
2072	enddo
2073	c
2074	ccalcul de ale_wake et alp_wake
2075	if (iflag_wake.eq.1) then
2076	if (itap .le. it_wape_prescr) then
2077	do i = 1,klon
2078	ale_wake(i) = wape_prescr
2079	alp_wake(i) = fip_prescr
2080	enddo
2081	else
2082	do i = 1,klon
2083	cjyg ALE=WAPE au lieu de ALE = 1/2 Cstar**2
2084	ccc ale_wake(i) = 0.5wake_cstar(i)*2
2085	ale_wake(i) = wake_pe(i)
2086	alp_wake(i) = wake_fip(i)
2087	enddo
2088	endif
2089	else
2090	do i = 1,klon
2091	ale_wake(i) = 0.
2092	alp_wake(i) = 0.
2093	enddo
2094	endif
2095	ccombinaison avec ale et alp de couche limite: constantes si pas de couplage, valeurs calculees
2096	cdans le thermique sinon
2097	if (iflag_coupl.eq.0) then
2098	if (debut.and.prt_level.gt.9)
2099	$ WRITE(lunout,*)'ALE et ALP imposes'
2100	do i = 1,klon
2101	con ne couple que ale
2102	c ALE(i) = max(ale_wake(i),Ale_bl(i))
2103	ALE(i) = max(ale_wake(i),ale_bl_prescr)
2104	con ne couple que alp
2105	c ALP(i) = alp_wake(i) + Alp_bl(i)
2106	ALP(i) = alp_wake(i) + alp_bl_prescr
2107	enddo
2108	else
2109	IF(prt_level>9)WRITE(lunout,*)'ALE et ALP couples au thermique'
2110	! do i = 1,klon
2111	! ALE(i) = max(ale_wake(i),Ale_bl(i))
2112	! avant ALP(i) = alp_wake(i) + Alp_bl(i)
2113	! ALP(i) = alp_wake(i) + Alp_bl(i) + alp_offset ! modif sb
2114	! write(20,*)'ALE',ALE(i),Ale_bl(i),ale_wake(i)
2115	! write(21,*)'ALP',ALP(i),Alp_bl(i),alp_wake(i)
2116	! enddo
2117
2118	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2119	! Modif FH 2010/04/27. Sans doute temporaire.
2120	! Deux options pour le alp_offset : constant si >Ã 0 ou proportionnel Ãa
2121	! w si <0
2122	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2123	do i = 1,klon
2124	ALE(i) = max(ale_wake(i),Ale_bl(i))
2125	if (alp_offset>=0.) then
2126	ALP(i) = alp_wake(i) + Alp_bl(i) + alp_offset ! modif sb
2127	else
2128	ALP(i)=alp_wake(i)+Alp_bl(i)+alp_offset*min(omega(i,6),0.)
2129	if (alp(i)<0.) then
2130	print*,'ALP ',alp(i),alp_wake(i)
2131	s ,Alp_bl(i),alp_offset*min(omega(i,6),0.)
2132	endif
2133	endif
2134	enddo
2135	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2136
2137	endif
2138	do i=1,klon
2139	if (alp(i)>alp_max) then
2140	IF(prt_level>9)WRITE(lunout,*) &
2141	& 'WARNING SUPER ALP (seuil=',alp_max,
2142	, '): i, alp, alp_wake,ale',i,alp(i),alp_wake(i),ale(i)
2143	alp(i)=alp_max
2144	endif
2145	if (ale(i)>ale_max) then
2146	IF(prt_level>9)WRITE(lunout,*) &
2147	& 'WARNING SUPER ALE (seuil=',ale_max,
2148	, '): i, alp, alp_wake,ale',i,ale(i),ale_wake(i),alp(i)
2149	ale(i)=ale_max
2150	endif
2151	enddo
2152
2153	cfin calcul ale et alp
2154	c=================================================================================================
2155
2156
2157	c sb, oct02:
2158	c Schema de convection modularise et vectorise:
2159	c (driver commun aux versions 3 et 4)
2160	c
2161	IF (ok_cvl) THEN ! new driver for convectL
2162
2163	CALL concvl (iflag_con,iflag_clos,
2164	. dtime,paprs,pplay,t_undi,q_undi,
2165	. t_wake,q_wake,wake_s,
2166	. u_seri,v_seri,tr_seri,nbtr,
2167	. ALE,ALP,
2168	. ema_work1,ema_work2,
2169	. d_t_con,d_q_con,d_u_con,d_v_con,d_tr,
2170	. rain_con, snow_con, ibas_con, itop_con, sigd,
2171	. upwd,dnwd,dnwd0,
2172	. Ma,mip,Vprecip,cape,cin,tvp,Tconv,iflagctrl,
2173	. pbase,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr,qcondc,wd,
2174	. pmflxr,pmflxs,da,phi,mp,
2175	. ftd,fqd,lalim_conv,wght_th)
2176
2177	cIM begin
2178	c print*,'physiq: cin pbase dnwd0 ftd fqd ',cin(1),pbase(1),
2179	c .dnwd0(1,1),ftd(1,1),fqd(1,1)
2180	cIM end
2181	cIM cf. FH
2182	clwcon0=qcondc
2183	pmfu(:,:)=upwd(:,:)+dnwd(:,:)
2184
2185	do i = 1, klon
2186	if (iflagctrl(i).le.1) itau_con(i)=itau_con(i)+1
2187	enddo
2188
2189	ELSE ! ok_cvl
2190	c MAF conema3 ne contient pas les traceurs
2191	CALL conema3 (dtime,
2192	. paprs,pplay,t_seri,q_seri,
2193	. u_seri,v_seri,tr_seri,ntra,
2194	. ema_work1,ema_work2,
2195	. d_t_con,d_q_con,d_u_con,d_v_con,d_tr,
2196	. rain_con, snow_con, ibas_con, itop_con,
2197	. upwd,dnwd,dnwd0,bas,top,
2198	. Ma,cape,tvp,rflag,
2199	. pbase
2200	. ,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr
2201	. ,clwcon0)
2202
2203	ENDIF ! ok_cvl
2204
2205	c
2206	c Correction precip
2207	rain_con = rain_con * cvl_corr
2208	snow_con = snow_con * cvl_corr
2209	c
2210
2211	IF (.NOT. ok_gust) THEN
2212	do i = 1, klon
2213	wd(i)=0.0
2214	enddo
2215	ENDIF
2216
2217	c =================================================================== c
2218	c Calcul des proprietes des nuages convectifs
2219	c
2220
2221	c calcul des proprietes des nuages convectifs
2222	clwcon0(:,:)=fact_cldcon*clwcon0(:,:)
2223	call clouds_gno
2224	s (klon,klev,q_seri,zqsat,clwcon0,ptconv,ratqsc,rnebcon0)
2225
2226	c =================================================================== c
2227
2228	DO i = 1, klon
2229	ema_pcb(i) = pbase(i)
2230	ENDDO
2231	DO i = 1, klon
2232
2233	! L'idicage de itop_con peut cacher un pb potentiel
2234	! FH sous la dictee de JYG, CR
2235	ema_pct(i) = paprs(i,itop_con(i)+1)
2236
2237	if (itop_con(i).gt.klev-3) then
2238	if(prt_level >= 9) then
2239	write(lunout,*)'La convection monte trop haut '
2240	write(lunout,*)'itop_con(,',i,',)=',itop_con(i)
2241	endif
2242	endif
2243	ENDDO
2244	DO i = 1, klon
2245	ema_cbmf(i) = ema_workcbmf(i)
2246	ENDDO
2247	ELSE IF (iflag_con.eq.0) THEN
2248	write(lunout,*) 'On n appelle pas la convection'
2249	clwcon0=0.
2250	rnebcon0=0.
2251	d_t_con=0.
2252	d_q_con=0.
2253	d_u_con=0.
2254	d_v_con=0.
2255	rain_con=0.
2256	snow_con=0.
2257	bas=1
2258	top=1
2259	ELSE
2260	WRITE(lunout,*) "iflag_con non-prevu", iflag_con
2261	CALL abort
2262	ENDIF
2263
2264	c CALL homogene(paprs, q_seri, d_q_con, u_seri,v_seri,
2265	c . d_u_con, d_v_con)
2266
2267	!-----------------------------------------------------------------------------------------
2268	! ajout des tendances de la diffusion turbulente
2269	CALL add_phys_tend(d_u_con,d_v_con,d_t_con,d_q_con,dql0,'con')
2270	!-----------------------------------------------------------------------------------------
2271
2272	if (mydebug) then
2273	call writefield_phy('u_seri',u_seri,llm)
2274	call writefield_phy('v_seri',v_seri,llm)
2275	call writefield_phy('t_seri',t_seri,llm)
2276	call writefield_phy('q_seri',q_seri,llm)
2277	endif
2278
2279	cIM
2280	IF (ip_ebil_phy.ge.2) THEN
2281	ztit='after convect'
2282	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
2283	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2284	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2285	call diagphy(airephy,ztit,ip_ebil_phy
2286	e , zero_v, zero_v, zero_v, zero_v, zero_v
2287	e , zero_v, rain_con, snow_con, ztsol
2288	e , d_h_vcol, d_qt, d_ec
2289	s , fs_bound, fq_bound )
2290	END IF
2291	C
2292	IF (check) THEN
2293	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
2294	WRITE(lunout,*)"aprescon=", za
2295	zx_t = 0.0
2296	za = 0.0
2297	DO i = 1, klon
2298	za = za + airephy(i)/REAL(klon)
2299	zx_t = zx_t + (rain_con(i)+
2300	. snow_con(i))*airephy(i)/REAL(klon)
2301	ENDDO
2302	zx_t = zx_t/za*dtime
2303	WRITE(lunout,*)"Precip=", zx_t
2304	ENDIF
2305	IF (zx_ajustq) THEN
2306	DO i = 1, klon
2307	z_apres(i) = 0.0
2308	ENDDO
2309	DO k = 1, klev
2310	DO i = 1, klon
2311	z_apres(i) = z_apres(i) + (q_seri(i,k)+ql_seri(i,k))
2312	. *(paprs(i,k)-paprs(i,k+1))/RG
2313	ENDDO
2314	ENDDO
2315	DO i = 1, klon
2316	z_factor(i) = (z_avant(i)-(rain_con(i)+snow_con(i))*dtime)
2317	. /z_apres(i)
2318	ENDDO
2319	DO k = 1, klev
2320	DO i = 1, klon
2321	IF (z_factor(i).GT.(1.0+1.0E-08) .OR.
2322	. z_factor(i).LT.(1.0-1.0E-08)) THEN
2323	q_seri(i,k) = q_seri(i,k) * z_factor(i)
2324	ENDIF
2325	ENDDO
2326	ENDDO
2327	ENDIF
2328	zx_ajustq=.FALSE.
2329
2330	c
2331	c=============================================================================
2332	cRR:Evolution de la poche froide: on ne fait pas de separation wake/env
2333	cpour la couche limite diffuse pour l instant
2334	c
2335	if (iflag_wake.eq.1) then
2336	DO k=1,klev
2337	DO i=1,klon
2338	dt_dwn(i,k) = ftd(i,k)
2339	wdt_PBL(i,k) = 0.
2340	dq_dwn(i,k) = fqd(i,k)
2341	wdq_PBL(i,k) = 0.
2342	M_dwn(i,k) = dnwd0(i,k)
2343	M_up(i,k) = upwd(i,k)
2344	dt_a(i,k) = d_t_con(i,k)/dtime - ftd(i,k)
2345	udt_PBL(i,k) = 0.
2346	dq_a(i,k) = d_q_con(i,k)/dtime - fqd(i,k)
2347	udq_PBL(i,k) = 0.
2348	ENDDO
2349	ENDDO
2350	c
2351	ccalcul caracteristiques de la poche froide
2352	call calWAKE (paprs,pplay,dtime
2353	: ,t_seri,q_seri,omega
2354	: ,dt_dwn,dq_dwn,M_dwn,M_up
2355	: ,dt_a,dq_a,sigd
2356	: ,wdt_PBL,wdq_PBL
2357	: ,udt_PBL,udq_PBL
2358	o ,wake_deltat,wake_deltaq,wake_dth
2359	o ,wake_h,wake_s,wake_dens
2360	o ,wake_pe,wake_fip,wake_gfl
2361	o ,dt_wake,dq_wake
2362	o ,wake_k, t_undi,q_undi
2363	o ,wake_omgbdth,wake_dp_omgb
2364	o ,wake_dtKE,wake_dqKE
2365	o ,wake_dtPBL,wake_dqPBL
2366	o ,wake_omg,wake_dp_deltomg
2367	o ,wake_spread,wake_Cstar,wake_d_deltat_gw
2368	o ,wake_ddeltat,wake_ddeltaq)
2369	c
2370	!-----------------------------------------------------------------------------------------
2371	! ajout des tendances des poches froides
2372	! Faire rapidement disparaitre l'ancien dt_wake pour garder un d_t_wake
2373	! coherent avec les autres d_t_...
2374	d_t_wake(:,:)=dt_wake(:,:)*dtime
2375	d_q_wake(:,:)=dq_wake(:,:)*dtime
2376	CALL add_phys_tend(du0,dv0,d_t_wake,d_q_wake,dql0,'wake')
2377	!-----------------------------------------------------------------------------------------
2378
2379	endif
2380	c
2381	c===================================================================
2382	cJYG
2383	IF (ip_ebil_phy.ge.2) THEN
2384	ztit='after wake'
2385	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
2386	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2387	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2388	call diagphy(airephy,ztit,ip_ebil_phy
2389	e , zero_v, zero_v, zero_v, zero_v, zero_v
2390	e , zero_v, zero_v, zero_v, ztsol
2391	e , d_h_vcol, d_qt, d_ec
2392	s , fs_bound, fq_bound )
2393	END IF
2394
2395	c print*,'apres callwake iflag_cldcon=', iflag_cldcon
2396	c
2397	c===================================================================
2398	c Convection seche (thermiques ou ajustement)
2399	c===================================================================
2400	c
2401	call stratocu_if(klon,klev,pctsrf,paprs, pplay,t_seri
2402	s ,seuil_inversion,weak_inversion,dthmin)
2403
2404
2405
2406	d_t_ajsb(:,:)=0.
2407	d_q_ajsb(:,:)=0.
2408	d_t_ajs(:,:)=0.
2409	d_u_ajs(:,:)=0.
2410	d_v_ajs(:,:)=0.
2411	d_q_ajs(:,:)=0.
2412	clwcon0th(:,:)=0.
2413	c
2414	c fm_therm(:,:)=0.
2415	c entr_therm(:,:)=0.
2416	c detr_therm(:,:)=0.
2417	c
2418	IF(prt_level>9)WRITE(lunout,*)
2419	. 'AVANT LA CONVECTION SECHE , iflag_thermals='
2420	s ,iflag_thermals,' nsplit_thermals=',nsplit_thermals
2421	if(iflag_thermals.lt.0) then
2422	c Rien
2423	c ====
2424	IF(prt_level>9)WRITE(lunout,*)'pas de convection'
2425
2426
2427	else
2428
2429	c Thermiques
2430	c ==========
2431	IF(prt_level>9)WRITE(lunout,*)'JUSTE AVANT , iflag_thermals='
2432	s ,iflag_thermals,' nsplit_thermals=',nsplit_thermals
2433
2434
2435	if (iflag_thermals.gt.1) then
2436	call calltherm(pdtphys
2437	s ,pplay,paprs,pphi,weak_inversion
2438	s ,u_seri,v_seri,t_seri,q_seri,zqsat,debut
2439	s ,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs
2440	s ,fm_therm,entr_therm,detr_therm
2441	s ,zqasc,clwcon0th,lmax_th,ratqscth
2442	s ,ratqsdiff,zqsatth
2443	con rajoute ale et alp, et les caracteristiques de la couche alim
2444	s ,Ale_bl,Alp_bl,lalim_conv,wght_th, zmax0, f0, zw2,fraca)
2445
2446	! ----------------------------------------------------------------------
2447	! Transport de la TKE par les panaches thermiques.
2448	! FH : 2010/02/01
2449	if (iflag_pbl.eq.10) then
2450	call thermcell_dtke(klon,klev,nbsrf,pdtphys,fm_therm,entr_therm,
2451	s rg,paprs,pbl_tke)
2452	endif
2453	! ----------------------------------------------------------------------
2454
2455	endif
2456
2457
2458
2459	c Ajustement sec
2460	c ==============
2461
2462	! Dans le cas où on active les thermiques, on fait partir l'ajustement
2463	! a partir du sommet des thermiques.
2464	! Dans le cas contraire, on demarre au niveau 1.
2465
2466	if (iflag_thermals.ge.13.or.iflag_thermals.eq.0) then
2467
2468	if(iflag_thermals.eq.0) then
2469	IF(prt_level>9)WRITE(lunout,*)'ajsec'
2470	limbas(:)=1
2471	else
2472	limbas(:)=lmax_th(:)
2473	endif
2474
2475	! Attention : le call ajsec_convV2 n'est maintenu que momentanneement
2476	! pour des test de convergence numerique.
2477	! Le nouveau ajsec est a priori mieux, meme pour le cas
2478	! iflag_thermals = 0 (l'ancienne version peut faire des tendances
2479	! non nulles numeriquement pour des mailles non concernees.
2480
2481	if (iflag_thermals.eq.0) then
2482	CALL ajsec_convV2(paprs, pplay, t_seri,q_seri
2483	s , d_t_ajsb, d_q_ajsb)
2484	else
2485	CALL ajsec(paprs, pplay, t_seri,q_seri,limbas
2486	s , d_t_ajsb, d_q_ajsb)
2487	endif
2488
2489	!-----------------------------------------------------------------------------------------
2490	! ajout des tendances de l'ajustement sec ou des thermiques
2491	CALL add_phys_tend(du0,dv0,d_t_ajsb,d_q_ajsb,dql0,'ajsb')
2492	d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_ajsb(:,:)
2493	d_q_ajs(:,:)=d_q_ajs(:,:)+d_q_ajsb(:,:)
2494
2495	!-----------------------------------------------------------------------------------------
2496
2497	endif
2498
2499	endif
2500	c
2501	c===================================================================
2502	cIM
2503	IF (ip_ebil_phy.ge.2) THEN
2504	ztit='after dry_adjust'
2505	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
2506	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2507	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2508	call diagphy(airephy,ztit,ip_ebil_phy
2509	e , zero_v, zero_v, zero_v, zero_v, zero_v
2510	e , zero_v, zero_v, zero_v, ztsol
2511	e , d_h_vcol, d_qt, d_ec
2512	s , fs_bound, fq_bound )
2513	END IF
2514
2515
2516	c-------------------------------------------------------------------------
2517	c Caclul des ratqs
2518	c-------------------------------------------------------------------------
2519
2520	c print*,'calcul des ratqs'
2521	c ratqs convectifs a l'ancienne en fonction de q(z=0)-q / q
2522	c ----------------
2523	c on ecrase le tableau ratqsc calcule par clouds_gno
2524	if (iflag_cldcon.eq.1) then
2525	do k=1,klev
2526	do i=1,klon
2527	if(ptconv(i,k)) then
2528	ratqsc(i,k)=ratqsbas
2529	s +fact_cldcon*(q_seri(i,1)-q_seri(i,k))/q_seri(i,k)
2530	else
2531	ratqsc(i,k)=0.
2532	endif
2533	enddo
2534	enddo
2535
2536	c-----------------------------------------------------------------------
2537	c par nversion de la fonction log normale
2538	c-----------------------------------------------------------------------
2539	else if (iflag_cldcon.eq.4) then
2540	ptconvth(:,:)=.false.
2541	ratqsc(:,:)=0.
2542	if(prt_level.ge.9) print*,'avant clouds_gno thermique'
2543	call clouds_gno
2544	s (klon,klev,q_seri,zqsat,clwcon0th,ptconvth,ratqsc,rnebcon0th)
2545	if(prt_level.ge.9) print*,' CLOUDS_GNO OK'
2546
2547	c-----------------------------------------------------------------------
2548	c par calcul direct de l'ecart-type
2549	c-----------------------------------------------------------------------
2550
2551	else if (iflag_cldcon>=5) then
2552	wmax_th(:)=0.
2553	zmax_th(:)=0.
2554	do k=1,klev
2555	do i=1,klon
2556	wmax_th(i)=max(wmax_th(i),zw2(i,k))
2557	if (detr_therm(i,k).gt.0.) zmax_th(i)=pphi(i,k)/rg
2558	enddo
2559	enddo
2560	tau_overturning_th(:)=zmax_th(:)/max(0.5*wmax_th(:),0.1)
2561	if(prt_level.ge.9)
2562	& write(lunout,*)'TAU TH OK ',
2563	& tau_overturning_th(1),detr_therm(1,3)
2564
2565	c On impose que l'air autour de la fraction couverte par le thermique
2566	c plus son air detraine durant tau_overturning_th soit superieur
2567	c a 0.1 q_seri
2568	zz=0.1
2569	do k=1,klev
2570	do i=1,klon
2571	lambda_th(i,k)=0.5*(fraca(i,k)+fraca(i,k+1))+
2572	s tau_overturning_th(i)*detr_therm(i,k)
2573	s *rg/(paprs(i,k)-paprs(i,k+1))
2574	znum=(1.-zz)*q_seri(i,k)
2575	zden=zqasc(i,k)-zz*q_seri(i,k)
2576	if (znum-lambda_th(i,k)*zden<0.) lambda_th(i,k)=znum/zden
2577	lambda_th(i,k)=min(lambda_th(i,k),0.9)
2578	enddo
2579	enddo
2580
2581	if(iflag_cldcon==5) then
2582	do k=1,klev
2583	do i=1,klon
2584	ratqsc(i,k)=sqrt(lambda_th(i,k)/(1.-lambda_th(i,k)))*
2585	s abs((zqasc(i,k)-q_seri(i,k))/q_seri(i,k))
2586	enddo
2587	enddo
2588	else if(iflag_cldcon==6) then
2589	do k=1,klev
2590	do i=1,klon
2591	ratqsc(i,k)=sqrt(lambda_th(i,k))*
2592	s (zqasc(i,k)-q_seri(i,k))/q_seri(i,k)
2593	enddo
2594	enddo
2595	endif
2596
2597	endif
2598
2599	c ratqs stables
2600	c -------------
2601
2602	if (iflag_ratqs.eq.0) then
2603
2604	! Le cas iflag_ratqs=0 correspond a la version IPCC 2005 du modele.
2605	do k=1,klev
2606	do i=1, klon
2607	ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*
2608	s min((paprs(i,1)-pplay(i,k))/(paprs(i,1)-30000.),1.)
2609	enddo
2610	enddo
2611
2612	! Pour iflag_ratqs=1 ou 2, le ratqs est constant au dessus de
2613	! 300 hPa (ratqshaut), varie lineariement en fonction de la pression
2614	! entre 600 et 300 hPa et est soit constant (ratqsbas) pour iflag_ratqs=1
2615	! soit lineaire (entre 0 a la surface et ratqsbas) pour iflag_ratqs=2
2616	! Il s'agit de differents tests dans la phase de reglage du modele
2617	! avec thermiques.
2618
2619	else if (iflag_ratqs.eq.1) then
2620
2621	do k=1,klev
2622	do i=1, klon
2623	if (pplay(i,k).ge.60000.) then
2624	ratqss(i,k)=ratqsbas
2625	else if ((pplay(i,k).ge.30000.).and.
2626	s (pplay(i,k).lt.60000.)) then
2627	ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*
2628	s (60000.-pplay(i,k))/(60000.-30000.)
2629	else
2630	ratqss(i,k)=ratqshaut
2631	endif
2632	enddo
2633	enddo
2634
2635	else
2636
2637	do k=1,klev
2638	do i=1, klon
2639	if (pplay(i,k).ge.60000.) then
2640	ratqss(i,k)=ratqsbas
2641	s *(paprs(i,1)-pplay(i,k))/(paprs(i,1)-60000.)
2642	else if ((pplay(i,k).ge.30000.).and.
2643	s (pplay(i,k).lt.60000.)) then
2644	ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*
2645	s (60000.-pplay(i,k))/(60000.-30000.)
2646	else
2647	ratqss(i,k)=ratqshaut
2648	endif
2649	enddo
2650	enddo
2651	endif
2652
2653
2654
2655
2656	c ratqs final
2657	c -----------
2658
2659	if (iflag_cldcon.eq.1 .or.iflag_cldcon.eq.2
2660	s .or.iflag_cldcon.ge.4) then
2661
2662	! On ajoute une constante au ratqsc*2 pour tenir compte de
2663	! fluctuations turbulentes de petite echelle
2664
2665	do k=1,klev
2666	do i=1,klon
2667	if ((fm_therm(i,k).gt.1.e-10)) then
2668	ratqsc(i,k)=sqrt(ratqsc(i,k)2+0.052)
2669	endif
2670	enddo
2671	enddo
2672
2673	! les ratqs sont une combinaison de ratqss et ratqsc
2674	if(prt_level.ge.9)
2675	$ write(lunout,*)'PHYLMD NOUVEAU TAU_RATQS ',tau_ratqs
2676
2677	if (tau_ratqs>1.e-10) then
2678	facteur=exp(-pdtphys/tau_ratqs)
2679	else
2680	facteur=0.
2681	endif
2682	ratqs(:,:)=ratqsc(:,:)(1.-facteur)+ratqs(:,:)facteur
2683	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2684	! FH 22/09/2009
2685	! La ligne ci-dessous faisait osciller le modele et donnait une solution
2686	! assymptotique bidon et dÃ©pendant fortement du pas de temps.
2687	! ratqs(:,:)=sqrt(ratqs(:,:)2+ratqss(:,:)2)
2688	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2689	ratqs(:,:)=max(ratqs(:,:),ratqss(:,:))
2690	else
2691	! on ne prend que le ratqs stable pour fisrtilp
2692	ratqs(:,:)=ratqss(:,:)
2693	endif
2694
2695
2696	c
2697	c Appeler le processus de condensation a grande echelle
2698	c et le processus de precipitation
2699	c-------------------------------------------------------------------------
2700	CALL fisrtilp(dtime,paprs,pplay,
2701	. t_seri, q_seri,ptconv,ratqs,
2702	. d_t_lsc, d_q_lsc, d_ql_lsc, rneb, cldliq,
2703	. rain_lsc, snow_lsc,
2704	. pfrac_impa, pfrac_nucl, pfrac_1nucl,
2705	. frac_impa, frac_nucl,
2706	. prfl, psfl, rhcl)
2707
2708	WHERE (rain_lsc < 0) rain_lsc = 0.
2709	WHERE (snow_lsc < 0) snow_lsc = 0.
2710	!-----------------------------------------------------------------------------------------
2711	! ajout des tendances de la diffusion turbulente
2712	CALL add_phys_tend(du0,dv0,d_t_lsc,d_q_lsc,d_ql_lsc,'lsc')
2713	!-----------------------------------------------------------------------------------------
2714	DO k = 1, klev
2715	DO i = 1, klon
2716	cldfra(i,k) = rneb(i,k)
2717	IF (.NOT.new_oliq) cldliq(i,k) = ql_seri(i,k)
2718	ENDDO
2719	ENDDO
2720	IF (check) THEN
2721	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
2722	WRITE(lunout,*)"apresilp=", za
2723	zx_t = 0.0
2724	za = 0.0
2725	DO i = 1, klon
2726	za = za + airephy(i)/REAL(klon)
2727	zx_t = zx_t + (rain_lsc(i)
2728	. + snow_lsc(i))*airephy(i)/REAL(klon)
2729	ENDDO
2730	zx_t = zx_t/za*dtime
2731	WRITE(lunout,*)"Precip=", zx_t
2732	ENDIF
2733	cIM
2734	IF (ip_ebil_phy.ge.2) THEN
2735	ztit='after fisrt'
2736	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
2737	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2738	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2739	call diagphy(airephy,ztit,ip_ebil_phy
2740	e , zero_v, zero_v, zero_v, zero_v, zero_v
2741	e , zero_v, rain_lsc, snow_lsc, ztsol
2742	e , d_h_vcol, d_qt, d_ec
2743	s , fs_bound, fq_bound )
2744	END IF
2745
2746	if (mydebug) then
2747	call writefield_phy('u_seri',u_seri,llm)
2748	call writefield_phy('v_seri',v_seri,llm)
2749	call writefield_phy('t_seri',t_seri,llm)
2750	call writefield_phy('q_seri',q_seri,llm)
2751	endif
2752
2753	c
2754	c-------------------------------------------------------------------
2755	c PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT
2756	c-------------------------------------------------------------------
2757
2758	c 1. NUAGES CONVECTIFS
2759	c
2760	cIM cf FH
2761	c IF (iflag_cldcon.eq.-1) THEN ! seulement pour Tiedtke
2762	IF (iflag_cldcon.le.-1) THEN ! seulement pour Tiedtke
2763	snow_tiedtke=0.
2764	c print*,'avant calcul de la pseudo precip '
2765	c print*,'iflag_cldcon',iflag_cldcon
2766	if (iflag_cldcon.eq.-1) then
2767	rain_tiedtke=rain_con
2768	else
2769	c print*,'calcul de la pseudo precip '
2770	rain_tiedtke=0.
2771	c print*,'calcul de la pseudo precip 0'
2772	do k=1,klev
2773	do i=1,klon
2774	if (d_q_con(i,k).lt.0.) then
2775	rain_tiedtke(i)=rain_tiedtke(i)-d_q_con(i,k)/pdtphys
2776	s *(paprs(i,k)-paprs(i,k+1))/rg
2777	endif
2778	enddo
2779	enddo
2780	endif
2781	c
2782	c call dump2d(iim,jjm,rain_tiedtke(2:klon-1),'PSEUDO PRECIP ')
2783	c
2784
2785	c Nuages diagnostiques pour Tiedtke
2786	CALL diagcld1(paprs,pplay,
2787	cIM cf FH . rain_con,snow_con,ibas_con,itop_con,
2788	. rain_tiedtke,snow_tiedtke,ibas_con,itop_con,
2789	. diafra,dialiq)
2790	DO k = 1, klev
2791	DO i = 1, klon
2792	IF (diafra(i,k).GT.cldfra(i,k)) THEN
2793	cldliq(i,k) = dialiq(i,k)
2794	cldfra(i,k) = diafra(i,k)
2795	ENDIF
2796	ENDDO
2797	ENDDO
2798
2799	ELSE IF (iflag_cldcon.ge.3) THEN
2800	c On prend pour les nuages convectifs le max du calcul de la
2801	c convection et du calcul du pas de temps precedent diminue d'un facteur
2802	c facttemps
2803	facteur = pdtphys *facttemps
2804	do k=1,klev
2805	do i=1,klon
2806	rnebcon(i,k)=rnebcon(i,k)*facteur
2807	if (rnebcon0(i,k)clwcon0(i,k).gt.rnebcon(i,k)clwcon(i,k))
2808	s then
2809	rnebcon(i,k)=rnebcon0(i,k)
2810	clwcon(i,k)=clwcon0(i,k)
2811	endif
2812	enddo
2813	enddo
2814
2815	c
2816	cjq - introduce the aerosol direct and first indirect radiative forcings
2817	cjq - Johannes Quaas, 27/11/2003 (quaas@lmd.jussieu.fr)
2818	IF (ok_ade.OR.ok_aie) THEN
2819	IF (.NOT. aerosol_couple)
2820	& CALL readaerosol_optic(
2821	& debut, new_aod, flag_aerosol, itap, jD_cur-jD_ref,
2822	& pdtphys, pplay, paprs, t_seri, rhcl, presnivs,
2823	& mass_solu_aero, mass_solu_aero_pi,
2824	& tau_aero, piz_aero, cg_aero,
2825	& tausum_aero, tau3d_aero)
2826	ELSE
2827	cIM 170310 BEG
2828	tausum_aero(:,:,:) = 0.
2829	cIM 170310 END
2830	tau_aero(:,:,:,:) = 0.
2831	piz_aero(:,:,:,:) = 0.
2832	cg_aero(:,:,:,:) = 0.
2833	ENDIF
2834
2835	cIM calcul nuages par le simulateur ISCCP
2836	c
2837	#ifdef histISCCP
2838	IF (ok_isccp) THEN
2839	c
2840	cIM lecture invtau, tautab des fichiers formattes
2841	c
2842	IF (debut) THEN
2843	c$OMP MASTER
2844	c
2845	open(99,file='tautab.formatted', FORM='FORMATTED')
2846	read(99,'(f30.20)') tautab_omp
2847	close(99)
2848	c
2849	open(99,file='invtau.formatted',form='FORMATTED')
2850	read(99,'(i10)') invtau_omp
2851
2852	c print*,'calcul_simulISCCP invtau_omp',invtau_omp
2853	c write(6,'(a,8i10)') 'invtau_omp',(invtau_omp(i),i=1,100)
2854
2855	close(99)
2856	c$OMP END MASTER
2857	c$OMP BARRIER
2858	tautab=tautab_omp
2859	invtau=invtau_omp
2860	c
2861	ENDIF !debut
2862	c
2863	cIM appel simulateur toutes les NINT(freq_ISCCP/dtime) heures
2864	IF (MOD(itap,NINT(freq_ISCCP/dtime)).EQ.0) THEN
2865	#include "calcul_simulISCCP.h"
2866	ENDIF !(MOD(itap,NINT(freq_ISCCP/dtime))
2867	ENDIF !ok_isccp
2868	#endif
2869
2870	c On prend la somme des fractions nuageuses et des contenus en eau
2871	cldfra(:,:)=min(max(cldfra(:,:),rnebcon(:,:)),1.)
2872	cldliq(:,:)=cldliq(:,:)+rnebcon(:,:)*clwcon(:,:)
2873
2874	ENDIF
2875	c
2876	c 2. NUAGES STARTIFORMES
2877	c
2878	IF (ok_stratus) THEN
2879	CALL diagcld2(paprs,pplay,t_seri,q_seri, diafra,dialiq)
2880	DO k = 1, klev
2881	DO i = 1, klon
2882	IF (diafra(i,k).GT.cldfra(i,k)) THEN
2883	cldliq(i,k) = dialiq(i,k)
2884	cldfra(i,k) = diafra(i,k)
2885	ENDIF
2886	ENDDO
2887	ENDDO
2888	ENDIF
2889	c
2890	c Precipitation totale
2891	c
2892	DO i = 1, klon
2893	rain_fall(i) = rain_con(i) + rain_lsc(i)
2894	snow_fall(i) = snow_con(i) + snow_lsc(i)
2895	ENDDO
2896	cIM
2897	IF (ip_ebil_phy.ge.2) THEN
2898	ztit="after diagcld"
2899	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
2900	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2901	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2902	call diagphy(airephy,ztit,ip_ebil_phy
2903	e , zero_v, zero_v, zero_v, zero_v, zero_v
2904	e , zero_v, zero_v, zero_v, ztsol
2905	e , d_h_vcol, d_qt, d_ec
2906	s , fs_bound, fq_bound )
2907	END IF
2908	c
2909	c Calculer l'humidite relative pour diagnostique
2910	c
2911	DO k = 1, klev
2912	DO i = 1, klon
2913	zx_t = t_seri(i,k)
2914	IF (thermcep) THEN
2915	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
2916	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
2917	zx_qs = MIN(0.5,zx_qs)
2918	zcor = 1./(1.-retv*zx_qs)
2919	zx_qs = zx_qs*zcor
2920	ELSE
2921	IF (zx_t.LT.t_coup) THEN
2922	zx_qs = qsats(zx_t)/pplay(i,k)
2923	ELSE
2924	zx_qs = qsatl(zx_t)/pplay(i,k)
2925	ENDIF
2926	ENDIF
2927	zx_rh(i,k) = q_seri(i,k)/zx_qs
2928	zqsat(i,k)=zx_qs
2929	ENDDO
2930	ENDDO
2931
2932	cIM Calcul temp.potentielle a 2m (tpot) et temp. potentielle
2933	c equivalente a 2m (tpote) pour diagnostique
2934	c
2935	DO i = 1, klon
2936	tpot(i)=zt2m(i)(100000./paprs(i,1))*RKAPPA
2937	IF (thermcep) THEN
2938	IF(zt2m(i).LT.RTT) then
2939	Lheat=RLSTT
2940	ELSE
2941	Lheat=RLVTT
2942	ENDIF
2943	ELSE
2944	IF (zt2m(i).LT.RTT) THEN
2945	Lheat=RLSTT
2946	ELSE
2947	Lheat=RLVTT
2948	ENDIF
2949	ENDIF
2950	tpote(i) = tpot(i)*
2951	. EXP((Lheat qsat2m(i))/(RCPDzt2m(i)))
2952	ENDDO
2953
2954	IF (config_inca /= 'none') THEN
2955	#ifdef INCA
2956	CALL VTe(VTphysiq)
2957	CALL VTb(VTinca)
2958	calday = REAL(days_elapsed + 1) + jH_cur
2959
2960	call chemtime(itap+itau_phy-1, date0, dtime)
2961	IF (config_inca == 'aero') THEN
2962	CALL AEROSOL_METEO_CALC(
2963	$ calday,pdtphys,pplay,paprs,t,pmflxr,pmflxs,
2964	$ prfl,psfl,pctsrf,airephy,rlat,rlon,u10m,v10m)
2965	END IF
2966
2967	zxsnow_dummy(:) = 0.0
2968
2969	CALL chemhook_begin (calday,
2970	$ days_elapsed+1,
2971	$ jH_cur,
2972	$ pctsrf(1,1),
2973	$ rlat,
2974	$ rlon,
2975	$ airephy,
2976	$ paprs,
2977	$ pplay,
2978	$ coefh,
2979	$ pphi,
2980	$ t_seri,
2981	$ u,
2982	$ v,
2983	$ wo(:, :, 1),
2984	$ q_seri,
2985	$ zxtsol,
2986	$ zxsnow_dummy,
2987	$ solsw,
2988	$ albsol1,
2989	$ rain_fall,
2990	$ snow_fall,
2991	$ itop_con,
2992	$ ibas_con,
2993	$ cldfra,
2994	$ iim,
2995	$ jjm,
2996	$ tr_seri,
2997	$ ftsol,
2998	$ paprs,
2999	$ cdragh,
3000	$ cdragm,
3001	$ pctsrf,
3002	$ pdtphys,
3003	$ itap)
3004
3005	CALL VTe(VTinca)
3006	CALL VTb(VTphysiq)
3007	#endif
3008	END IF !config_inca /= 'none'
3009	c
3010	c Calculer les parametres optiques des nuages et quelques
3011	c parametres pour diagnostiques:
3012	c
3013
3014	IF (aerosol_couple) THEN
3015	mass_solu_aero(:,:) = ccm(:,:,1)
3016	mass_solu_aero_pi(:,:) = ccm(:,:,2)
3017	END IF
3018
3019	if (ok_newmicro) then
3020	CALL newmicro (paprs, pplay,ok_newmicro,
3021	. t_seri, cldliq, cldfra, cldtau, cldemi,
3022	. cldh, cldl, cldm, cldt, cldq,
3023	. flwp, fiwp, flwc, fiwc,
3024	e ok_aie,
3025	e mass_solu_aero, mass_solu_aero_pi,
3026	e bl95_b0, bl95_b1,
3027	s cldtaupi, re, fl, ref_liq, ref_ice)
3028	else
3029	CALL nuage (paprs, pplay,
3030	. t_seri, cldliq, cldfra, cldtau, cldemi,
3031	. cldh, cldl, cldm, cldt, cldq,
3032	e ok_aie,
3033	e mass_solu_aero, mass_solu_aero_pi,
3034	e bl95_b0, bl95_b1,
3035	s cldtaupi, re, fl)
3036
3037	endif
3038	c
3039	c Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
3040	c
3041	IF (MOD(itaprad,radpas).EQ.0) THEN
3042
3043	DO i = 1, klon
3044	albsol1(i) = falb1(i,is_oce) * pctsrf(i,is_oce)
3045	. + falb1(i,is_lic) * pctsrf(i,is_lic)
3046	. + falb1(i,is_ter) * pctsrf(i,is_ter)
3047	. + falb1(i,is_sic) * pctsrf(i,is_sic)
3048	albsol2(i) = falb2(i,is_oce) * pctsrf(i,is_oce)
3049	. + falb2(i,is_lic) * pctsrf(i,is_lic)
3050	. + falb2(i,is_ter) * pctsrf(i,is_ter)
3051	. + falb2(i,is_sic) * pctsrf(i,is_sic)
3052	ENDDO
3053
3054	if (mydebug) then
3055	call writefield_phy('u_seri',u_seri,llm)
3056	call writefield_phy('v_seri',v_seri,llm)
3057	call writefield_phy('t_seri',t_seri,llm)
3058	call writefield_phy('q_seri',q_seri,llm)
3059	endif
3060
3061	IF (aerosol_couple) THEN
3062	#ifdef INCA
3063	CALL radlwsw_inca
3064	e (kdlon,kflev,dist, rmu0, fract, solaire,
3065	e paprs, pplay,zxtsol,albsol1, albsol2, t_seri,q_seri,
3066	e wo(:, :, 1),
3067	e cldfra, cldemi, cldtau,
3068	s heat,heat0,cool,cool0,radsol,albpla,
3069	s topsw,toplw,solsw,sollw,
3070	s sollwdown,
3071	s topsw0,toplw0,solsw0,sollw0,
3072	s lwdn0, lwdn, lwup0, lwup,
3073	s swdn0, swdn, swup0, swup,
3074	e ok_ade, ok_aie,
3075	e tau_aero, piz_aero, cg_aero,
3076	s topswad_aero, solswad_aero,
3077	s topswad0_aero, solswad0_aero,
3078	s topsw_aero, topsw0_aero,
3079	s solsw_aero, solsw0_aero,
3080	e cldtaupi,
3081	s topswai_aero, solswai_aero)
3082
3083	#endif
3084	ELSE
3085
3086	CALL radlwsw
3087	e (dist, rmu0, fract,
3088	e paprs, pplay,zxtsol,albsol1, albsol2,
3089	e t_seri,q_seri,wo,
3090	e cldfra, cldemi, cldtau,
3091	e ok_ade, ok_aie,
3092	e tau_aero, piz_aero, cg_aero,
3093	e cldtaupi,new_aod,
3094	e zqsat, flwc, fiwc,
3095	s heat,heat0,cool,cool0,radsol,albpla,
3096	s topsw,toplw,solsw,sollw,
3097	s sollwdown,
3098	s topsw0,toplw0,solsw0,sollw0,
3099	s lwdn0, lwdn, lwup0, lwup,
3100	s swdn0, swdn, swup0, swup,
3101	s topswad_aero, solswad_aero,
3102	s topswai_aero, solswai_aero,
3103	o topswad0_aero, solswad0_aero,
3104	o topsw_aero, topsw0_aero,
3105	o solsw_aero, solsw0_aero,
3106	o topswcf_aero, solswcf_aero)
3107
3108
3109	ENDIF ! aerosol_couple
3110	itaprad = 0
3111	ENDIF ! MOD(itaprad,radpas)
3112	itaprad = itaprad + 1
3113
3114	if (iflag_radia.eq.0 .and. prt_level.ge.9) then
3115	print *,'--------------------------------------------------'
3116	print *,'>>>> ATTENTION rayonnement desactive pour ce cas'
3117	print *,'>>>> heat et cool mis a zero '
3118	print *,'--------------------------------------------------'
3119	heat=0.
3120	cool=0.
3121	endif
3122
3123	c
3124	c Ajouter la tendance des rayonnements (tous les pas)
3125	c
3126	DO k = 1, klev
3127	DO i = 1, klon
3128	t_seri(i,k) = t_seri(i,k)
3129	. + (heat(i,k)-cool(i,k)) * dtime/RDAY
3130	ENDDO
3131	ENDDO
3132	c
3133	if (mydebug) then
3134	call writefield_phy('u_seri',u_seri,llm)
3135	call writefield_phy('v_seri',v_seri,llm)
3136	call writefield_phy('t_seri',t_seri,llm)
3137	call writefield_phy('q_seri',q_seri,llm)
3138	endif
3139
3140	cIM
3141	IF (ip_ebil_phy.ge.2) THEN
3142	ztit='after rad'
3143	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
3144	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
3145	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3146	call diagphy(airephy,ztit,ip_ebil_phy
3147	e , topsw, toplw, solsw, sollw, zero_v
3148	e , zero_v, zero_v, zero_v, ztsol
3149	e , d_h_vcol, d_qt, d_ec
3150	s , fs_bound, fq_bound )
3151	END IF
3152	c
3153	c
3154	c Calculer l'hydrologie de la surface
3155	c
3156	c CALL hydrol(dtime,pctsrf,rain_fall, snow_fall, zxevap,
3157	c . agesno, ftsol,fqsurf,fsnow, ruis)
3158	c
3159
3160	c
3161	c Calculer le bilan du sol et la derive de temperature (couplage)
3162	c
3163	DO i = 1, klon
3164	c bils(i) = radsol(i) - sens(i) - evap(i)*RLVTT
3165	c a la demande de JLD
3166	bils(i) = radsol(i) - sens(i) + zxfluxlat(i)
3167	ENDDO
3168	c
3169	cmoddeblott(jan95)
3170	c Appeler le programme de parametrisation de l'orographie
3171	c a l'echelle sous-maille:
3172	c
3173	IF (ok_orodr) THEN
3174	c
3175	c selection des points pour lesquels le shema est actif:
3176	igwd=0
3177	DO i=1,klon
3178	itest(i)=0
3179	c IF ((zstd(i).gt.10.0)) THEN
3180	IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
3181	itest(i)=1
3182	igwd=igwd+1
3183	idx(igwd)=i
3184	ENDIF
3185	ENDDO
3186	c igwdim=MAX(1,igwd)
3187	c
3188	IF (ok_strato) THEN
3189
3190	CALL drag_noro_strato(klon,klev,dtime,paprs,pplay,
3191	e zmea,zstd, zsig, zgam, zthe,zpic,zval,
3192	e igwd,idx,itest,
3193	e t_seri, u_seri, v_seri,
3194	s zulow, zvlow, zustrdr, zvstrdr,
3195	s d_t_oro, d_u_oro, d_v_oro)
3196
3197	ELSE
3198	CALL drag_noro(klon,klev,dtime,paprs,pplay,
3199	e zmea,zstd, zsig, zgam, zthe,zpic,zval,
3200	e igwd,idx,itest,
3201	e t_seri, u_seri, v_seri,
3202	s zulow, zvlow, zustrdr, zvstrdr,
3203	s d_t_oro, d_u_oro, d_v_oro)
3204	ENDIF
3205	c
3206	c ajout des tendances
3207	!-----------------------------------------------------------------------------------------
3208	! ajout des tendances de la trainee de l'orographie
3209	CALL add_phys_tend(d_u_oro,d_v_oro,d_t_oro,dq0,dql0,'oro')
3210	!-----------------------------------------------------------------------------------------
3211	c
3212	ENDIF ! fin de test sur ok_orodr
3213	c
3214	if (mydebug) then
3215	call writefield_phy('u_seri',u_seri,llm)
3216	call writefield_phy('v_seri',v_seri,llm)
3217	call writefield_phy('t_seri',t_seri,llm)
3218	call writefield_phy('q_seri',q_seri,llm)
3219	endif
3220
3221	IF (ok_orolf) THEN
3222	c
3223	c selection des points pour lesquels le shema est actif:
3224	igwd=0
3225	DO i=1,klon
3226	itest(i)=0
3227	IF ((zpic(i)-zmea(i)).GT.100.) THEN
3228	itest(i)=1
3229	igwd=igwd+1
3230	idx(igwd)=i
3231	ENDIF
3232	ENDDO
3233	c igwdim=MAX(1,igwd)
3234	c
3235	IF (ok_strato) THEN
3236
3237	CALL lift_noro_strato(klon,klev,dtime,paprs,pplay,
3238	e rlat,zmea,zstd,zpic,zgam,zthe,zpic,zval,
3239	e igwd,idx,itest,
3240	e t_seri, u_seri, v_seri,
3241	s zulow, zvlow, zustrli, zvstrli,
3242	s d_t_lif, d_u_lif, d_v_lif )
3243
3244	ELSE
3245	CALL lift_noro(klon,klev,dtime,paprs,pplay,
3246	e rlat,zmea,zstd,zpic,
3247	e itest,
3248	e t_seri, u_seri, v_seri,
3249	s zulow, zvlow, zustrli, zvstrli,
3250	s d_t_lif, d_u_lif, d_v_lif)
3251	ENDIF
3252	c
3253	!-----------------------------------------------------------------------------------------
3254	! ajout des tendances de la portance de l'orographie
3255	CALL add_phys_tend(d_u_lif,d_v_lif,d_t_lif,dq0,dql0,'lif')
3256	!-----------------------------------------------------------------------------------------
3257	c
3258	ENDIF ! fin de test sur ok_orolf
3259	C HINES GWD PARAMETRIZATION
3260
3261	IF (ok_hines) then
3262
3263	CALL hines_gwd(klon,klev,dtime,paprs,pplay,
3264	i rlat,t_seri,u_seri,v_seri,
3265	o zustrhi,zvstrhi,
3266	o d_t_hin, d_u_hin, d_v_hin)
3267	c
3268	c ajout des tendances
3269	CALL add_phys_tend(d_u_hin,d_v_hin,d_t_hin,dq0,dql0,'lif')
3270
3271	ENDIF
3272	c
3273
3274	c
3275	cIM cf. FLott BEG
3276	C STRESS NECESSAIRES: TOUTE LA PHYSIQUE
3277
3278	if (mydebug) then
3279	call writefield_phy('u_seri',u_seri,llm)
3280	call writefield_phy('v_seri',v_seri,llm)
3281	call writefield_phy('t_seri',t_seri,llm)
3282	call writefield_phy('q_seri',q_seri,llm)
3283	endif
3284
3285	DO i = 1, klon
3286	zustrph(i)=0.
3287	zvstrph(i)=0.
3288	ENDDO
3289	DO k = 1, klev
3290	DO i = 1, klon
3291	zustrph(i)=zustrph(i)+(u_seri(i,k)-u(i,k))/dtime*
3292	c (paprs(i,k)-paprs(i,k+1))/rg
3293	zvstrph(i)=zvstrph(i)+(v_seri(i,k)-v(i,k))/dtime*
3294	c (paprs(i,k)-paprs(i,k+1))/rg
3295	ENDDO
3296	ENDDO
3297	c
3298	cIM calcul composantes axiales du moment angulaire et couple des montagnes
3299	c
3300	IF (is_sequential) THEN
3301
3302	CALL aaam_bud (27,klon,klev,jD_cur-jD_ref,jH_cur,
3303	C ra,rg,romega,
3304	C rlat,rlon,pphis,
3305	C zustrdr,zustrli,zustrph,
3306	C zvstrdr,zvstrli,zvstrph,
3307	C paprs,u,v,
3308	C aam, torsfc)
3309	ENDIF
3310	cIM cf. FLott END
3311	cIM
3312	IF (ip_ebil_phy.ge.2) THEN
3313	ztit='after orography'
3314	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
3315	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
3316	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3317	call diagphy(airephy,ztit,ip_ebil_phy
3318	e , zero_v, zero_v, zero_v, zero_v, zero_v
3319	e , zero_v, zero_v, zero_v, ztsol
3320	e , d_h_vcol, d_qt, d_ec
3321	s , fs_bound, fq_bound )
3322	END IF
3323	c
3324	c
3325	!====================================================================
3326	! Interface Simulateur COSP (Calipso, ISCCP, MISR, ..)
3327	!====================================================================
3328	! Abderrahmane 24.08.09
3329
3330	IF (ok_cosp) THEN
3331	! adeclarer
3332	#ifdef CPP_COSP
3333	IF (MOD(itap,NINT(freq_cosp/dtime)).EQ.0) THEN
3334
3335	print*,'freq_cosp',freq_cosp
3336	mr_ozone=wo(:, :, 1) * dobson_u * 1e3 / zmasse
3337	! print*,'Dans physiq.F avant appel cosp ref_liq,ref_ice=',
3338	! s ref_liq,ref_ice
3339	call phys_cosp(itap,dtime,freq_cosp,
3340	$ ecrit_mth,ecrit_day,ecrit_hf,overlap,
3341	$ klon,klev,rlon,rlat,presnivs,
3342	$ ref_liq,ref_ice,
3343	$ pctsrf(:,is_ter)+pctsrf(:,is_lic),
3344	$ zu10m,zv10m,
3345	$ zphi,paprs(:,1:klev),pplay,zxtsol,t_seri,
3346	$ qx(:,:,ivap),zx_rh,cldfra,rnebcon,flwc,fiwc,
3347	$ prfl(:,1:klev),psfl(:,1:klev),
3348	$ pmflxr(:,1:klev),pmflxs(:,1:klev),
3349	$ mr_ozone,cldtau, cldemi)
3350	! L calipso2D,calipso3D,cfadlidar,parasolrefl,atb,betamol,
3351	! L cfaddbze,clcalipso2,dbze,cltlidarradar,
3352	! M clMISR,
3353	! R clisccp2,boxtauisccp,boxptopisccp,tclisccp,ctpisccp,
3354	! I tauisccp,albisccp,meantbisccp,meantbclrisccp)
3355
3356	ENDIF
3357
3358	#endif
3359	ENDIF !ok_cosp
3360	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3361	cAA
3362	cAA Installation de l'interface online-offline pour traceurs
3363	cAA
3364	c====================================================================
3365	c Calcul des tendances traceurs
3366	c====================================================================
3367	C
3368
3369	call phytrac (
3370	I itap, days_elapsed+1, jH_cur, debut,
3371	I lafin, dtime, u, v, t,
3372	I paprs, pplay, pmfu, pmfd,
3373	I pen_u, pde_u, pen_d, pde_d,
3374	I cdragh, coefh, fm_therm, entr_therm,
3375	I u1, v1, ftsol, pctsrf,
3376	I rlat, frac_impa, frac_nucl,rlon,
3377	I presnivs, pphis, pphi, albsol1,
3378	I qx(:,:,ivap),rhcl, cldfra, rneb,
3379	I diafra, cldliq, itop_con, ibas_con,
3380	I pmflxr, pmflxs, prfl, psfl,
3381	I da, phi, mp, upwd,
3382	I dnwd, aerosol_couple, flxmass_w,
3383	I tau_aero, piz_aero, cg_aero, ccm,
3384	I rfname,
3385	O tr_seri)
3386
3387	IF (offline) THEN
3388
3389	IF (prt_level.ge.9)
3390	$ print*,'Attention on met a 0 les thermiques pour phystoke'
3391	call phystokenc (
3392	I nlon,klev,pdtphys,rlon,rlat,
3393	I t,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
3394	I fm_therm,entr_therm,
3395	I cdragh,coefh,u1,v1,ftsol,pctsrf,
3396	I frac_impa, frac_nucl,
3397	I pphis,airephy,dtime,itap)
3398
3399
3400	ENDIF
3401
3402	c
3403	c Calculer le transport de l'eau et de l'energie (diagnostique)
3404	c
3405	CALL transp (paprs,zxtsol,
3406	e t_seri, q_seri, u_seri, v_seri, zphi,
3407	s ve, vq, ue, uq)
3408	c
3409	cIM global posePB BEG
3410	IF(1.EQ.0) THEN
3411	c
3412	CALL transp_lay (paprs,zxtsol,
3413	e t_seri, q_seri, u_seri, v_seri, zphi,
3414	s ve_lay, vq_lay, ue_lay, uq_lay)
3415	c
3416	ENDIF !(1.EQ.0) THEN
3417	cIM global posePB END
3418	c Accumuler les variables a stocker dans les fichiers histoire:
3419	c
3420	c+jld ec_conser
3421	DO k = 1, klev
3422	DO i = 1, klon
3423	ZRCPD = RCPD(1.0+RVTMP2q_seri(i,k))
3424	d_t_ec(i,k)=0.5/ZRCPD
3425	$ (u(i,k)2+v(i,k)2-u_seri(i,k)2-v_seri(i,k)*2)
3426	ENDDO
3427	ENDDO
3428
3429	DO k = 1, klev
3430	DO i = 1, klon
3431	t_seri(i,k)=t_seri(i,k)+d_t_ec(i,k)
3432	d_t_ec(i,k) = d_t_ec(i,k)/dtime
3433	END DO
3434	END DO
3435	c-jld ec_conser
3436	cIM
3437	IF (ip_ebil_phy.ge.1) THEN
3438	ztit='after physic'
3439	CALL diagetpq(airephy,ztit,ip_ebil_phy,1,1,dtime
3440	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
3441	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3442	C Comme les tendances de la physique sont ajoute dans la dynamique,
3443	C on devrait avoir que la variation d'entalpie par la dynamique
3444	C est egale a la variation de la physique au pas de temps precedent.
3445	C Donc la somme de ces 2 variations devrait etre nulle.
3446
3447	call diagphy(airephy,ztit,ip_ebil_phy
3448	e , topsw, toplw, solsw, sollw, sens
3449	e , evap, rain_fall, snow_fall, ztsol
3450	e , d_h_vcol, d_qt, d_ec
3451	s , fs_bound, fq_bound )
3452	C
3453	d_h_vcol_phy=d_h_vcol
3454	C
3455	END IF
3456	C
3457	c=======================================================================
3458	c SORTIES
3459	c=======================================================================
3460
3461	cIM Interpolation sur les niveaux de pression du NMC
3462	c -------------------------------------------------
3463	c
3464	#include "calcul_STDlev.h"
3465	twriteSTD(:,:,1)=tsumSTD(:,:,1)
3466	qwriteSTD(:,:,1)=qsumSTD(:,:,1)
3467	rhwriteSTD(:,:,1)=rhsumSTD(:,:,1)
3468	phiwriteSTD(:,:,1)=phisumSTD(:,:,1)
3469	uwriteSTD(:,:,1)=usumSTD(:,:,1)
3470	vwriteSTD(:,:,1)=vsumSTD(:,:,1)
3471	wwriteSTD(:,:,1)=wsumSTD(:,:,1)
3472
3473	twriteSTD(:,:,2)=tsumSTD(:,:,2)
3474	qwriteSTD(:,:,2)=qsumSTD(:,:,2)
3475	rhwriteSTD(:,:,2)=rhsumSTD(:,:,2)
3476	phiwriteSTD(:,:,2)=phisumSTD(:,:,2)
3477	uwriteSTD(:,:,2)=usumSTD(:,:,2)
3478	vwriteSTD(:,:,2)=vsumSTD(:,:,2)
3479	wwriteSTD(:,:,2)=wsumSTD(:,:,2)
3480
3481	twriteSTD(:,:,3)=tlevSTD(:,:)
3482	qwriteSTD(:,:,3)=qlevSTD(:,:)
3483	rhwriteSTD(:,:,3)=rhlevSTD(:,:)
3484	phiwriteSTD(:,:,3)=philevSTD(:,:)
3485	uwriteSTD(:,:,3)=ulevSTD(:,:)
3486	vwriteSTD(:,:,3)=vlevSTD(:,:)
3487	wwriteSTD(:,:,3)=wlevSTD(:,:)
3488
3489	twriteSTD(:,:,4)=tlevSTD(:,:)
3490	qwriteSTD(:,:,4)=qlevSTD(:,:)
3491	rhwriteSTD(:,:,4)=rhlevSTD(:,:)
3492	phiwriteSTD(:,:,4)=philevSTD(:,:)
3493	uwriteSTD(:,:,4)=ulevSTD(:,:)
3494	vwriteSTD(:,:,4)=vlevSTD(:,:)
3495	wwriteSTD(:,:,4)=wlevSTD(:,:)
3496	c
3497	cIM ajoute 5eme niveau 170310 BEG
3498	twriteSTD(:,:,5)=tlevSTD(:,:)
3499	qwriteSTD(:,:,5)=qlevSTD(:,:)
3500	rhwriteSTD(:,:,5)=rhlevSTD(:,:)
3501	phiwriteSTD(:,:,5)=philevSTD(:,:)
3502	uwriteSTD(:,:,5)=ulevSTD(:,:)
3503	vwriteSTD(:,:,5)=vlevSTD(:,:)
3504	wwriteSTD(:,:,5)=wlevSTD(:,:)
3505	cIM for NMC files
3506	DO n=1, nlevSTD3
3507	DO k=1, nlevSTD
3508	if(rlevSTD3(n).EQ.rlevSTD(k)) THEN
3509	twriteSTD3(:,n)=tlevSTD(:,k)
3510	qwriteSTD3(:,n)=qlevSTD(:,k)
3511	rhwriteSTD3(:,n)=rhlevSTD(:,k)
3512	phiwriteSTD3(:,n)=philevSTD(:,k)
3513	uwriteSTD3(:,n)=ulevSTD(:,k)
3514	vwriteSTD3(:,n)=vlevSTD(:,k)
3515	wwriteSTD3(:,n)=wlevSTD(:,k)
3516	endif !rlevSTD3(n).EQ.rlevSTD(k)
3517	ENDDO
3518	ENDDO
3519	c
3520	DO n=1, nlevSTD8
3521	DO k=1, nlevSTD
3522	if(rlevSTD8(n).EQ.rlevSTD(k)) THEN
3523	tnondefSTD8(:,n)=tnondef(:,k,2)
3524	twriteSTD8(:,n)=tsumSTD(:,k,2)
3525	qwriteSTD8(:,n)=qsumSTD(:,k,2)
3526	rhwriteSTD8(:,n)=rhsumSTD(:,k,2)
3527	phiwriteSTD8(:,n)=phisumSTD(:,k,2)
3528	uwriteSTD8(:,n)=usumSTD(:,k,2)
3529	vwriteSTD8(:,n)=vsumSTD(:,k,2)
3530	wwriteSTD8(:,n)=wsumSTD(:,k,2)
3531	endif !rlevSTD8(n).EQ.rlevSTD(k)
3532	ENDDO
3533	ENDDO
3534	c
3535	c slp sea level pressure
3536	slp(:) = paprs(:,1)exp(pphis(:)/(RDt_seri(:,1)))
3537	c
3538	ccc prw = eau precipitable
3539	DO i = 1, klon
3540	prw(i) = 0.
3541	DO k = 1, klev
3542	prw(i) = prw(i) +
3543	. q_seri(i,k)*(paprs(i,k)-paprs(i,k+1))/RG
3544	ENDDO
3545	ENDDO
3546	c
3547	cIM initialisation + calculs divers diag AMIP2
3548	c
3549	#include "calcul_divers.h"
3550	c
3551	IF (config_inca /= 'none') THEN
3552	#ifdef INCA
3553	CALL VTe(VTphysiq)
3554	CALL VTb(VTinca)
3555
3556	CALL chemhook_end (
3557	$ dtime,
3558	$ pplay,
3559	$ t_seri,
3560	$ tr_seri,
3561	$ nbtr,
3562	$ paprs,
3563	$ q_seri,
3564	$ airephy,
3565	$ pphi,
3566	$ pphis,
3567	$ zx_rh)
3568
3569	CALL VTe(VTinca)
3570	CALL VTb(VTphysiq)
3571	#endif
3572	END IF
3573
3574	c=============================================================
3575	c
3576	c Convertir les incrementations en tendances
3577	c
3578	if (mydebug) then
3579	call writefield_phy('u_seri',u_seri,llm)
3580	call writefield_phy('v_seri',v_seri,llm)
3581	call writefield_phy('t_seri',t_seri,llm)
3582	call writefield_phy('q_seri',q_seri,llm)
3583	endif
3584
3585	DO k = 1, klev
3586	DO i = 1, klon
3587	d_u(i,k) = ( u_seri(i,k) - u(i,k) ) / dtime
3588	d_v(i,k) = ( v_seri(i,k) - v(i,k) ) / dtime
3589	d_t(i,k) = ( t_seri(i,k)-t(i,k) ) / dtime
3590	d_qx(i,k,ivap) = ( q_seri(i,k) - qx(i,k,ivap) ) / dtime
3591	d_qx(i,k,iliq) = ( ql_seri(i,k) - qx(i,k,iliq) ) / dtime
3592	ENDDO
3593	ENDDO
3594	c
3595	IF (nqtot.GE.3) THEN
3596	DO iq = 3, nqtot
3597	DO k = 1, klev
3598	DO i = 1, klon
3599	d_qx(i,k,iq) = ( tr_seri(i,k,iq-2) - qx(i,k,iq) ) / dtime
3600	ENDDO
3601	ENDDO
3602	ENDDO
3603	ENDIF
3604	c
3605	cIM rajout diagnostiques bilan KP pour analyse MJO par Jun-Ichi Yano
3606	cIM global posePB#include "write_bilKP_ins.h"
3607	cIM global posePB#include "write_bilKP_ave.h"
3608	c
3609	c Sauvegarder les valeurs de t et q a la fin de la physique:
3610	c
3611	DO k = 1, klev
3612	DO i = 1, klon
3613	u_ancien(i,k) = u_seri(i,k)
3614	v_ancien(i,k) = v_seri(i,k)
3615	t_ancien(i,k) = t_seri(i,k)
3616	q_ancien(i,k) = q_seri(i,k)
3617	ENDDO
3618	ENDDO
3619	c
3620	!==========================================================================
3621	! Sorties des tendances pour un point particulier
3622	! a utiliser en 1D, avec igout=1 ou en 3D sur un point particulier
3623	! pour le debug
3624	! La valeur de igout est attribuee plus haut dans le programme
3625	!==========================================================================
3626
3627	if (prt_level.ge.1) then
3628	write(lunout,*) 'FIN DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
3629	write(lunout,*)
3630	s 'nlon,klev,nqtot,debut,lafin,jD_cur, jH_cur, pdtphys pct tlos'
3631	write(lunout,*)
3632	s nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur ,pdtphys,
3633	s pctsrf(igout,is_ter), pctsrf(igout,is_lic),pctsrf(igout,is_oce),
3634	s pctsrf(igout,is_sic)
3635	write(lunout,*) 'd_t_dyn,d_t_con,d_t_lsc,d_t_ajsb,d_t_ajs,d_t_eva'
3636	do k=1,klev
3637	write(lunout,*) d_t_dyn(igout,k),d_t_con(igout,k),
3638	s d_t_lsc(igout,k),d_t_ajsb(igout,k),d_t_ajs(igout,k),
3639	s d_t_eva(igout,k)
3640	enddo
3641	write(lunout,*) 'cool,heat'
3642	do k=1,klev
3643	write(lunout,*) cool(igout,k),heat(igout,k)
3644	enddo
3645
3646	write(lunout,*) 'd_t_oli,d_t_vdf,d_t_oro,d_t_lif,d_t_ec'
3647	do k=1,klev
3648	write(lunout,*) d_t_oli(igout,k),d_t_vdf(igout,k),
3649	s d_t_oro(igout,k),d_t_lif(igout,k),d_t_ec(igout,k)
3650	enddo
3651
3652	write(lunout,*) 'd_ps ',d_ps(igout)
3653	write(lunout,*) 'd_u, d_v, d_t, d_qx1, d_qx2 '
3654	do k=1,klev
3655	write(lunout,*) d_u(igout,k),d_v(igout,k),d_t(igout,k),
3656	s d_qx(igout,k,1),d_qx(igout,k,2)
3657	enddo
3658	endif
3659
3660	!==========================================================================
3661
3662	c============================================================
3663	c Calcul de la temperature potentielle
3664	c============================================================
3665	DO k = 1, klev
3666	DO i = 1, klon
3667	theta(i,k)=t(i,k)(100000./pplay(i,k))*(RD/RCPD)
3668	ENDDO
3669	ENDDO
3670	c
3671
3672	c 22.03.04 BEG
3673	c=============================================================
3674	c Ecriture des sorties
3675	c=============================================================
3676	#ifdef CPP_IOIPSL
3677
3678	c Recupere des varibles calcule dans differents modules
3679	c pour ecriture dans histxxx.nc
3680
3681	! Get some variables from module fonte_neige_mod
3682	CALL fonte_neige_get_vars(pctsrf,
3683	. zxfqcalving, zxfqfonte, zxffonte)
3684
3685
3686	#include "phys_output_write.h"
3687
3688	#ifdef histISCCP
3689	#include "write_histISCCP.h"
3690	#endif
3691
3692	#ifdef histNMC
3693	#include "write_histhfNMC.h"
3694	#include "write_histdayNMC.h"
3695	#include "write_histmthNMC.h"
3696	#endif
3697
3698	#include "write_histday_seri.h"
3699
3700	#include "write_paramLMDZ_phy.h"
3701
3702	#endif
3703
3704	c 22.03.04 END
3705	c
3706	c====================================================================
3707	c Si c'est la fin, il faut conserver l'etat de redemarrage
3708	c====================================================================
3709	c
3710
3711
3712	IF (lafin) THEN
3713	itau_phy = itau_phy + itap
3714	CALL phyredem ("restartphy.nc")
3715	! open(97,form="unformatted",file="finbin")
3716	! write(97) u_seri,v_seri,t_seri,q_seri
3717	! close(97)
3718	C$OMP MASTER
3719	if (read_climoz >= 1) then
3720	if (is_mpi_root) then
3721	call nf95_close(ncid_climoz)
3722	end if
3723	deallocate(press_climoz) ! pointer
3724	end if
3725	C$OMP END MASTER
3726	ENDIF
3727
3728	! first=.false.
3729
3730	RETURN
3731	END
3732	FUNCTION qcheck(klon,klev,paprs,q,ql,aire)
3733	IMPLICIT none
3734	c
3735	c Calculer et imprimer l'eau totale. A utiliser pour verifier
3736	c la conservation de l'eau
3737	c
3738	#include "YOMCST.h"
3739	INTEGER klon,klev
3740	REAL paprs(klon,klev+1), q(klon,klev), ql(klon,klev)
3741	REAL aire(klon)
3742	REAL qtotal, zx, qcheck
3743	INTEGER i, k
3744	c
3745	zx = 0.0
3746	DO i = 1, klon
3747	zx = zx + aire(i)
3748	ENDDO
3749	qtotal = 0.0
3750	DO k = 1, klev
3751	DO i = 1, klon
3752	qtotal = qtotal + (q(i,k)+ql(i,k)) * aire(i)
3753	. *(paprs(i,k)-paprs(i,k+1))/RG
3754	ENDDO
3755	ENDDO
3756	c
3757	qcheck = qtotal/zx
3758	c
3759	RETURN
3760	END
3761	SUBROUTINE gr_fi_ecrit(nfield,nlon,iim,jjmp1,fi,ecrit)
3762	IMPLICIT none
3763	c
3764	c Tranformer une variable de la grille physique a
3765	c la grille d'ecriture
3766	c
3767	INTEGER nfield,nlon,iim,jjmp1, jjm
3768	REAL fi(nlon,nfield), ecrit(iim*jjmp1,nfield)
3769	c
3770	INTEGER i, n, ig
3771	c
3772	jjm = jjmp1 - 1
3773	DO n = 1, nfield
3774	DO i=1,iim
3775	ecrit(i,n) = fi(1,n)
3776	ecrit(i+jjm*iim,n) = fi(nlon,n)
3777	ENDDO
3778	DO ig = 1, nlon - 2
3779	ecrit(iim+ig,n) = fi(1+ig,n)
3780	ENDDO
3781	ENDDO
3782	RETURN
3783	END

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format