Context Navigation

physiq_mod.F @ 3556

Last change on this file since 3556 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

Rev	Line
[3]	1	!
[1549]	2	! $Id: $
[3]	3	!
[1549]	4	MODULE physiq_mod
	5
	6	IMPLICIT NONE
	7
	8	CONTAINS
	9
[3]	10	SUBROUTINE physiq (nlon,nlev,nqmax,
	11	. debut,lafin,rjourvrai,gmtime,pdtphys,
[97]	12	. paprs,pplay,ppk,pphi,pphis,presnivs,
[3]	13	. u,v,t,qx,
[2125]	14	. flxmw,
[3]	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
[953]	40	c gmtime--input-R-fraction de la journee (0 a 1)
[3]	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)
[1301]	53	c flxmw---input-R-flux de masse vertical en kg/s
[3]	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======================================================================
[101]	61	USE ioipsl
[849]	62	! USE histcom ! not needed; histcom is included in ioipsl
[101]	63	use dimphy
[1545]	64	USE geometry_mod,only: longitude, latitude, ! in radians
	65	& longitude_deg,latitude_deg, ! in degrees
	66	& cell_area,dx,dy
[892]	67	USE phys_state_var_mod ! Variables sauvegardees de la physique
[1621]	68	USE cpdet_phy_mod, only: cpdet, t2tpot
[1310]	69	USE chemparam_mod
[3451]	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
[1310]	72	USE conc
[2464]	73	USE param_v4_h
[1310]	74	USE compo_hedin83_mod2
[2486]	75	use radlwsw_newtoncool_mod, only: radlwsw_newtoncool
[1506]	76	! use ieee_arithmetic
[1524]	77	use time_phylmdz_mod, only: annee_ref, day_ref, itau_phy
[1530]	78	use mod_grid_phy_lmdz, only: nbp_lon
[1621]	79	use infotrac_phy, only: iflag_trac, tname, ttext
[1642]	80	use vertical_layers_mod, only: pseudoalt
[1723]	81	use turb_mod, only : sens, turb_resolved
[2580]	82	use nonoro_gwd_ran_mod, only: nonoro_gwd_ran
[2259]	83	use sed_and_prod_mad, only: aer_sedimentation, drop_sedimentation
[2836]	84	use iono_h, only: temp_elect, temp_ion
[1642]	85	#ifdef CPP_XIOS
	86	use xios_output_mod, only: initialize_xios_output,
	87	& update_xios_timestep,
	88	& send_xios_field
[1682]	89	use wxios, only: wxios_context_init, xios_context_finalize
[1642]	90	#endif
[1723]	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,
[2193]	99	& is_south_pole_phy,
	100	& is_master
[1723]	101	#endif
[101]	102	IMPLICIT none
	103	c======================================================================
	104	c CLEFS CPP POUR LES IO
	105	c =====================
[1723]	106	#ifndef MESOSCALE
[101]	107	c#define histhf
	108	#define histday
	109	#define histmth
	110	#define histins
[1723]	111	#endif
[101]	112	c======================================================================
[3]	113	#include "dimsoil.h"
	114	#include "clesphys.h"
	115	#include "iniprint.h"
	116	#include "timerad.h"
[815]	117	#include "tabcontrol.h"
[1310]	118	#include "nirdata.h"
	119	#include "hedin.h"
[3]	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
[1301]	161	REAL flxmw(klon,klev)
[1665]	162
[3]	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
[2193]	181
	182	integer,save :: itap ! physics counter
[3]	183	REAL delp(klon,klev) ! epaisseur d'une couche
[2523]	184	REAL omega(klon,klev) ! vitesse verticale en Pa/s (+ downward)
	185	REAL vertwind(klon,klev) ! vitesse verticale en m/s (+ upward)
[3]	186
	187	INTEGER igwd,idx(klon),itest(klon)
[2193]	188
[3]	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)
[2047]	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)
[3]	199
	200	c Pour calcul GW drag oro et nonoro: CALCUL de N2:
[1665]	201	real zdtlev(klon,klev),zdzlev(klon,klev)
	202	real ztlev(klon,klev)
[3]	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
[105]	225	REAL,save,allocatable :: source(:,:)
[3]	226	REAL ycoefh(klon,klev) ! coef d'echange pour phytrac
[3535]	227	REAL :: kzz_p(klev) ! coef d'echange pour phytrac pour le 1D
[3]	228	REAL yu1(klon) ! vents dans la premiere couche U
	229	REAL yv1(klon) ! vents dans la premiere couche V
	230
[1723]	231	REAL dsens(klon) ! derivee chaleur sensible
[3]	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
[1305]	237	REAL Fsedim(klon,klev+1) ! Flux de sedimentation (kg.m-2)
[3]	238
	239	c======================================================================
	240	c
	241	c Declaration des procedures appelees
	242	c
[2836]	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
[3]	253	EXTERNAL printflag
	254	EXTERNAL zenang
	255	EXTERNAL diagetpq
	256	EXTERNAL conf_phys
	257	EXTERNAL diagphy
	258	EXTERNAL mucorr
