Context Navigation

physiq.F90 @ 1973

Last change on this file since 1973 was 1973, checked in by fhourdin, 11 years ago
Le add_pbl_tend est passé de phy1d à phylmd et appelé dans phylmd/physiq.F90 si klon_glo=1 add_pbl_tend is moved from phy1d to phylmd and called in phylmd/physiq.F90 if klon_glo=1
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: 120.8 KB

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

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

Download in other formats:

Original Format