Context Navigation

physiq.F @ 1368

Last change on this file since 1368 was 1368, checked in by idelkadi, 14 years ago
Rajout des cles logiques pour activer les sorties COSP par frequence Reglages des niveaux de sorties pour les variables CFMIP
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 120.3 KB

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

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

Download in other formats:

Original Format