[1310]	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
[1442]	270	EXTERNAL moldiff_red
	271
[3]	272	c
	273	c Variables locales
	274	c
	275	CXXX PB
[2836]	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
[3]	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
[2836]	284	REAL tmpout(klon,klev) ! [K/s]
[3]	285	c
	286	REAL dist, rmu0(klon), fract(klon)
	287	REAL zdtime, zlongi
	288	c
[1718]	289	INTEGER i, k, iq, ig, j, ll, ilon, ilat, ilev, isoil
[3]	290	c
	291	REAL zphi(klon,klev)
[1301]	292	REAL zzlev(klon,klev+1),zzlay(klon,klev),z1,z2
[2686]	293	REAL tlaymean ! valeur temporaire pour calculer zzlay
[1310]	294	real tsurf(klon)
[1301]	295
[1310]	296	c va avec nlte_tcool
	297	INTEGER ierr_nlte
	298	REAL varerr
	299
[2193]	300	! photochemistry
	301
	302	integer :: chempas
[3035]	303	real :: zctime
[2193]	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
[3]	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
[1310]	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]
[2836]	340	real d_u_molvis(klon,klev) ! [m/s] /s
	341	real d_v_molvis(klon,klev) ! [m/s] /s
[1310]	342
[1442]	343	c Tendencies due to molecular diffusion
	344	real d_q_moldif(klon,klev,nqmax)
	345
[3035]	346	c Tendencies due to ambipolar ion diffusion
	347	real d_q_iondif(klon,klev,nqmax)
	348
[3]	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
[1305]	368	REAL :: tr_seri(klon,klev,nqmax)
[2464]	369	REAL :: tr_hedin(klon,klev,nqmax)
[1305]	370	REAL :: d_tr(klon,klev,nqmax)
[2523]	371	c pour sorties
	372	REAL :: col_dens_tr(klon,nqmax)
[2580]	373	REAL,allocatable,save :: prod_tr(:,:,:)
	374	REAL,allocatable,save :: loss_tr(:,:,:)
[1305]	375
[3]	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
[97]	380	REAL zsto, zout, zsto1, zsto2, zero
[3]	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
[2580]	415	c ALBEDO VARIATIONS (VCD)
	416	REAL factAlb
[3]	417	c TEST VENUS...
	418	REAL mang(klon,klev) ! moment cinetique
	419	REAL mangtot ! moment cinetique total
	420
[1572]	421	c cell_area for outputs in hist*
	422	REAL cell_area_out(klon)
[1723]	423	#ifdef MESOSCALE
	424	REAL :: dt_dyn(klev)
	425	#endif
[3]	426	c Declaration des constantes et des fonctions thermodynamiques
	427	c
	428	#include "YOMCST.h"
	429
	430	c======================================================================
[2464]	431	c======================================================================
[3]	432	c INITIALISATIONS
[2464]	433	c======================================================================
[3]	434
	435	modname = 'physiq'
	436	ok_sync=.TRUE.
	437
[119]	438	bilansmc = 0
[3]	439	ballons = 0
[892]	440	! NE FONCTIONNENT PAS ENCORE EN PARALLELE !!!
[1723]	441	#ifndef MESOSCALE
[892]	442	if (is_parallel) then
	443	bilansmc = 0
	444	ballons = 0
	445	endif
[1723]	446	#endif
[3]	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
[2464]	458	c======================================================================
[3451]	459	c DONE ONLY AT FIRST CALL
[3]	460	c========================
[3451]	461	IF (debut) THEN
[105]	462	allocate(source(klon,nqmax))
[2580]	463	allocate(prod_tr(klon,klev,nqmax))
	464	allocate(loss_tr(klon,klev,nqmax))
[2795]	465	allocate(no_emission(klon,klev))
	466	allocate(o2_emission(klon,klev))
[3]	467
[1682]	468	#ifdef CPP_XIOS
	469	! Initialize XIOS context
	470	write(,) "physiq: call wxios_context_init"
	471	CALL wxios_context_init
	472	#endif
	473
[2135]	474	! The call to suphec is now done in iniphysiq_mod (interface)
	475	! CALL suphec ! initialiser constantes et parametres phys.
[3]	476
	477	IF (if_ebil.ge.1) d_h_vcol_phy=0.
[3451]	478	!
	479	! load flag and parameter values from physiq.def
	480	!
[3]	481	call conf_phys(ok_journe, ok_mensuel,
	482	. ok_instan,
	483	. if_ebil)
	484
[2193]	485	call phys_state_var_init(nqmax)
[3]	486	c
[1310]	487	c Initialising Hedin model for upper atm
	488	c (to be revised when coupled to chemistry) :
	489	call conc_init
[1723]	490
[2193]	491	! initialise physics counter
	492
	493	itap = 0
	494
[1723]	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
[3]	519	c
	520	c Lecture startphy.nc :
	521	c
