Context Navigation

physiq.F @ 1356

Last change on this file since 1356 was 1356, checked in by musat, 14 years ago
use of phys_cal_mod in conf_phys (and change call order in physiq.F) to automatically calculate the output frequency of standard pressure monthly file (histmthNMC.nc) IM
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 122.5 KB

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

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

Download in other formats:

Original Format