Context Navigation

physiq_mod.F @ 3567

Last change on this file since 3567 was 3535, checked in by mlefevre, 2 weeks ago
Addition of vertical eddy coeficient profile for the 1D model. ML
File size: 74.4 KB

Line
1	!
2	! $Id: $
3	!
4	MODULE physiq_mod
5
6	IMPLICIT NONE
7
8	CONTAINS
9
10	SUBROUTINE physiq (nlon,nlev,nqmax,
11	. debut,lafin,rjourvrai,gmtime,pdtphys,
12	. paprs,pplay,ppk,pphi,pphis,presnivs,
13	. u,v,t,qx,
14	. flxmw,
15	. d_u, d_v, d_t, d_qx, d_ps)
16
17	c======================================================================
18	c
19	c Modifications pour la physique de Venus
20	c S. Lebonnois (LMD/CNRS) Septembre 2005
21	c
22	c ---------------------------------------------------------------------
23	c Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
24	c
25	c Objet: Moniteur general de la physique du modele
26	cAA Modifications quant aux traceurs :
27	cAA - uniformisation des parametrisations ds phytrac
28	cAA - stockage des moyennes des champs necessaires
29	cAA en mode traceur off-line
30	c modif ( P. Le Van , 12/10/98 )
31	c
32	c Arguments:
33	c
34	c nlon----input-I-nombre de points horizontaux
35	c nlev----input-I-nombre de couches verticales
36	c nqmax---input-I-nombre de traceurs
37	c debut---input-L-variable logique indiquant le premier passage
38	c lafin---input-L-variable logique indiquant le dernier passage
39	c rjour---input-R-numero du jour de l'experience
40	c gmtime--input-R-fraction de la journee (0 a 1)
41	c pdtphys-input-R-pas d'integration pour la physique (seconde)
42	c paprs---input-R-pression pour chaque inter-couche (en Pa)
43	c pplay---input-R-pression pour le mileu de chaque couche (en Pa)
44	c ppk ---input-R-fonction d'Exner au milieu de couche
45	c pphi----input-R-geopotentiel de chaque couche (g z) (reference sol)
46	c pphis---input-R-geopotentiel du sol
47	c presnivs-input_R_pressions approximat. des milieux couches ( en PA)
48	c u-------input-R-vitesse dans la direction X (de O a E) en m/s
49	c v-------input-R-vitesse Y (de S a N) en m/s
50	c t-------input-R-temperature (K)
51	c qx------input-R-mass mixing ratio traceurs (kg/kg)
52	c d_t_dyn-input-R-tendance dynamique pour "t" (K/s)
53	c flxmw---input-R-flux de masse vertical en kg/s
54	c
55	c d_u-----output-R-tendance physique de "u" (m/s/s)
56	c d_v-----output-R-tendance physique de "v" (m/s/s)
57	c d_t-----output-R-tendance physique de "t" (K/s)
58	c d_qx----output-R-tendance physique de "qx" (kg/kg/s)
59	c d_ps----output-R-tendance physique de la pression au sol
60	c======================================================================
61	USE ioipsl
62	! USE histcom ! not needed; histcom is included in ioipsl
63	use dimphy
64	USE geometry_mod,only: longitude, latitude, ! in radians
65	& longitude_deg,latitude_deg, ! in degrees
66	& cell_area,dx,dy
67	USE phys_state_var_mod ! Variables sauvegardees de la physique
68	USE cpdet_phy_mod, only: cpdet, t2tpot
69	USE chemparam_mod
70	USE age_of_air_mod, ONLY: ok_aoa, reinit_aoa, i_aoa, init_age
71	USE age_of_air_mod, ONLY: aoa_ini, age_of_air
72	USE conc
73	USE param_v4_h
74	USE compo_hedin83_mod2
75	use radlwsw_newtoncool_mod, only: radlwsw_newtoncool
76	! use ieee_arithmetic
77	use time_phylmdz_mod, only: annee_ref, day_ref, itau_phy
78	use mod_grid_phy_lmdz, only: nbp_lon
79	use infotrac_phy, only: iflag_trac, tname, ttext
80	use vertical_layers_mod, only: pseudoalt
81	use turb_mod, only : sens, turb_resolved
82	use nonoro_gwd_ran_mod, only: nonoro_gwd_ran
83	use sed_and_prod_mad, only: aer_sedimentation, drop_sedimentation
84	use iono_h, only: temp_elect, temp_ion
85	#ifdef CPP_XIOS
86	use xios_output_mod, only: initialize_xios_output,
87	& update_xios_timestep,
88	& send_xios_field
89	use wxios, only: wxios_context_init, xios_context_finalize
90	#endif
91	#ifdef MESOSCALE
92	use comm_wrf
93	#else
94	use iophy
95	use write_field_phy
96	use mod_phys_lmdz_omp_data, ONLY: is_omp_master
97	USE mod_phys_lmdz_para, only : is_parallel,jj_nb,
98	& is_north_pole_phy,
99	& is_south_pole_phy,
100	& is_master
101	#endif
102	IMPLICIT none
103	c======================================================================
104	c CLEFS CPP POUR LES IO
105	c =====================
106	#ifndef MESOSCALE
107	c#define histhf
108	#define histday
109	#define histmth
110	#define histins
111	#endif
112	c======================================================================
113	#include "dimsoil.h"
114	#include "clesphys.h"
115	#include "iniprint.h"
116	#include "timerad.h"
117	#include "tabcontrol.h"
118	#include "nirdata.h"
119	#include "hedin.h"
120	c======================================================================
121	LOGICAL ok_journe ! sortir le fichier journalier
122	save ok_journe
123	c PARAMETER (ok_journe=.true.)
124	c
125	LOGICAL ok_mensuel ! sortir le fichier mensuel
126	save ok_mensuel
127	c PARAMETER (ok_mensuel=.true.)
128	c
129	LOGICAL ok_instan ! sortir le fichier instantane
130	save ok_instan
131	c PARAMETER (ok_instan=.true.)
132	c
133	c======================================================================
134	c
135	c Variables argument:
136	c
137	INTEGER nlon
138	INTEGER nlev
139	INTEGER nqmax
140	REAL rjourvrai
141	REAL gmtime
142	REAL pdtphys
143	LOGICAL debut, lafin
144	REAL paprs(klon,klev+1)
145	REAL pplay(klon,klev)
146	REAL pphi(klon,klev)
147	REAL pphis(klon)
148	REAL presnivs(klev)
149
150	! ADAPTATION GCM POUR CP(T)
151	REAL ppk(klon,klev)
152
153	REAL u(klon,klev)
154	REAL v(klon,klev)
155	REAL t(klon,klev)
156	REAL qx(klon,klev,nqmax)
157
158	REAL d_u_dyn(klon,klev)
159	REAL d_t_dyn(klon,klev)
160
161	REAL flxmw(klon,klev)
162
163	REAL d_u(klon,klev)
164	REAL d_v(klon,klev)
165	REAL d_t(klon,klev)
166	REAL d_qx(klon,klev,nqmax)
167	REAL d_ps(klon)
168
169	logical ok_hf
170	real ecrit_hf
171	integer nid_hf
172	save ok_hf, ecrit_hf, nid_hf
173
174	#ifdef histhf
175	data ok_hf,ecrit_hf/.true.,0.25/
176	#else
177	data ok_hf/.false./
178	#endif
179
180	c Variables propres a la physique
181
182	integer,save :: itap ! physics counter
183	REAL delp(klon,klev) ! epaisseur d'une couche
184	REAL omega(klon,klev) ! vitesse verticale en Pa/s (+ downward)
185	REAL vertwind(klon,klev) ! vitesse verticale en m/s (+ upward)
186
187	INTEGER igwd,idx(klon),itest(klon)
188
189	c Diagnostiques 2D de drag_noro, lift_noro et gw_nonoro
190
191	REAL zulow(klon),zvlow(klon)
192	REAL zustrdr(klon), zvstrdr(klon)
193	REAL zustrli(klon), zvstrli(klon)
194	REAL zustrhi(klon), zvstrhi(klon)
195	REAL zublstrdr(klon), zvblstrdr(klon)
196	REAL znlow(klon), zeff(klon)
197	REAL zbl(klon), knu2(klon),kbreak(nlon)
198	REAL tau0(klon), ztau(klon,klev)
199
200	c Pour calcul GW drag oro et nonoro: CALCUL de N2:
201	real zdtlev(klon,klev),zdzlev(klon,klev)
202	real ztlev(klon,klev)
203	real zn2(klon,klev) ! BV^2 at plev
204
205	c Pour les bilans de moment angulaire,
206	integer bilansmc
207	c Pour le transport de ballons
208	integer ballons
209	c j'ai aussi besoin
210	c du stress de couche limite a la surface:
211
212	REAL zustrcl(klon),zvstrcl(klon)
213
214	c et du stress total c de la physique:
215
216	REAL zustrph(klon),zvstrph(klon)
217
218	c Variables locales:
219	c
220	REAL cdragh(klon) ! drag coefficient pour T and Q
221	REAL cdragm(klon) ! drag coefficient pour vent
222	c
223	cAA Pour TRACEURS
224	cAA
225	REAL,save,allocatable :: source(:,:)
226	REAL ycoefh(klon,klev) ! coef d'echange pour phytrac
227	REAL :: kzz_p(klev) ! coef d'echange pour phytrac pour le 1D
228	REAL yu1(klon) ! vents dans la premiere couche U
229	REAL yv1(klon) ! vents dans la premiere couche V
230
231	REAL dsens(klon) ! derivee chaleur sensible
232	REAL ve(klon) ! integr. verticale du transport meri. de l'energie
233	REAL vq(klon) ! integr. verticale du transport meri. de l'eau
234	REAL ue(klon) ! integr. verticale du transport zonal de l'energie
235	REAL uq(klon) ! integr. verticale du transport zonal de l'eau
236	c
237	REAL Fsedim(klon,klev+1) ! Flux de sedimentation (kg.m-2)
238
239	c======================================================================
240	c
241	c Declaration des procedures appelees
242	c
243	EXTERNAL ajsec ! ajustement sec
244	EXTERNAL clmain ! couche limite
245	EXTERNAL clmain_ideal ! couche limite simple
246	EXTERNAL hgardfou ! verifier les temperatures
247	c EXTERNAL orbite ! calculer l'orbite
248	EXTERNAL phyetat0 ! lire l'etat initial de la physique
249	EXTERNAL phyredem ! ecrire l'etat de redemarrage de la physique
250	EXTERNAL radlwsw ! rayonnements solaire et infrarouge
251	! EXTERNAL suphec ! initialiser certaines constantes
252	EXTERNAL transp ! transport total de l'eau et de l'energie
253	EXTERNAL printflag
254	EXTERNAL zenang
255	EXTERNAL diagetpq
256	EXTERNAL conf_phys
257	EXTERNAL diagphy
258	EXTERNAL mucorr
259	EXTERNAL nirco2abs
260	EXTERNAL nir_leedat
261	EXTERNAL nltecool
262	EXTERNAL nlte_tcool
263	EXTERNAL nlte_setup
264	EXTERNAL blendrad
265	EXTERNAL nlthermeq
266	EXTERNAL euvheat
267	EXTERNAL param_read_e107
268	EXTERNAL conduction
269	EXTERNAL molvis
270	EXTERNAL moldiff_red
271
272	c
273	c Variables locales
274	c
275	CXXX PB
276	REAL fluxt(klon,klev) ! flux turbulent de chaleur
277	REAL fluxu(klon,klev) ! flux turbulent de vitesse u
278	REAL fluxv(klon,klev) ! flux turbulent de vitesse v
279	c
280	REAL flux_dyn(klon,klev) ! flux de chaleur produit par la dynamique
281	REAL flux_ajs(klon,klev) ! flux de chaleur ajustement sec
282	REAL flux_ec(klon,klev) ! flux de chaleur Ec
283	c
284	REAL tmpout(klon,klev) ! [K/s]
285	c
286	REAL dist, rmu0(klon), fract(klon)
287	REAL zdtime, zlongi
288	c
289	INTEGER i, k, iq, ig, j, ll, ilon, ilat, ilev, isoil
290	c
291	REAL zphi(klon,klev)
292	REAL zzlev(klon,klev+1),zzlay(klon,klev),z1,z2
293	REAL tlaymean ! valeur temporaire pour calculer zzlay
294	real tsurf(klon)
295
296	c va avec nlte_tcool
297	INTEGER ierr_nlte
298	REAL varerr
299
300	! photochemistry
301
302	integer :: chempas
303	real :: zctime
304
305	! sedimentation
306
307	REAL :: m0_mode1(klev,2),m0_mode2(klev,2)
308	REAL :: m3_mode1(klev,3),m3_mode2 (klev,3)
309	REAL :: d_drop_sed(klev),d_ccn_sed(klev,2),d_liq_sed(klev,2)
310	REAL :: aer_flux(klev)
311	c
312	c Variables du changement
313	c
314	c ajs: ajustement sec
315	c vdf: couche limite (Vertical DiFfusion)
316	REAL d_t_ajs(klon,klev), d_tr_ajs(klon,klev,nqmax)
317	REAL d_u_ajs(klon,klev), d_v_ajs(klon,klev)
318	c
319	REAL d_ts(klon)
320	c
321	REAL d_u_vdf(klon,klev), d_v_vdf(klon,klev)
322	REAL d_t_vdf(klon,klev), d_tr_vdf(klon,klev,nqmax)
323	c
324	CMOD LOTT: Tendances Orography Sous-maille
325	REAL d_u_oro(klon,klev), d_v_oro(klon,klev)
326	REAL d_t_oro(klon,klev)
327	REAL d_u_lif(klon,klev), d_v_lif(klon,klev)
328	REAL d_t_lif(klon,klev)
329	C Tendances Ondes de G non oro (runs strato).
330	REAL d_u_hin(klon,klev), d_v_hin(klon,klev)
331	REAL d_t_hin(klon,klev)
332
333	c Tendencies due to radiative scheme [K/s]
334	c d_t_rad,dtsw,dtlw,d_t_nirco2,d_t_nlte,d_t_euv
335	c are not computed at each physical timestep
336	c therefore, they are defined and saved in phys_state_var_mod
337
338	c Tendencies due to molecular viscosity and conduction
339	real d_t_conduc(klon,klev) ! [K/s]
340	real d_u_molvis(klon,klev) ! [m/s] /s
341	real d_v_molvis(klon,klev) ! [m/s] /s
342
343	c Tendencies due to molecular diffusion
344	real d_q_moldif(klon,klev,nqmax)
345
346	c Tendencies due to ambipolar ion diffusion
347	real d_q_iondif(klon,klev,nqmax)
348
349	c
350	c Variables liees a l'ecriture de la bande histoire physique
351	c
352	INTEGER ecrit_mth
353	SAVE ecrit_mth ! frequence d'ecriture (fichier mensuel)
354	c
355	INTEGER ecrit_day
356	SAVE ecrit_day ! frequence d'ecriture (fichier journalier)
357	c
358	INTEGER ecrit_ins
359	SAVE ecrit_ins ! frequence d'ecriture (fichier instantane)
360	c
361	integer itau_w ! pas de temps ecriture = itap + itau_phy
362
363	c Variables locales pour effectuer les appels en serie
364	c
365	REAL t_seri(klon,klev)
366	REAL u_seri(klon,klev), v_seri(klon,klev)
367	c
368	REAL :: tr_seri(klon,klev,nqmax)
369	REAL :: tr_hedin(klon,klev,nqmax)
370	REAL :: d_tr(klon,klev,nqmax)
371	c pour sorties
372	REAL :: col_dens_tr(klon,nqmax)
373	REAL,allocatable,save :: prod_tr(:,:,:)
374	REAL,allocatable,save :: loss_tr(:,:,:)
375
376	c pour ioipsl
377	INTEGER nid_day, nid_mth, nid_ins
378	SAVE nid_day, nid_mth, nid_ins
379	INTEGER nhori, nvert, idayref
380	REAL zsto, zout, zsto1, zsto2, zero
381	parameter (zero=0.0e0)
382	real zjulian
383	save zjulian
384
385	CHARACTER*2 str2
386	character*20 modname
387	character*80 abort_message
388	logical ok_sync
389
390	character*30 nom_fichier
391	character*10 varname
392	character*40 vartitle
393	character*20 varunits
394	C Variables liees au bilan d'energie et d'enthalpi
395	REAL ztsol(klon)
396	REAL h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
397	$ , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
398	SAVE h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
399	$ , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
400	REAL d_h_vcol, d_h_dair, d_qt, d_qw, d_ql, d_qs, d_ec
401	REAL d_h_vcol_phy
402	REAL fs_bound, fq_bound
403	SAVE d_h_vcol_phy
404	REAL zero_v(klon),zero_v2(klon,klev)
405	CHARACTER*15 ztit
406	INTEGER ip_ebil ! PRINT level for energy conserv. diag.
407	SAVE ip_ebil
408	DATA ip_ebil/2/
409	INTEGER if_ebil ! level for energy conserv. dignostics
410	SAVE if_ebil
411	c+jld ec_conser
412	REAL d_t_ec(klon,klev) ! tendance du a la conversion Ec -> E thermique
413	c-jld ec_conser
414
415	c ALBEDO VARIATIONS (VCD)
416	REAL factAlb
417	c TEST VENUS...
418	REAL mang(klon,klev) ! moment cinetique
419	REAL mangtot ! moment cinetique total
420
421	c cell_area for outputs in hist*
422	REAL cell_area_out(klon)
423	#ifdef MESOSCALE
424	REAL :: dt_dyn(klev)
425	#endif
426	c Declaration des constantes et des fonctions thermodynamiques
427	c
428	#include "YOMCST.h"
429
430	c======================================================================
431	c======================================================================
432	c INITIALISATIONS
433	c======================================================================
434
435	modname = 'physiq'
436	ok_sync=.TRUE.
437
438	bilansmc = 0
439	ballons = 0
440	! NE FONCTIONNENT PAS ENCORE EN PARALLELE !!!
441	#ifndef MESOSCALE
442	if (is_parallel) then
443	bilansmc = 0
444	ballons = 0
445	endif
446	#endif
447	IF (if_ebil.ge.1) THEN
448	DO i=1,klon
449	zero_v(i)=0.
450	END DO
451	DO i=1,klon
452	DO j=1,klev
453	zero_v2(i,j)=0.
454	END DO
455	END DO
456	END IF
457
458	c======================================================================
459	c DONE ONLY AT FIRST CALL
460	c========================
461	IF (debut) THEN
462	allocate(source(klon,nqmax))
463	allocate(prod_tr(klon,klev,nqmax))
464	allocate(loss_tr(klon,klev,nqmax))
465	allocate(no_emission(klon,klev))
466	allocate(o2_emission(klon,klev))
467
468	#ifdef CPP_XIOS
469	! Initialize XIOS context
470	write(,) "physiq: call wxios_context_init"
471	CALL wxios_context_init
472	#endif
473
474	! The call to suphec is now done in iniphysiq_mod (interface)
475	! CALL suphec ! initialiser constantes et parametres phys.
476
477	IF (if_ebil.ge.1) d_h_vcol_phy=0.
478	!
479	! load flag and parameter values from physiq.def
480	!
481	call conf_phys(ok_journe, ok_mensuel,
482	. ok_instan,
483	. if_ebil)
484
485	call phys_state_var_init(nqmax)
486	c
487	c Initialising Hedin model for upper atm
488	c (to be revised when coupled to chemistry) :
489	call conc_init
490
491	! initialise physics counter
492
493	itap = 0
494
495	#ifdef MESOSCALE
496	print*,'check pdtphys',pdtphys
497	PRINT*,'check phisfi ',pphis(1),pphis(klon)
498	PRINT*,'check geop',pphi(1,1),pphi(klon,klev)
499	PRINT*,'check radsol',radsol(1),radsol(klon)
500	print*,'check ppk',ppk(1,1),ppk(klon,klev)
501	print*,'check ftsoil',ftsoil(1,1),ftsoil(klon,nsoilmx)
502	print*,'check ftsol',ftsol(1),ftsol(klon)
503	print*, "check temp", t(1,1),t(klon,klev)
504	print*, "check pres",paprs(1,1),paprs(klon,klev),pplay(1,1),
505	. pplay(klon,klev)
506	print*, "check u", u(1,1),u(klon,klev)
507	print*, "check v", v(1,1),v(klon,klev)
508	print*,'check falbe',falbe(1),falbe(klon)
509	!nqtot=nqmax
510	!ALLOCATE(tname(nqtot))
511	!tname=noms
512	zmea=0.
513	zstd=0.
514	zsig=0.
515	zgam=0.
516	zthe=0.
517	dtime=pdtphys
518	#else
519	c
520	c Lecture startphy.nc :
521	c
522	CALL phyetat0 ("startphy.nc")
523	IF (.not.startphy_file) THEN
524	! Additionnal academic initializations
525	ftsol(:)=t(:,1) ! surface temperature as in first atm. layer
526	DO isoil=1, nsoilmx
527	! subsurface temperatures equal to surface temperature
528	ftsoil(:,isoil)=ftsol(:)
529	ENDDO
530	ENDIF
531	#endif
532
533	c dtime est defini dans tabcontrol.h et lu dans startphy
534	c pdtphys est calcule a partir des nouvelles conditions:
535	c Reinitialisation du pas de temps physique quand changement
536	IF (ABS(dtime-pdtphys).GT.0.001) THEN
537	WRITE(lunout,*) 'Pas physique a change',dtime,
538	. pdtphys
539	c abort_message='Pas physique n est pas correct '
540	c call abort_physic(modname,abort_message,1)
541	c----------------
542	c pour initialiser convenablement le time_counter, il faut tenir compte
543	c du changement de dtime en changeant itau_phy (point de depart)
544	itau_phy = NINT(itau_phy*dtime/pdtphys)
545	c----------------
546	dtime=pdtphys
547	ENDIF
548
549	radpas = NINT(RDAY/pdtphys/nbapp_rad)
550
551	CALL printflag( ok_journe,ok_instan )
552
553	#ifdef CPP_XIOS
554	write(,) "physiq: call initialize_xios_output"
555	call initialize_xios_output(rjourvrai,gmtime,pdtphys,RDAY,
556	& presnivs,pseudoalt)
557	#endif
558
559	c
560	c---------
561	c FLOTT
562	IF (ok_orodr) THEN
563	DO i=1,klon
564	rugoro(i) = MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
565	ENDDO
566	CALL SUGWD(klon,klev,paprs,pplay)
567	DO i=1,klon
568	zuthe(i)=0.
569	zvthe(i)=0.
570	if(zstd(i).gt.10.)then
571	zuthe(i)=(1.-zgam(i))*cos(zthe(i))
572	zvthe(i)=(1.-zgam(i))*sin(zthe(i))
573	endif
574	ENDDO
575	ENDIF
576
577	if (bilansmc.eq.1) then
578	C OUVERTURE D'UN FICHIER FORMATTE POUR STOCKER LES COMPOSANTES
579	C DU BILAN DE MOMENT ANGULAIRE.
580	open(27,file='aaam_bud.out',form='formatted')
581	open(28,file='fields_2d.out',form='formatted')
582	write(,)'Ouverture de aaam_bud.out (FL Vous parle)'
583	write(,)'Ouverture de fields_2d.out (FL Vous parle)'
584	endif !bilansmc
585
586	c--------------SLEBONNOIS
587	C OUVERTURE DES FICHIERS FORMATTES CONTENANT LES POSITIONS ET VITESSES
588	C DES BALLONS
589	if (ballons.eq.1) then
590	open(30,file='ballons-lat.out',form='formatted')
591	open(31,file='ballons-lon.out',form='formatted')
592	open(32,file='ballons-u.out',form='formatted')
593	open(33,file='ballons-v.out',form='formatted')
594	open(34,file='ballons-alt.out',form='formatted')
595	write(,)'Ouverture des ballons*.out'
596	endif !ballons
597	c-------------
598
599	c---------
600	C TRACEURS
601	C source dans couche limite
602	source(:,:) = 0.0 ! pas de source, pour l'instant
603	c---------
604
605	c---------
606	c INITIALIZE THERMOSPHERIC PARAMETERS
607	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
608
609	if (callthermos) then
610	call fill_data_thermos
611	call allocate_param_thermos(klev)
612	call param_read_e107
613	endif
614
615	c Initialisation (recomputed in concentration2)
616	do ig=1,klon
617	do j=1,klev
618	rnew(ig,j)=RD
619	cpnew(ig,j)=cpdet(t(ig,j))
620	mmean(ig,j)=RMD
621	akknew(ig,j)=1.e-4
622	rho(ig,j)=pplay(ig,j)/(rnew(ig,j)*t(ig,j))
623	enddo
624
625	enddo
626
627	IF(callthermos.or.callnlte.or.callnirco2) THEN
628	call compo_hedin83_init2
629	ENDIF
630	if (callnlte.and.nltemodel.eq.2) call nlte_setup
631	if (callnirco2.and.nircorr.eq.1) call nir_leedat
632	c---------
633
634	c
635	c Verifications:
636	c
637	IF (nlon .NE. klon) THEN
638	WRITE(lunout,*)'nlon et klon ne sont pas coherents', nlon,
639	. klon
640	abort_message='nlon et klon ne sont pas coherents'
641	call abort_physic(modname,abort_message,1)
642	ENDIF
643	IF (nlev .NE. klev) THEN
644	WRITE(lunout,*)'nlev et klev ne sont pas coherents', nlev,
645	. klev
646	abort_message='nlev et klev ne sont pas coherents'
647	call abort_physic(modname,abort_message,1)
648	ENDIF
649	c
650	IF (dtimeREAL(radpas).GT.(RDAY0.25).AND.cycle_diurne)
651	$ THEN
652	WRITE(lunout,*)'Nbre d appels au rayonnement insuffisant'
653	WRITE(lunout,*)"Au minimum 4 appels par jour si cycle diurne"
654	abort_message='Nbre d appels au rayonnement insuffisant'
655	call abort_physic(modname,abort_message,1)
656	ENDIF
657	c
658	WRITE(lunout,*)"Clef pour la convection seche, iflag_ajs=",
659	. iflag_ajs
660	c
661	ecrit_mth = NINT(RDAY/dtime*ecriphy) ! tous les ecritphy jours
662	IF (ok_mensuel) THEN
663	WRITE(lunout,*)'La frequence de sortie mensuelle est de ',
664	. ecrit_mth
665	ENDIF
666
667	ecrit_day = NINT(RDAY/dtime *1.0) ! tous les jours
668	IF (ok_journe) THEN
669	WRITE(lunout,*)'La frequence de sortie journaliere est de ',
670	. ecrit_day
671	ENDIF
672
673	ecrit_ins = NINT(RDAY/dtime*ecriphy) ! Fraction de jour reglable
674	IF (ok_instan) THEN
675	WRITE(lunout,*)'La frequence de sortie instant. est de ',
676	. ecrit_ins
677	ENDIF
678
679	c Verification synchronize AMBIPOLAR DIFFUSION & CHEMISTRY
680
681	IF ((ok_iondiff) .and. (NINT(RDAY/dtime).ne.nbapp_chem)) THEN
682	WRITE(lunout,*)'nbapp_chem .NE. day_step'
683	WRITE(lunout,*)'nbapp_chem ', nbapp_chem
684	WRITE(lunout,*)'day_step ', NINT(RDAY/dtime)
685	WRITE(lunout,*)'nbapp_chem must be equal to day_step'
686	STOP
687	ENDIF
688
689	c Initialisation des sorties
690	c===========================
691
692	#ifdef CPP_IOIPSL
693
694	#ifdef histhf
695	#include "ini_histhf.h"
696	#endif
697
698	#ifdef histday
699	#include "ini_histday.h"
700	#endif
701
702	#ifdef histmth
703	#include "ini_histmth.h"
704	#endif
705
706	#ifdef histins
707	#include "ini_histins.h"
708	#endif
709
710	#endif
711
712	c
713	c Initialiser les valeurs de u pour calculs tendances
714	c (pour T, c'est fait dans phyetat0)
715	c
716	DO k = 1, klev
717	DO i = 1, klon
718	u_ancien(i,k) = u(i,k)
719	ENDDO
720	ENDDO
721
722	! initialisation of microphysical and chemical parameters
723
724	if (ok_chem .and. .not. ok_cloud) then
725	print*,"chemistry requires clouds"
726	print*,"ok_cloud must be .true."
727	stop
728	end if
729
730	if (.not. ok_chem .and. ok_cloud .and. (cl_scheme == 1)) then
731	print*,"cl_scheme = 1 requires chemistry"
732	print*,"ok_chem must be .true."
733	stop
734	end if
735
736	! number of microphysical tracers
737
738	nmicro = 0
739	if (ok_cloud .and. (cl_scheme == 1)) nmicro = 2
740	if (ok_cloud .and. (cl_scheme == 2)) nmicro = 12
741
742	! initialise chemical parameters. includes the indexation of microphys tracers
743
744	if (ok_chem .or. cl_scheme == 2) then
745	call chemparam_ini()
746	end if
747
748	if ((iflag_trac.eq.1) .and. ok_aoa) then
749	write(lunout,*) 'Initialising age of air subroutine'
750	allocate(init_age(klon,klev))
751	call aoa_ini(ok_chem, tr_scheme) ! get index of age of air tracer in the tracer array
752	if (reinit_aoa) then
753	age(:,:) = 0. ! Set initial value of age of air to 0 everywhere
754	tr_seri(:,:,i_aoa) = 1e-30 ! Set initial value of tracer to tiny everywhere
755	endif ! reinitialisation loop
756	init_age(:,:) = age(:,:) ! save initial value of age, either read in from restartphy or set to 0
757	endif ! age of air initialisation
758
759	! initialise cloud parameters (for cl_scheme = 1)
760
761	if (ok_cloud .and. (cl_scheme == 1)) then
762	call cloud_ini(nlon,nlev,nb_mode)
763	end if
764
765	! initialise mmean
766
767	if (callthermos) then
768	call concentrations2(pplay,t,qx,nqmax)
769	end if
770
771	c========================
772	ENDIF ! debut
773	c========================
774
775	c======================================================================
776	! ------------------------------------------------------
777	! Initializations done at every physical timestep:
778	! ------------------------------------------------------
779
780	c Mettre a zero des variables de sortie (pour securite)
781	c
782	DO i = 1, klon
783	d_ps(i) = 0.0
784	ENDDO
785	DO k = 1, klev
786	DO i = 1, klon
787	d_t(i,k) = 0.0
788	d_u(i,k) = 0.0
789	d_v(i,k) = 0.0
790	ENDDO
791	ENDDO
792	DO iq = 1, nqmax
793	DO k = 1, klev
794	DO i = 1, klon
795	d_qx(i,k,iq) = 0.0
796	ENDDO
797	ENDDO
798	ENDDO
799	c
800	c Ne pas affecter les valeurs entrees de u, v, h, et q
801	c
802	DO k = 1, klev
803	DO i = 1, klon
804	t_seri(i,k) = t(i,k)
805	u_seri(i,k) = u(i,k)
806	v_seri(i,k) = v(i,k)
807	ENDDO
808	ENDDO
809	DO iq = 1, nqmax
810	DO k = 1, klev
811	DO i = 1, klon
812	tr_seri(i,k,iq) = qx(i,k,iq)
813	ENDDO
814	ENDDO
815	ENDDO
816	C
817	DO i = 1, klon
818	ztsol(i) = ftsol(i)
819	ENDDO
820	C
821	IF (if_ebil.ge.1) THEN
822	ztit='after dynamic'
823	CALL diagetpq(cell_area,ztit,ip_ebil,1,1,dtime
824	e , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
825	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
826	C Comme les tendances de la physique sont ajoute dans la dynamique,
827	C on devrait avoir que la variation d'entalpie par la dynamique
828	C est egale a la variation de la physique au pas de temps precedent.
829	C Donc la somme de ces 2 variations devrait etre nulle.
830	call diagphy(cell_area,ztit,ip_ebil
831	e , zero_v, zero_v, zero_v, zero_v, zero_v
832	e , zero_v, zero_v, zero_v, ztsol
833	e , d_h_vcol+d_h_vcol_phy, d_qt, 0.
834	s , fs_bound, fq_bound )
835	END IF
836
837	c====================================================================
838	c XIOS outputs
839
840	#ifdef CPP_XIOS
841	! update XIOS time/calendar
842	call update_xios_timestep
843	#endif
844
845	c====================================================================
846	c Diagnostiquer la tendance dynamique
847	c
848	IF (ancien_ok) THEN
849	DO k = 1, klev
850	DO i = 1, klon
851	d_u_dyn(i,k) = (u_seri(i,k)-u_ancien(i,k))/dtime
852	d_t_dyn(i,k) = (t_seri(i,k)-t_ancien(i,k))/dtime
853	ENDDO
854	ENDDO
855
856	! ADAPTATION GCM POUR CP(T)
857	do i=1,klon
858	flux_dyn(i,1) = 0.0
859	do j=2,klev
860	flux_dyn(i,j) = flux_dyn(i,j-1)
861	. +cpnew(i,j-1)/RGd_t_dyn(i,j-1)(paprs(i,j-1)-paprs(i,j))
862	enddo
863	enddo
864
865	ELSE
866	DO k = 1, klev
867	DO i = 1, klon
868	d_u_dyn(i,k) = 0.0
869	d_t_dyn(i,k) = 0.0
870	ENDDO
871	ENDDO
872	ancien_ok = .TRUE.
873	ENDIF
874	c====================================================================
875
876	! Compute vertical velocity (Pa/s) from vertical mass flux
877	! Need to linearly interpolate mass flux to mid-layers
878	do k=1,klev-1
879	omega(1:klon,k) = 0.5RG(flxmw(1:klon,k)+flxmw(1:klon,k+1))
880	. / cell_area(1:klon)
881	enddo
882	omega(1:klon,klev) = 0.5RGflxmw(1:klon,klev) / cell_area(1:klon)
883
884	c======
885	c GEOP CORRECTION
886	c
887	c Ajouter le geopotentiel du sol:
888	c
889	DO k = 1, klev
890	DO i = 1, klon
891	zphi(i,k) = pphi(i,k) + pphis(i)
892	ENDDO
893	ENDDO
894
895	c............................
896	c CETTE CORRECTION VA DE PAIR AVEC DES MODIFS DE LEAPFROG(_p)
897	c ELLE MARCHE A 50 NIVEAUX (si mmean constante...)
898	c MAIS PAS A 78 NIVEAUX (quand mmean varie...)
899	c A ANALYSER PLUS EN DETAIL AVANT D'UTILISER
900	c............................
901	c zphi is recomputed (pphi is not ok if mean molecular mass varies)
902	c with dphi = RT/mmean d(ln p) [evaluated at interface]
903
904	c DO i = 1, klon
905	c zphi(i,1) = pphis(i) + Rt_seri(i,1)/mmean(i,1)1000.
906	c * *( log(paprs(i,1)) - log(pplay(i,1)) )
907	c DO k = 2, klev
908	c zphi(i,k) = zphi(i,k-1)
909	c * + R500.(t_seri(i,k)/mmean(i,k)+t_seri(i,k-1)/mmean(i,k-1))
910	c * * (log(pplay(i,k-1)) - log(pplay(i,k)))
911	c ENDDO
912	c ENDDO
913	c............................
914	c=====
915
916	c calcul de l'altitude aux niveaux intercouches
917	c ponderation des altitudes au niveau des couches en dp/p
918
919	DO i=1,klon
920	zzlay(i,1)=zphi(i,1)/RG ! [m]
921	zzlev(i,1)=pphis(i)/RG ! [m]
922	ENDDO
923	DO k=2,klev
924	DO i=1,klon
925	tlaymean=t_seri(i,k)
926	IF (t_seri(i,k).ne.t_seri(i,k-1))
927	& tlaymean=(t_seri(i,k)-t_seri(i,k-1))
928	& /log(t_seri(i,k)/t_seri(i,k-1))
929
930	zzlay(i,k)=zzlay(i,k-1)
931	& -(log(pplay(i,k)/pplay(i,k-1))rnew(i,k-1)tlaymean
932	& /(RG(RA/(RA+zzlay(i,k-1)))*2))
933	ENDDO
934	ENDDO
935	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
936	! Old : from geopotential. Problem when mu varies in the upper atm...
937	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
938	! DO k=1,klev
939	! DO i=1,klon
940	! zzlay(i,k)=zphi(i,k)/RG ! [m]
941	! ENDDO
942	! ENDDO
943	! DO i=1,klon
944	! zzlev(i,1)=pphis(i)/RG ! [m]
945	! ENDDO
946	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
947	DO k=2,klev
948	DO i=1,klon
949	z1=(pplay(i,k-1)+paprs(i,k))/(pplay(i,k-1)-paprs(i,k))
950	z2=(paprs(i,k) +pplay(i,k))/(paprs(i,k) -pplay(i,k))
951	zzlev(i,k)=(z1zzlay(i,k-1)+z2zzlay(i,k))/(z1+z2)
952	ENDDO
953	ENDDO
954	DO i=1,klon
955	zzlev(i,klev+1)=zzlay(i,klev)+(zzlay(i,klev)-zzlev(i,klev))
956	ENDDO
957
958	c====================================================================
959	c
960	c Verifier les temperatures
961	c
962	CALL hgardfou(t_seri,ftsol,'debutphy')
963
964	c====================================================================
965	c Orbite et eclairement
966	c=======================
967
968	c Pour VENUS, on fixe l'obliquite a 0 et l'eccentricite a 0.
969	c donc pas de variations de Ls, ni de dist.
970	c La seule chose qui compte, c'est la rotation de la planete devant
971	c le Soleil...
972
973	zlongi = 0.0
974	dist = 0.72333 ! en UA
975
976	c Si on veut remettre l'obliquite a 3 degres et/ou l'eccentricite
977	c a sa valeur, et prendre en compte leur evolution,
978	c il faudra refaire un orbite.F...
979	c CALL orbite(zlongi,dist)
980
981	IF (cycle_diurne) THEN
982	zdtime=dtime*REAL(radpas) ! pas de temps du rayonnement (s)
983	CALL zenang(zlongi,gmtime,zdtime,latitude_deg,longitude_deg,
984	& rmu0,fract)
985	ELSE
986	call mucorr(klon,zlongi,latitude_deg,rmu0,fract)
987	ENDIF
988
989	! print fraction of venus day
990
991	if (is_master) then
992	print*, 'gmtime = ', gmtime
993	end if
994
995	c======================================================================
996	c FIN INITIALISATIONS
997	c======================================================================
998	c======================================================================
999
1000	c=======================================================
1001	! CONDUCTION and MOLECULAR VISCOSITY
1002	c=======================================================
1003
1004	d_t_conduc(:,:)=0.
1005	d_u_molvis(:,:)=0.
1006	d_v_molvis(:,:)=0.
1007
1008	IF (callthermos) THEN
1009
1010	tsurf(:)=t_seri(:,1)
1011	call conduction(klon, klev,pdtphys,
1012	$ pplay,paprs,t_seri,
1013	$ tsurf,zzlev,zzlay,d_t_conduc)
1014
1015	call molvis(klon, klev,pdtphys,
1016	$ pplay,paprs,t_seri,
1017	$ u,tsurf,zzlev,zzlay,d_u_molvis)
1018
1019	call molvis(klon, klev, pdtphys,
1020	$ pplay,paprs,t_seri,
1021	$ v,tsurf,zzlev,zzlay,d_v_molvis)
1022
1023	DO k=1,klev
1024	DO ig=1,klon
1025	t_seri(ig,k)= t_seri(ig,k)+ d_t_conduc(ig,k)*dtime ! [K]
1026	u_seri(ig,k)= u_seri(ig,k)+ d_u_molvis(ig,k)*dtime ! [m/s]
1027	v_seri(ig,k)= v_seri(ig,k)+ d_v_molvis(ig,k)*dtime ! [m/s]
1028	ENDDO
1029	ENDDO
1030	ENDIF
1031
1032	c====================================================================
1033	c Compute mean mass, cp and R :
1034	c------------------------------------
1035
1036	if(callthermos) then
1037	call concentrations2(pplay,t_seri,tr_seri, nqmax)
1038	endif
1039
1040
1041	c=======================================================
1042	! CHEMISTRY AND MICROPHYSICS
1043	c=======================================================
1044
1045	if (iflag_trac.eq.1) then
1046
1047	if ( ok_aoa ) then
1048	call age_of_air(tr_seri(:,:,i_aoa),klon,klev,
1049	I itap,pdtphys,init_age,
1050	O age)
1051	end if
1052
1053	!====================================================================
1054	! Case 1: pseudo-chemistry with relaxation toward fixed profile
1055	!=========
1056	if (tr_scheme.eq.1) then
1057
1058	call phytrac_relax (debut,lafin,nqmax,
1059	I nlon,nlev,dtime,pplay,
1060	O tr_seri)
1061
1062	elseif (tr_scheme.eq.2) then
1063	!====================================================================
1064	! Case 2: surface emission
1065	! For the moment, inspired from Mars version
1066	! However, the variable 'source' could be used in physiq
1067	! so the call to phytrac_emiss could be to initialise it.
1068	!=========
1069
1070	call phytrac_emiss (debut,lafin,nqmax,
1071	I nlon,nlev,dtime,paprs,
1072	I latitude_deg,longitude_deg,
1073	O tr_seri)
1074
1075	else if (tr_scheme.eq.3) then
1076	!====================================================================
1077	! Case 3: Full chemistry and/or clouds.
1078	! routines are called every "chempas" physical timestep.
1079	!
1080	! if the physics is called 96000 times per venus day:
1081	!
1082	! nbapp_chem = 24000 => chempas = 4 => zctime = 420 s
1083	! nbapp_chem = 12000 => chempas = 8 => zctime = 840 s
1084	!=========
1085
1086
1087	chempas = nint(rday/pdtphys/nbapp_chem)
1088	zctime = dtime*real(chempas) ! chemical timestep
1089
1090	if (mod(itap,chempas) == 0) then ! <------- start of chemistry supercycling
1091
1092	! photochemistry and microphysics
1093
1094	call phytrac_chimie(debut,
1095	$ gmtime,
1096	$ nqmax,
1097	$ klon,
1098	$ latitude_deg,
1099	$ longitude_deg,
1100	$ zzlay,
1101	$ nlev,
1102	$ zctime,
1103	$ t_seri,
1104	$ pplay,
1105	$ tr_seri,
1106	$ d_tr_chem,
1107	$ iter,
1108	$ prod_tr,
1109	$ loss_tr,
1110	$ no_emission,
1111	$ o2_emission)
1112
1113	if (ok_sedim) then
1114	if (cl_scheme == 1) then
1115
1116	! sedimentation for simplified microphysics
1117
1118	#ifndef MESOSCALE
1119	call new_cloud_sedim(nb_mode,
1120	$ klon,
1121	$ nlev,
1122	$ zctime,
1123	$ pplay,
1124	$ paprs,
1125	$ t_seri,
1126	$ tr_seri,
1127	$ d_tr_chem,
1128	$ d_tr_sed(:,:,1:2),
1129	$ nqmax,
1130	$ Fsedim)
1131
1132	! test to avoid nans
1133
1134	do k = 1, klev
1135	do i = 1, klon
1136	if ((d_tr_sed(i,k,1) /= d_tr_sed(i,k,1)) .or.
1137	$ (d_tr_sed(i,k,2) /= d_tr_sed(i,k,2))) then
1138	print*,'sedim NaN PROBLEM'
1139	print*,'d_tr_sed Nan?',d_tr_sed(i,k,:)
1140	print*,'Temp',t_seri(i,k)
1141	print*,'lat-lon',i,'level',k,'zctime',zctime
1142	print*,'F_sed',Fsedim(i,k)
1143	d_tr_sed(i,k,:) = 0.
1144	end if
1145	end do
1146	end do
1147
1148	! tendency due to condensation and sedimentation
1149
1150	d_tr_sed(:,:,1:2) = d_tr_sed(:,:,1:2)/zctime
1151	Fsedim(:,1:klev) = Fsedim(:,1:klev)/zctime
1152	Fsedim(:,klev+1) = 0.
1153
1154	else if (cl_scheme == 2) then
1155
1156	! sedimentation for detailed microphysics
1157
1158	d_tr_sed(:,:,:) = 0.
1159
1160	do i = 1, klon
1161
1162	! mode 1
1163
1164	m0_mode1(:,1) = tr_seri(i,:,i_m0_mode1drop)
1165	m0_mode1(:,2) = tr_seri(i,:,i_m0_mode1ccn)
1166	m3_mode1(:,1) = tr_seri(i,:,i_m3_mode1sa)
1167	m3_mode1(:,2) = tr_seri(i,:,i_m3_mode1w)
1168	m3_mode1(:,3) = tr_seri(i,:,i_m3_mode1ccn)
1169
1170	call drop_sedimentation(zctime,klev,m0_mode1,m3_mode1,
1171	$ paprs(i,:),zzlev(i,:),
1172	$ zzlay(i,:),t_seri(i,:),1,
1173	$ d_ccn_sed(:,1),d_drop_sed,
1174	$ d_ccn_sed(:,2),d_liq_sed)
1175
1176	d_tr_sed(i,:,i_m0_mode1drop)= d_tr_sed(i,:,i_m0_mode1drop)
1177	$ + d_drop_sed
1178	d_tr_sed(i,:,i_m0_mode1ccn) = d_tr_sed(i,:,i_m0_mode1ccn)
1179	$ + d_ccn_sed(:,1)
1180	d_tr_sed(i,:,i_m3_mode1ccn) = d_tr_sed(i,:,i_m3_mode1ccn)
1181	$ + d_ccn_sed(:,2)
1182	d_tr_sed(i,:,i_m3_mode1sa) = d_tr_sed(i,:,i_m3_mode1sa)
1183	$ + d_liq_sed(:,1)
1184	d_tr_sed(i,:,i_m3_mode1w) = d_tr_sed(i,:,i_m3_mode1w)
1185	$ + d_liq_sed(:,2)
1186
1187	! mode 2
1188
1189	m0_mode2(:,1) = tr_seri(i,:,i_m0_mode2drop)
1190	m0_mode2(:,2) = tr_seri(i,:,i_m0_mode2ccn)
1191	m3_mode2(:,1) = tr_seri(i,:,i_m3_mode2sa)
1192	m3_mode2(:,2) = tr_seri(i,:,i_m3_mode2w)
1193	m3_mode2(:,3) = tr_seri(i,:,i_m3_mode2ccn)
1194
1195	call drop_sedimentation(zctime,klev,m0_mode2,m3_mode2,
1196	$ paprs(i,:),zzlev(i,:),
1197	$ zzlay(i,:),t_seri(i,:),2,
1198	$ d_ccn_sed(:,1),d_drop_sed,
1199	$ d_ccn_sed(:,2),d_liq_sed)
1200
1201	d_tr_sed(i,:,i_m0_mode2drop)= d_tr_sed(i,:,i_m0_mode2drop)
1202	$ + d_drop_sed
1203	d_tr_sed(i,:,i_m0_mode2ccn) = d_tr_sed(i,:,i_m0_mode2ccn)
1204	$ + d_ccn_sed(:,1)
1205	d_tr_sed(i,:,i_m3_mode2ccn) = d_tr_sed(i,:,i_m3_mode2ccn)
1206	$ + d_ccn_sed(:,2)
1207	d_tr_sed(i,:,i_m3_mode2sa) = d_tr_sed(i,:,i_m3_mode2sa)
1208	$ + d_liq_sed(:,1)
1209	d_tr_sed(i,:,i_m3_mode2w) = d_tr_sed(i,:,i_m3_mode2w)
1210	$ + d_liq_sed(:,2)
1211
1212	! aer
1213
1214	call aer_sedimentation(zctime,klev,
1215	$ tr_seri(i,:,i_m0_aer),
1216	$ tr_seri(i,:,i_m3_aer),
1217	$ paprs(i,:),zzlev(i,:),
1218	$ zzlay(i,:),t_seri(i,:),
1219	$ d_tr_sed(i,:,i_m0_aer),
1220	$ d_tr_sed(i,:,i_m3_aer),
1221	$ aer_flux)
1222
1223	end do
1224
1225	! tendency due to sedimentation
1226
1227	do iq = nqmax-nmicro+1,nqmax
1228	d_tr_sed(:,:,iq) = d_tr_sed(:,:,iq)/zctime
1229	end do
1230	#endif
1231	end if ! cl_scheme
1232
1233	! update gaseous tracers (chemistry)
1234
1235	do iq = 1, nqmax - nmicro
1236	tr_seri(:,:,iq) = max(tr_seri(:,:,iq)
1237	$ + d_tr_chem(:,:,iq)*zctime,1.e-30)
1238	end do
1239
1240	! update condensed tracers (condensation + sedimentation)
1241
1242	if (cl_scheme == 1) then
1243	tr_seri(:,:,i_h2so4liq) = max(tr_seri(:,:,i_h2so4liq)
1244	$ + d_tr_sed(:,:,1)*zctime, 1.e-30)
1245	tr_seri(:,:,i_h2oliq) = max(tr_seri(:,:,i_h2oliq)
1246	$ + d_tr_sed(:,:,2)*zctime, 1.e-30)
1247	else if (cl_scheme == 2) then
1248	do iq = nqmax-nmicro+1,nqmax
1249	tr_seri(:,:,iq) = max(tr_seri(:,:,iq)
1250	$ + d_tr_sed(:,:,iq)*zctime,1.e-30)
1251	end do
1252	end if ! cl_scheme
1253
1254	end if ! ok_sedim
1255	end if ! mod(itap,chempas) <------- end of chemistry supercycling
1256
1257	!=========
1258	! End Case 3: Full chemistry and/or clouds.
1259	!====================================================================
1260
1261	end if ! tr_scheme
1262	end if ! iflag_trac
1263
1264	c====================================================================
1265	c Appeler la diffusion verticale (programme de couche limite)
1266	c====================================================================
1267
1268	c-------------------------------
1269	c VENUS TEST: on ne tient pas compte des calculs de clmain mais on force
1270	c l'equilibre radiatif du sol
1271	if (.not. ok_clmain) then
1272	if (debut) then
1273	print*,"ATTENTION, CLMAIN SHUNTEE..."
1274	endif
1275
1276	DO i = 1, klon
1277	sens(i) = 0.0e0 ! flux de chaleur sensible au sol
1278	fder(i) = 0.0e0
1279	dlw(i) = 0.0e0
1280	ENDDO
1281
1282	c Incrementer la temperature du sol
1283	c
1284	DO i = 1, klon
1285	d_ts(i) = dtime * radsol(i)/22000. !valeur calculee par GCM pour I=200
1286	ftsol(i) = ftsol(i) + d_ts(i)
1287	do j=1,nsoilmx
1288	ftsoil(i,j)=ftsol(i)
1289	enddo
1290	ENDDO
1291
1292	c-------------------------------
1293	else
1294	c-------------------------------
1295
1296	fder = dlw
1297
1298	! ADAPTATION GCM POUR CP(T)
1299
1300	if (physideal) then
1301	CALL clmain_ideal(dtime,itap,
1302	e t_seri,u_seri,v_seri,
1303	e rmu0,
1304	e ftsol,
1305	$ ftsoil,
1306	$ paprs,pplay,ppk,radsol,falbe,
1307	e solsw, sollw, sollwdown, fder,
1308	e longitude_deg, latitude_deg, dx, dy,
1309	e debut, lafin,
1310	s d_t_vdf,d_u_vdf,d_v_vdf,d_ts,
1311	s fluxt,fluxu,fluxv,cdragh,cdragm,
1312	s dsens,
1313	s ycoefh,yu1,yv1)
1314	else
1315	CALL clmain(dtime,itap,
1316	e t_seri,u_seri,v_seri,
1317	e rmu0,
1318	e ftsol,
1319	$ ftsoil,
1320	$ paprs,pplay,ppk,radsol,falbe,
1321	e solsw, sollw, sollwdown, fder,
1322	e longitude_deg, latitude_deg, dx, dy,
1323	& q2,
1324	e debut, lafin,
1325	s d_t_vdf,d_u_vdf,d_v_vdf,d_ts,
1326	s fluxt,fluxu,fluxv,cdragh,cdragm,
1327	s dsens,
1328	s ycoefh,yu1,yv1)
1329	endif
1330
1331	CXXX Incrementation des flux
1332	DO i = 1, klon
1333	sens(i) = - fluxt(i,1) ! flux de chaleur sensible au sol
1334	fder(i) = dlw(i) + dsens(i)
1335	ENDDO
1336	CXXX
1337	IF (.not. turb_resolved) then !True only for LES
1338	DO k = 1, klev
1339	DO i = 1, klon
1340	t_seri(i,k) = t_seri(i,k) + d_t_vdf(i,k)
1341	d_t_vdf(i,k)= d_t_vdf(i,k)/dtime ! K/s
1342	u_seri(i,k) = u_seri(i,k) + d_u_vdf(i,k)
1343	d_u_vdf(i,k)= d_u_vdf(i,k)/dtime ! (m/s)/s
1344	v_seri(i,k) = v_seri(i,k) + d_v_vdf(i,k)
1345	d_v_vdf(i,k)= d_v_vdf(i,k)/dtime ! (m/s)/s
1346	ENDDO
1347	ENDDO
1348	ENDIF
1349	C TRACEURS
1350
1351	if (iflag_trac.eq.1) then
1352	DO k = 1, klev
1353	DO i = 1, klon
1354	delp(i,k) = paprs(i,k)-paprs(i,k+1)
1355	ENDDO
1356	ENDDO
1357
1358	if (klon.eq.1) then !For the 1D model
1359	!Reading the prescribed profile from deftank
1360	open(35,file='kzz_p.txt',form='formatted',status='old')
1361	rewind(35)
1362	DO k=1,klev
1363	read(35,*) kzz_p(k)
1364	ENDDO
1365	close(35)
1366
1367	!Implementing the new profile in m2/s
1368	DO k = 1, klev
1369	ycoefh(1,k) = kzz_p(k)
1370	ENDDO
1371	endif
1372
1373	DO iq=1, nqmax
1374
1375	CALL cltrac(dtime,ycoefh,t_seri,
1376	s tr_seri(:,:,iq),source(:,iq),
1377	e paprs, pplay,delp,
1378	s d_tr_vdf(:,:,iq))
1379
1380	tr_seri(:,:,iq) = tr_seri(:,:,iq) + d_tr_vdf(:,:,iq)
1381	d_tr_vdf(:,:,iq)= d_tr_vdf(:,:,iq)/dtime ! /s
1382
1383	ENDDO !nqmax
1384
1385	endif
1386
1387	IF (if_ebil.ge.2) THEN
1388	ztit='after clmain'
1389	CALL diagetpq(cell_area,ztit,ip_ebil,2,1,dtime
1390	e , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
1391	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1392	call diagphy(cell_area,ztit,ip_ebil
1393	e , zero_v, zero_v, zero_v, zero_v, sens
1394	e , zero_v, zero_v, zero_v, ztsol
1395	e , d_h_vcol, d_qt, d_ec
1396	s , fs_bound, fq_bound )
1397	END IF
1398	C
1399	c
1400	c Incrementer la temperature du sol
1401	c
1402	DO i = 1, klon
1403	ftsol(i) = ftsol(i) + d_ts(i)
1404	ENDDO
1405
1406	c Calculer la derive du flux infrarouge
1407	c
1408	DO i = 1, klon
1409	dlw(i) = - 4.0RSIGMAftsol(i)**3
1410	ENDDO
1411
1412	c-------------------------------
1413	endif ! fin du VENUS TEST
1414
1415	! tests: output tendencies
1416	! call writefield_phy('physiq_d_t_vdf',d_t_vdf,klev)
1417	! call writefield_phy('physiq_d_u_vdf',d_u_vdf,klev)
1418	! call writefield_phy('physiq_d_v_vdf',d_v_vdf,klev)
1419	! call writefield_phy('physiq_d_ts',d_ts,1)
1420
1421	c===================================================================
1422	c Convection seche
1423	c===================================================================
1424	c
1425	d_t_ajs(:,:)=0.
1426	d_u_ajs(:,:)=0.
1427	d_v_ajs(:,:)=0.
1428	d_tr_ajs(:,:,:)=0.
1429	c
1430	IF(prt_level>9)WRITE(lunout,*)
1431	. 'AVANT LA CONVECTION SECHE , iflag_ajs='
1432	s ,iflag_ajs
1433
1434	if(iflag_ajs.eq.0) then
1435	c Rien
1436	c ====
1437	IF(prt_level>9)WRITE(lunout,*)'pas de convection'
1438
1439	else if(iflag_ajs.eq.1) then
1440
1441	c Ajustement sec
1442	c ==============
1443	IF(prt_level>9)WRITE(lunout,*)'ajsec'
1444
1445	! ADAPTATION GCM POUR CP(T)
1446	CALL ajsec(paprs, pplay, ppk, t_seri, u_seri, v_seri, nqmax,
1447	. tr_seri, d_t_ajs, d_u_ajs, d_v_ajs, d_tr_ajs)
1448
1449	! ADAPTATION GCM POUR CP(T)
1450	do i=1,klon
1451	flux_ajs(i,1) = 0.0
1452	do j=2,klev
1453	flux_ajs(i,j) = flux_ajs(i,j-1)
1454	. + cpnew(i,j-1)/RG*d_t_ajs(i,j-1)/dtime
1455	. *(paprs(i,j-1)-paprs(i,j))
1456	enddo
1457	enddo
1458
1459	t_seri(:,:) = t_seri(:,:) + d_t_ajs(:,:)
1460	d_t_ajs(:,:)= d_t_ajs(:,:)/dtime ! K/s
1461	u_seri(:,:) = u_seri(:,:) + d_u_ajs(:,:)
1462	d_u_ajs(:,:)= d_u_ajs(:,:)/dtime ! (m/s)/s
1463	v_seri(:,:) = v_seri(:,:) + d_v_ajs(:,:)
1464	d_v_ajs(:,:)= d_v_ajs(:,:)/dtime ! (m/s)/s
1465
1466	if (iflag_trac.eq.1) then
1467	tr_seri(:,:,:) = tr_seri(:,:,:) + d_tr_ajs(:,:,:)
1468	d_tr_ajs(:,:,:)= d_tr_ajs(:,:,:)/dtime ! /s
1469	endif
1470	endif
1471
1472	! tests: output tendencies
1473	! call writefield_phy('physiq_d_t_ajs',d_t_ajs,klev)
1474	! call writefield_phy('physiq_d_u_ajs',d_u_ajs,klev)
1475	! call writefield_phy('physiq_d_v_ajs',d_v_ajs,klev)
1476	c
1477	IF (if_ebil.ge.2) THEN
1478	ztit='after dry_adjust'
1479	CALL diagetpq(cell_area,ztit,ip_ebil,2,2,dtime
1480	e , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
1481	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1482	call diagphy(cell_area,ztit,ip_ebil
1483	e , zero_v, zero_v, zero_v, zero_v, sens
1484	e , zero_v, zero_v, zero_v, ztsol
1485	e , d_h_vcol, d_qt, d_ec
1486	s , fs_bound, fq_bound )
1487	END IF
1488
1489	c====================================================================
1490	c RAYONNEMENT
1491	c====================================================================
1492	if (mod(itap,radpas) == 0) then
1493
1494	dtimerad = dtime*REAL(radpas) ! pas de temps du rayonnement (s)
1495
1496	! update mmean
1497	if (ok_chem) then
1498	mmean(:,:) = 0.
1499	do iq = 1,nqmax - nmicro
1500	mmean(:,:) = mmean(:,:)+tr_seri(:,:,iq)/m_tr(iq)
1501	enddo
1502	mmean(:,:) = 1./mmean(:,:)
1503	rnew(:,:) = 8.314/mmean(:,:)*1.e3 ! J/kg K
1504	endif
1505
1506	cc---------------------------------------------
1507	if (callnlte .or. callthermos) then
1508	if (ok_chem) then
1509
1510	! nlte : use computed chemical species
1511
1512	co2vmr_gcm(:,:) = tr_seri(:,:,i_co2)*mmean(:,:)/m_tr(i_co2)
1513	covmr_gcm(:,:) = tr_seri(:,:,i_co)*mmean(:,:)/m_tr(i_co)
1514	ovmr_gcm(:,:) = tr_seri(:,:,i_o)*mmean(:,:)/m_tr(i_o)
1515	n2vmr_gcm(:,:) = tr_seri(:,:,i_n2)*mmean(:,:)/m_tr(i_n2)
1516
1517	else
1518
1519	! nlte : use hedin climatology
1520
1521	call compo_hedin83_mod(pplay,rmu0,
1522	$ co2vmr_gcm,covmr_gcm,ovmr_gcm,n2vmr_gcm,nvmr_gcm)
1523	end if
1524	end if
1525
1526	c NLTE cooling from CO2 emission
1527	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1528
1529	IF(callnlte) THEN
1530	if(nltemodel.eq.0.or.nltemodel.eq.1) then
1531	c nltecool call not correct...
1532	c CALL nltecool(klon, klev, nqmax, pplay*9.869e-6, t_seri,
1533	c $ tr_seri, d_t_nlte)
1534	abort_message='nltemodel=0 or 1 should not be used...'
1535	call abort_physic(modname,abort_message,1)
1536	else if(nltemodel.eq.2) then
1537	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1538	!! HEDIN instead of compo for this calculation
1539	! call compo_hedin83_mod(pplay,rmu0,
1540	! $ co2vmr_gcm,covmr_gcm,ovmr_gcm,n2vmr_gcm,nvmr_gcm)
1541	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1542	CALL nlte_tcool(klon,klev,pplay*9.869e-6,
1543	$ t_seri,zzlay,co2vmr_gcm, n2vmr_gcm, covmr_gcm,
1544	$ ovmr_gcm,d_t_nlte,ierr_nlte,varerr )
1545	if(ierr_nlte.gt.0) then
1546	write(,)
1547	$ 'WARNING: nlte_tcool output with error message',
1548	$ 'ierr_nlte=',ierr_nlte,'varerr=',varerr
1549	write(,)'I will continue anyway'
1550	endif
1551	endif
1552
1553	ELSE
1554
1555	d_t_nlte(:,:)=0.
1556
1557	ENDIF
1558
1559	c Find number of layers for LTE radiation calculations
1560
1561	IF(callnlte .or. callnirco2)
1562	$ CALL nlthermeq(klon, klev, paprs, pplay)
1563
1564	cc---------------------------------------------
1565	c LTE radiative transfert / solar / IR matrix
1566	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1567	if (physideal) then
1568	CALL radlwsw_newtoncool(presnivs,t_seri)
1569	else
1570	CALL radlwsw
1571	e (dist, rmu0, fract, zzlev,
1572	e paprs, pplay,ftsol, t_seri)
1573	endif
1574
1575	c ALBEDO VARIATIONS: test for Yeon Joo Lee
1576	c increment to increase it for 20 Vd => +80%
1577	c heat(:,:)=heat(:,:)(1.+0.80((rjourvrai-356)+gmtime)/20.)
1578	c or to decrease it for 20 Vd => 1/1.8
1579	c heat(:,:)=heat(:,:)/(1.+0.80*((rjourvrai-356)+gmtime)/20.)
1580
1581	c ------------ ALBEDO VARIATIONS: scenarios for VCD
1582	c shape of relative variation from Lee et al 2019 (Fig 13b)
1583	c between 57 km (4e4 Pa) and 72 km (2.5e3 Pa), peak at 67 km (6e3 Pa)
1584	c do j=1,klev
1585	c factAlb = 0.
1586	c if ((presnivs(j).gt.6.e3).and.(presnivs(j).lt.4.e4)) then
1587	c factAlb = (log(presnivs(j))-log(4.e4))/(log(6.e3)-log(4.e4))
1588	c elseif ((presnivs(j).lt.6.e3).and.(presnivs(j).gt.2.5e3)) then
1589	c factAlb = (log(presnivs(j))-log(2.5e3))/(log(6.e3)-log(2.5e3))
1590	c endif
1591	c Increase by 50% (Minimum albedo)
1592	c heat(:,j)=heat(:,j)(1+factAlb0.5)
1593	c Decrease by 30% (Maximum albedo)
1594	c heat(:,j)=heat(:,j)(1-factAlb0.3)
1595	c enddo
1596	c ------------ END ALBEDO VARIATIONS
1597
1598	cc---------------------------------------------
1599	c CO2 near infrared absorption
1600	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1601
1602	d_t_nirco2(:,:)=0.
1603	if (callnirco2) then
1604	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1605	!! HEDIN instead of compo for this calculation
1606	! call compo_hedin83_mod(pplay,rmu0,
1607	! $ co2vmr_gcm,covmr_gcm,ovmr_gcm,n2vmr_gcm,nvmr_gcm)
1608	! tr_hedin(:,:,:)=tr_seri(:,:,:)
1609	! tr_hedin(:,:,i_co2)=co2vmr_gcm(:,:)/mmean(:,:)*m_tr(i_co2)
1610	! tr_hedin(:,:,i_co) = covmr_gcm(:,:)/mmean(:,:)*m_tr(i_co)
1611	! tr_hedin(:,:,i_o) = ovmr_gcm(:,:)/mmean(:,:)*m_tr(i_o)
1612	! tr_hedin(:,:,i_n2) = n2vmr_gcm(:,:)/mmean(:,:)*m_tr(i_n2)
1613	! call nirco2abs (klon, klev, pplay, dist, nqmax, tr_hedin,
1614	! . rmu0, fract, d_t_nirco2)
1615	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1616	call nirco2abs (klon, klev, pplay, dist, nqmax, tr_seri,
1617	. rmu0, fract, d_t_nirco2)
1618	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1619	endif
1620
1621
1622	cc---------------------------------------------
1623	c Net atmospheric radiative heating rate (K.s-1)
1624	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1625
1626	IF(callnlte.or.callnirco2) THEN
1627	CALL blendrad(klon, klev, pplay,heat,
1628	& cool, d_t_nirco2,d_t_nlte, dtsw, dtlw)
1629	ELSE
1630	dtsw(:,:)=heat(:,:)
1631	dtlw(:,:)=-1*cool(:,:)
1632	ENDIF
1633
1634	DO k=1,klev
1635	DO i=1,klon
1636	d_t_rad(i,k) = dtsw(i,k) + dtlw(i,k) ! K/s
1637	ENDDO
1638	ENDDO
1639
1640
1641	cc---------------------------------------------
1642	c EUV heating rate (K.s-1)
1643	c ~~~~~~~~~~~~~~~~~~~~~~~~
1644
1645	d_t_euv(:,:)=0.
1646
1647	IF (callthermos) THEN
1648
1649	call euvheat(klon, klev, nqmax, t_seri,paprs,pplay,zzlay,
1650	$ rmu0,dtimerad,gmtime,
1651	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1652	!! HEDIN instead of compo for this calculation
1653	!! cf nlte_tcool and nirco2abs above !!
1654	! $ tr_hedin, d_tr, d_t_euv )
1655	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1656	$ tr_seri, d_tr, d_t_euv )
1657	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1658
1659	DO k=1,klev
1660	DO ig=1,klon
1661	d_t_rad(ig,k)=d_t_rad(ig,k)+d_t_euv(ig,k)
1662	ENDDO
1663	ENDDO
1664
1665	ENDIF ! callthermos
1666
1667	ENDIF ! radpas
1668	c====================================================================
1669	c
1670	c Ajouter la tendance des rayonnements (tous les pas)
1671	c
1672	DO k = 1, klev
1673	DO i = 1, klon
1674	t_seri(i,k) = t_seri(i,k) + d_t_rad(i,k) * dtime
1675	ENDDO
1676	ENDDO
1677
1678	! increment physics counter
1679
1680	itap = itap + 1
1681	c====================================================================
1682
1683
1684	c==============================
1685	! -- MOLECULAR DIFFUSION ---
1686	c==============================
1687
1688	d_q_moldif(:,:,:)=0
1689
1690	IF (callthermos .and. ok_chem) THEN
1691
1692	call moldiff_red(klon, klev, nqmax,
1693	& pplay,paprs,t_seri, tr_seri, pdtphys,
1694	& d_t_euv,d_t_conduc,d_q_moldif)
1695
1696
1697	! --- update tendencies tracers ---
1698
1699	DO iq = 1, nqmax
1700	DO k=1,klev
1701	DO ig=1,klon
1702	tr_seri(ig,k,iq)= max(tr_seri(ig,k,iq)+
1703	& d_q_moldif(ig,k,iq)*dtime,1.e-30)
1704	ENDDO
1705	ENDDO
1706	ENDDO
1707
1708	ENDIF ! callthermos & ok_chem
1709
1710	c====================================================================
1711
1712	c==================================
1713	! -- ION AMBIPOLAR DIFFUSION ---
1714	c==================================
1715
1716	d_q_iondif(:,:,:)=0
1717
1718	IF (callthermos .and. ok_chem .and. ok_ionchem) THEN
1719	IF (ok_iondiff) THEN
1720
1721	call iondiff_red(klon,klev,nqmax,pplay,paprs,
1722	& t_seri,tr_seri,pphis,
1723	& gmtime,latitude_deg,longitude_deg,
1724	& pdtphys,d_t_euv,d_t_conduc,d_q_moldif,
1725	& d_q_iondif)
1726
1727	!write (,) 'TITI EST PASSE PAR LA'
1728	! --- update tendencies tracers ---
1729
1730	DO iq = 1, nqmax
1731	IF (type_tr(iq) .eq. 2) THEN
1732	DO k=1,klev
1733	DO ig=1,klon
1734	tr_seri(ig,k,iq)= max(tr_seri(ig,k,iq)+
1735	& d_q_iondif(ig,k,iq)*dtime,1.e-30)
1736	ENDDO
1737	ENDDO
1738	ENDIF
1739	ENDDO
1740	ENDIF ! ok_iondiff
1741	ENDIF ! callthermos & ok_chem & ok_ionchem
1742
1743	c====================================================================
1744	! tests: output tendencies
1745	! call writefield_phy('physiq_dtrad',dtrad,klev)
1746
1747	IF (if_ebil.ge.2) THEN
1748	ztit='after rad'
1749	CALL diagetpq(cell_area,ztit,ip_ebil,2,2,dtime
1750	e , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
1751	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1752	call diagphy(cell_area,ztit,ip_ebil
1753	e , topsw, toplw, solsw, sollw, zero_v
1754	e , zero_v, zero_v, zero_v, ztsol
1755	e , d_h_vcol, d_qt, d_ec
1756	s , fs_bound, fq_bound )
1757	END IF
1758	c
1759
1760	c====================================================================
1761	c Calcul des gravity waves FLOTT
1762	c====================================================================
1763	c
1764	c if (ok_orodr.or.ok_gw_nonoro) then
1765
1766	c CALCUL DE N2
1767	c UTILISE LA RELATION ENTRE N2 ET STABILITE
1768	c N2 = RG/T (dT/dz+RG/cp(T))
1769	c ET DONC EN N'UTILISE QUE LA TEMPERATURE, PAS teta.
1770
1771	do i=1,klon
1772	do k=2,klev
1773	ztlev(i,k) = (t_seri(i,k)+t_seri(i,k-1))/2.
1774	enddo
1775	enddo
1776	do i=1,klon
1777	do k=2,klev
1778	ztlev(i,k) = (t_seri(i,k)+t_seri(i,k-1))/2.
1779	zdtlev(i,k) = t_seri(i,k)-t_seri(i,k-1)
1780	zdzlev(i,k) = (zzlay(i,k)-zzlay(i,k-1))
1781	zn2(i,k) = RG/ztlev(i,k) * ( zdtlev(i,k)/zdzlev(i,k)
1782	. + RG/cpnew(i,k) )
1783	zn2(i,k) = max(zn2(i,k),1.e-12) ! securite
1784	enddo
1785	zn2(i,1) = 1.e-12 ! securite
1786	enddo
1787
1788	c endif
1789
1790	c ----------------------------ORODRAG
1791	IF (ok_orodr) THEN
1792	c
1793	c selection des points pour lesquels le shema est actif:
1794	igwd=0
1795	DO i=1,klon
1796	itest(i)=0
1797	c IF ((zstd(i).gt.10.0)) THEN
1798	IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
1799	itest(i)=1
1800	igwd=igwd+1
1801	idx(igwd)=i
1802	ENDIF
1803	ENDDO
1804	c igwdim=MAX(1,igwd)
1805	c
1806	c A ADAPTER POUR VENUS!!! [ TN: c'est fait ! ]
1807	CALL drag_noro(klon,klev,dtime,paprs,pplay,pphi,zn2,
1808	e zmea,zstd, zsig, zgam, zthe,zpic,zval,
1809	e igwd,idx,itest,
1810	e t_seri, u_seri, v_seri,
1811	s zulow, zvlow, zustrdr, zvstrdr,
1812	s d_t_oro, d_u_oro, d_v_oro,
1813	s zublstrdr,zvblstrdr,znlow,zeff,zbl,
1814	s ztau,tau0,knu2,kbreak)
1815
1816	c print*,"d_u_oro=",d_u_oro(klon/2,:)
1817	c ajout des tendances
1818	t_seri(:,:) = t_seri(:,:) + d_t_oro(:,:)
1819	d_t_oro(:,:)= d_t_oro(:,:)/dtime ! K/s
1820	u_seri(:,:) = u_seri(:,:) + d_u_oro(:,:)
1821	d_u_oro(:,:)= d_u_oro(:,:)/dtime ! (m/s)/s
1822	v_seri(:,:) = v_seri(:,:) + d_v_oro(:,:)
1823	d_v_oro(:,:)= d_v_oro(:,:)/dtime ! (m/s)/s
1824	c
1825	ELSE
1826	d_t_oro = 0.
1827	d_u_oro = 0.
1828	d_v_oro = 0.
1829	zustrdr = 0.
1830	zvstrdr = 0.
1831	zublstrdr = 0.
1832	zvblstrdr = 0.
1833	znlow = 0.
1834	zeff = 0.
1835	zbl = 0
1836	knu2 = 0
1837	kbreak = 0
1838	ztau = 0
1839	tau0 = 0.
1840	c
1841	ENDIF ! fin de test sur ok_orodr
1842	c
1843	! tests: output tendencies
1844	! call writefield_phy('physiq_d_t_oro',d_t_oro,klev)
1845	! call writefield_phy('physiq_d_u_oro',d_u_oro,klev)
1846	! call writefield_phy('physiq_d_v_oro',d_v_oro,klev)
1847
1848	c ----------------------------OROLIFT
1849	IF (ok_orolf) THEN
1850	print*,"ok_orolf NOT IMPLEMENTED !"
1851	stop
1852	c
1853	c selection des points pour lesquels le shema est actif:
1854	igwd=0
1855	DO i=1,klon
1856	itest(i)=0
1857	IF ((zpic(i)-zmea(i)).GT.100.) THEN
1858	itest(i)=1
1859	igwd=igwd+1
1860	idx(igwd)=i
1861	ENDIF
1862	ENDDO
1863	c igwdim=MAX(1,igwd)
1864	c
1865	c A ADAPTER POUR VENUS!!!
1866	c CALL lift_noro(klon,klev,dtime,paprs,pplay,
1867	c e latitude_deg,zmea,zstd,zpic,zgam,zthe,zpic,zval,
1868	c e igwd,idx,itest,
1869	c e t_seri, u_seri, v_seri,
1870	c s zulow, zvlow, zustrli, zvstrli,
1871	c s d_t_lif, d_u_lif, d_v_lif )
1872
1873	c
1874	c ajout des tendances
1875	t_seri(:,:) = t_seri(:,:) + d_t_lif(:,:)
1876	d_t_lif(:,:)= d_t_lif(:,:)/dtime ! K/s
1877	u_seri(:,:) = u_seri(:,:) + d_u_lif(:,:)
1878	d_u_lif(:,:)= d_u_lif(:,:)/dtime ! (m/s)/s
1879	v_seri(:,:) = v_seri(:,:) + d_v_lif(:,:)
1880	d_v_lif(:,:)= d_v_lif(:,:)/dtime ! (m/s)/s
1881	c
1882	ELSE
1883	d_t_lif = 0.
1884	d_u_lif = 0.
1885	d_v_lif = 0.
1886	zustrli = 0.
1887	zvstrli = 0.
1888	c
1889	ENDIF ! fin de test sur ok_orolf
1890
1891	c ---------------------------- NON-ORO GRAVITY WAVES
1892	IF(ok_gw_nonoro) then
1893
1894	! Obsolete
1895	! but used for VCD 1.1
1896	! call flott_gwd_ran(klon,klev,dtime,pplay,zn2,
1897	! e t_seri, u_seri, v_seri, paprs(klon/2+1,:),
1898	! o zustrhi,zvstrhi,
1899	! o d_t_hin, d_u_hin, d_v_hin)
1900
1901	! New routine based on Generic
1902	! used after VCD 1.1, for VCD 2.0
1903	call nonoro_gwd_ran(klon,klev,dtime,pplay,zn2,presnivs,
1904	e t_seri, u_seri, v_seri,
1905	o zustrhi,zvstrhi,
1906	o d_t_hin, d_u_hin, d_v_hin)
1907
1908	c ajout des tendances
1909
1910	t_seri(:,:) = t_seri(:,:) + d_t_hin(:,:)
1911	d_t_hin(:,:)= d_t_hin(:,:)/dtime ! K/s
1912	u_seri(:,:) = u_seri(:,:) + d_u_hin(:,:)
1913	d_u_hin(:,:)= d_u_hin(:,:)/dtime ! (m/s)/s
1914	v_seri(:,:) = v_seri(:,:) + d_v_hin(:,:)
1915	d_v_hin(:,:)= d_v_hin(:,:)/dtime ! (m/s)/s
1916
1917	ELSE
1918	d_t_hin = 0.
1919	d_u_hin = 0.
1920	d_v_hin = 0.
1921	zustrhi = 0.
1922	zvstrhi = 0.
1923
1924	ENDIF ! fin de test sur ok_gw_nonoro
1925
1926	! tests: output tendencies
1927	! call writefield_phy('physiq_d_t_hin',d_t_hin,klev)
1928	! call writefield_phy('physiq_d_u_hin',d_u_hin,klev)
1929	! call writefield_phy('physiq_d_v_hin',d_v_hin,klev)
1930
1931	c====================================================================
1932	c Transport de ballons
1933	c====================================================================
1934	if (ballons.eq.1) then
1935	CALL ballon(30,pdtphys,rjourvrai,gmtime*RDAY,
1936	& latitude_deg,longitude_deg,
1937	c C t,pplay,u,v,pphi) ! alt above surface (smoothed for GCM)
1938	C t,pplay,u,v,zphi) ! alt above planet average radius
1939	endif !ballons
1940
1941	c====================================================================
1942	c Bilan de mmt angulaire
1943	c====================================================================
1944	if (bilansmc.eq.1) then
1945	CMODDEB FLOTT
1946	C CALCULER LE BILAN DE MOMENT ANGULAIRE (DIAGNOSTIQUE)
1947	C STRESS NECESSAIRES: COUCHE LIMITE ET TOUTE LA PHYSIQUE
1948
1949	DO i = 1, klon
1950	zustrph(i)=0.
1951	zvstrph(i)=0.
1952	zustrcl(i)=0.
1953	zvstrcl(i)=0.
1954	ENDDO
1955	DO k = 1, klev
1956	DO i = 1, klon
1957	zustrph(i)=zustrph(i)+(u_seri(i,k)-u(i,k))/dtime*
1958	c (paprs(i,k)-paprs(i,k+1))/rg
1959	zvstrph(i)=zvstrph(i)+(v_seri(i,k)-v(i,k))/dtime*
1960	c (paprs(i,k)-paprs(i,k+1))/rg
1961	zustrcl(i)=zustrcl(i)+d_u_vdf(i,k)*
1962	c (paprs(i,k)-paprs(i,k+1))/rg
1963	zvstrcl(i)=zvstrcl(i)+d_v_vdf(i,k)*
1964	c (paprs(i,k)-paprs(i,k+1))/rg
1965	ENDDO
1966	ENDDO
1967
1968	CALL aaam_bud (27,klon,klev,rjourvrai,gmtime*RDAY,
1969	C ra,rg,romega,
1970	C latitude_deg,longitude_deg,pphis,
1971	C zustrdr,zustrli,zustrcl,
1972	C zvstrdr,zvstrli,zvstrcl,
1973	C paprs,u,v)
1974
1975	CCMODFIN FLOTT
1976	endif !bilansmc
1977
1978	c====================================================================
1979	c====================================================================
1980	c Calculer le transport de l'eau et de l'energie (diagnostique)
1981	c
1982	c A REVOIR POUR VENUS...
1983	c
1984	c CALL transp (paprs,ftsol,
1985	c e t_seri, q_seri, u_seri, v_seri, zphi,
1986	c s ve, vq, ue, uq)
1987	c
1988	c====================================================================
1989	c+jld ec_conser
1990	DO k = 1, klev
1991	DO i = 1, klon
1992	d_t_ec(i,k)=0.5/cpnew(i,k)
1993	$ (u(i,k)2+v(i,k)2-u_seri(i,k)2-v_seri(i,k)*2)
1994	t_seri(i,k)=t_seri(i,k)+d_t_ec(i,k)
1995	d_t_ec(i,k) = d_t_ec(i,k)/dtime
1996	END DO
1997	END DO
1998	do i=1,klon
1999	flux_ec(i,1) = 0.0
2000	do j=2,klev
2001	flux_ec(i,j) = flux_ec(i,j-1)
2002	. +cpnew(i,j-1)/RGd_t_ec(i,j-1)(paprs(i,j-1)-paprs(i,j))
2003	enddo
2004	enddo
2005
2006	c-jld ec_conser
2007	c====================================================================
2008	IF (if_ebil.ge.1) THEN
2009	ztit='after physic'
2010	CALL diagetpq(cell_area,ztit,ip_ebil,1,1,dtime
2011	e , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
2012	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2013	C Comme les tendances de la physique sont ajoute dans la dynamique,
2014	C on devrait avoir que la variation d'entalpie par la dynamique
2015	C est egale a la variation de la physique au pas de temps precedent.
2016	C Donc la somme de ces 2 variations devrait etre nulle.
2017	call diagphy(cell_area,ztit,ip_ebil
2018	e , topsw, toplw, solsw, sollw, sens
2019	e , zero_v, zero_v, zero_v, ztsol
2020	e , d_h_vcol, d_qt, d_ec
2021	s , fs_bound, fq_bound )
2022	C
2023	d_h_vcol_phy=d_h_vcol
2024	C
2025	END IF
2026	C
2027	c=======================================================================
2028	c SORTIES
2029	c=======================================================================
2030
2031	c Convertir les incrementations en tendances
2032	c
2033	DO k = 1, klev
2034	DO i = 1, klon
2035	d_u(i,k) = ( u_seri(i,k) - u(i,k) ) / dtime
2036	d_v(i,k) = ( v_seri(i,k) - v(i,k) ) / dtime
2037	d_t(i,k) = ( t_seri(i,k) - t(i,k) ) / dtime
2038	ENDDO
2039	ENDDO
2040	c
2041	DO iq = 1, nqmax
2042	DO k = 1, klev
2043	DO i = 1, klon
2044	d_qx(i,k,iq) = ( tr_seri(i,k,iq) - qx(i,k,iq) ) / dtime
2045	ENDDO
2046	ENDDO
2047	ENDDO
2048
2049	c mise à jour rho,mmean pour sorties
2050	if(callthermos) then
2051	call concentrations2(pplay,t_seri,tr_seri, nqmax)
2052	endif
2053
2054	c calcul vitesse verticale en m/s
2055	DO k = 1, klev
2056	DO i = 1, klon
2057	vertwind(i,k) = -omega(i,k)/(rho(i,k)*RG)
2058	END DO
2059	END DO
2060
2061	c------------------------
2062	c Calcul moment cinetique
2063	c------------------------
2064	c TEST VENUS...
2065	c mangtot = 0.0
2066	c DO k = 1, klev
2067	c DO i = 1, klon
2068	c mang(i,k) = RA*cos(latitude(i))
2069	c . (u_seri(i,k)+RAcos(latitude(i))*ROMEGA)
2070	c . cell_area(i)(paprs(i,k)-paprs(i,k+1))/RG
2071	c mangtot=mangtot+mang(i,k)
2072	c ENDDO
2073	c ENDDO
2074	c print*,"Moment cinetique total = ",mangtot
2075	c
2076	c------------------------
2077	c
2078	c Sauvegarder les valeurs de t et u a la fin de la physique:
2079	c
2080	DO k = 1, klev
2081	DO i = 1, klon
2082	u_ancien(i,k) = u_seri(i,k)
2083	t_ancien(i,k) = t_seri(i,k)
2084	ENDDO
2085	ENDDO
2086	c
2087	c=============================================================
2088	c Ecriture des sorties
2089	c=============================================================
2090	#ifndef MESOSCALE
2091	#ifdef CPP_IOIPSL
2092
2093	#ifdef histhf
2094	#include "write_histhf.h"
2095	#endif
2096
2097	#ifdef histday
2098	#include "write_histday.h"
2099	#endif
2100
2101	#ifdef histmth
2102	#include "write_histmth.h"
2103	#endif
2104
2105	#ifdef histins
2106	#include "write_histins.h"
2107	#endif
2108
2109	#endif
2110
2111	! XIOS outputs
2112	! This can be done ANYWHERE in the physics routines !
2113
2114	#ifdef CPP_XIOS
2115	! Send fields to XIOS: (NB these fields must also be defined as
2116	! <field id="..." /> in context_lmdz_physics.xml to be correctly used)
2117
2118	! 2D fields
2119
2120	CALL send_xios_field("phis",pphis)
2121	cell_area_out(:)=cell_area(:)
2122	if (is_north_pole_phy) cell_area_out(1)=cell_area(1)/nbp_lon
2123	if (is_south_pole_phy) cell_area_out(klon)=cell_area(klon)/nbp_lon
2124	CALL send_xios_field("aire",cell_area_out)
2125	CALL send_xios_field("tsol",ftsol)
2126	CALL send_xios_field("psol",paprs(:,1))
2127	CALL send_xios_field("cdragh",cdragh)
2128	CALL send_xios_field("cdragm",cdragm)
2129
2130	CALL send_xios_field("tops",topsw)
2131	CALL send_xios_field("topl",toplw)
2132	CALL send_xios_field("sols",solsw)
2133	CALL send_xios_field("soll",sollw)
2134
2135	! 3D fields
2136
2137	CALL send_xios_field("temp",t_seri)
2138	CALL send_xios_field("pres",pplay)
2139	CALL send_xios_field("geop",zphi)
2140	CALL send_xios_field("vitu",u_seri)
2141	c VENUS: regardee a l envers!!!!!!!!!!!!!!!
2142	CALL send_xios_field("vitv",-1.*v_seri)
2143	CALL send_xios_field("vitw",omega)
2144	CALL send_xios_field("vitwz",vertwind)
2145	CALL send_xios_field("Kz",ycoefh)
2146	CALL send_xios_field("mmean",mmean)
2147	CALL send_xios_field("rho",rho)
2148	CALL send_xios_field("BV2",zn2)
2149
2150	CALL send_xios_field("dudyn",d_u_dyn)
2151	CALL send_xios_field("duvdf",d_u_vdf)
2152	c VENUS: regardee a l envers!!!!!!!!!!!!!!!
2153	CALL send_xios_field("dvvdf",-1.*d_v_vdf)
2154	CALL send_xios_field("duajs",d_u_ajs)
2155	CALL send_xios_field("dugwo",d_u_oro)
2156	CALL send_xios_field("dugwno",d_u_hin)
2157	CALL send_xios_field("dvgwno",-1.*d_v_hin)
2158	CALL send_xios_field("dumolvis",d_u_molvis)
2159	c VENUS: regardee a l envers!!!!!!!!!!!!!!!
2160	CALL send_xios_field("dvmolvis",-1.*d_v_molvis)
2161	CALL send_xios_field("dtdyn",d_t_dyn)
2162	CALL send_xios_field("dtphy",d_t)
2163	CALL send_xios_field("dtvdf",d_t_vdf)
2164	CALL send_xios_field("dtajs",d_t_ajs)
2165	CALL send_xios_field("dtswr",dtsw)
2166	CALL send_xios_field("dtswrNLTE",d_t_nirco2)
2167	CALL send_xios_field("dtswrLTE",heat)
2168	CALL send_xios_field("dtlwr",dtlw)
2169	CALL send_xios_field("dtlwrNLTE",d_t_nlte)
2170	CALL send_xios_field("dtlwrLTE",-1.*cool)
2171	CALL send_xios_field("dteuv",d_t_euv)
2172	CALL send_xios_field("dtcond",d_t_conduc)
2173	CALL send_xios_field("dtec",d_t_ec)
2174
2175	CALL send_xios_field("SWnet",swnet(:,1:klev))
2176	CALL send_xios_field("LWnet",lwnet(:,1:klev))
2177	CALL send_xios_field("fluxvdf",fluxt)
2178	CALL send_xios_field("fluxdyn",flux_dyn)
2179	CALL send_xios_field("fluxajs",flux_ajs)
2180	CALL send_xios_field("fluxec",flux_ec)
2181
2182	! when using tracers
2183
2184	if (iflag_trac == 1) then
2185
2186	if (ok_aoa) then
2187	call send_xios_field("age",age)
2188	call send_xios_field("aoa",tr_seri(:,:,i_aoa))
2189	endif
2190
2191	! production and destruction rate, cm-3.s-1
2192	! Beware of the context*.xml file !!
2193	if ((tr_scheme == 3) .and. (ok_chem)) THEN
2194	do iq = 1,nqmax - nmicro
2195	if ((iq.eq.i_o).or.(iq.eq.i_co)) THEN
2196	call send_xios_field("prod_"//tname(iq),
2197	$ prod_tr(:,:,iq))
2198	call send_xios_field("loss_"//tname(iq),
2199	$ loss_tr(:,:,iq))
2200	end if
2201	!if (iq.eq.i_o) then
2202	! call send_xios_field('prod_o', prod_tr(:,:,iq))
2203	! call send_xios_field('loss_o', loss_tr(:,:,iq))
2204	!end if
2205	!if (iq.eq.i_co) then
2206	! call send_xios_field('prod_co', prod_tr(:,:,iq))
2207	! call send_xios_field('loss_co', loss_tr(:,:,iq))
2208	!end if
2209	end do
2210	end if
2211
2212	! tracers in gas phase, volume mixing ratio
2213
2214	do iq = 1,nqmax - nmicro
2215	call send_xios_field(tname(iq),
2216	$ tr_seri(:,:,iq)*mmean(:,:)/m_tr(iq))
2217	end do
2218
2219	! tracers in gas phase, column densities
2220
2221	do iq = 1,nqmax - nmicro
2222	col_dens_tr(:,iq)=0.
2223	if (type_tr(iq).eq.1) THEN
2224	do k = 1, klev
2225	col_dens_tr(:,iq) = col_dens_tr(:,iq) +
2226	$ tr_seri(:,k,iq)* (paprs(:,k)-paprs(:,k+1)) / RG
2227	end do
2228	call send_xios_field("col_"//tname(iq),col_dens_tr(:,iq))
2229	end if
2230	end do
2231
2232	! tracers in liquid phase, volume mixing ratio
2233
2234	if ((tr_scheme == 3) .and. (cl_scheme == 1)) THEN
2235	call send_xios_field(tname(i_h2oliq),
2236	$ tr_seri(:,:,i_h2oliq)*mmean(:,:)/m_tr(i_h2oliq))
2237	call send_xios_field(tname(i_h2so4liq),
2238	$ tr_seri(:,:,i_h2so4liq)*mmean(:,:)/m_tr(i_h2so4liq))
2239	if (ok_sedim) then
2240	call send_xios_field("Fsedim",fsedim(:,1:klev))
2241	end if
2242	end if
2243
2244	! aeronomical emissions
2245
2246	call send_xios_field("no_emis",no_emission)
2247	call send_xios_field("o2_emis",o2_emission)
2248
2249	! chemical iterations
2250
2251	if (tr_scheme.eq.3) call send_xios_field("iter",real(iter))
2252
2253	end if
2254
2255	IF (callthermos .and. ok_chem) THEN
2256	CALL send_xios_field("d_qmoldifCO2",d_q_moldif(:,:,i_co2))
2257	CALL send_xios_field("d_qmoldifO3p",d_q_moldif(:,:,i_o))
2258	CALL send_xios_field("d_qmoldifN2",d_q_moldif(:,:,i_n2))
2259	ENDIF
2260
2261	!! DEBUG
2262	! if (is_master) print*,"itauphy=",itap
2263	! if (itap.eq.10) lafin=.true.
2264
2265	if (lafin.and.is_omp_master) then
2266	write(,) "physiq: call xios_context_finalize"
2267	call xios_context_finalize
2268	endif
2269
2270	#endif
2271	#else
2272	! Outputs MESOSCALE
2273	CALL allocate_comm_wrf(klon,klev)
2274	comm_HR_SW(1:klon,1:klev) = dtsw(1:klon,1:klev)
2275	comm_HR_LW(1:klon,1:klev) = dtlw(1:klon,1:klev)
2276	comm_DT_RAD(1:klon,1:klev) = d_t_rad(1:klon,1:klev)
2277	IF (turb_resolved) THEN
2278	open(17,file='hrdyn.txt',form='formatted',status='old')
2279	rewind(17)
2280	DO k=1,klev
2281	read(17,*) dt_dyn(k)
2282	ENDDO
2283	close(17)
2284
2285	do i=1,klon
2286	d_t(i,:)=d_t(i,:)+dt_dyn(:)
2287	comm_HR_DYN(i,:) = dt_dyn(:)
2288	enddo
2289	ELSE
2290	comm_HR_DYN(1:klon,1:klev) = d_t_dyn(1:klon,1:klev)
2291	comm_DT_VDF(1:klon,1:klev) = d_t_vdf(1:klon,1:klev)
2292	comm_DT_AJS(1:klon,1:klev) = d_t_ajs(1:klon,1:klev)
2293	ENDIF
2294	comm_DT(1:klon,1:klev)=d_t(1:klon,1:klev)
2295	#endif
2296
2297
2298	c====================================================================
2299	c Si c'est la fin, il faut conserver l'etat de redemarrage
2300	c====================================================================
2301	c
2302	IF (lafin) THEN
2303	itau_phy = itau_phy + itap
2304	CALL phyredem ("restartphy.nc")
2305
2306	c--------------FLOTT
2307	CMODEB LOTT
2308	C FERMETURE DU FICHIER FORMATTE CONTENANT LES COMPOSANTES
2309	C DU BILAN DE MOMENT ANGULAIRE.
2310	if (bilansmc.eq.1) then
2311	write(,)'Fermeture de aaam_bud.out (FL Vous parle)'
2312	close(27)
2313	close(28)
2314	endif !bilansmc
2315	CMODFIN
2316	c-------------
2317	c--------------SLEBONNOIS
2318	C FERMETURE DES FICHIERS FORMATTES CONTENANT LES POSITIONS ET VITESSES
2319	C DES BALLONS
2320	if (ballons.eq.1) then
2321	write(,)'Fermeture des ballons*.out'
2322	close(30)
2323	close(31)
2324	close(32)
2325	close(33)
2326	close(34)
2327	endif !ballons
2328	c-------------
2329	ENDIF
2330
2331	END SUBROUTINE physiq
2332
2333	END MODULE physiq_mod
2334

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

Download in other formats:

Original Format