[892]	522	CALL phyetat0 ("startphy.nc")
[1718]	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
[1723]	531	#endif
[3]	532
[815]	533	c dtime est defini dans tabcontrol.h et lu dans startphy
[150]	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 '
[1719]	540	c call abort_physic(modname,abort_message,1)
[152]	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----------------
[150]	546	dtime=pdtphys
	547	ENDIF
	548
[2193]	549	radpas = NINT(RDAY/pdtphys/nbapp_rad)
[3]	550
[815]	551	CALL printflag( ok_journe,ok_instan )
[1642]	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
[3]	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
[1442]	602	source(:,:) = 0.0 ! pas de source, pour l'instant
[3]	603	c---------
[1310]	604
	605	c---------
	606	c INITIALIZE THERMOSPHERIC PARAMETERS
	607	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	608
	609	if (callthermos) then
[2464]	610	call fill_data_thermos
	611	call allocate_param_thermos(klev)
	612	call param_read_e107
[1310]	613	endif
	614
	615	c Initialisation (recomputed in concentration2)
	616	do ig=1,klon
	617	do j=1,klev
[2788]	618	rnew(ig,j)=RD
[1665]	619	cpnew(ig,j)=cpdet(t(ig,j))
[1442]	620	mmean(ig,j)=RMD
[1310]	621	akknew(ig,j)=1.e-4
[2788]	622	rho(ig,j)=pplay(ig,j)/(rnew(ig,j)*t(ig,j))
[1310]	623	enddo
	624
	625	enddo
	626
	627	IF(callthermos.or.callnlte.or.callnirco2) THEN
	628	call compo_hedin83_init2
	629	ENDIF
[1442]	630	if (callnlte.and.nltemodel.eq.2) call nlte_setup
	631	if (callnirco2.and.nircorr.eq.1) call nir_leedat
[1310]	632	c---------
	633
[3]	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'
[1719]	641	call abort_physic(modname,abort_message,1)
[3]	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'
[1719]	647	call abort_physic(modname,abort_message,1)
[3]	648	ENDIF
	649	c
[808]	650	IF (dtimeREAL(radpas).GT.(RDAY0.25).AND.cycle_diurne)
[3]	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'
[1719]	655	call abort_physic(modname,abort_message,1)
[3]	656	ENDIF
	657	c
	658	WRITE(lunout,*)"Clef pour la convection seche, iflag_ajs=",
	659	. iflag_ajs
	660	c
[1120]	661	ecrit_mth = NINT(RDAY/dtime*ecriphy) ! tous les ecritphy jours
[3]	662	IF (ok_mensuel) THEN
	663	WRITE(lunout,*)'La frequence de sortie mensuelle est de ',
	664	. ecrit_mth
	665	ENDIF
[97]	666
[3]	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
[97]	672
[1120]	673	ecrit_ins = NINT(RDAY/dtime*ecriphy) ! Fraction de jour reglable
[3]	674	IF (ok_instan) THEN
	675	WRITE(lunout,*)'La frequence de sortie instant. est de ',
	676	. ecrit_ins
	677	ENDIF
	678
[3035]	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
[3]	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
[2683]	722	! initialisation of microphysical and chemical parameters
[1661]	723
[2683]	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
[1661]	735
[2683]	736	! number of microphysical tracers
[2193]	737
	738	nmicro = 0
	739	if (ok_cloud .and. (cl_scheme == 1)) nmicro = 2
	740	if (ok_cloud .and. (cl_scheme == 2)) nmicro = 12
[1661]	741
[2683]	742	! initialise chemical parameters. includes the indexation of microphys tracers
[1305]	743
[2683]	744	if (ok_chem .or. cl_scheme == 2) then
	745	call chemparam_ini()
	746	end if
[3]	747
[3451]	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
[2683]	759	! initialise cloud parameters (for cl_scheme = 1)
[1442]	760
[2683]	761	if (ok_cloud .and. (cl_scheme == 1)) then
[3323]	762	call cloud_ini(nlon,nlev,nb_mode)
[2683]	763	end if
	764
[2286]	765	! initialise mmean
	766
[2683]	767	if (callthermos) then
[2464]	768	call concentrations2(pplay,t,qx,nqmax)
[2683]	769	end if
[2286]	770
[2464]	771	c========================
[3]	772	ENDIF ! debut
[2464]	773	c========================
[2193]	774
[2464]	775	c======================================================================
[1642]	776	! ------------------------------------------------------
	777	! Initializations done at every physical timestep:
	778	! ------------------------------------------------------
[3]	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'
[1543]	823	CALL diagetpq(cell_area,ztit,ip_ebil,1,1,dtime
[3]	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.
[1543]	830	call diagphy(cell_area,ztit,ip_ebil
[3]	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====================================================================
[1642]	838	c XIOS outputs
	839
	840	#ifdef CPP_XIOS
	841	! update XIOS time/calendar
	842	call update_xios_timestep
	843	#endif
	844
	845	c====================================================================
[3]	846	c Diagnostiquer la tendance dynamique
	847	c
	848	IF (ancien_ok) THEN
	849	DO k = 1, klev
[2836]	850	DO i = 1, klon
[3]	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
[2836]	853	ENDDO
[3]	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)
[2686]	861	. +cpnew(i,j-1)/RGd_t_dyn(i,j-1)(paprs(i,j-1)-paprs(i,j))
[3]	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====================================================================
[1301]	875
[2523]	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)
[1301]	883
[2047]	884	c======
	885	c GEOP CORRECTION
