Context Navigation

physiq_mod.F90 @ 2518

Last change on this file since 2518 was 2517, checked in by oboucher, 9 years ago
Adding runofflic diagnostic in outputs This is the freshwater flux going into the coupler
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: 147.6 KB

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

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

Download in other formats:

Original Format