Context Navigation

physiq.F90 @ 1999

Last change on this file since 1999 was 1999, checked in by Laurent Fairhead, 10 years ago
Merged trunk changes r1920:1997 into testing branch
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 121.2 KB

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

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

Download in other formats:

Original Format