[3]	886	c
	887	c Ajouter le geopotentiel du sol:
	888	c
	889	DO k = 1, klev
[2836]	890	DO i = 1, klon
[3]	891	zphi(i,k) = pphi(i,k) + pphis(i)
[2836]	892	ENDDO
[3]	893	ENDDO
[1301]	894
[2047]	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
[2686]	916	c calcul de l'altitude aux niveaux intercouches
[1301]	917	c ponderation des altitudes au niveau des couches en dp/p
	918
	919	DO i=1,klon
[2686]	920	zzlay(i,1)=zphi(i,1)/RG ! [m]
[1310]	921	zzlev(i,1)=pphis(i)/RG ! [m]
[1301]	922	ENDDO
	923	DO k=2,klev
[2686]	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
[1301]	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
[3]	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
[2464]	966	c=======================
[3]	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
[2464]	974	dist = 0.72333 ! en UA
[3]	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
[808]	982	zdtime=dtime*REAL(radpas) ! pas de temps du rayonnement (s)
[1545]	983	CALL zenang(zlongi,gmtime,zdtime,latitude_deg,longitude_deg,
	984	& rmu0,fract)
[3]	985	ELSE
[1545]	986	call mucorr(klon,zlongi,latitude_deg,rmu0,fract)
[3]	987	ENDIF
	988
[2193]	989	! print fraction of venus day
[3]	990
[2193]	991	if (is_master) then
	992	print*, 'gmtime = ', gmtime
	993	end if
	994
[2464]	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]
[2836]	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]
[2464]	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
[3]	1045	if (iflag_trac.eq.1) then
[3451]	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
[2193]	1053	!====================================================================
	1054	! Case 1: pseudo-chemistry with relaxation toward fixed profile
[2464]	1055	!=========
[1160]	1056	if (tr_scheme.eq.1) then
[1305]	1057
[1160]	1058	call phytrac_relax (debut,lafin,nqmax,
	1059	I nlon,nlev,dtime,pplay,
[3]	1060	O tr_seri)
[1160]	1061
	1062	elseif (tr_scheme.eq.2) then
[2193]	1063	!====================================================================
[1160]	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.
[2464]	1068	!=========
[1305]	1069
[2464]	1070	call phytrac_emiss (debut,lafin,nqmax,
[1160]	1071	I nlon,nlev,dtime,paprs,
[1545]	1072	I latitude_deg,longitude_deg,
[1160]	1073	O tr_seri)
[3]	1074
[2193]	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.
[2262]	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
[2464]	1084	!=========
[1305]	1085
[3035]	1086
[2193]	1087	chempas = nint(rday/pdtphys/nbapp_chem)
[2262]	1088	zctime = dtime*real(chempas) ! chemical timestep
[1305]	1089
[2193]	1090	if (mod(itap,chempas) == 0) then ! <------- start of chemistry supercycling
[1305]	1091
[2193]	1092	! photochemistry and microphysics
[1661]	1093
[2193]	1094	call phytrac_chimie(debut,
	1095	$ gmtime,
	1096	$ nqmax,
	1097	$ klon,
	1098	$ latitude_deg,
	1099	$ longitude_deg,
[2780]	1100	$ zzlay,
[2193]	1101	$ nlev,
	1102	$ zctime,
	1103	$ t_seri,
	1104	$ pplay,
	1105	$ tr_seri,
	1106	$ d_tr_chem,
[2580]	1107	$ iter,
	1108	$ prod_tr,
[2795]	1109	$ loss_tr,
	1110	$ no_emission,
	1111	$ o2_emission)
[2193]	1112
[1305]	1113	if (ok_sedim) then
[2193]	1114	if (cl_scheme == 1) then
[1305]	1115
[2193]	1116	! sedimentation for simplified microphysics
[1661]	1117
[1723]	1118	#ifndef MESOSCALE
[3323]	1119	call new_cloud_sedim(nb_mode,
	1120	$ klon,
[2193]	1121	$ nlev,
	1122	$ zctime,
	1123	$ pplay,
	1124	$ paprs,
	1125	$ t_seri,
	1126	$ tr_seri,
[2200]	1127	$ d_tr_chem,
[2193]	1128	$ d_tr_sed(:,:,1:2),
	1129	$ nqmax,
	1130	$ Fsedim)
[1305]	1131
[2193]	1132	! test to avoid nans
[1442]	1133
[2193]	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
[1305]	1147
[2203]	1148	! tendency due to condensation and sedimentation
[1305]	1149
[2203]	1150	d_tr_sed(:,:,1:2) = d_tr_sed(:,:,1:2)/zctime
[2193]	1151	Fsedim(:,1:klev) = Fsedim(:,1:klev)/zctime
	1152	Fsedim(:,klev+1) = 0.
