Context Navigation

physiq.F @ 1398

Last change on this file since 1398 was 1398, checked in by musat, 14 years ago

Last corrections for CMIP5:

Add O3 at standard level files histmthNMC.nc
Add positive attribute "down" for vertical axes for all output files
Replace "inst" by "ave" for hist*NMC.nc files to have time_counter and bounds for time axis (Marie-Alice's hint)
Correct units for vertical axes : mb instead of hPa
Add mass flux at the bottom of clouds
Comment non initialized variables (s_capCL, s_oliqCL, s_cteiCL, s_trmb1, s_trmb2, s_trmb3) for the output files
Geopotential field phy850, phi700, phi500, etc are modified to "geopotential height and are called z850, z700, z500, etc
Meaning of specific humidity outputs - ovapinit and ovap - were interchanged
Fields albs, albslw become alb1, alb2 in output files
Correct title for rugs_* fields
Correct units for pbase and ptop are Pa (not mb)
Correct ndayrain field

FH/JLD/JYG/MAF/IM

Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 120.3 KB

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

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

Download in other formats:

Original Format