Context Navigation

physiq_mod.F @ 1704

Last change on this file since 1704 was 1691, checked in by slebonnois, 8 years ago
SL: corrections after testing of cloud microphysics in 1D
File size: 59.1 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,
[1621]	14	. flxmw, plevmoy, tmoy,
[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
[1572]	67	USE mod_phys_lmdz_para, only : is_parallel,jj_nb,
	68	& is_north_pole_phy,
	69	& is_south_pole_phy
[892]	70	USE phys_state_var_mod ! Variables sauvegardees de la physique
[973]	71	USE write_field_phy
[892]	72	USE iophy
[1621]	73	USE cpdet_phy_mod, only: cpdet, t2tpot
[1310]	74	USE chemparam_mod
	75	USE conc
	76	USE compo_hedin83_mod2
[1506]	77	! use ieee_arithmetic
[1524]	78	use time_phylmdz_mod, only: annee_ref, day_ref, itau_phy
[1530]	79	use mod_grid_phy_lmdz, only: nbp_lon
[1621]	80	use infotrac_phy, only: iflag_trac, tname, ttext
[1642]	81	use vertical_layers_mod, only: pseudoalt
[1682]	82	use mod_phys_lmdz_omp_data, ONLY: is_omp_master
[1642]	83	#ifdef CPP_XIOS
	84	use xios_output_mod, only: initialize_xios_output,
	85	& update_xios_timestep,
	86	& send_xios_field
[1682]	87	use wxios, only: wxios_context_init, xios_context_finalize
[1642]	88	#endif
[101]	89	IMPLICIT none
	90	c======================================================================
	91	c CLEFS CPP POUR LES IO
	92	c =====================
	93	c#define histhf
	94	#define histday
	95	#define histmth
	96	#define histins
	97	c======================================================================
[3]	98	#include "dimsoil.h"
	99	#include "clesphys.h"
	100	#include "iniprint.h"
	101	#include "timerad.h"
[815]	102	#include "tabcontrol.h"
[1310]	103	#include "nirdata.h"
	104	#include "hedin.h"
[3]	105	c======================================================================
	106	LOGICAL ok_journe ! sortir le fichier journalier
	107	save ok_journe
	108	c PARAMETER (ok_journe=.true.)
	109	c
	110	LOGICAL ok_mensuel ! sortir le fichier mensuel
	111	save ok_mensuel
	112	c PARAMETER (ok_mensuel=.true.)
	113	c
	114	LOGICAL ok_instan ! sortir le fichier instantane
	115	save ok_instan
	116	c PARAMETER (ok_instan=.true.)
	117	c
	118	c======================================================================
	119	c
	120	c Variables argument:
	121	c
	122	INTEGER nlon
	123	INTEGER nlev
	124	INTEGER nqmax
	125	REAL rjourvrai
	126	REAL gmtime
	127	REAL pdtphys
	128	LOGICAL debut, lafin
	129	REAL paprs(klon,klev+1)
	130	REAL pplay(klon,klev)
	131	REAL pphi(klon,klev)
	132	REAL pphis(klon)
	133	REAL presnivs(klev)
	134
	135	! ADAPTATION GCM POUR CP(T)
	136	REAL ppk(klon,klev)
	137
	138	REAL u(klon,klev)
	139	REAL v(klon,klev)
	140	REAL t(klon,klev)
	141	REAL qx(klon,klev,nqmax)
	142
	143	REAL d_u_dyn(klon,klev)
	144	REAL d_t_dyn(klon,klev)
	145
[1301]	146	REAL flxmw(klon,klev)
[1665]	147
	148	! A VIRER POUR DYNAMICO
[1621]	149	REAL,INTENT(IN) :: plevmoy(klev+1) ! planet-averaged mean pressure (Pa) at interfaces
	150	REAL,INTENT(IN) :: tmoy(klev) ! planet-averaged mean temperature (Pa) at mid-layers
[3]	151
	152	REAL d_u(klon,klev)
	153	REAL d_v(klon,klev)
	154	REAL d_t(klon,klev)
	155	REAL d_qx(klon,klev,nqmax)
	156	REAL d_ps(klon)
	157
	158	logical ok_hf
	159	real ecrit_hf
	160	integer nid_hf
	161	save ok_hf, ecrit_hf, nid_hf
	162
	163	#ifdef histhf
	164	data ok_hf,ecrit_hf/.true.,0.25/
	165	#else
	166	data ok_hf/.false./
	167	#endif
	168
	169	c Variables propres a la physique
	170	c
[101]	171	INTEGER,save :: itap ! compteur pour la physique
[3]	172	REAL delp(klon,klev) ! epaisseur d'une couche
[1301]	173	REAL omega(klon,klev) ! vitesse verticale en Pa/s
	174
[3]	175
	176	INTEGER igwd,idx(klon),itest(klon)
	177	c
	178	c Diagnostiques 2D de drag_noro, lift_noro et gw_nonoro
	179
	180	REAL zulow(klon),zvlow(klon)
	181	REAL zustrdr(klon), zvstrdr(klon)
	182	REAL zustrli(klon), zvstrli(klon)
	183	REAL zustrhi(klon), zvstrhi(klon)
	184
	185	c Pour calcul GW drag oro et nonoro: CALCUL de N2:
[1665]	186	real zdtlev(klon,klev),zdzlev(klon,klev)
	187	real ztlev(klon,klev)
[3]	188	real zn2(klon,klev) ! BV^2 at plev
	189
	190	c Pour les bilans de moment angulaire,
	191	integer bilansmc
	192	c Pour le transport de ballons
	193	integer ballons
	194	c j'ai aussi besoin
	195	c du stress de couche limite a la surface:
	196
	197	REAL zustrcl(klon),zvstrcl(klon)
	198
	199	c et du stress total c de la physique:
	200
	201	REAL zustrph(klon),zvstrph(klon)
	202
	203	c Variables locales:
	204	c
	205	REAL cdragh(klon) ! drag coefficient pour T and Q
	206	REAL cdragm(klon) ! drag coefficient pour vent
	207	c
	208	cAA Pour TRACEURS
	209	cAA
[105]	210	REAL,save,allocatable :: source(:,:)
[3]	211	REAL ycoefh(klon,klev) ! coef d'echange pour phytrac
	212	REAL yu1(klon) ! vents dans la premiere couche U
	213	REAL yv1(klon) ! vents dans la premiere couche V
	214
	215	REAL sens(klon), dsens(klon) ! chaleur sensible et sa derivee
	216	REAL ve(klon) ! integr. verticale du transport meri. de l'energie
	217	REAL vq(klon) ! integr. verticale du transport meri. de l'eau
	218	REAL ue(klon) ! integr. verticale du transport zonal de l'energie
	219	REAL uq(klon) ! integr. verticale du transport zonal de l'eau
	220	c
[1305]	221	REAL Fsedim(klon,klev+1) ! Flux de sedimentation (kg.m-2)
[3]	222
	223	c======================================================================
	224	c
	225	c Declaration des procedures appelees
	226	c
	227	EXTERNAL ajsec ! ajustement sec
	228	EXTERNAL clmain ! couche limite
	229	EXTERNAL hgardfou ! verifier les temperatures
	230	c EXTERNAL orbite ! calculer l'orbite
	231	EXTERNAL phyetat0 ! lire l'etat initial de la physique
	232	EXTERNAL phyredem ! ecrire l'etat de redemarrage de la physique
	233	EXTERNAL radlwsw ! rayonnements solaire et infrarouge
	234	EXTERNAL suphec ! initialiser certaines constantes
	235	EXTERNAL transp ! transport total de l'eau et de l'energie
	236	EXTERNAL abort_gcm
	237	EXTERNAL printflag
	238	EXTERNAL zenang
	239	EXTERNAL diagetpq
	240	EXTERNAL conf_phys
	241	EXTERNAL diagphy
	242	EXTERNAL mucorr
[1310]	243	EXTERNAL nirco2abs
	244	EXTERNAL nir_leedat
	245	EXTERNAL nltecool
	246	EXTERNAL nlte_tcool
	247	EXTERNAL nlte_setup
	248	EXTERNAL blendrad
	249	EXTERNAL nlthermeq
	250	EXTERNAL euvheat
	251	EXTERNAL param_read
	252	EXTERNAL param_read_e107
	253	EXTERNAL conduction
	254	EXTERNAL molvis
[1442]	255	EXTERNAL moldiff_red
	256
[3]	257	c
	258	c Variables locales
	259	c
	260	CXXX PB
	261	REAL fluxt(klon,klev) ! flux turbulent de chaleur
	262	REAL fluxu(klon,klev) ! flux turbulent de vitesse u
	263	REAL fluxv(klon,klev) ! flux turbulent de vitesse v
	264	c
	265	REAL flux_dyn(klon,klev) ! flux de chaleur produit par la dynamique
	266	REAL flux_ajs(klon,klev) ! flux de chaleur ajustement sec
	267	REAL flux_ec(klon,klev) ! flux de chaleur Ec
	268	c
	269	REAL tmpout(klon,klev) ! K s-1
	270
	271	INTEGER itaprad
	272	SAVE itaprad
	273	c
	274	REAL dist, rmu0(klon), fract(klon)
	275	REAL zdtime, zlongi
	276	c
[1442]	277	INTEGER i, k, iq, ig, j, ll, ilon, ilat, ilev
[3]	278	c
	279	REAL zphi(klon,klev)
[1301]	280	REAL zzlev(klon,klev+1),zzlay(klon,klev),z1,z2
[1310]	281	real tsurf(klon)
[1301]	282
[1310]	283	c va avec nlte_tcool
	284	INTEGER ierr_nlte
	285	REAL varerr
	286
[3]	287	c Variables du changement
	288	c
	289	c ajs: ajustement sec
	290	c vdf: couche limite (Vertical DiFfusion)
	291	REAL d_t_ajs(klon,klev), d_tr_ajs(klon,klev,nqmax)
	292	REAL d_u_ajs(klon,klev), d_v_ajs(klon,klev)
	293	c
	294	REAL d_ts(klon)
	295	c
	296	REAL d_u_vdf(klon,klev), d_v_vdf(klon,klev)
	297	REAL d_t_vdf(klon,klev), d_tr_vdf(klon,klev,nqmax)
	298	c
	299	CMOD LOTT: Tendances Orography Sous-maille
	300	REAL d_u_oro(klon,klev), d_v_oro(klon,klev)
	301	REAL d_t_oro(klon,klev)
	302	REAL d_u_lif(klon,klev), d_v_lif(klon,klev)
	303	REAL d_t_lif(klon,klev)
	304	C Tendances Ondes de G non oro (runs strato).
	305	REAL d_u_hin(klon,klev), d_v_hin(klon,klev)
	306	REAL d_t_hin(klon,klev)
	307
[1310]	308	c Tendencies due to radiative scheme [K/s]
	309	c d_t_rad,dtsw,dtlw,d_t_nirco2,d_t_nlte,d_t_euv
	310	c are not computed at each physical timestep
	311	c therefore, they are defined and saved in phys_state_var_mod
	312
	313	c Tendencies due to molecular viscosity and conduction
	314	real d_t_conduc(klon,klev) ! [K/s]
	315	real d_u_molvis(klon,klev) ! (m/s) /s
	316	real d_v_molvis(klon,klev) ! (m/s) /s
	317
[1442]	318	c Tendencies due to molecular diffusion
	319	real d_q_moldif(klon,klev,nqmax)
	320
[3]	321	c
	322	c Variables liees a l'ecriture de la bande histoire physique
	323	c
	324	INTEGER ecrit_mth
	325	SAVE ecrit_mth ! frequence d'ecriture (fichier mensuel)
	326	c
	327	INTEGER ecrit_day
	328	SAVE ecrit_day ! frequence d'ecriture (fichier journalier)
	329	c
	330	INTEGER ecrit_ins
	331	SAVE ecrit_ins ! frequence d'ecriture (fichier instantane)
	332	c
	333	integer itau_w ! pas de temps ecriture = itap + itau_phy
	334
	335	c Variables locales pour effectuer les appels en serie
	336	c
	337	REAL t_seri(klon,klev)
	338	REAL u_seri(klon,klev), v_seri(klon,klev)
	339	c
[1305]	340	REAL :: tr_seri(klon,klev,nqmax)
	341	REAL :: d_tr(klon,klev,nqmax)
	342
	343	c Variables tendance sedimentation
	344
[1661]	345	REAL :: m0_mode1(klev,2),m0_mode2(klev,2)
	346	REAL :: m3_mode1(klev,3),m3_mode2 (klev,3)
	347	REAL :: d_drop_sed(klev),d_ccn_sed(klev,2),d_liq_sed(klev,2)
	348	REAL :: aer_flux(klev)
	349	REAL :: d_tr_sed(klon,klev,nqmax)
[1305]	350	REAL :: d_tr_ssed(klon)
[3]	351	c
	352	c pour ioipsl
	353	INTEGER nid_day, nid_mth, nid_ins
	354	SAVE nid_day, nid_mth, nid_ins
	355	INTEGER nhori, nvert, idayref
[97]	356	REAL zsto, zout, zsto1, zsto2, zero
[3]	357	parameter (zero=0.0e0)
	358	real zjulian
	359	save zjulian
	360
	361	CHARACTER*2 str2
	362	character*20 modname
	363	character*80 abort_message
	364	logical ok_sync
	365
	366	character*30 nom_fichier
	367	character*10 varname
	368	character*40 vartitle
	369	character*20 varunits
	370	C Variables liees au bilan d'energie et d'enthalpi
	371	REAL ztsol(klon)
	372	REAL h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
	373	$ , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
	374	SAVE h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
	375	$ , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
	376	REAL d_h_vcol, d_h_dair, d_qt, d_qw, d_ql, d_qs, d_ec
	377	REAL d_h_vcol_phy
	378	REAL fs_bound, fq_bound
	379	SAVE d_h_vcol_phy
	380	REAL zero_v(klon),zero_v2(klon,klev)
	381	CHARACTER*15 ztit
	382	INTEGER ip_ebil ! PRINT level for energy conserv. diag.
	383	SAVE ip_ebil
	384	DATA ip_ebil/2/
	385	INTEGER if_ebil ! level for energy conserv. dignostics
	386	SAVE if_ebil
	387	c+jld ec_conser
	388	REAL d_t_ec(klon,klev) ! tendance du a la conversion Ec -> E thermique
	389	c-jld ec_conser
	390
	391	c TEST VENUS...
	392	REAL mang(klon,klev) ! moment cinetique
	393	REAL mangtot ! moment cinetique total
	394
[1572]	395	c cell_area for outputs in hist*
	396	REAL cell_area_out(klon)
	397
[3]	398	c Declaration des constantes et des fonctions thermodynamiques
	399	c
	400	#include "YOMCST.h"
	401
	402	c======================================================================
	403	c INITIALISATIONS
	404	c================
	405
	406	modname = 'physiq'
	407	ok_sync=.TRUE.
	408
[119]	409	bilansmc = 0
[3]	410	ballons = 0
[892]	411	! NE FONCTIONNENT PAS ENCORE EN PARALLELE !!!
	412	if (is_parallel) then
	413	bilansmc = 0
	414	ballons = 0
	415	endif
[3]	416
	417	IF (if_ebil.ge.1) THEN
	418	DO i=1,klon
	419	zero_v(i)=0.
	420	END DO
	421	DO i=1,klon
	422	DO j=1,klev
	423	zero_v2(i,j)=0.
	424	END DO
	425	END DO
	426	END IF
	427
	428	c PREMIER APPEL SEULEMENT
	429	c========================
	430	IF (debut) THEN
[105]	431	allocate(source(klon,nqmax))
[3]	432
[1682]	433	#ifdef CPP_XIOS
	434	! Initialize XIOS context
	435	write(,) "physiq: call wxios_context_init"
	436	CALL wxios_context_init
	437	#endif
	438
[3]	439	CALL suphec ! initialiser constantes et parametres phys.
	440
	441	IF (if_ebil.ge.1) d_h_vcol_phy=0.
	442	c
	443	c appel a la lecture du physiq.def
	444	c
	445	call conf_phys(ok_journe, ok_mensuel,
	446	. ok_instan,
	447	. if_ebil)
	448
[892]	449	call phys_state_var_init
[3]	450	c
[1310]	451	c Initialising Hedin model for upper atm
	452	c (to be revised when coupled to chemistry) :
	453	call conc_init
	454	c
[3]	455	c Initialiser les compteurs:
	456	c
	457	itap = 0
	458	itaprad = 0
	459	c
	460	c Lecture startphy.nc :
	461	c
[892]	462	CALL phyetat0 ("startphy.nc")
[3]	463
[815]	464	c dtime est defini dans tabcontrol.h et lu dans startphy
[150]	465	c pdtphys est calcule a partir des nouvelles conditions:
	466	c Reinitialisation du pas de temps physique quand changement
	467	IF (ABS(dtime-pdtphys).GT.0.001) THEN
	468	WRITE(lunout,*) 'Pas physique a change',dtime,
	469	. pdtphys
	470	c abort_message='Pas physique n est pas correct '
	471	c call abort_gcm(modname,abort_message,1)
[152]	472	c----------------
	473	c pour initialiser convenablement le time_counter, il faut tenir compte
	474	c du changement de dtime en changeant itau_phy (point de depart)
	475	itau_phy = NINT(itau_phy*dtime/pdtphys)
	476	c----------------
[150]	477	dtime=pdtphys
	478	ENDIF
	479
[3]	480	radpas = NINT( RDAY/pdtphys/nbapp_rad)
	481
[815]	482	CALL printflag( ok_journe,ok_instan )
[1642]	483
	484	#ifdef CPP_XIOS
	485
	486	write(,) "physiq: call initialize_xios_output"
	487	call initialize_xios_output(rjourvrai,gmtime,pdtphys,RDAY,
	488	& presnivs,pseudoalt)
	489	#endif
	490
[3]	491	c
	492	c---------
	493	c FLOTT
	494	IF (ok_orodr) THEN
	495	DO i=1,klon
	496	rugoro(i) = MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
	497	ENDDO
	498	CALL SUGWD(klon,klev,paprs,pplay)
	499	DO i=1,klon
	500	zuthe(i)=0.
	501	zvthe(i)=0.
	502	if(zstd(i).gt.10.)then
	503	zuthe(i)=(1.-zgam(i))*cos(zthe(i))
	504	zvthe(i)=(1.-zgam(i))*sin(zthe(i))
	505	endif
	506	ENDDO
	507	ENDIF
	508
	509	if (bilansmc.eq.1) then
	510	C OUVERTURE D'UN FICHIER FORMATTE POUR STOCKER LES COMPOSANTES
	511	C DU BILAN DE MOMENT ANGULAIRE.
	512	open(27,file='aaam_bud.out',form='formatted')
	513	open(28,file='fields_2d.out',form='formatted')
	514	write(,)'Ouverture de aaam_bud.out (FL Vous parle)'
	515	write(,)'Ouverture de fields_2d.out (FL Vous parle)'
	516	endif !bilansmc
	517
	518	c--------------SLEBONNOIS
	519	C OUVERTURE DES FICHIERS FORMATTES CONTENANT LES POSITIONS ET VITESSES
	520	C DES BALLONS
	521	if (ballons.eq.1) then
	522	open(30,file='ballons-lat.out',form='formatted')
	523	open(31,file='ballons-lon.out',form='formatted')
	524	open(32,file='ballons-u.out',form='formatted')
	525	open(33,file='ballons-v.out',form='formatted')
	526	open(34,file='ballons-alt.out',form='formatted')
	527	write(,)'Ouverture des ballons*.out'
	528	endif !ballons
	529	c-------------
	530
	531	c---------
	532	C TRACEURS
	533	C source dans couche limite
[1442]	534	source(:,:) = 0.0 ! pas de source, pour l'instant
[3]	535	c---------
[1310]	536
	537	c---------
	538	c INITIALIZE THERMOSPHERIC PARAMETERS
	539	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	540
	541	if (callthermos) then
	542	if(solvarmod.eq.0) call param_read
	543	if(solvarmod.eq.1) call param_read_e107
	544	endif
	545
	546	c Initialisation (recomputed in concentration2)
	547	do ig=1,klon
	548	do j=1,klev
	549	rnew(ig,j)=R
[1665]	550	cpnew(ig,j)=cpdet(t(ig,j))
[1442]	551	mmean(ig,j)=RMD
[1310]	552	akknew(ig,j)=1.e-4
[1665]	553	rho(ig,j)=pplay(ig,j)mmean(ig,j)1e-3/(rnew(ig,j)*t(ig,j))
[1310]	554	enddo
	555	c stop
	556
	557	enddo
	558
	559	IF(callthermos.or.callnlte.or.callnirco2) THEN
	560	call compo_hedin83_init2
	561	ENDIF
[1442]	562	if (callnlte.and.nltemodel.eq.2) call nlte_setup
	563	if (callnirco2.and.nircorr.eq.1) call nir_leedat
[1310]	564	c---------
	565
[3]	566	c
	567	c Verifications:
	568	c
	569	IF (nlon .NE. klon) THEN
	570	WRITE(lunout,*)'nlon et klon ne sont pas coherents', nlon,
	571	. klon
	572	abort_message='nlon et klon ne sont pas coherents'
	573	call abort_gcm(modname,abort_message,1)
	574	ENDIF
	575	IF (nlev .NE. klev) THEN
	576	WRITE(lunout,*)'nlev et klev ne sont pas coherents', nlev,
	577	. klev
	578	abort_message='nlev et klev ne sont pas coherents'
	579	call abort_gcm(modname,abort_message,1)
	580	ENDIF
	581	c
[808]	582	IF (dtimeREAL(radpas).GT.(RDAY0.25).AND.cycle_diurne)
[3]	583	$ THEN
	584	WRITE(lunout,*)'Nbre d appels au rayonnement insuffisant'
	585	WRITE(lunout,*)"Au minimum 4 appels par jour si cycle diurne"
	586	abort_message='Nbre d appels au rayonnement insuffisant'
	587	call abort_gcm(modname,abort_message,1)
	588	ENDIF
	589	c
	590	WRITE(lunout,*)"Clef pour la convection seche, iflag_ajs=",
	591	. iflag_ajs
	592	c
[1120]	593	ecrit_mth = NINT(RDAY/dtime*ecriphy) ! tous les ecritphy jours
[3]	594	IF (ok_mensuel) THEN
	595	WRITE(lunout,*)'La frequence de sortie mensuelle est de ',
	596	. ecrit_mth
	597	ENDIF
[97]	598
[3]	599	ecrit_day = NINT(RDAY/dtime *1.0) ! tous les jours
	600	IF (ok_journe) THEN
	601	WRITE(lunout,*)'La frequence de sortie journaliere est de ',
	602	. ecrit_day
	603	ENDIF
[97]	604
[1120]	605	ecrit_ins = NINT(RDAY/dtime*ecriphy) ! Fraction de jour reglable
[3]	606	IF (ok_instan) THEN
	607	WRITE(lunout,*)'La frequence de sortie instant. est de ',
	608	. ecrit_ins
	609	ENDIF
	610
	611	c Initialisation des sorties
	612	c===========================
	613
	614	#ifdef CPP_IOIPSL
	615
	616	#ifdef histhf
	617	#include "ini_histhf.h"
	618	#endif
	619
	620	#ifdef histday
	621	#include "ini_histday.h"
	622	#endif
	623
	624	#ifdef histmth
	625	#include "ini_histmth.h"
	626	#endif
	627
	628	#ifdef histins
	629	#include "ini_histins.h"
	630	#endif
	631
	632	#endif
	633
	634	c
	635	c Initialiser les valeurs de u pour calculs tendances
	636	c (pour T, c'est fait dans phyetat0)
	637	c
	638	DO k = 1, klev
	639	DO i = 1, klon
	640	u_ancien(i,k) = u(i,k)
	641	ENDDO
	642	ENDDO
	643
[1305]	644	c---------
	645	c Ecriture fichier initialisation
	646	c PRINT*,'Ecriture Initial_State.csv'
	647	c OPEN(88,file='Trac_Point.csv',
	648	c & form='formatted')
	649	c---------
	650
	651	c---------
	652	c Initialisation des parametres des nuages
	653	c===============================================
	654
[1661]	655	c MICROPHY SANS CHIMIE: seulement si full microphy (cl_scheme=2)
	656
	657	if (ok_chem.and..not.ok_cloud) then
	658	print*,"LA CHIMIE A BESOIN DE LA MICROPHYSIQUE"
	659	print*,"ok_cloud doit etre = a ok_chem"
	660	stop
	661	endif
	662	if (.not.ok_chem.and.ok_cloud.and.(cl_scheme.eq.1)) then
	663	print*,"cl_scheme=1 doesnot work without chemistry"
	664	stop
	665	endif
	666	if (.not.ok_chem.and.ok_cloud.and.(cl_scheme.eq.2)) then
	667	print*,"Full microphysics without chemistry"
	668	c indexation of microphys tracers
	669	CALL chemparam_ini()
	670	endif
	671
	672	c number of microphysical tracers
	673	nmicro=0
	674	if (ok_cloud .AND. (cl_scheme.eq.1)) nmicro=2
[1665]	675	if (ok_cloud .AND. (cl_scheme.eq.2)) nmicro=12
[1661]	676
	677	c CAS 1D POUR MICROPHYS Aurelien
	678	if ((nlon .EQ. 1) .AND. ok_cloud .AND. (cl_scheme.eq.1)) then
[1305]	679	PRINT*,'Open profile_cloud_parameters.csv'
	680	OPEN(66,file='profile_cloud_parameters.csv',
	681	& form='formatted')
	682	endif
	683
[1661]	684	if ((nlon .EQ. 1) .AND. ok_sedim .AND. (cl_scheme.eq.1)) then
[1305]	685	PRINT*,'Open profile_cloud_sedim.csv'
	686	OPEN(77,file='profile_cloud_sedim.csv',
	687	& form='formatted')
	688	endif
	689
[1661]	690	c INIT PHOTOCHEMISTRY ! includes the indexation of microphys tracers
	691	if ((nlon .GT. 1) .AND. ok_chem) then
[1305]	692	c !!! DONC 3D !!!
	693	CALL chemparam_ini()
	694	endif
[3]	695
[1661]	696	c INIT MICROPHYS SCHEME 1 (AURELIEN)
	697	if ((nlon .GT. 1) .AND. ok_cloud .AND. (cl_scheme.eq.1)) then
[1305]	698	c !!! DONC 3D !!!
	699	CALL cloud_ini(nlon,nlev)
	700	endif
[1442]	701
[3]	702	ENDIF ! debut
	703	c======================================================================
[1442]	704	c======================================================================
[1642]	705	! ------------------------------------------------------
	706	! Initializations done at every physical timestep:
	707	! ------------------------------------------------------
[3]	708
	709	c Mettre a zero des variables de sortie (pour securite)
	710	c
	711	DO i = 1, klon
	712	d_ps(i) = 0.0
	713	ENDDO
	714	DO k = 1, klev
	715	DO i = 1, klon
	716	d_t(i,k) = 0.0
	717	d_u(i,k) = 0.0
	718	d_v(i,k) = 0.0
	719	ENDDO
	720	ENDDO
	721	DO iq = 1, nqmax
	722	DO k = 1, klev
	723	DO i = 1, klon
	724	d_qx(i,k,iq) = 0.0
	725	ENDDO
	726	ENDDO
	727	ENDDO
	728	c
	729	c Ne pas affecter les valeurs entrees de u, v, h, et q
	730	c
	731	DO k = 1, klev
	732	DO i = 1, klon
	733	t_seri(i,k) = t(i,k)
	734	u_seri(i,k) = u(i,k)
	735	v_seri(i,k) = v(i,k)
	736	ENDDO
	737	ENDDO
	738	DO iq = 1, nqmax
	739	DO k = 1, klev
	740	DO i = 1, klon
	741	tr_seri(i,k,iq) = qx(i,k,iq)
	742	ENDDO
	743	ENDDO
	744	ENDDO
	745	C
	746	DO i = 1, klon
	747	ztsol(i) = ftsol(i)
	748	ENDDO
	749	C
	750	IF (if_ebil.ge.1) THEN
	751	ztit='after dynamic'
[1543]	752	CALL diagetpq(cell_area,ztit,ip_ebil,1,1,dtime
[3]	753	e , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
	754	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
	755	C Comme les tendances de la physique sont ajoute dans la dynamique,
	756	C on devrait avoir que la variation d'entalpie par la dynamique
	757	C est egale a la variation de la physique au pas de temps precedent.
	758	C Donc la somme de ces 2 variations devrait etre nulle.
[1543]	759	call diagphy(cell_area,ztit,ip_ebil
[3]	760	e , zero_v, zero_v, zero_v, zero_v, zero_v
	761	e , zero_v, zero_v, zero_v, ztsol
	762	e , d_h_vcol+d_h_vcol_phy, d_qt, 0.
	763	s , fs_bound, fq_bound )
	764	END IF
	765
	766	c====================================================================
[1642]	767	c XIOS outputs
	768
	769	#ifdef CPP_XIOS
	770	! update XIOS time/calendar
	771	call update_xios_timestep
	772	#endif
	773
	774	c====================================================================
[3]	775	c Diagnostiquer la tendance dynamique
	776	c
	777	IF (ancien_ok) THEN
	778	DO k = 1, klev
	779	DO i = 1, klon
	780	d_u_dyn(i,k) = (u_seri(i,k)-u_ancien(i,k))/dtime
	781	d_t_dyn(i,k) = (t_seri(i,k)-t_ancien(i,k))/dtime
	782	ENDDO
	783	ENDDO
	784
	785	! ADAPTATION GCM POUR CP(T)
	786	do i=1,klon
	787	flux_dyn(i,1) = 0.0
	788	do j=2,klev
	789	flux_dyn(i,j) = flux_dyn(i,j-1)
	790	. +cpdet(t_seri(i,j-1))/RGd_t_dyn(i,j-1)(paprs(i,j-1)-paprs(i,j))
	791	enddo
	792	enddo
	793
	794	ELSE
	795	DO k = 1, klev
	796	DO i = 1, klon
	797	d_u_dyn(i,k) = 0.0
	798	d_t_dyn(i,k) = 0.0
	799	ENDDO
	800	ENDDO
	801	ancien_ok = .TRUE.
	802	ENDIF
	803	c====================================================================
[1301]	804
	805	c Calcule de vitesse verticale a partir de flux de masse verticale
	806	DO k = 1, klev
	807	DO i = 1, klon
[1543]	808	omega(i,k) = RG*flxmw(i,k) / cell_area(i)
[1301]	809	END DO
	810	END DO
	811
[3]	812	c
	813	c Ajouter le geopotentiel du sol:
	814	c
	815	DO k = 1, klev
	816	DO i = 1, klon
	817	zphi(i,k) = pphi(i,k) + pphis(i)
	818	ENDDO
	819	ENDDO
[1301]	820
	821	c calcul du geopotentiel aux niveaux intercouches
	822	c ponderation des altitudes au niveau des couches en dp/p
	823
	824	DO k=1,klev
	825	DO i=1,klon
[1310]	826	zzlay(i,k)=zphi(i,k)/RG ! [m]
[1301]	827	ENDDO
	828	ENDDO
	829	DO i=1,klon
[1310]	830	zzlev(i,1)=pphis(i)/RG ! [m]
[1301]	831	ENDDO
	832	DO k=2,klev
	833	DO i=1,klon
	834	z1=(pplay(i,k-1)+paprs(i,k))/(pplay(i,k-1)-paprs(i,k))
	835	z2=(paprs(i,k) +pplay(i,k))/(paprs(i,k) -pplay(i,k))
	836	zzlev(i,k)=(z1zzlay(i,k-1)+z2zzlay(i,k))/(z1+z2)
	837	ENDDO
	838	ENDDO
	839	DO i=1,klon
	840	zzlev(i,klev+1)=zzlay(i,klev)+(zzlay(i,klev)-zzlev(i,klev))
	841	ENDDO
	842
[3]	843	c====================================================================
	844	c
	845	c Verifier les temperatures
	846	c
	847	CALL hgardfou(t_seri,ftsol,'debutphy')
	848	c====================================================================
	849	c
	850	c Incrementer le compteur de la physique
	851	c
	852	itap = itap + 1
	853
	854	c====================================================================
	855	c Orbite et eclairement
	856	c====================================================================
	857
	858	c Pour VENUS, on fixe l'obliquite a 0 et l'eccentricite a 0.
	859	c donc pas de variations de Ls, ni de dist.
	860	c La seule chose qui compte, c'est la rotation de la planete devant
	861	c le Soleil...
	862
	863	zlongi = 0.0
	864	dist = 0.72 ! en UA
	865
	866	c Si on veut remettre l'obliquite a 3 degres et/ou l'eccentricite
	867	c a sa valeur, et prendre en compte leur evolution,
	868	c il faudra refaire un orbite.F...
	869	c CALL orbite(zlongi,dist)
	870
	871	IF (cycle_diurne) THEN
[808]	872	zdtime=dtime*REAL(radpas) ! pas de temps du rayonnement (s)
[1545]	873	CALL zenang(zlongi,gmtime,zdtime,latitude_deg,longitude_deg,
	874	& rmu0,fract)
[3]	875	ELSE
[1545]	876	call mucorr(klon,zlongi,latitude_deg,rmu0,fract)
[3]	877	ENDIF
	878
	879	c====================================================================
	880	c Calcul des tendances traceurs
	881	c====================================================================
	882
	883	if (iflag_trac.eq.1) then
[1160]	884
	885	if (tr_scheme.eq.1) then
	886	! Case 1: pseudo-chemistry with relaxation toward fixed profile
[1305]	887
[1160]	888	call phytrac_relax (debut,lafin,nqmax,
	889	I nlon,nlev,dtime,pplay,
[3]	890	O tr_seri)
[1160]	891
	892	elseif (tr_scheme.eq.2) then
	893	! Case 2: surface emission
	894	! For the moment, inspired from Mars version
	895	! However, the variable 'source' could be used in physiq
	896	! so the call to phytrac_emiss could be to initialise it.
[1305]	897
[1160]	898	call phytrac_emiss ( (rjourvrai+gmtime)*RDAY,
	899	I debut,lafin,nqmax,
	900	I nlon,nlev,dtime,paprs,
[1545]	901	I latitude_deg,longitude_deg,
[1160]	902	O tr_seri)
[3]	903
[1305]	904	elseif (tr_scheme.eq.3) then ! identical to ok_chem.or.ok_cloud
	905	! Case 3: Full chemistry and/or clouds
	906
[1442]	907	call phytrac_chimie(
[1305]	908	I debut,
	909	I gmtime,
	910	I nqmax,
[1442]	911	I klon,
[1545]	912	I latitude_deg,
	913	I longitude_deg,
[1305]	914	I nlev,
	915	I dtime,
[1442]	916	I t_seri,
	917	I pplay,
	918	O tr_seri)
[1305]	919
	920	c CALL WriteField_phy('Pression',pplay,nlev)
	921	c CALL WriteField_phy('PressionBnd',paprs,nlev+1)
	922	c CALL WriteField_phy('Temp',t_seri,nlev)
	923	c IF (ok_cloud) THEN
	924	c CALL WriteField_phy('NBRTOT',NBRTOT,nlev)
	925	c ENDIF
	926	c CALL WriteField_phy('SAl',tr_seri(:,:,i_h2so4liq),nlev)
	927	c CALL WriteField_phy('SAg',tr_seri(:,:,i_h2so4),nlev)
	928
[1661]	929	C === SEDIMENTATION ===
	930
[1305]	931	if (ok_sedim) then
	932
[1661]	933	c !! Depend on cl_scheme !!
	934
	935	if (cl_scheme.eq.1) then
	936	c ================
	937
[1305]	938	CALL new_cloud_sedim(
[1442]	939	I klon,
[1543]	940	I nlev,
	941	I dtime,
[1442]	942	I pplay,
[1543]	943	I paprs,
	944	I t_seri,
[1442]	945	I tr_seri,
[1661]	946	O d_tr_sed(:,:,1:2),
[1543]	947	O d_tr_ssed,
	948	I nqmax,
[1442]	949	O Fsedim)
[1305]	950
[1442]	951	DO k = 1, klev
	952	DO i = 1, klon
	953
[1305]	954	c--------------------
	955	c Ce test est necessaire pour eviter Xliq=NaN
[1506]	956	! IF (ieee_is_nan(d_tr_sed(i,k,1)).OR.
	957	! & ieee_is_nan(d_tr_sed(i,k,2))) THEN
	958	IF ((d_tr_sed(i,k,1).ne.d_tr_sed(i,k,1)).OR.
	959	& (d_tr_sed(i,k,2).ne.d_tr_sed(i,k,2))) THEN
[1305]	960	PRINT*,'sedim NaN PROBLEM'
	961	PRINT*,'d_tr_sed Nan?',d_tr_sed(i,k,:),'Temp',t_seri(i,k)
	962	PRINT*,'lat-lon',i,'level',k,'dtime',dtime
[1442]	963	PRINT*,'F_sed',Fsedim(i,k)
[1305]	964	PRINT*,'==============================================='
	965	d_tr_sed(i,k,:)=0.
	966	ENDIF
	967	c--------------------
	968
	969	tr_seri(i,k,i_h2so4liq) = tr_seri(i,k,i_h2so4liq)+
	970	& d_tr_sed(i,k,1)
	971	tr_seri(i,k,i_h2oliq) = tr_seri(i,k,i_h2oliq)+
	972	& d_tr_sed(i,k,2)
	973	d_tr_sed(i,k,:) = d_tr_sed(i,k,:) / dtime
	974	Fsedim(i,k) = Fsedim(i,k) / dtime
	975
	976	ENDDO
[1442]	977	ENDDO
[1305]	978
	979	Fsedim(:,klev+1) = 0.
[1661]	980
	981	elseif (cl_scheme.eq.2) then
	982	c ====================
[1305]	983
[1661]	984	d_tr_sed(:,:,:) = 0.D0
	985
	986	DO i=1, klon
	987
	988	c Mode 1
	989	m0_mode1(:,1)=tr_seri(i,:,i_m0_mode1drop)
	990	m0_mode1(:,2)=tr_seri(i,:,i_m0_mode1ccn)
	991	m3_mode1(:,1)=tr_seri(i,:,i_m3_mode1sa)
	992	m3_mode1(:,2)=tr_seri(i,:,i_m3_mode1w)
	993	m3_mode1(:,3)=tr_seri(i,:,i_m3_mode1ccn)
	994
	995	call drop_sedimentation(dtime,klev,m0_mode1,m3_mode1,
	996	& paprs(i,:),zzlev(i,:),zzlay(i,:),t_seri(i,:),
	997	& 1,
	998	& d_ccn_sed(:,1),d_drop_sed,d_ccn_sed(:,2),d_liq_sed)
	999
	1000	d_tr_sed(i,:,i_m0_mode1drop)= d_tr_sed(i,:,i_m0_mode1drop)
	1001	& + d_drop_sed
	1002	d_tr_sed(i,:,i_m0_mode1ccn) = d_tr_sed(i,:,i_m0_mode1ccn)
	1003	& + d_ccn_sed(:,1)
	1004	d_tr_sed(i,:,i_m3_mode1ccn) = d_tr_sed(i,:,i_m3_mode1ccn)
	1005	& + d_ccn_sed(:,2)
	1006	d_tr_sed(i,:,i_m3_mode1sa) = d_tr_sed(i,:,i_m3_mode1sa)
	1007	& + d_liq_sed(:,1)
	1008	d_tr_sed(i,:,i_m3_mode1w) = d_tr_sed(i,:,i_m3_mode1w)
	1009	& + d_liq_sed(:,2)
	1010
	1011	c Mode 2
	1012	m0_mode2(:,1)=tr_seri(i,:,i_m0_mode2drop)
	1013	m0_mode2(:,2)=tr_seri(i,:,i_m0_mode2ccn)
	1014	m3_mode2(:,1)=tr_seri(i,:,i_m3_mode2sa)
	1015	m3_mode2(:,2)=tr_seri(i,:,i_m3_mode2w)
	1016	m3_mode2(:,3)=tr_seri(i,:,i_m3_mode2ccn)
	1017
	1018	call drop_sedimentation(dtime,klev,m0_mode2,m3_mode2,
	1019	& paprs(i,:),zzlev(i,:),zzlay(i,:),t_seri(i,:),
	1020	& 2,
	1021	& d_ccn_sed(:,1),d_drop_sed,d_ccn_sed(:,2),d_liq_sed)
	1022
	1023	d_tr_sed(i,:,i_m0_mode2drop)= d_tr_sed(i,:,i_m0_mode2drop)
	1024	& + d_drop_sed
	1025	d_tr_sed(i,:,i_m0_mode2ccn) = d_tr_sed(i,:,i_m0_mode2ccn)
	1026	& + d_ccn_sed(:,1)
	1027	d_tr_sed(i,:,i_m3_mode2ccn) = d_tr_sed(i,:,i_m3_mode2ccn)
	1028	& + d_ccn_sed(:,2)
	1029	d_tr_sed(i,:,i_m3_mode2sa) = d_tr_sed(i,:,i_m3_mode2sa)
	1030	& + d_liq_sed(:,1)
	1031	d_tr_sed(i,:,i_m3_mode2w) = d_tr_sed(i,:,i_m3_mode2w)
	1032	& + d_liq_sed(:,2)
	1033
	1034	c Aer
	1035	call aer_sedimentation(dtime,klev,tr_seri(i,:,i_m0_aer),
	1036	& tr_seri(i,:,i_m3_aer),paprs(i,:),
	1037	& zzlev(i,:),zzlay(i,:),t_seri(i,:),
	1038	& d_tr_sed(i,:,i_m0_aer),d_tr_sed(i,:,i_m3_aer),aer_flux)
	1039
	1040	END DO
	1041
	1042	DO iq=nqmax-nmicro+1,nqmax
	1043	tr_seri(:,:,iq) = tr_seri(:,:,iq) + d_tr_sed(:,:,iq)
	1044	d_tr_sed(:,:,iq) = d_tr_sed(:,:,iq) / dtime
	1045	END DO
	1046
	1047	endif
	1048	c ====================
	1049
	1050	C === FIN SEDIMENTATION ===
	1051
[1305]	1052	endif ! ok_sedim
	1053
	1054	endif ! tr_scheme
	1055	endif ! iflag_trac
	1056
[3]	1057	c
	1058	c====================================================================
	1059	c Appeler la diffusion verticale (programme de couche limite)
	1060	c====================================================================
	1061
	1062	c-------------------------------
	1063	c VENUS TEST: on ne tient pas compte des calculs de clmain mais on force
	1064	c l'equilibre radiatif du sol
	1065	if (1.eq.0) then
	1066	if (debut) then
	1067	print*,"ATTENTION, CLMAIN SHUNTEE..."
	1068	endif
	1069
	1070	DO i = 1, klon
	1071	sens(i) = 0.0e0 ! flux de chaleur sensible au sol
	1072	fder(i) = 0.0e0
	1073	dlw(i) = 0.0e0
	1074	ENDDO
	1075
	1076	c Incrementer la temperature du sol
	1077	c
	1078	DO i = 1, klon
	1079	d_ts(i) = dtime * radsol(i)/22000. !valeur calculee par GCM pour I=200
	1080	ftsol(i) = ftsol(i) + d_ts(i)
	1081	do j=1,nsoilmx
	1082	ftsoil(i,j)=ftsol(i)
	1083	enddo
	1084	ENDDO
	1085
	1086	c-------------------------------
	1087	else
	1088	c-------------------------------
	1089
	1090	fder = dlw
	1091
	1092	! ADAPTATION GCM POUR CP(T)
[808]	1093
[3]	1094	CALL clmain(dtime,itap,
	1095	e t_seri,u_seri,v_seri,
	1096	e rmu0,
	1097	e ftsol,
	1098	$ ftsoil,
	1099	$ paprs,pplay,ppk,radsol,falbe,
	1100	e solsw, sollw, sollwdown, fder,
[1545]	1101	e longitude_deg, latitude_deg, dx, dy,
[3]	1102	e debut, lafin,
	1103	s d_t_vdf,d_u_vdf,d_v_vdf,d_ts,
	1104	s fluxt,fluxu,fluxv,cdragh,cdragm,
	1105	s dsens,
	1106	s ycoefh,yu1,yv1)
	1107
	1108	CXXX Incrementation des flux
	1109	DO i = 1, klon
	1110	sens(i) = - fluxt(i,1) ! flux de chaleur sensible au sol
	1111	fder(i) = dlw(i) + dsens(i)
	1112	ENDDO
	1113	CXXX
	1114
	1115	DO k = 1, klev
	1116	DO i = 1, klon
	1117	t_seri(i,k) = t_seri(i,k) + d_t_vdf(i,k)
	1118	d_t_vdf(i,k)= d_t_vdf(i,k)/dtime ! K/s
	1119	u_seri(i,k) = u_seri(i,k) + d_u_vdf(i,k)
	1120	d_u_vdf(i,k)= d_u_vdf(i,k)/dtime ! (m/s)/s
	1121	v_seri(i,k) = v_seri(i,k) + d_v_vdf(i,k)
	1122	d_v_vdf(i,k)= d_v_vdf(i,k)/dtime ! (m/s)/s
	1123	ENDDO
	1124	ENDDO
	1125
	1126	C TRACEURS
	1127
	1128	if (iflag_trac.eq.1) then
[101]	1129	DO k = 1, klev
	1130	DO i = 1, klon
	1131	delp(i,k) = paprs(i,k)-paprs(i,k+1)
	1132	ENDDO
	1133	ENDDO
[1442]	1134
[3]	1135	DO iq=1, nqmax
[1442]	1136
[3]	1137	CALL cltrac(dtime,ycoefh,t_seri,
[1442]	1138	s tr_seri(:,:,iq),source(:,iq),
[3]	1139	e paprs, pplay,delp,
[1442]	1140	s d_tr_vdf(:,:,iq))
	1141
[3]	1142	tr_seri(:,:,iq) = tr_seri(:,:,iq) + d_tr_vdf(:,:,iq)
	1143	d_tr_vdf(:,:,iq)= d_tr_vdf(:,:,iq)/dtime ! /s
[1442]	1144
	1145	ENDDO !nqmax
[3]	1146
[1442]	1147	endif
	1148
[3]	1149	IF (if_ebil.ge.2) THEN
	1150	ztit='after clmain'
[1543]	1151	CALL diagetpq(cell_area,ztit,ip_ebil,2,1,dtime
[3]	1152	e , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
	1153	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
[1543]	1154	call diagphy(cell_area,ztit,ip_ebil
[3]	1155	e , zero_v, zero_v, zero_v, zero_v, sens
	1156	e , zero_v, zero_v, zero_v, ztsol
	1157	e , d_h_vcol, d_qt, d_ec
	1158	s , fs_bound, fq_bound )
	1159	END IF
	1160	C
	1161	c
	1162	c Incrementer la temperature du sol
	1163	c
	1164	DO i = 1, klon
	1165	ftsol(i) = ftsol(i) + d_ts(i)
	1166	ENDDO
	1167
	1168	c Calculer la derive du flux infrarouge
	1169	c
	1170	DO i = 1, klon
	1171	dlw(i) = - 4.0RSIGMAftsol(i)**3
	1172	ENDDO
	1173
	1174	c-------------------------------
	1175	endif ! fin du VENUS TEST
	1176
[973]	1177	! tests: output tendencies
	1178	! call writefield_phy('physiq_d_t_vdf',d_t_vdf,klev)
	1179	! call writefield_phy('physiq_d_u_vdf',d_u_vdf,klev)
	1180	! call writefield_phy('physiq_d_v_vdf',d_v_vdf,klev)
	1181	! call writefield_phy('physiq_d_ts',d_ts,1)
	1182
[3]	1183	c
	1184	c Appeler l'ajustement sec
	1185	c
	1186	c===================================================================
	1187	c Convection seche
	1188	c===================================================================
	1189	c
	1190	d_t_ajs(:,:)=0.
	1191	d_u_ajs(:,:)=0.
	1192	d_v_ajs(:,:)=0.
	1193	d_tr_ajs(:,:,:)=0.
	1194	c
	1195	IF(prt_level>9)WRITE(lunout,*)
	1196	. 'AVANT LA CONVECTION SECHE , iflag_ajs='
	1197	s ,iflag_ajs
	1198
	1199	if(iflag_ajs.eq.0) then
	1200	c Rien
	1201	c ====
	1202	IF(prt_level>9)WRITE(lunout,*)'pas de convection'
	1203
	1204	else if(iflag_ajs.eq.1) then
	1205
	1206	c Ajustement sec
	1207	c ==============
	1208	IF(prt_level>9)WRITE(lunout,*)'ajsec'
	1209
	1210	! ADAPTATION GCM POUR CP(T)
	1211	CALL ajsec(paprs, pplay, ppk, t_seri, u_seri, v_seri, nqmax,
	1212	. tr_seri, d_t_ajs, d_u_ajs, d_v_ajs, d_tr_ajs)
	1213
	1214	! ADAPTATION GCM POUR CP(T)
	1215	do i=1,klon
	1216	flux_ajs(i,1) = 0.0
	1217	do j=2,klev
	1218	flux_ajs(i,j) = flux_ajs(i,j-1)
	1219	. + cpdet(t_seri(i,j-1))/RG*d_t_ajs(i,j-1)/dtime
	1220	. *(paprs(i,j-1)-paprs(i,j))
	1221	enddo
	1222	enddo
	1223
	1224	t_seri(:,:) = t_seri(:,:) + d_t_ajs(:,:)
	1225	d_t_ajs(:,:)= d_t_ajs(:,:)/dtime ! K/s
	1226	u_seri(:,:) = u_seri(:,:) + d_u_ajs(:,:)
	1227	d_u_ajs(:,:)= d_u_ajs(:,:)/dtime ! (m/s)/s
	1228	v_seri(:,:) = v_seri(:,:) + d_v_ajs(:,:)
	1229	d_v_ajs(:,:)= d_v_ajs(:,:)/dtime ! (m/s)/s
[1301]	1230
	1231	if (iflag_trac.eq.1) then
[3]	1232	tr_seri(:,:,:) = tr_seri(:,:,:) + d_tr_ajs(:,:,:)
	1233	d_tr_ajs(:,:,:)= d_tr_ajs(:,:,:)/dtime ! /s
[1301]	1234	endif
[3]	1235	endif
[973]	1236
	1237	! tests: output tendencies
	1238	! call writefield_phy('physiq_d_t_ajs',d_t_ajs,klev)
	1239	! call writefield_phy('physiq_d_u_ajs',d_u_ajs,klev)
	1240	! call writefield_phy('physiq_d_v_ajs',d_v_ajs,klev)
[3]	1241	c
	1242	IF (if_ebil.ge.2) THEN
	1243	ztit='after dry_adjust'
[1543]	1244	CALL diagetpq(cell_area,ztit,ip_ebil,2,2,dtime
[3]	1245	e , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
	1246	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
[1543]	1247	call diagphy(cell_area,ztit,ip_ebil
[3]	1248	e , zero_v, zero_v, zero_v, zero_v, sens
	1249	e , zero_v, zero_v, zero_v, ztsol
	1250	e , d_h_vcol, d_qt, d_ec
	1251	s , fs_bound, fq_bound )
	1252	END IF
	1253
[1442]	1254
[3]	1255	c====================================================================
	1256	c RAYONNEMENT
	1257	c====================================================================
	1258
[1442]	1259	c------------------------------------
[1310]	1260	c . Compute radiative tendencies :
	1261	c------------------------------------
	1262	c====================================================================
[3]	1263	IF (MOD(itaprad,radpas).EQ.0) THEN
[1310]	1264	c====================================================================
[3]	1265
[808]	1266	dtimerad = dtime*REAL(radpas) ! pas de temps du rayonnement (s)
[3]	1267	c PRINT*,'dtimerad,dtime,radpas',dtimerad,dtime,radpas
	1268
[1310]	1269
[1442]	1270	c------------------------------------
	1271	c . Compute mean mass, cp and R :
	1272	c------------------------------------
	1273
[1310]	1274	if(callthermos) then
[1452]	1275	call concentrations2(pplay,t_seri,d_t,tr_seri, nqmax)
[1442]	1276
[1310]	1277	endif
	1278
	1279
[1442]	1280	cc!!! ADD key callhedin
[1310]	1281
[1442]	1282	IF(callnlte.or.callthermos) THEN
	1283	call compo_hedin83_mod(pplay,rmu0,
[1543]	1284	& co2vmr_gcm,covmr_gcm,ovmr_gcm,n2vmr_gcm,nvmr_gcm)
[1442]	1285
	1286	IF(ok_chem) then
	1287
	1288	CC !! GG : Using only mayor species tracers abundances to compute NLTE heating/cooling
	1289
	1290	CC Conversion [mmr] ---> [vmr]
	1291
	1292	co2vmr_gcm(:,:) = tr_seri(1:nlon,1:nlev,i_co2)*
	1293	& mmean(1:nlon,1:nlev)/M_tr(i_co2)
	1294	covmr_gcm(:,:) = tr_seri(1:nlon,1:nlev,i_co)*
	1295	& mmean(1:nlon,1:nlev)/M_tr(i_co)
	1296	ovmr_gcm(:,:) = tr_seri(1:nlon,1:nlev,i_o)*
	1297	& mmean(1:nlon,1:nlev)/M_tr(i_o)
	1298	n2vmr_gcm(:,:) = tr_seri(1:nlon,1:nlev,i_n2)*
	1299	& mmean(1:nlon,1:nlev)/M_tr(i_n2)
	1300
	1301	ENDIF
	1302
[1543]	1303	ENDIF
[1442]	1304
	1305	c
	1306	c NLTE cooling from CO2 emission
	1307	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1308
[1310]	1309	IF(callnlte) THEN
	1310	if(nltemodel.eq.0.or.nltemodel.eq.1) then
[1442]	1311	CALL nltecool(klon, klev, nqmax, pplay*9.869e-6, t_seri,
	1312	$ tr_seri, d_t_nlte)
[1310]	1313	else if(nltemodel.eq.2) then
[1442]	1314	CALL nlte_tcool(klon,klev,pplay*9.869e-6,
[1310]	1315	$ t_seri,zzlay,co2vmr_gcm, n2vmr_gcm, covmr_gcm,
	1316	$ ovmr_gcm,d_t_nlte,ierr_nlte,varerr )
	1317	if(ierr_nlte.gt.0) then
	1318	write(,)
	1319	$ 'WARNING: nlte_tcool output with error message',
	1320	$ 'ierr_nlte=',ierr_nlte,'varerr=',varerr
	1321	write(,)'I will continue anyway'
	1322	endif
	1323
	1324	endif
	1325
	1326	ELSE
	1327
	1328	d_t_nlte(:,:)=0.
	1329
[1543]	1330	ENDIF
[1310]	1331
	1332	c Find number of layers for LTE radiation calculations
	1333
	1334	IF(callnlte .or. callnirco2)
	1335	$ CALL nlthermeq(klon, klev, paprs, pplay)
	1336
[1442]	1337	c
	1338	c LTE radiative transfert / solar / IR matrix
	1339	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[3]	1340	CALL radlwsw
[1301]	1341	e (dist, rmu0, fract, zzlev,
[1310]	1342	e paprs, pplay,ftsol, t_seri)
[808]	1343
[1301]	1344
[1442]	1345	c CO2 near infrared absorption
	1346	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[1310]	1347
	1348	d_t_nirco2(:,:)=0.
	1349	if (callnirco2) then
[1442]	1350	call nirco2abs (klon, klev, pplay, dist, nqmax, tr_seri,
	1351	. rmu0, fract, d_t_nirco2)
[1310]	1352	endif
	1353
	1354
	1355	c Net atmospheric radiative heating rate (K.s-1)
	1356	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1357
	1358	IF(callnlte.or.callnirco2) THEN
	1359	CALL blendrad(klon, klev, pplay,heat,
	1360	& cool, d_t_nirco2,d_t_nlte, dtsw, dtlw)
	1361	ELSE
[1543]	1362	dtsw(:,:)=heat(:,:)
	1363	dtlw(:,:)=-1*cool(:,:)
	1364	ENDIF
[1310]	1365
	1366	DO k=1,klev
	1367	DO i=1,klon
	1368	d_t_rad(i,k) = dtsw(i,k) + dtlw(i,k) ! K/s
	1369	ENDDO
	1370	ENDDO
	1371
	1372
	1373	cc---------------------------------------------
	1374
	1375	c EUV heating rate (K.s-1)
	1376	c ~~~~~~~~~~~~~~~~~~~~~~~~
	1377
	1378	d_t_euv(:,:)=0.
	1379
	1380	IF (callthermos) THEN
	1381
[1442]	1382	c call euvheat(klon, klev,t_seri,paprs,pplay,zzlay,
	1383	c $ rmu0,pdtphys,gmtime,rjourvrai, co2vmr_gcm, n2vmr_gcm,
	1384	c $ covmr_gcm, ovmr_gcm,d_t_euv )
	1385	call euvheat(klon, klev, nqmax, t_seri,paprs,pplay,zzlay,
	1386	$ rmu0,pdtphys,gmtime,rjourvrai,
	1387	$ tr_seri, d_tr, d_t_euv )
[1543]	1388
[1310]	1389	DO k=1,klev
	1390	DO ig=1,klon
	1391	d_t_rad(ig,k)=d_t_rad(ig,k)+d_t_euv(ig,k)
	1392
	1393	ENDDO
	1394	ENDDO
	1395
	1396	ENDIF ! callthermos
	1397
[3]	1398	c====================================================================
[1310]	1399	itaprad = 0
	1400	ENDIF ! radpas
	1401	c====================================================================
[3]	1402	c
	1403	c Ajouter la tendance des rayonnements (tous les pas)
	1404	c
	1405	DO k = 1, klev
	1406	DO i = 1, klon
[1310]	1407	t_seri(i,k) = t_seri(i,k) + d_t_rad(i,k) * dtime
[3]	1408	ENDDO
	1409	ENDDO
[973]	1410
[1310]	1411	! CONDUCTION and MOLECULAR VISCOSITY
	1412	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1413
	1414	d_t_conduc(:,:)=0.
	1415	d_u_molvis(:,:)=0.
	1416	d_v_molvis(:,:)=0.
	1417
	1418	IF (callthermos) THEN
	1419
	1420	tsurf(:)=t_seri(:,1)
	1421	call conduction(klon, klev,pdtphys,
	1422	$ pplay,paprs,t_seri,
	1423	$ tsurf,zzlev,zzlay,d_t_conduc)
	1424
	1425	call molvis(klon, klev,pdtphys,
	1426	$ pplay,paprs,t_seri,
	1427	$ u,tsurf,zzlev,zzlay,d_u_molvis)
	1428
	1429	call molvis(klon, klev, pdtphys,
	1430	$ pplay,paprs,t_seri,
[1442]	1431	$ v,tsurf,zzlev,zzlay,d_v_molvis)
[1310]	1432
	1433	DO k=1,klev
	1434	DO ig=1,klon
	1435	t_seri(ig,k)= t_seri(ig,k)+ d_t_conduc(ig,k)*dtime ! [K]
	1436	u_seri(ig,k)= u_seri(ig,k)+ d_u_molvis(ig,k)*dtime ! m/s
	1437	v_seri(ig,k)= v_seri(ig,k)+ d_v_molvis(ig,k)*dtime ! m/s
	1438	ENDDO
	1439	ENDDO
[1442]	1440	ENDIF
[1310]	1441
	1442
[1442]	1443	! -- MOLECULAR DIFFUSION ---
	1444
	1445	d_q_moldif(:,:,:)=0
	1446
	1447	IF (callthermos .and. ok_chem) THEN
	1448
	1449	call moldiff_red(klon, klev, nqmax,
	1450	& pplay,paprs,t_seri, tr_seri, pdtphys,
	1451	& zzlay,d_t_euv,d_t_conduc,d_q_moldif)
	1452
	1453
	1454	! --- update tendencies tracers ---
	1455
	1456	DO iq = 1, nqmax
	1457	DO k=1,klev
	1458	DO ig=1,klon
	1459	tr_seri(ig,k,iq)= tr_seri(ig,k,iq)+
	1460	& d_q_moldif(ig,k,iq)*dtime ! [Kg/kg]?
	1461	ENDDO
	1462	ENDDO
	1463	ENDDO
	1464
	1465
	1466	ENDIF ! callthermos & ok_chem
	1467
[1310]	1468	c====================================================================
[973]	1469	! tests: output tendencies
	1470	! call writefield_phy('physiq_dtrad',dtrad,klev)
[3]	1471
	1472	IF (if_ebil.ge.2) THEN
	1473	ztit='after rad'
[1543]	1474	CALL diagetpq(cell_area,ztit,ip_ebil,2,2,dtime
[3]	1475	e , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
	1476	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
[1543]	1477	call diagphy(cell_area,ztit,ip_ebil
[3]	1478	e , topsw, toplw, solsw, sollw, zero_v
	1479	e , zero_v, zero_v, zero_v, ztsol
	1480	e , d_h_vcol, d_qt, d_ec
	1481	s , fs_bound, fq_bound )
	1482	END IF
	1483	c
	1484
	1485	c====================================================================
	1486	c Calcul des gravity waves FLOTT
	1487	c====================================================================
	1488	c
	1489	if (ok_orodr.or.ok_gw_nonoro) then
[1665]	1490
	1491	c CALCUL DE N2
	1492	c UTILISE LA RELATION ENTRE N2 ET STABILITE
	1493	c N2 = RG/T (dT/dz+RG/cp(T))
	1494	c ET DONC EN N'UTILISE QUE LA TEMPERATURE, PAS teta.
	1495
[3]	1496	do i=1,klon
	1497	do k=2,klev
[1543]	1498	ztlev(i,k) = (t_seri(i,k)+t_seri(i,k-1))/2.
	1499	enddo
[3]	1500	enddo
	1501	do i=1,klon
	1502	do k=2,klev
[1665]	1503	ztlev(i,k) = (t_seri(i,k)+t_seri(i,k-1))/2.
	1504	zdtlev(i,k) = t_seri(i,k)-t_seri(i,k-1)
	1505	zdzlev(i,k) = (zphi(i,k)-zphi(i,k-1))/RG
	1506	zn2(i,k) = RG/ztlev(i,k) * ( zdtlev(i,k)/zdzlev(i,k)
	1507	. + RG/cpdet(ztlev(i,k)) )
[3]	1508	zn2(i,k) = max(zn2(i,k),1.e-12) ! securite
[1543]	1509	enddo
[808]	1510	zn2(i,1) = 1.e-12 ! securite
[3]	1511	enddo
	1512
	1513	endif
	1514
	1515	c ----------------------------ORODRAG
	1516	IF (ok_orodr) THEN
	1517	c
	1518	c selection des points pour lesquels le shema est actif:
	1519	igwd=0
	1520	DO i=1,klon
	1521	itest(i)=0
	1522	c IF ((zstd(i).gt.10.0)) THEN
	1523	IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
	1524	itest(i)=1
	1525	igwd=igwd+1
	1526	idx(igwd)=i
	1527	ENDIF
	1528	ENDDO
	1529	c igwdim=MAX(1,igwd)
	1530	c
	1531	c A ADAPTER POUR VENUS!!!
[1306]	1532	CALL drag_noro(klon,klev,dtime,paprs,pplay,pphi,zn2,
[3]	1533	e zmea,zstd, zsig, zgam, zthe,zpic,zval,
	1534	e igwd,idx,itest,
	1535	e t_seri, u_seri, v_seri,
	1536	s zulow, zvlow, zustrdr, zvstrdr,
	1537	s d_t_oro, d_u_oro, d_v_oro)
	1538
[808]	1539	c print*,"d_u_oro=",d_u_oro(klon/2,:)
[3]	1540	c ajout des tendances
	1541	t_seri(:,:) = t_seri(:,:) + d_t_oro(:,:)
	1542	d_t_oro(:,:)= d_t_oro(:,:)/dtime ! K/s
	1543	u_seri(:,:) = u_seri(:,:) + d_u_oro(:,:)
	1544	d_u_oro(:,:)= d_u_oro(:,:)/dtime ! (m/s)/s
	1545	v_seri(:,:) = v_seri(:,:) + d_v_oro(:,:)
	1546	d_v_oro(:,:)= d_v_oro(:,:)/dtime ! (m/s)/s
	1547	c
	1548	ELSE
	1549	d_t_oro = 0.
	1550	d_u_oro = 0.
	1551	d_v_oro = 0.
[1543]	1552	zustrdr = 0.
	1553	zvstrdr = 0.
[3]	1554	c
	1555	ENDIF ! fin de test sur ok_orodr
	1556	c
[973]	1557	! tests: output tendencies
	1558	! call writefield_phy('physiq_d_t_oro',d_t_oro,klev)
	1559	! call writefield_phy('physiq_d_u_oro',d_u_oro,klev)
	1560	! call writefield_phy('physiq_d_v_oro',d_v_oro,klev)
	1561
[3]	1562	c ----------------------------OROLIFT
	1563	IF (ok_orolf) THEN
[808]	1564	print*,"ok_orolf NOT IMPLEMENTED !"
	1565	stop
[3]	1566	c
	1567	c selection des points pour lesquels le shema est actif:
	1568	igwd=0
	1569	DO i=1,klon
	1570	itest(i)=0
	1571	IF ((zpic(i)-zmea(i)).GT.100.) THEN
	1572	itest(i)=1
	1573	igwd=igwd+1
	1574	idx(igwd)=i
	1575	ENDIF
	1576	ENDDO
	1577	c igwdim=MAX(1,igwd)
	1578	c
	1579	c A ADAPTER POUR VENUS!!!
	1580	c CALL lift_noro(klon,klev,dtime,paprs,pplay,
[1545]	1581	c e latitude_deg,zmea,zstd,zpic,zgam,zthe,zpic,zval,
[3]	1582	c e igwd,idx,itest,
	1583	c e t_seri, u_seri, v_seri,
	1584	c s zulow, zvlow, zustrli, zvstrli,
	1585	c s d_t_lif, d_u_lif, d_v_lif )
	1586
	1587	c
	1588	c ajout des tendances
	1589	t_seri(:,:) = t_seri(:,:) + d_t_lif(:,:)
	1590	d_t_lif(:,:)= d_t_lif(:,:)/dtime ! K/s
	1591	u_seri(:,:) = u_seri(:,:) + d_u_lif(:,:)
	1592	d_u_lif(:,:)= d_u_lif(:,:)/dtime ! (m/s)/s
	1593	v_seri(:,:) = v_seri(:,:) + d_v_lif(:,:)
	1594	d_v_lif(:,:)= d_v_lif(:,:)/dtime ! (m/s)/s
	1595	c
	1596	ELSE
	1597	d_t_lif = 0.
	1598	d_u_lif = 0.
	1599	d_v_lif = 0.
	1600	zustrli = 0.
	1601	zvstrli = 0.
	1602	c
	1603	ENDIF ! fin de test sur ok_orolf
	1604
	1605	c ---------------------------- NON-ORO GRAVITY WAVES
	1606	IF(ok_gw_nonoro) then
	1607
	1608	call flott_gwd_ran(klon,klev,dtime,pplay,zn2,
[1621]	1609	e t_seri, u_seri, v_seri, plevmoy,
[3]	1610	o zustrhi,zvstrhi,
	1611	o d_t_hin, d_u_hin, d_v_hin)
	1612
	1613	c ajout des tendances
	1614
	1615	t_seri(:,:) = t_seri(:,:) + d_t_hin(:,:)
	1616	d_t_hin(:,:)= d_t_hin(:,:)/dtime ! K/s
	1617	u_seri(:,:) = u_seri(:,:) + d_u_hin(:,:)
	1618	d_u_hin(:,:)= d_u_hin(:,:)/dtime ! (m/s)/s
	1619	v_seri(:,:) = v_seri(:,:) + d_v_hin(:,:)
	1620	d_v_hin(:,:)= d_v_hin(:,:)/dtime ! (m/s)/s
	1621
	1622	ELSE
	1623	d_t_hin = 0.
	1624	d_u_hin = 0.
	1625	d_v_hin = 0.
	1626	zustrhi = 0.
	1627	zvstrhi = 0.
	1628
	1629	ENDIF ! fin de test sur ok_gw_nonoro
	1630
[973]	1631	! tests: output tendencies
	1632	! call writefield_phy('physiq_d_t_hin',d_t_hin,klev)
	1633	! call writefield_phy('physiq_d_u_hin',d_u_hin,klev)
	1634	! call writefield_phy('physiq_d_v_hin',d_v_hin,klev)
	1635
[3]	1636	c====================================================================
	1637	c Transport de ballons
	1638	c====================================================================
	1639	if (ballons.eq.1) then
[1545]	1640	CALL ballon(30,pdtphys,rjourvrai,gmtime*RDAY,
	1641	& latitude_deg,longitude_deg,
[1543]	1642	c C t,pplay,u,v,pphi) ! alt above surface (smoothed for GCM)
	1643	C t,pplay,u,v,zphi) ! alt above planet average radius
[3]	1644	endif !ballons
	1645
	1646	c====================================================================
	1647	c Bilan de mmt angulaire
	1648	c====================================================================
	1649	if (bilansmc.eq.1) then
	1650	CMODDEB FLOTT
	1651	C CALCULER LE BILAN DE MOMENT ANGULAIRE (DIAGNOSTIQUE)
	1652	C STRESS NECESSAIRES: COUCHE LIMITE ET TOUTE LA PHYSIQUE
	1653
	1654	DO i = 1, klon
	1655	zustrph(i)=0.
	1656	zvstrph(i)=0.
	1657	zustrcl(i)=0.
	1658	zvstrcl(i)=0.
	1659	ENDDO
	1660	DO k = 1, klev
	1661	DO i = 1, klon
	1662	zustrph(i)=zustrph(i)+(u_seri(i,k)-u(i,k))/dtime*
	1663	c (paprs(i,k)-paprs(i,k+1))/rg
	1664	zvstrph(i)=zvstrph(i)+(v_seri(i,k)-v(i,k))/dtime*
	1665	c (paprs(i,k)-paprs(i,k+1))/rg
	1666	zustrcl(i)=zustrcl(i)+d_u_vdf(i,k)*
	1667	c (paprs(i,k)-paprs(i,k+1))/rg
	1668	zvstrcl(i)=zvstrcl(i)+d_v_vdf(i,k)*
	1669	c (paprs(i,k)-paprs(i,k+1))/rg
	1670	ENDDO
	1671	ENDDO
	1672
[953]	1673	CALL aaam_bud (27,klon,klev,rjourvrai,gmtime*RDAY,
[3]	1674	C ra,rg,romega,
[1545]	1675	C latitude_deg,longitude_deg,pphis,
[3]	1676	C zustrdr,zustrli,zustrcl,
	1677	C zvstrdr,zvstrli,zvstrcl,
	1678	C paprs,u,v)
	1679
	1680	CCMODFIN FLOTT
	1681	endif !bilansmc
	1682
	1683	c====================================================================
	1684	c====================================================================
	1685	c Calculer le transport de l'eau et de l'energie (diagnostique)
	1686	c
	1687	c A REVOIR POUR VENUS...
	1688	c
	1689	c CALL transp (paprs,ftsol,
	1690	c e t_seri, q_seri, u_seri, v_seri, zphi,
	1691	c s ve, vq, ue, uq)
	1692	c
	1693	c====================================================================
	1694	c+jld ec_conser
	1695	DO k = 1, klev
	1696	DO i = 1, klon
	1697	d_t_ec(i,k)=0.5/cpdet(t_seri(i,k))
	1698	$ (u(i,k)2+v(i,k)2-u_seri(i,k)2-v_seri(i,k)*2)
	1699	t_seri(i,k)=t_seri(i,k)+d_t_ec(i,k)
	1700	d_t_ec(i,k) = d_t_ec(i,k)/dtime
	1701	END DO
	1702	END DO
	1703	do i=1,klon
	1704	flux_ec(i,1) = 0.0
	1705	do j=2,klev
	1706	flux_ec(i,j) = flux_ec(i,j-1)
	1707	. +cpdet(t_seri(i,j-1))/RGd_t_ec(i,j-1)(paprs(i,j-1)-paprs(i,j))
	1708	enddo
	1709	enddo
	1710
	1711	c-jld ec_conser
	1712	c====================================================================
	1713	IF (if_ebil.ge.1) THEN
	1714	ztit='after physic'
[1543]	1715	CALL diagetpq(cell_area,ztit,ip_ebil,1,1,dtime
[3]	1716	e , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
	1717	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
	1718	C Comme les tendances de la physique sont ajoute dans la dynamique,
	1719	C on devrait avoir que la variation d'entalpie par la dynamique
	1720	C est egale a la variation de la physique au pas de temps precedent.
	1721	C Donc la somme de ces 2 variations devrait etre nulle.
[1543]	1722	call diagphy(cell_area,ztit,ip_ebil
[3]	1723	e , topsw, toplw, solsw, sollw, sens
	1724	e , zero_v, zero_v, zero_v, ztsol
	1725	e , d_h_vcol, d_qt, d_ec
	1726	s , fs_bound, fq_bound )
	1727	C
	1728	d_h_vcol_phy=d_h_vcol
	1729	C
	1730	END IF
	1731	C
	1732	c=======================================================================
	1733	c SORTIES
	1734	c=======================================================================
	1735
	1736	c Convertir les incrementations en tendances
	1737	c
	1738	DO k = 1, klev
	1739	DO i = 1, klon
	1740	d_u(i,k) = ( u_seri(i,k) - u(i,k) ) / dtime
	1741	d_v(i,k) = ( v_seri(i,k) - v(i,k) ) / dtime
	1742	d_t(i,k) = ( t_seri(i,k) - t(i,k) ) / dtime
	1743	ENDDO
	1744	ENDDO
	1745	c
	1746	DO iq = 1, nqmax
	1747	DO k = 1, klev
	1748	DO i = 1, klon
	1749	d_qx(i,k,iq) = ( tr_seri(i,k,iq) - qx(i,k,iq) ) / dtime
	1750	ENDDO
	1751	ENDDO
	1752	ENDDO
	1753
	1754	c------------------------
	1755	c Calcul moment cinetique
	1756	c------------------------
	1757	c TEST VENUS...
	1758	c mangtot = 0.0
	1759	c DO k = 1, klev
	1760	c DO i = 1, klon
[1545]	1761	c mang(i,k) = RA*cos(latitude(i))
	1762	c . (u_seri(i,k)+RAcos(latitude(i))*ROMEGA)
[1543]	1763	c . cell_area(i)(paprs(i,k)-paprs(i,k+1))/RG
[3]	1764	c mangtot=mangtot+mang(i,k)
	1765	c ENDDO
	1766	c ENDDO
	1767	c print*,"Moment cinetique total = ",mangtot
	1768	c
	1769	c------------------------
	1770	c
	1771	c Sauvegarder les valeurs de t et u a la fin de la physique:
	1772	c
	1773	DO k = 1, klev
	1774	DO i = 1, klon
	1775	u_ancien(i,k) = u_seri(i,k)
	1776	t_ancien(i,k) = t_seri(i,k)
	1777	ENDDO
	1778	ENDDO
	1779	c
	1780	c=============================================================
	1781	c Ecriture des sorties
	1782	c=============================================================
	1783
	1784	#ifdef CPP_IOIPSL
	1785
	1786	#ifdef histhf
	1787	#include "write_histhf.h"
	1788	#endif
	1789
	1790	#ifdef histday
	1791	#include "write_histday.h"
	1792	#endif
	1793
	1794	#ifdef histmth
	1795	#include "write_histmth.h"
	1796	#endif
	1797
	1798	#ifdef histins
	1799	#include "write_histins.h"
	1800	#endif
	1801
	1802	#endif
	1803
[1642]	1804	! XIOS outputs
	1805	! This can be done ANYWHERE in the physics routines !
	1806
	1807	#ifdef CPP_XIOS
	1808	! Send fields to XIOS: (NB these fields must also be defined as
	1809	! <field id="..." /> in context_lmdz_physics.xml to be correctly used)
	1810
[1665]	1811	! 2D fields
	1812
	1813	CALL send_xios_field("phis",pphis)
[1642]	1814	cell_area_out(:)=cell_area(:)
	1815	if (is_north_pole_phy) cell_area_out(1)=cell_area(1)/nbp_lon
	1816	if (is_south_pole_phy) cell_area_out(klon)=cell_area(klon)/nbp_lon
	1817	CALL send_xios_field("aire",cell_area_out)
	1818	CALL send_xios_field("tsol",ftsol)
	1819	CALL send_xios_field("psol",paprs(:,1))
[1665]	1820	CALL send_xios_field("cdragh",cdragh)
	1821	CALL send_xios_field("cdragm",cdragm)
[1642]	1822
[1665]	1823	CALL send_xios_field("tops",topsw)
	1824	CALL send_xios_field("topl",toplw)
	1825	CALL send_xios_field("sols",solsw)
	1826	CALL send_xios_field("soll",sollw)
	1827
	1828	! 3D fields
	1829
[1642]	1830	CALL send_xios_field("temp",t_seri)
[1665]	1831	CALL send_xios_field("pres",pplay)
	1832	CALL send_xios_field("geop",zphi)
[1642]	1833	CALL send_xios_field("vitu",u_seri)
	1834	c VENUS: regardee a l envers!!!!!!!!!!!!!!!
	1835	CALL send_xios_field("vitv",-1.*v_seri)
[1665]	1836	CALL send_xios_field("vitw",omega)
	1837	CALL send_xios_field("Kz",ycoefh)
	1838	CALL send_xios_field("mmean",mmean)
	1839	CALL send_xios_field("rho",rho)
	1840
	1841	CALL send_xios_field("dudyn",d_u_dyn)
	1842	CALL send_xios_field("duvdf",d_u_vdf)
	1843	c VENUS: regardee a l envers!!!!!!!!!!!!!!!
	1844	CALL send_xios_field("dvvdf",-1.*d_v_vdf)
	1845	CALL send_xios_field("duajs",d_u_ajs)
	1846	CALL send_xios_field("dugwo",d_u_oro)
	1847	CALL send_xios_field("dugwno",d_u_hin)
	1848	CALL send_xios_field("dumolvis",d_u_molvis)
	1849	c VENUS: regardee a l envers!!!!!!!!!!!!!!!
	1850	CALL send_xios_field("dvmolvis",-1.*d_v_molvis)
	1851	CALL send_xios_field("dtdyn",d_t_dyn)
	1852	CALL send_xios_field("dtphy",d_t)
	1853	CALL send_xios_field("dtvdf",d_t_vdf)
	1854	CALL send_xios_field("dtajs",d_t_ajs)
	1855	CALL send_xios_field("dtswr",dtsw)
	1856	CALL send_xios_field("dtswrNLTE",d_t_nirco2)
	1857	CALL send_xios_field("dtswrLTE",heat)
	1858	CALL send_xios_field("dtlwr",dtlw)
	1859	CALL send_xios_field("dtlwrNLTE",d_t_nlte)
	1860	CALL send_xios_field("dtlwrLTE",-1.*cool)
	1861	CALL send_xios_field("dteuv",d_t_euv)
	1862	CALL send_xios_field("dtcond",d_t_conduc)
	1863	CALL send_xios_field("dtec",d_t_ec)
	1864
[1687]	1865	CALL send_xios_field("SWnet",swnet(:,1:klev))
	1866	CALL send_xios_field("LWnet",lwnet(:,1:klev))
[1665]	1867	CALL send_xios_field("fluxvdf",fluxt)
	1868	CALL send_xios_field("fluxdyn",flux_dyn)
	1869	CALL send_xios_field("fluxajs",flux_ajs)
	1870	CALL send_xios_field("fluxec",flux_ec)
	1871
[1691]	1872	c When using tracers
[1665]	1873	IF (iflag_trac.eq.1) THEN
[1691]	1874	c photochemical compounds !!!outputs in [vmr]
	1875	DO iq=1,nqmax-nmicro
[1665]	1876	CALL send_xios_field(tname(iq),qx(:,:,iq)*mmean(:,:)/M_tr(iq))
	1877	ENDDO
[1687]	1878
	1879	IF ((tr_scheme.eq.3).and.(ok_sedim).and.(cl_scheme.eq.1)) THEN
	1880	CALL send_xios_field("Fsedim",Fsedim(:,1:klev))
	1881	ENDIF
[1665]	1882	ENDIF
	1883
	1884	IF (callthermos .and. ok_chem) THEN
	1885	CALL send_xios_field("d_qmoldifCO2",d_q_moldif(:,:,i_co2))
	1886	CALL send_xios_field("d_qmoldifO3p",d_q_moldif(:,:,i_o))
	1887	CALL send_xios_field("d_qmoldifN2",d_q_moldif(:,:,i_n2))
	1888	ENDIF
[1682]	1889
	1890	if (lafin.and.is_omp_master) then
	1891	write(,) "physiq: call xios_context_finalize"
	1892	call xios_context_finalize
	1893	endif
	1894
[1642]	1895	#endif
	1896
[3]	1897	c====================================================================
	1898	c Si c'est la fin, il faut conserver l'etat de redemarrage
	1899	c====================================================================
	1900	c
	1901	IF (lafin) THEN
	1902	itau_phy = itau_phy + itap
[892]	1903	CALL phyredem ("restartphy.nc")
[3]	1904
	1905	c--------------FLOTT
	1906	CMODEB LOTT
	1907	C FERMETURE DU FICHIER FORMATTE CONTENANT LES COMPOSANTES
	1908	C DU BILAN DE MOMENT ANGULAIRE.
	1909	if (bilansmc.eq.1) then
	1910	write(,)'Fermeture de aaam_bud.out (FL Vous parle)'
	1911	close(27)
	1912	close(28)
	1913	endif !bilansmc
	1914	CMODFIN
	1915	c-------------
	1916	c--------------SLEBONNOIS
	1917	C FERMETURE DES FICHIERS FORMATTES CONTENANT LES POSITIONS ET VITESSES
	1918	C DES BALLONS
	1919	if (ballons.eq.1) then
	1920	write(,)'Fermeture des ballons*.out'
	1921	close(30)
	1922	close(31)
	1923	close(32)
	1924	close(33)
	1925	close(34)
	1926	endif !ballons
	1927	c-------------
	1928	ENDIF
	1929
[1549]	1930	END SUBROUTINE physiq
[3]	1931
[1549]	1932	END MODULE physiq_mod
[3]	1933

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

Download in other formats:

Original Format