[1661]	1153
[2193]	1154	else if (cl_scheme == 2) then
[1661]	1155
[2193]	1156	! sedimentation for detailed microphysics
[1661]	1157
[2193]	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
[1661]	1176	d_tr_sed(i,:,i_m0_mode1drop)= d_tr_sed(i,:,i_m0_mode1drop)
[2193]	1177	$ + d_drop_sed
[1661]	1178	d_tr_sed(i,:,i_m0_mode1ccn) = d_tr_sed(i,:,i_m0_mode1ccn)
[2193]	1179	$ + d_ccn_sed(:,1)
[1661]	1180	d_tr_sed(i,:,i_m3_mode1ccn) = d_tr_sed(i,:,i_m3_mode1ccn)
[2193]	1181	$ + d_ccn_sed(:,2)
[1661]	1182	d_tr_sed(i,:,i_m3_mode1sa) = d_tr_sed(i,:,i_m3_mode1sa)
[2193]	1183	$ + d_liq_sed(:,1)
[1661]	1184	d_tr_sed(i,:,i_m3_mode1w) = d_tr_sed(i,:,i_m3_mode1w)
[2193]	1185	$ + d_liq_sed(:,2)
[1661]	1186
[2193]	1187	! mode 2
[1661]	1188
[2193]	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
[1661]	1201	d_tr_sed(i,:,i_m0_mode2drop)= d_tr_sed(i,:,i_m0_mode2drop)
[2193]	1202	$ + d_drop_sed
[1661]	1203	d_tr_sed(i,:,i_m0_mode2ccn) = d_tr_sed(i,:,i_m0_mode2ccn)
[2193]	1204	$ + d_ccn_sed(:,1)
[1661]	1205	d_tr_sed(i,:,i_m3_mode2ccn) = d_tr_sed(i,:,i_m3_mode2ccn)
[2193]	1206	$ + d_ccn_sed(:,2)
[1661]	1207	d_tr_sed(i,:,i_m3_mode2sa) = d_tr_sed(i,:,i_m3_mode2sa)
[2193]	1208	$ + d_liq_sed(:,1)
[1661]	1209	d_tr_sed(i,:,i_m3_mode2w) = d_tr_sed(i,:,i_m3_mode2w)
[2193]	1210	$ + d_liq_sed(:,2)
[1661]	1211
[2193]	1212	! aer
[1661]	1213
[2193]	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
[1661]	1224
[2193]	1225	! tendency due to sedimentation
	1226
[2194]	1227	do iq = nqmax-nmicro+1,nqmax
[2193]	1228	d_tr_sed(:,:,iq) = d_tr_sed(:,:,iq)/zctime
	1229	end do
[1723]	1230	#endif
[2193]	1231	end if ! cl_scheme
[1661]	1232
[2203]	1233	! update gaseous tracers (chemistry)
[1661]	1234
[2203]	1235	do iq = 1, nqmax - nmicro
[2464]	1236	tr_seri(:,:,iq) = max(tr_seri(:,:,iq)
	1237	$ + d_tr_chem(:,:,iq)*zctime,1.e-30)
[2203]	1238	end do
[1305]	1239
[2203]	1240	! update condensed tracers (condensation + sedimentation)
[1305]	1241
[2194]	1242	if (cl_scheme == 1) then
[2203]	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)
[2194]	1247	else if (cl_scheme == 2) then
[2203]	1248	do iq = nqmax-nmicro+1,nqmax
[2464]	1249	tr_seri(:,:,iq) = max(tr_seri(:,:,iq)
	1250	$ + d_tr_sed(:,:,iq)*zctime,1.e-30)
[2203]	1251	end do
[2194]	1252	end if ! cl_scheme
[2203]	1253
[2194]	1254	end if ! ok_sedim
[2203]	1255	end if ! mod(itap,chempas) <------- end of chemistry supercycling
	1256
[2464]	1257	!=========
[2194]	1258	! End Case 3: Full chemistry and/or clouds.
	1259	!====================================================================
[2193]	1260
	1261	end if ! tr_scheme
	1262	end if ! iflag_trac
	1263
[3]	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
[2135]	1271	if (.not. ok_clmain) then
[3]	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)
[808]	1299
[2135]	1300	if (physideal) then
	1301	CALL clmain_ideal(dtime,itap,
[3]	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,
[1545]	1308	e longitude_deg, latitude_deg, dx, dy,
[3]	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)
[2135]	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,
[2535]	1322	e longitude_deg, latitude_deg, dx, dy,
	1323	& q2,
[2135]	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
[3]	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
[1723]	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
[3]	1349	C TRACEURS
	1350
	1351	if (iflag_trac.eq.1) then
[101]	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
[3535]	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
[1442]	1372
[3]	1373	DO iq=1, nqmax
[1442]	1374
[3]	1375	CALL cltrac(dtime,ycoefh,t_seri,
[1442]	1376	s tr_seri(:,:,iq),source(:,iq),
[3]	1377	e paprs, pplay,delp,
[1442]	1378	s d_tr_vdf(:,:,iq))
	1379
[3]	1380	tr_seri(:,:,iq) = tr_seri(:,:,iq) + d_tr_vdf(:,:,iq)
	1381	d_tr_vdf(:,:,iq)= d_tr_vdf(:,:,iq)/dtime ! /s
[1442]	1382
	1383	ENDDO !nqmax
[3]	1384
[1442]	1385	endif
	1386
[3]	1387	IF (if_ebil.ge.2) THEN
	1388	ztit='after clmain'
[1543]	1389	CALL diagetpq(cell_area,ztit,ip_ebil,2,1,dtime
[3]	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)
[1543]	1392	call diagphy(cell_area,ztit,ip_ebil
[3]	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
[973]	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
[3]	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)
[2686]	1454	. + cpnew(i,j-1)/RG*d_t_ajs(i,j-1)/dtime
[3]	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
[1301]	1465
	1466	if (iflag_trac.eq.1) then
[3]	1467	tr_seri(:,:,:) = tr_seri(:,:,:) + d_tr_ajs(:,:,:)
	1468	d_tr_ajs(:,:,:)= d_tr_ajs(:,:,:)/dtime ! /s
[1301]	1469	endif
[3]	1470	endif
[973]	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)
[3]	1476	c
	1477	IF (if_ebil.ge.2) THEN
	1478	ztit='after dry_adjust'
[1543]	1479	CALL diagetpq(cell_area,ztit,ip_ebil,2,2,dtime
[3]	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)
[1543]	1482	call diagphy(cell_area,ztit,ip_ebil
[3]	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
[2193]	1489	c====================================================================
	1490	c RAYONNEMENT
	1491	c====================================================================
	1492	if (mod(itap,radpas) == 0) then
[1310]	1493
[2193]	1494	dtimerad = dtime*REAL(radpas) ! pas de temps du rayonnement (s)
[1310]	1495
[2464]	1496	! update mmean
	1497	if (ok_chem) then
	1498	mmean(:,:) = 0.
	1499	do iq = 1,nqmax - nmicro
[2622]	1500	mmean(:,:) = mmean(:,:)+tr_seri(:,:,iq)/m_tr(iq)
[2464]	1501	enddo
[2622]	1502	mmean(:,:) = 1./mmean(:,:)
	1503	rnew(:,:) = 8.314/mmean(:,:)*1.e3 ! J/kg K
[2464]	1504	endif
	1505
	1506	cc---------------------------------------------
[2193]	1507	if (callnlte .or. callthermos) then
	1508	if (ok_chem) then
[1442]	1509
[2193]	1510	! nlte : use computed chemical species
[1442]	1511
[2193]	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)
[1442]	1516
[2193]	1517	else
[1442]	1518
[2193]	1519	! nlte : use hedin climatology
[1442]	1520
[2193]	1521	call compo_hedin83_mod(pplay,rmu0,
	1522	$ co2vmr_gcm,covmr_gcm,ovmr_gcm,n2vmr_gcm,nvmr_gcm)
	1523	end if
	1524	end if
[1442]	1525
	1526	c NLTE cooling from CO2 emission
	1527	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1528
[1310]	1529	IF(callnlte) THEN
	1530	if(nltemodel.eq.0.or.nltemodel.eq.1) then
[2464]	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)
[1310]	1536	else if(nltemodel.eq.2) then
[2464]	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	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
[1442]	1542	CALL nlte_tcool(klon,klev,pplay*9.869e-6,
[1310]	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
[1543]	1557	ENDIF
[1310]	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
[2464]	1564	cc---------------------------------------------
[1442]	1565	c LTE radiative transfert / solar / IR matrix
	1566	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[2135]	1567	if (physideal) then
[2486]	1568	CALL radlwsw_newtoncool(presnivs,t_seri)
[2135]	1569	else
	1570	CALL radlwsw
	1571	e (dist, rmu0, fract, zzlev,
	1572	e paprs, pplay,ftsol, t_seri)
	1573	endif
[808]	1574
[2580]	1575	c ALBEDO VARIATIONS: test for Yeon Joo Lee
[2135]	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.)
[1301]	1580
[2580]	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
[2464]	1598	cc---------------------------------------------
[1442]	1599	c CO2 near infrared absorption
	1600	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[1310]	1601
	1602	d_t_nirco2(:,:)=0.
	1603	if (callnirco2) then
[2464]	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	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
[1442]	1616	call nirco2abs (klon, klev, pplay, dist, nqmax, tr_seri,
	1617	. rmu0, fract, d_t_nirco2)
[2464]	1618	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
[1310]	1619	endif
	1620
	1621
[2464]	1622	cc---------------------------------------------
[1310]	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
[1543]	1630	dtsw(:,:)=heat(:,:)
	1631	dtlw(:,:)=-1*cool(:,:)
	1632	ENDIF
[1310]	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
[1442]	1649	call euvheat(klon, klev, nqmax, t_seri,paprs,pplay,zzlay,
[2464]	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	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
[1442]	1656	$ tr_seri, d_tr, d_t_euv )
[2464]	1657	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
[1543]	1658
[1310]	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====================================================================
[3]	1669	c
	1670	c Ajouter la tendance des rayonnements (tous les pas)
	1671	c
	1672	DO k = 1, klev
	1673	DO i = 1, klon
[1310]	1674	t_seri(i,k) = t_seri(i,k) + d_t_rad(i,k) * dtime
[3]	1675	ENDDO
	1676	ENDDO
[973]	1677
[2193]	1678	! increment physics counter
	1679
	1680	itap = itap + 1
[2464]	1681	c====================================================================
[2193]	1682
[1310]	1683
[2464]	1684	c==============================
[1442]	1685	! -- MOLECULAR DIFFUSION ---
[2464]	1686	c==============================
[1442]	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,
[2464]	1694	& d_t_euv,d_t_conduc,d_q_moldif)
[1442]	1695
	1696
	1697	! --- update tendencies tracers ---
	1698
	1699	DO iq = 1, nqmax
	1700	DO k=1,klev
	1701	DO ig=1,klon
[2464]	1702	tr_seri(ig,k,iq)= max(tr_seri(ig,k,iq)+
	1703	& d_q_moldif(ig,k,iq)*dtime,1.e-30)
[1442]	1704	ENDDO
	1705	ENDDO
	1706	ENDDO
	1707
	1708	ENDIF ! callthermos & ok_chem
	1709
[1310]	1710	c====================================================================
[3035]	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====================================================================
[973]	1744	! tests: output tendencies
	1745	! call writefield_phy('physiq_dtrad',dtrad,klev)
[3]	1746
	1747	IF (if_ebil.ge.2) THEN
	1748	ztit='after rad'
[1543]	1749	CALL diagetpq(cell_area,ztit,ip_ebil,2,2,dtime
[3]	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)
[1543]	1752	call diagphy(cell_area,ztit,ip_ebil
[3]	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
[2047]	1764	c if (ok_orodr.or.ok_gw_nonoro) then
[1665]	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
[3]	1771	do i=1,klon
	1772	do k=2,klev
[1543]	1773	ztlev(i,k) = (t_seri(i,k)+t_seri(i,k-1))/2.
	1774	enddo
[3]	1775	enddo
	1776	do i=1,klon
	1777	do k=2,klev
[1665]	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)
[2686]	1780	zdzlev(i,k) = (zzlay(i,k)-zzlay(i,k-1))
[1665]	1781	zn2(i,k) = RG/ztlev(i,k) * ( zdtlev(i,k)/zdzlev(i,k)
[2686]	1782	. + RG/cpnew(i,k) )
[3]	1783	zn2(i,k) = max(zn2(i,k),1.e-12) ! securite
[1543]	1784	enddo
[808]	1785	zn2(i,1) = 1.e-12 ! securite
[3]	1786	enddo
	1787
[2047]	1788	c endif
[3]	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
[2047]	1806	c A ADAPTER POUR VENUS!!! [ TN: c'est fait ! ]
[1306]	1807	CALL drag_noro(klon,klev,dtime,paprs,pplay,pphi,zn2,
[3]	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,
[2047]	1812	s d_t_oro, d_u_oro, d_v_oro,
	1813	s zublstrdr,zvblstrdr,znlow,zeff,zbl,
	1814	s ztau,tau0,knu2,kbreak)
[3]	1815
[808]	1816	c print*,"d_u_oro=",d_u_oro(klon/2,:)
[3]	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.
[1543]	1829	zustrdr = 0.
	1830	zvstrdr = 0.
[2047]	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.
[3]	1840	c
	1841	ENDIF ! fin de test sur ok_orodr
	1842	c
[973]	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
[3]	1848	c ----------------------------OROLIFT
	1849	IF (ok_orolf) THEN
[808]	1850	print*,"ok_orolf NOT IMPLEMENTED !"
	1851	stop
[3]	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,
[1545]	1867	c e latitude_deg,zmea,zstd,zpic,zgam,zthe,zpic,zval,
[3]	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
[2580]	1894	! Obsolete
	1895	! but used for VCD 1.1
[2622]	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,
[3]	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
[973]	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
[3]	1931	c====================================================================
	1932	c Transport de ballons
	1933	c====================================================================
	1934	if (ballons.eq.1) then
[1545]	1935	CALL ballon(30,pdtphys,rjourvrai,gmtime*RDAY,
	1936	& latitude_deg,longitude_deg,
[1543]	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
[3]	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
[953]	1968	CALL aaam_bud (27,klon,klev,rjourvrai,gmtime*RDAY,
[3]	1969	C ra,rg,romega,
[1545]	1970	C latitude_deg,longitude_deg,pphis,
[3]	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
[2686]	1992	d_t_ec(i,k)=0.5/cpnew(i,k)
[3]	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)
[2686]	2002	. +cpnew(i,j-1)/RGd_t_ec(i,j-1)(paprs(i,j-1)-paprs(i,j))
[3]	2003	enddo
	2004	enddo
	2005
	2006	c-jld ec_conser
	2007	c====================================================================
	2008	IF (if_ebil.ge.1) THEN
	2009	ztit='after physic'
[1543]	2010	CALL diagetpq(cell_area,ztit,ip_ebil,1,1,dtime
[3]	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.
[1543]	2017	call diagphy(cell_area,ztit,ip_ebil
[3]	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
[2464]	2049	c mise à jour rho,mmean pour sorties
	2050	if(callthermos) then
	2051	call concentrations2(pplay,t_seri,tr_seri, nqmax)
	2052	endif
	2053
[2523]	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
[3]	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
[1545]	2068	c mang(i,k) = RA*cos(latitude(i))
	2069	c . (u_seri(i,k)+RAcos(latitude(i))*ROMEGA)
[1543]	2070	c . cell_area(i)(paprs(i,k)-paprs(i,k+1))/RG
[3]	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=============================================================
[1723]	2090	#ifndef MESOSCALE
[3]	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
[1642]	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
[1665]	2118	! 2D fields
	2119
	2120	CALL send_xios_field("phis",pphis)
[1642]	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))
[1665]	2127	CALL send_xios_field("cdragh",cdragh)
	2128	CALL send_xios_field("cdragm",cdragm)
[1642]	2129
[1665]	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
[1642]	2137	CALL send_xios_field("temp",t_seri)
[1665]	2138	CALL send_xios_field("pres",pplay)
	2139	CALL send_xios_field("geop",zphi)
[1642]	2140	CALL send_xios_field("vitu",u_seri)
	2141	c VENUS: regardee a l envers!!!!!!!!!!!!!!!
	2142	CALL send_xios_field("vitv",-1.*v_seri)
[1665]	2143	CALL send_xios_field("vitw",omega)
[2523]	2144	CALL send_xios_field("vitwz",vertwind)
[1665]	2145	CALL send_xios_field("Kz",ycoefh)
	2146	CALL send_xios_field("mmean",mmean)
	2147	CALL send_xios_field("rho",rho)
[2047]	2148	CALL send_xios_field("BV2",zn2)
[1665]	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)
[2836]	2157	CALL send_xios_field("dvgwno",-1.*d_v_hin)
[1665]	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
[1687]	2175	CALL send_xios_field("SWnet",swnet(:,1:klev))
	2176	CALL send_xios_field("LWnet",lwnet(:,1:klev))
[1665]	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
[2193]	2182	! when using tracers
[1687]	2183
[2193]	2184	if (iflag_trac == 1) then
[1665]	2185
[3451]	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
[2580]	2191	! production and destruction rate, cm-3.s-1
	2192	! Beware of the context*.xml file !!
[2836]	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
[2580]	2211
[2193]	2212	! tracers in gas phase, volume mixing ratio
	2213
	2214	do iq = 1,nqmax - nmicro
	2215	call send_xios_field(tname(iq),
[2464]	2216	$ tr_seri(:,:,iq)*mmean(:,:)/m_tr(iq))
[2193]	2217	end do
	2218
[2523]	2219	! tracers in gas phase, column densities
[2795]	2220
[2523]	2221	do iq = 1,nqmax - nmicro
	2222	col_dens_tr(:,iq)=0.
[2836]	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
[2523]	2230	end do
	2231
[2193]	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),
[2464]	2236	$ tr_seri(:,:,i_h2oliq)*mmean(:,:)/m_tr(i_h2oliq))
[2193]	2237	call send_xios_field(tname(i_h2so4liq),
[2464]	2238	$ tr_seri(:,:,i_h2so4liq)*mmean(:,:)/m_tr(i_h2so4liq))
[2193]	2239	if (ok_sedim) then
	2240	call send_xios_field("Fsedim",fsedim(:,1:klev))
	2241	end if
	2242	end if
	2243
[2795]	2244	! aeronomical emissions
	2245
[2836]	2246	call send_xios_field("no_emis",no_emission)
	2247	call send_xios_field("o2_emis",o2_emission)
[2795]	2248
[2193]	2249	! chemical iterations
	2250
[2464]	2251	if (tr_scheme.eq.3) call send_xios_field("iter",real(iter))
[2193]	2252
	2253	end if
	2254
[1665]	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
[1682]	2260
[2464]	2261	!! DEBUG
	2262	! if (is_master) print*,"itauphy=",itap
	2263	! if (itap.eq.10) lafin=.true.
	2264
[1682]	2265	if (lafin.and.is_omp_master) then
	2266	write(,) "physiq: call xios_context_finalize"
	2267	call xios_context_finalize
	2268	endif
	2269
[1642]	2270	#endif
[1723]	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)
[1642]	2284
[1723]	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
[3]	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
[892]	2304	CALL phyredem ("restartphy.nc")
[3]	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
[1549]	2331	END SUBROUTINE physiq
[3]	2332
[1549]	2333	END MODULE physiq_mod
[3]	2334

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

Download in other formats:

Original Format