Context Navigation

physiq.F @ 1310

Last change on this file since 1310 was 1310, checked in by slebonnois, 11 years ago
SL: VENUS VERTICAL EXTENSION. NLTE and thermospheric processes, to be run with 78 levels and specific inputs.
File size: 50.1 KB

Rev	Line
[3]	1	!
	2	! $Header: /home/cvsroot/LMDZ4/libf/phylmd/physiq.F,v 1.8 2005/02/24 09:58:18 fairhead Exp $
	3	!
	4	c
	5	SUBROUTINE physiq (nlon,nlev,nqmax,
	6	. debut,lafin,rjourvrai,gmtime,pdtphys,
[97]	7	. paprs,pplay,ppk,pphi,pphis,presnivs,
[3]	8	. u,v,t,qx,
[1301]	9	. flxmw,
[3]	10	. d_u, d_v, d_t, d_qx, d_ps)
	11
	12	c======================================================================
	13	c
	14	c Modifications pour la physique de Venus
	15	c S. Lebonnois (LMD/CNRS) Septembre 2005
	16	c
	17	c ---------------------------------------------------------------------
	18	c Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
	19	c
	20	c Objet: Moniteur general de la physique du modele
	21	cAA Modifications quant aux traceurs :
	22	cAA - uniformisation des parametrisations ds phytrac
	23	cAA - stockage des moyennes des champs necessaires
	24	cAA en mode traceur off-line
	25	c modif ( P. Le Van , 12/10/98 )
	26	c
	27	c Arguments:
	28	c
	29	c nlon----input-I-nombre de points horizontaux
	30	c nlev----input-I-nombre de couches verticales
	31	c nqmax---input-I-nombre de traceurs
	32	c debut---input-L-variable logique indiquant le premier passage
	33	c lafin---input-L-variable logique indiquant le dernier passage
	34	c rjour---input-R-numero du jour de l'experience
[953]	35	c gmtime--input-R-fraction de la journee (0 a 1)
[3]	36	c pdtphys-input-R-pas d'integration pour la physique (seconde)
	37	c paprs---input-R-pression pour chaque inter-couche (en Pa)
	38	c pplay---input-R-pression pour le mileu de chaque couche (en Pa)
	39	c ppk ---input-R-fonction d'Exner au milieu de couche
	40	c pphi----input-R-geopotentiel de chaque couche (g z) (reference sol)
	41	c pphis---input-R-geopotentiel du sol
	42	c presnivs-input_R_pressions approximat. des milieux couches ( en PA)
	43	c u-------input-R-vitesse dans la direction X (de O a E) en m/s
	44	c v-------input-R-vitesse Y (de S a N) en m/s
	45	c t-------input-R-temperature (K)
	46	c qx------input-R-mass mixing ratio traceurs (kg/kg)
	47	c d_t_dyn-input-R-tendance dynamique pour "t" (K/s)
[1301]	48	c flxmw---input-R-flux de masse vertical en kg/s
[3]	49	c
	50	c d_u-----output-R-tendance physique de "u" (m/s/s)
	51	c d_v-----output-R-tendance physique de "v" (m/s/s)
	52	c d_t-----output-R-tendance physique de "t" (K/s)
	53	c d_qx----output-R-tendance physique de "qx" (kg/kg/s)
	54	c d_ps----output-R-tendance physique de la pression au sol
	55	c======================================================================
[101]	56	USE ioipsl
[849]	57	! USE histcom ! not needed; histcom is included in ioipsl
[101]	58	USE infotrac
	59	USE control_mod
	60	use dimphy
	61	USE comgeomphy
[892]	62	USE mod_phys_lmdz_para, only : is_parallel,jj_nb
	63	USE phys_state_var_mod ! Variables sauvegardees de la physique
[973]	64	USE write_field_phy
[892]	65	USE iophy
[1017]	66	use cpdet_mod, only: cpdet, t2tpot
[1310]	67	USE chemparam_mod
	68	USE conc
	69	USE compo_hedin83_mod2
	70	use moyzon_mod, only: tmoy
[1305]	71	use ieee_arithmetic
[101]	72	IMPLICIT none
	73	c======================================================================
	74	c CLEFS CPP POUR LES IO
	75	c =====================
	76	c#define histhf
	77	#define histday
	78	#define histmth
	79	#define histins
	80	c======================================================================
[3]	81	#include "dimensions.h"
	82	integer jjmp1
	83	parameter (jjmp1=jjm+1-1/jjm)
	84	#include "dimsoil.h"
	85	#include "clesphys.h"
	86	#include "temps.h"
	87	#include "iniprint.h"
	88	#include "timerad.h"
	89	#include "logic.h"
[815]	90	#include "tabcontrol.h"
[1310]	91	#include "nirdata.h"
	92	#include "hedin.h"
[3]	93	c======================================================================
	94	LOGICAL ok_journe ! sortir le fichier journalier
	95	save ok_journe
	96	c PARAMETER (ok_journe=.true.)
	97	c
	98	LOGICAL ok_mensuel ! sortir le fichier mensuel
	99	save ok_mensuel
	100	c PARAMETER (ok_mensuel=.true.)
	101	c
	102	LOGICAL ok_instan ! sortir le fichier instantane
	103	save ok_instan
	104	c PARAMETER (ok_instan=.true.)
	105	c
	106	c======================================================================
	107	c
	108	c Variables argument:
	109	c
	110	INTEGER nlon
	111	INTEGER nlev
	112	INTEGER nqmax
	113	REAL rjourvrai
	114	REAL gmtime
	115	REAL pdtphys
	116	LOGICAL debut, lafin
	117	REAL paprs(klon,klev+1)
	118	REAL pplay(klon,klev)
	119	REAL pphi(klon,klev)
	120	REAL pphis(klon)
	121	REAL presnivs(klev)
	122
	123	! ADAPTATION GCM POUR CP(T)
	124	REAL ppk(klon,klev)
	125
	126	REAL u(klon,klev)
	127	REAL v(klon,klev)
	128	REAL t(klon,klev)
	129	REAL qx(klon,klev,nqmax)
	130
	131	REAL d_u_dyn(klon,klev)
	132	REAL d_t_dyn(klon,klev)
	133
[1301]	134	REAL flxmw(klon,klev)
[3]	135
	136	REAL d_u(klon,klev)
	137	REAL d_v(klon,klev)
	138	REAL d_t(klon,klev)
	139	REAL d_qx(klon,klev,nqmax)
	140	REAL d_ps(klon)
	141
	142	logical ok_hf
	143	real ecrit_hf
	144	integer nid_hf
	145	save ok_hf, ecrit_hf, nid_hf
	146
	147	#ifdef histhf
	148	data ok_hf,ecrit_hf/.true.,0.25/
	149	#else
	150	data ok_hf/.false./
	151	#endif
	152
	153	c Variables propres a la physique
	154	c
[101]	155	INTEGER,save :: itap ! compteur pour la physique
[3]	156	REAL delp(klon,klev) ! epaisseur d'une couche
[1301]	157	REAL omega(klon,klev) ! vitesse verticale en Pa/s
	158
[3]	159
	160	INTEGER igwd,idx(klon),itest(klon)
	161	c
	162	c Diagnostiques 2D de drag_noro, lift_noro et gw_nonoro
	163
	164	REAL zulow(klon),zvlow(klon)
	165	REAL zustrdr(klon), zvstrdr(klon)
	166	REAL zustrli(klon), zvstrli(klon)
	167	REAL zustrhi(klon), zvstrhi(klon)
	168
	169	c Pour calcul GW drag oro et nonoro: CALCUL de N2:
	170	real zdzlev(klon,klev)
	171	real ztlev(klon,klev),zpklev(klon,klev)
	172	real ztetalay(klon,klev),ztetalev(klon,klev)
	173	real zdtetalev(klon,klev)
	174	real zn2(klon,klev) ! BV^2 at plev
	175
	176	c Pour les bilans de moment angulaire,
	177	integer bilansmc
	178	c Pour le transport de ballons
	179	integer ballons
	180	c j'ai aussi besoin
	181	c du stress de couche limite a la surface:
	182
	183	REAL zustrcl(klon),zvstrcl(klon)
	184
	185	c et du stress total c de la physique:
	186
	187	REAL zustrph(klon),zvstrph(klon)
	188
	189	c Variables locales:
	190	c
	191	REAL cdragh(klon) ! drag coefficient pour T and Q
	192	REAL cdragm(klon) ! drag coefficient pour vent
	193	c
	194	cAA Pour TRACEURS
	195	cAA
[105]	196	REAL,save,allocatable :: source(:,:)
[3]	197	REAL ycoefh(klon,klev) ! coef d'echange pour phytrac
	198	REAL yu1(klon) ! vents dans la premiere couche U
	199	REAL yv1(klon) ! vents dans la premiere couche V
	200
	201	REAL sens(klon), dsens(klon) ! chaleur sensible et sa derivee
	202	REAL ve(klon) ! integr. verticale du transport meri. de l'energie
	203	REAL vq(klon) ! integr. verticale du transport meri. de l'eau
	204	REAL ue(klon) ! integr. verticale du transport zonal de l'energie
	205	REAL uq(klon) ! integr. verticale du transport zonal de l'eau
	206	c
[1305]	207	REAL Fsedim(klon,klev+1) ! Flux de sedimentation (kg.m-2)
[3]	208
	209	c======================================================================
	210	c
	211	c Declaration des procedures appelees
	212	c
	213	EXTERNAL ajsec ! ajustement sec
	214	EXTERNAL clmain ! couche limite
	215	EXTERNAL hgardfou ! verifier les temperatures
	216	c EXTERNAL orbite ! calculer l'orbite
	217	EXTERNAL phyetat0 ! lire l'etat initial de la physique
	218	EXTERNAL phyredem ! ecrire l'etat de redemarrage de la physique
	219	EXTERNAL radlwsw ! rayonnements solaire et infrarouge
	220	EXTERNAL suphec ! initialiser certaines constantes
	221	EXTERNAL transp ! transport total de l'eau et de l'energie
	222	EXTERNAL abort_gcm
	223	EXTERNAL printflag
	224	EXTERNAL zenang
	225	EXTERNAL diagetpq
	226	EXTERNAL conf_phys
	227	EXTERNAL diagphy
	228	EXTERNAL mucorr
	229	EXTERNAL phytrac
[1310]	230	EXTERNAL nirco2abs
	231	EXTERNAL nir_leedat
	232	EXTERNAL nltecool
	233	EXTERNAL nlte_tcool
	234	EXTERNAL nlte_setup
	235	EXTERNAL blendrad
	236	EXTERNAL nlthermeq
	237	EXTERNAL euvheat
	238	EXTERNAL param_read
	239	EXTERNAL param_read_e107
	240	EXTERNAL conduction
	241	EXTERNAL molvis
[3]	242	c
	243	c Variables locales
	244	c
	245	CXXX PB
	246	REAL fluxt(klon,klev) ! flux turbulent de chaleur
	247	REAL fluxu(klon,klev) ! flux turbulent de vitesse u
	248	REAL fluxv(klon,klev) ! flux turbulent de vitesse v
	249	c
	250	REAL flux_dyn(klon,klev) ! flux de chaleur produit par la dynamique
	251	REAL flux_ajs(klon,klev) ! flux de chaleur ajustement sec
	252	REAL flux_ec(klon,klev) ! flux de chaleur Ec
	253	c
	254	REAL tmpout(klon,klev) ! K s-1
	255
	256	INTEGER itaprad
	257	SAVE itaprad
	258	c
	259	REAL dist, rmu0(klon), fract(klon)
	260	REAL zdtime, zlongi
	261	c
	262	INTEGER i, k, iq, ig, j, ll
	263	c
	264	REAL zphi(klon,klev)
[1301]	265	REAL zzlev(klon,klev+1),zzlay(klon,klev),z1,z2
[1310]	266	real tsurf(klon)
[1301]	267
[1310]	268	c va avec nlte_tcool
	269	INTEGER ierr_nlte
	270	REAL varerr
	271
[3]	272	c Variables du changement
	273	c
	274	c ajs: ajustement sec
	275	c vdf: couche limite (Vertical DiFfusion)
	276	REAL d_t_ajs(klon,klev), d_tr_ajs(klon,klev,nqmax)
	277	REAL d_u_ajs(klon,klev), d_v_ajs(klon,klev)
	278	c
	279	REAL d_ts(klon)
	280	c
	281	REAL d_u_vdf(klon,klev), d_v_vdf(klon,klev)
	282	REAL d_t_vdf(klon,klev), d_tr_vdf(klon,klev,nqmax)
	283	c
	284	CMOD LOTT: Tendances Orography Sous-maille
	285	REAL d_u_oro(klon,klev), d_v_oro(klon,klev)
	286	REAL d_t_oro(klon,klev)
	287	REAL d_u_lif(klon,klev), d_v_lif(klon,klev)
	288	REAL d_t_lif(klon,klev)
	289	C Tendances Ondes de G non oro (runs strato).
	290	REAL d_u_hin(klon,klev), d_v_hin(klon,klev)
	291	REAL d_t_hin(klon,klev)
	292
[1310]	293	c Tendencies due to radiative scheme [K/s]
	294	c d_t_rad,dtsw,dtlw,d_t_nirco2,d_t_nlte,d_t_euv
	295	c are not computed at each physical timestep
	296	c therefore, they are defined and saved in phys_state_var_mod
	297
	298	c Tendencies due to molecular viscosity and conduction
	299	real d_t_conduc(klon,klev) ! [K/s]
	300	real d_u_molvis(klon,klev) ! (m/s) /s
	301	real d_v_molvis(klon,klev) ! (m/s) /s
	302
[3]	303	c
	304	c Variables liees a l'ecriture de la bande histoire physique
	305	c
	306	INTEGER ecrit_mth
	307	SAVE ecrit_mth ! frequence d'ecriture (fichier mensuel)
	308	c
	309	INTEGER ecrit_day
	310	SAVE ecrit_day ! frequence d'ecriture (fichier journalier)
	311	c
	312	INTEGER ecrit_ins
	313	SAVE ecrit_ins ! frequence d'ecriture (fichier instantane)
	314	c
	315	integer itau_w ! pas de temps ecriture = itap + itau_phy
	316
	317	c Variables locales pour effectuer les appels en serie
	318	c
	319	REAL t_seri(klon,klev)
	320	REAL u_seri(klon,klev), v_seri(klon,klev)
	321	c
[1305]	322	REAL :: tr_seri(klon,klev,nqmax)
	323	REAL :: d_tr(klon,klev,nqmax)
	324
	325	c Variables tendance sedimentation
	326
	327	REAL :: d_tr_sed(klon,klev,2)
	328	REAL :: d_tr_ssed(klon)
[3]	329	c
	330	c pour ioipsl
	331	INTEGER nid_day, nid_mth, nid_ins
	332	SAVE nid_day, nid_mth, nid_ins
	333	INTEGER nhori, nvert, idayref
[97]	334	REAL zsto, zout, zsto1, zsto2, zero
[3]	335	parameter (zero=0.0e0)
	336	real zjulian
	337	save zjulian
	338
	339	CHARACTER*2 str2
	340	character*20 modname
	341	character*80 abort_message
	342	logical ok_sync
	343
	344	character*30 nom_fichier
	345	character*10 varname
	346	character*40 vartitle
	347	character*20 varunits
	348	C Variables liees au bilan d'energie et d'enthalpi
	349	REAL ztsol(klon)
	350	REAL h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
	351	$ , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
	352	SAVE h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
	353	$ , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
	354	REAL d_h_vcol, d_h_dair, d_qt, d_qw, d_ql, d_qs, d_ec
	355	REAL d_h_vcol_phy
	356	REAL fs_bound, fq_bound
	357	SAVE d_h_vcol_phy
	358	REAL zero_v(klon),zero_v2(klon,klev)
	359	CHARACTER*15 ztit
	360	INTEGER ip_ebil ! PRINT level for energy conserv. diag.
	361	SAVE ip_ebil
	362	DATA ip_ebil/2/
	363	INTEGER if_ebil ! level for energy conserv. dignostics
	364	SAVE if_ebil
	365	c+jld ec_conser
	366	REAL d_t_ec(klon,klev) ! tendance du a la conversion Ec -> E thermique
	367	c-jld ec_conser
	368
	369	c TEST VENUS...
	370	REAL mang(klon,klev) ! moment cinetique
	371	REAL mangtot ! moment cinetique total
	372
	373	c Declaration des constantes et des fonctions thermodynamiques
	374	c
	375	#include "YOMCST.h"
	376
	377	c======================================================================
	378	c INITIALISATIONS
	379	c================
	380
	381	modname = 'physiq'
	382	ok_sync=.TRUE.
	383
[119]	384	bilansmc = 0
[3]	385	ballons = 0
[892]	386	! NE FONCTIONNENT PAS ENCORE EN PARALLELE !!!
	387	if (is_parallel) then
	388	bilansmc = 0
	389	ballons = 0
	390	endif
[3]	391
	392	IF (if_ebil.ge.1) THEN
	393	DO i=1,klon
	394	zero_v(i)=0.
	395	END DO
	396	DO i=1,klon
	397	DO j=1,klev
	398	zero_v2(i,j)=0.
	399	END DO
	400	END DO
	401	END IF
	402
	403	c PREMIER APPEL SEULEMENT
	404	c========================
	405	IF (debut) THEN
[105]	406	allocate(source(klon,nqmax))
[3]	407
	408	CALL suphec ! initialiser constantes et parametres phys.
	409
	410	IF (if_ebil.ge.1) d_h_vcol_phy=0.
	411	c
	412	c appel a la lecture du physiq.def
	413	c
	414	call conf_phys(ok_journe, ok_mensuel,
	415	. ok_instan,
	416	. if_ebil)
	417
[892]	418	call phys_state_var_init
[3]	419	c
[1310]	420	c Initialising Hedin model for upper atm
	421	c (to be revised when coupled to chemistry) :
	422	call conc_init
	423	c
[3]	424	c Initialiser les compteurs:
	425	c
	426	itap = 0
	427	itaprad = 0
	428	c
	429	c Lecture startphy.nc :
	430	c
[892]	431	CALL phyetat0 ("startphy.nc")
[3]	432
[815]	433	c dtime est defini dans tabcontrol.h et lu dans startphy
[150]	434	c pdtphys est calcule a partir des nouvelles conditions:
	435	c Reinitialisation du pas de temps physique quand changement
	436	IF (ABS(dtime-pdtphys).GT.0.001) THEN
	437	WRITE(lunout,*) 'Pas physique a change',dtime,
	438	. pdtphys
	439	c abort_message='Pas physique n est pas correct '
	440	c call abort_gcm(modname,abort_message,1)
[152]	441	c----------------
	442	c pour initialiser convenablement le time_counter, il faut tenir compte
	443	c du changement de dtime en changeant itau_phy (point de depart)
	444	itau_phy = NINT(itau_phy*dtime/pdtphys)
	445	c----------------
[150]	446	dtime=pdtphys
	447	ENDIF
	448
[3]	449	radpas = NINT( RDAY/pdtphys/nbapp_rad)
	450
[815]	451	CALL printflag( ok_journe,ok_instan )
[3]	452	c
	453	c---------
	454	c FLOTT
	455	IF (ok_orodr) THEN
	456	DO i=1,klon
	457	rugoro(i) = MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
	458	ENDDO
	459	CALL SUGWD(klon,klev,paprs,pplay)
	460	DO i=1,klon
	461	zuthe(i)=0.
	462	zvthe(i)=0.
	463	if(zstd(i).gt.10.)then
	464	zuthe(i)=(1.-zgam(i))*cos(zthe(i))
	465	zvthe(i)=(1.-zgam(i))*sin(zthe(i))
	466	endif
	467	ENDDO
	468	ENDIF
	469
	470	if (bilansmc.eq.1) then
	471	C OUVERTURE D'UN FICHIER FORMATTE POUR STOCKER LES COMPOSANTES
	472	C DU BILAN DE MOMENT ANGULAIRE.
	473	open(27,file='aaam_bud.out',form='formatted')
	474	open(28,file='fields_2d.out',form='formatted')
	475	write(,)'Ouverture de aaam_bud.out (FL Vous parle)'
	476	write(,)'Ouverture de fields_2d.out (FL Vous parle)'
	477	endif !bilansmc
	478
	479	c--------------SLEBONNOIS
	480	C OUVERTURE DES FICHIERS FORMATTES CONTENANT LES POSITIONS ET VITESSES
	481	C DES BALLONS
	482	if (ballons.eq.1) then
	483	open(30,file='ballons-lat.out',form='formatted')
	484	open(31,file='ballons-lon.out',form='formatted')
	485	open(32,file='ballons-u.out',form='formatted')
	486	open(33,file='ballons-v.out',form='formatted')
	487	open(34,file='ballons-alt.out',form='formatted')
	488	write(,)'Ouverture des ballons*.out'
	489	endif !ballons
	490	c-------------
	491
	492	c---------
	493	C TRACEURS
	494	C source dans couche limite
	495	source = 0.0 ! pas de source, pour l'instant
	496	c---------
[1310]	497
	498	c---------
	499	c INITIALIZE THERMOSPHERIC PARAMETERS
	500	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	501
	502	if (callthermos) then
	503	if(solvarmod.eq.0) call param_read
	504	if(solvarmod.eq.1) call param_read_e107
	505	endif
	506
	507	c Initialisation (recomputed in concentration2)
	508	do ig=1,klon
	509	do j=1,klev
	510	rnew(ig,j)=R
	511	cpnew(ig,j)=cpdet(tmoy(j))
	512	c print*, ' physique l503'
	513	c print*, j, cpdet(tmoy(j))
	514	mmean(ig,j)=RMD
	515	akknew(ig,j)=1.e-4
	516	enddo
	517	c stop
	518
	519	enddo
	520
	521	IF(callthermos.or.callnlte.or.callnirco2) THEN
	522	call compo_hedin83_init2
	523	ENDIF
	524	if (callnlte) call nlte_setup
	525	if(callnirco2.and.(nircorr.eq.1)) call nir_leedat
	526	c---------
	527
[3]	528	c
	529	c Verifications:
	530	c
	531	IF (nlon .NE. klon) THEN
	532	WRITE(lunout,*)'nlon et klon ne sont pas coherents', nlon,
	533	. klon
	534	abort_message='nlon et klon ne sont pas coherents'
	535	call abort_gcm(modname,abort_message,1)
	536	ENDIF
	537	IF (nlev .NE. klev) THEN
	538	WRITE(lunout,*)'nlev et klev ne sont pas coherents', nlev,
	539	. klev
	540	abort_message='nlev et klev ne sont pas coherents'
	541	call abort_gcm(modname,abort_message,1)
	542	ENDIF
	543	c
[808]	544	IF (dtimeREAL(radpas).GT.(RDAY0.25).AND.cycle_diurne)
[3]	545	$ THEN
	546	WRITE(lunout,*)'Nbre d appels au rayonnement insuffisant'
	547	WRITE(lunout,*)"Au minimum 4 appels par jour si cycle diurne"
	548	abort_message='Nbre d appels au rayonnement insuffisant'
	549	call abort_gcm(modname,abort_message,1)
	550	ENDIF
	551	c
	552	WRITE(lunout,*)"Clef pour la convection seche, iflag_ajs=",
	553	. iflag_ajs
	554	c
[1120]	555	ecrit_mth = NINT(RDAY/dtime*ecriphy) ! tous les ecritphy jours
[3]	556	IF (ok_mensuel) THEN
	557	WRITE(lunout,*)'La frequence de sortie mensuelle est de ',
	558	. ecrit_mth
	559	ENDIF
[97]	560
[3]	561	ecrit_day = NINT(RDAY/dtime *1.0) ! tous les jours
	562	IF (ok_journe) THEN
	563	WRITE(lunout,*)'La frequence de sortie journaliere est de ',
	564	. ecrit_day
	565	ENDIF
[97]	566
[1120]	567	ecrit_ins = NINT(RDAY/dtime*ecriphy) ! Fraction de jour reglable
[3]	568	IF (ok_instan) THEN
	569	WRITE(lunout,*)'La frequence de sortie instant. est de ',
	570	. ecrit_ins
	571	ENDIF
	572
	573	c Initialisation des sorties
	574	c===========================
	575
	576	#ifdef CPP_IOIPSL
	577
	578	#ifdef histhf
	579	#include "ini_histhf.h"
	580	#endif
	581
	582	#ifdef histday
	583	#include "ini_histday.h"
	584	#endif
	585
	586	#ifdef histmth
	587	#include "ini_histmth.h"
	588	#endif
	589
	590	#ifdef histins
	591	#include "ini_histins.h"
	592	#endif
	593
	594	#endif
	595
	596	c
	597	c Initialiser les valeurs de u pour calculs tendances
	598	c (pour T, c'est fait dans phyetat0)
	599	c
	600	DO k = 1, klev
	601	DO i = 1, klon
	602	u_ancien(i,k) = u(i,k)
	603	ENDDO
	604	ENDDO
	605
[1305]	606	c---------
	607	c Ecriture fichier initialisation
	608	c PRINT*,'Ecriture Initial_State.csv'
	609	c OPEN(88,file='Trac_Point.csv',
	610	c & form='formatted')
	611	c---------
	612
	613	c---------
	614	c Initialisation des parametres des nuages
	615	c===============================================
	616
	617	if ((nlon .EQ. 1) .AND. ok_cloud) then
	618	PRINT*,'Open profile_cloud_parameters.csv'
	619	OPEN(66,file='profile_cloud_parameters.csv',
	620	& form='formatted')
	621	endif
	622
	623	if ((nlon .EQ. 1) .AND. ok_sedim) then
	624	PRINT*,'Open profile_cloud_sedim.csv'
	625	OPEN(77,file='profile_cloud_sedim.csv',
	626	& form='formatted')
	627	endif
	628
	629	if ((nlon .GT. 1) .AND. ok_chem) then
	630	c !!! DONC 3D !!!
	631	CALL chemparam_ini()
	632	endif
[3]	633
[1305]	634	if ((nlon .GT. 1) .AND. ok_cloud) then
	635	c !!! DONC 3D !!!
	636	CALL cloud_ini(nlon,nlev)
	637	endif
	638
[3]	639	ENDIF ! debut
	640	c====================================================================
	641	c======================================================================
	642
	643	c Mettre a zero des variables de sortie (pour securite)
	644	c
	645	DO i = 1, klon
	646	d_ps(i) = 0.0
	647	ENDDO
	648	DO k = 1, klev
	649	DO i = 1, klon
	650	d_t(i,k) = 0.0
	651	d_u(i,k) = 0.0
	652	d_v(i,k) = 0.0
	653	ENDDO
	654	ENDDO
	655	DO iq = 1, nqmax
	656	DO k = 1, klev
	657	DO i = 1, klon
	658	d_qx(i,k,iq) = 0.0
	659	ENDDO
	660	ENDDO
	661	ENDDO
	662	c
	663	c Ne pas affecter les valeurs entrees de u, v, h, et q
	664	c
	665	DO k = 1, klev
	666	DO i = 1, klon
	667	t_seri(i,k) = t(i,k)
	668	u_seri(i,k) = u(i,k)
	669	v_seri(i,k) = v(i,k)
	670	ENDDO
	671	ENDDO
	672	DO iq = 1, nqmax
	673	DO k = 1, klev
	674	DO i = 1, klon
	675	tr_seri(i,k,iq) = qx(i,k,iq)
	676	ENDDO
	677	ENDDO
	678	ENDDO
	679	C
	680	DO i = 1, klon
	681	ztsol(i) = ftsol(i)
	682	ENDDO
	683	C
	684	IF (if_ebil.ge.1) THEN
	685	ztit='after dynamic'
	686	CALL diagetpq(airephy,ztit,ip_ebil,1,1,dtime
	687	e , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
	688	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
	689	C Comme les tendances de la physique sont ajoute dans la dynamique,
	690	C on devrait avoir que la variation d'entalpie par la dynamique
	691	C est egale a la variation de la physique au pas de temps precedent.
	692	C Donc la somme de ces 2 variations devrait etre nulle.
	693	call diagphy(airephy,ztit,ip_ebil
	694	e , zero_v, zero_v, zero_v, zero_v, zero_v
	695	e , zero_v, zero_v, zero_v, ztsol
	696	e , d_h_vcol+d_h_vcol_phy, d_qt, 0.
	697	s , fs_bound, fq_bound )
	698	END IF
	699
	700	c====================================================================
	701	c Diagnostiquer la tendance dynamique
	702	c
	703	IF (ancien_ok) THEN
	704	DO k = 1, klev
	705	DO i = 1, klon
	706	d_u_dyn(i,k) = (u_seri(i,k)-u_ancien(i,k))/dtime
	707	d_t_dyn(i,k) = (t_seri(i,k)-t_ancien(i,k))/dtime
	708	ENDDO
	709	ENDDO
	710
	711	! ADAPTATION GCM POUR CP(T)
	712	do i=1,klon
	713	flux_dyn(i,1) = 0.0
	714	do j=2,klev
	715	flux_dyn(i,j) = flux_dyn(i,j-1)
	716	. +cpdet(t_seri(i,j-1))/RGd_t_dyn(i,j-1)(paprs(i,j-1)-paprs(i,j))
	717	enddo
	718	enddo
	719
	720	ELSE
	721	DO k = 1, klev
	722	DO i = 1, klon
	723	d_u_dyn(i,k) = 0.0
	724	d_t_dyn(i,k) = 0.0
	725	ENDDO
	726	ENDDO
	727	ancien_ok = .TRUE.
	728	ENDIF
	729	c====================================================================
[1301]	730
	731	c Calcule de vitesse verticale a partir de flux de masse verticale
	732	DO k = 1, klev
	733	DO i = 1, klon
	734	omega(i,k) = RG*flxmw(i,k) / airephy(i)
	735	END DO
	736	END DO
	737
[3]	738	c
	739	c Ajouter le geopotentiel du sol:
	740	c
	741	DO k = 1, klev
	742	DO i = 1, klon
	743	zphi(i,k) = pphi(i,k) + pphis(i)
	744	ENDDO
	745	ENDDO
[1301]	746
	747	c calcul du geopotentiel aux niveaux intercouches
	748	c ponderation des altitudes au niveau des couches en dp/p
	749
	750	DO k=1,klev
	751	DO i=1,klon
[1310]	752	zzlay(i,k)=zphi(i,k)/RG ! [m]
[1301]	753	ENDDO
	754	ENDDO
	755	DO i=1,klon
[1310]	756	zzlev(i,1)=pphis(i)/RG ! [m]
[1301]	757	ENDDO
	758	DO k=2,klev
	759	DO i=1,klon
	760	z1=(pplay(i,k-1)+paprs(i,k))/(pplay(i,k-1)-paprs(i,k))
	761	z2=(paprs(i,k) +pplay(i,k))/(paprs(i,k) -pplay(i,k))
	762	zzlev(i,k)=(z1zzlay(i,k-1)+z2zzlay(i,k))/(z1+z2)
	763	ENDDO
	764	ENDDO
	765	DO i=1,klon
	766	zzlev(i,klev+1)=zzlay(i,klev)+(zzlay(i,klev)-zzlev(i,klev))
	767	ENDDO
	768
[3]	769	c====================================================================
	770	c
	771	c Verifier les temperatures
	772	c
	773	CALL hgardfou(t_seri,ftsol,'debutphy')
	774	c====================================================================
	775	c
	776	c Incrementer le compteur de la physique
	777	c
	778	itap = itap + 1
	779
	780	c====================================================================
	781	c Orbite et eclairement
	782	c====================================================================
	783
	784	c Pour VENUS, on fixe l'obliquite a 0 et l'eccentricite a 0.
	785	c donc pas de variations de Ls, ni de dist.
	786	c La seule chose qui compte, c'est la rotation de la planete devant
	787	c le Soleil...
	788
	789	zlongi = 0.0
	790	dist = 0.72 ! en UA
	791
	792	c Si on veut remettre l'obliquite a 3 degres et/ou l'eccentricite
	793	c a sa valeur, et prendre en compte leur evolution,
	794	c il faudra refaire un orbite.F...
	795	c CALL orbite(zlongi,dist)
	796
	797	IF (cycle_diurne) THEN
[808]	798	zdtime=dtime*REAL(radpas) ! pas de temps du rayonnement (s)
[3]	799	CALL zenang(zlongi,gmtime,zdtime,rlatd,rlond,rmu0,fract)
	800	ELSE
	801	call mucorr(klon,zlongi,rlatd,rmu0,fract)
	802	ENDIF
	803
	804	c====================================================================
	805	c Calcul des tendances traceurs
	806	c====================================================================
	807
	808	if (iflag_trac.eq.1) then
[1160]	809
	810	if (tr_scheme.eq.1) then
	811	! Case 1: pseudo-chemistry with relaxation toward fixed profile
[1305]	812
[1160]	813	call phytrac_relax (debut,lafin,nqmax,
	814	I nlon,nlev,dtime,pplay,
[3]	815	O tr_seri)
[1160]	816
	817	elseif (tr_scheme.eq.2) then
	818	! Case 2: surface emission
	819	! For the moment, inspired from Mars version
	820	! However, the variable 'source' could be used in physiq
	821	! so the call to phytrac_emiss could be to initialise it.
[1305]	822
[1160]	823	call phytrac_emiss ( (rjourvrai+gmtime)*RDAY,
	824	I debut,lafin,nqmax,
	825	I nlon,nlev,dtime,paprs,
	826	I rlatd,rlond,
	827	O tr_seri)
[3]	828
[1305]	829	elseif (tr_scheme.eq.3) then ! identical to ok_chem.or.ok_cloud
	830	! Case 3: Full chemistry and/or clouds
	831
	832	call phytrac_chimie(
	833	I debut,
	834	I gmtime,
	835	I nqmax,
	836	I nlon,
	837	I rlatd,
	838	I rlond,
	839	I nlev,
	840	I dtime,
	841	I t_seri,pplay,
	842	O tr_seri,
	843	O NBRTOT,
	844	O WH2SO4,
	845	O rho_droplet)
	846
	847	c CALL WriteField_phy('Pression',pplay,nlev)
	848	c CALL WriteField_phy('PressionBnd',paprs,nlev+1)
	849	c CALL WriteField_phy('Temp',t_seri,nlev)
	850	c IF (ok_cloud) THEN
	851	c CALL WriteField_phy('NBRTOT',NBRTOT,nlev)
	852	c ENDIF
	853	c CALL WriteField_phy('SAl',tr_seri(:,:,i_h2so4liq),nlev)
	854	c CALL WriteField_phy('SAg',tr_seri(:,:,i_h2so4),nlev)
	855
	856	if (ok_sedim) then
	857
	858	CALL new_cloud_sedim(
	859	I klon,
	860	I nlev,
	861	I dtime,
	862	I pplay,
	863	I paprs,
	864	I t_seri,
	865	I WH2SO4,
	866	I tr_seri,
	867	I nqmax,
	868	I NBRTOT,
	869	I rho_droplet,
	870	O Fsedim,
	871	O d_tr_sed,
	872	O d_tr_ssed)
	873
	874	DO k = 1, klev
	875	DO i = 1, klon
	876
	877	c WRITE(88,"(11(e15.8,','))") pplay(5,25),
	878	c & t_seri(5,25),tr_seri(5,25,i_h2oliq),
	879	c & tr_seri(5,25,i_h2o),tr_seri(5,25,i_h2so4liq),
	880	c & tr_seri(5,25,i_h2so4),NBRTOT(5,25),WH2SO4(5,25),
	881	c & Fsedim(5,25),d_tr_sed(5,25,1),d_tr_sed(5,25,2)
	882
	883	c--------------------
	884	c Ce test est necessaire pour eviter Xliq=NaN
	885	IF (ieee_is_nan(d_tr_sed(i,k,1)).OR.
	886	& ieee_is_nan(d_tr_sed(i,k,2))) THEN
	887	PRINT*,'sedim NaN PROBLEM'
	888	PRINT*,'d_tr_sed Nan?',d_tr_sed(i,k,:),'Temp',t_seri(i,k)
	889	PRINT*,'lat-lon',i,'level',k,'dtime',dtime
	890	PRINT*,'NBRTOT',NBRTOT(i,k),'F_sed',Fsedim(i,k)
	891	PRINT*,'==============================================='
	892	d_tr_sed(i,k,:)=0.
	893	ENDIF
	894	c--------------------
	895
	896	tr_seri(i,k,i_h2so4liq) = tr_seri(i,k,i_h2so4liq)+
	897	& d_tr_sed(i,k,1)
	898	tr_seri(i,k,i_h2oliq) = tr_seri(i,k,i_h2oliq)+
	899	& d_tr_sed(i,k,2)
	900	d_tr_sed(i,k,:) = d_tr_sed(i,k,:) / dtime
	901	Fsedim(i,k) = Fsedim(i,k) / dtime
	902
	903	ENDDO
	904	ENDDO
	905
	906	Fsedim(:,klev+1) = 0.
	907
	908	endif ! ok_sedim
	909
	910	endif ! tr_scheme
	911	endif ! iflag_trac
	912
[3]	913	c
	914	c====================================================================
	915	c Appeler la diffusion verticale (programme de couche limite)
	916	c====================================================================
	917
	918	c-------------------------------
	919	c VENUS TEST: on ne tient pas compte des calculs de clmain mais on force
	920	c l'equilibre radiatif du sol
	921	if (1.eq.0) then
	922	if (debut) then
	923	print*,"ATTENTION, CLMAIN SHUNTEE..."
	924	endif
	925
	926	DO i = 1, klon
	927	sens(i) = 0.0e0 ! flux de chaleur sensible au sol
	928	fder(i) = 0.0e0
	929	dlw(i) = 0.0e0
	930	ENDDO
	931
	932	c Incrementer la temperature du sol
	933	c
	934	DO i = 1, klon
	935	d_ts(i) = dtime * radsol(i)/22000. !valeur calculee par GCM pour I=200
	936	ftsol(i) = ftsol(i) + d_ts(i)
	937	do j=1,nsoilmx
	938	ftsoil(i,j)=ftsol(i)
	939	enddo
	940	ENDDO
	941
	942	c-------------------------------
	943	else
	944	c-------------------------------
	945
	946	fder = dlw
	947
	948	! ADAPTATION GCM POUR CP(T)
[808]	949
[3]	950	CALL clmain(dtime,itap,
	951	e t_seri,u_seri,v_seri,
	952	e rmu0,
	953	e ftsol,
	954	$ ftsoil,
	955	$ paprs,pplay,ppk,radsol,falbe,
	956	e solsw, sollw, sollwdown, fder,
	957	e rlond, rlatd, cuphy, cvphy,
	958	e debut, lafin,
	959	s d_t_vdf,d_u_vdf,d_v_vdf,d_ts,
	960	s fluxt,fluxu,fluxv,cdragh,cdragm,
	961	s dsens,
	962	s ycoefh,yu1,yv1)
	963
	964	CXXX Incrementation des flux
	965	DO i = 1, klon
	966	sens(i) = - fluxt(i,1) ! flux de chaleur sensible au sol
	967	fder(i) = dlw(i) + dsens(i)
	968	ENDDO
	969	CXXX
	970
	971	DO k = 1, klev
	972	DO i = 1, klon
	973	t_seri(i,k) = t_seri(i,k) + d_t_vdf(i,k)
	974	d_t_vdf(i,k)= d_t_vdf(i,k)/dtime ! K/s
	975	u_seri(i,k) = u_seri(i,k) + d_u_vdf(i,k)
	976	d_u_vdf(i,k)= d_u_vdf(i,k)/dtime ! (m/s)/s
	977	v_seri(i,k) = v_seri(i,k) + d_v_vdf(i,k)
	978	d_v_vdf(i,k)= d_v_vdf(i,k)/dtime ! (m/s)/s
	979	ENDDO
	980	ENDDO
	981
	982	C TRACEURS
	983
	984	if (iflag_trac.eq.1) then
[101]	985	DO k = 1, klev
	986	DO i = 1, klon
	987	delp(i,k) = paprs(i,k)-paprs(i,k+1)
	988	ENDDO
	989	ENDDO
[3]	990	DO iq=1, nqmax
	991	CALL cltrac(dtime,ycoefh,t_seri,
	992	s tr_seri(1,1,iq),source,
	993	e paprs, pplay,delp,
	994	s d_tr_vdf(1,1,iq))
	995	tr_seri(:,:,iq) = tr_seri(:,:,iq) + d_tr_vdf(:,:,iq)
	996	d_tr_vdf(:,:,iq)= d_tr_vdf(:,:,iq)/dtime ! /s
	997	ENDDO
	998	endif
	999
	1000	IF (if_ebil.ge.2) THEN
	1001	ztit='after clmain'
	1002	CALL diagetpq(airephy,ztit,ip_ebil,2,1,dtime
	1003	e , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
	1004	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
	1005	call diagphy(airephy,ztit,ip_ebil
	1006	e , zero_v, zero_v, zero_v, zero_v, sens
	1007	e , zero_v, zero_v, zero_v, ztsol
	1008	e , d_h_vcol, d_qt, d_ec
	1009	s , fs_bound, fq_bound )
	1010	END IF
	1011	C
	1012	c
	1013	c Incrementer la temperature du sol
	1014	c
	1015	DO i = 1, klon
	1016	ftsol(i) = ftsol(i) + d_ts(i)
	1017	ENDDO
	1018
	1019	c Calculer la derive du flux infrarouge
	1020	c
	1021	DO i = 1, klon
	1022	dlw(i) = - 4.0RSIGMAftsol(i)**3
	1023	ENDDO
	1024
	1025	c-------------------------------
	1026	endif ! fin du VENUS TEST
	1027
[973]	1028	! tests: output tendencies
	1029	! call writefield_phy('physiq_d_t_vdf',d_t_vdf,klev)
	1030	! call writefield_phy('physiq_d_u_vdf',d_u_vdf,klev)
	1031	! call writefield_phy('physiq_d_v_vdf',d_v_vdf,klev)
	1032	! call writefield_phy('physiq_d_ts',d_ts,1)
	1033
[3]	1034	c
	1035	c Appeler l'ajustement sec
	1036	c
	1037	c===================================================================
	1038	c Convection seche
	1039	c===================================================================
	1040	c
	1041	d_t_ajs(:,:)=0.
	1042	d_u_ajs(:,:)=0.
	1043	d_v_ajs(:,:)=0.
	1044	d_tr_ajs(:,:,:)=0.
	1045	c
	1046	IF(prt_level>9)WRITE(lunout,*)
	1047	. 'AVANT LA CONVECTION SECHE , iflag_ajs='
	1048	s ,iflag_ajs
	1049
	1050	if(iflag_ajs.eq.0) then
	1051	c Rien
	1052	c ====
	1053	IF(prt_level>9)WRITE(lunout,*)'pas de convection'
	1054
	1055	else if(iflag_ajs.eq.1) then
	1056
	1057	c Ajustement sec
	1058	c ==============
	1059	IF(prt_level>9)WRITE(lunout,*)'ajsec'
	1060
	1061	! ADAPTATION GCM POUR CP(T)
	1062	CALL ajsec(paprs, pplay, ppk, t_seri, u_seri, v_seri, nqmax,
	1063	. tr_seri, d_t_ajs, d_u_ajs, d_v_ajs, d_tr_ajs)
	1064
	1065	! ADAPTATION GCM POUR CP(T)
	1066	do i=1,klon
	1067	flux_ajs(i,1) = 0.0
	1068	do j=2,klev
	1069	flux_ajs(i,j) = flux_ajs(i,j-1)
	1070	. + cpdet(t_seri(i,j-1))/RG*d_t_ajs(i,j-1)/dtime
	1071	. *(paprs(i,j-1)-paprs(i,j))
	1072	enddo
	1073	enddo
	1074
	1075	t_seri(:,:) = t_seri(:,:) + d_t_ajs(:,:)
	1076	d_t_ajs(:,:)= d_t_ajs(:,:)/dtime ! K/s
	1077	u_seri(:,:) = u_seri(:,:) + d_u_ajs(:,:)
	1078	d_u_ajs(:,:)= d_u_ajs(:,:)/dtime ! (m/s)/s
	1079	v_seri(:,:) = v_seri(:,:) + d_v_ajs(:,:)
	1080	d_v_ajs(:,:)= d_v_ajs(:,:)/dtime ! (m/s)/s
[1301]	1081
	1082	if (iflag_trac.eq.1) then
[3]	1083	tr_seri(:,:,:) = tr_seri(:,:,:) + d_tr_ajs(:,:,:)
	1084	d_tr_ajs(:,:,:)= d_tr_ajs(:,:,:)/dtime ! /s
[1301]	1085	endif
[3]	1086
	1087	endif
[973]	1088
	1089	! tests: output tendencies
	1090	! call writefield_phy('physiq_d_t_ajs',d_t_ajs,klev)
	1091	! call writefield_phy('physiq_d_u_ajs',d_u_ajs,klev)
	1092	! call writefield_phy('physiq_d_v_ajs',d_v_ajs,klev)
[3]	1093	c
	1094	IF (if_ebil.ge.2) THEN
	1095	ztit='after dry_adjust'
	1096	CALL diagetpq(airephy,ztit,ip_ebil,2,2,dtime
	1097	e , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
	1098	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
	1099	call diagphy(airephy,ztit,ip_ebil
	1100	e , zero_v, zero_v, zero_v, zero_v, sens
	1101	e , zero_v, zero_v, zero_v, ztsol
	1102	e , d_h_vcol, d_qt, d_ec
	1103	s , fs_bound, fq_bound )
	1104	END IF
	1105
	1106	c====================================================================
	1107	c RAYONNEMENT
	1108	c====================================================================
	1109
[1310]	1110	c-----------------------------------------------------------------------
	1111	c . Compute radiative tendencies :
	1112	c------------------------------------
	1113	c====================================================================
[3]	1114	IF (MOD(itaprad,radpas).EQ.0) THEN
[1310]	1115	c====================================================================
[3]	1116
[808]	1117	dtimerad = dtime*REAL(radpas) ! pas de temps du rayonnement (s)
[3]	1118	c PRINT*,'dtimerad,dtime,radpas',dtimerad,dtime,radpas
	1119
[1310]	1120	c Calcul pour Cp rnew et mmean avec traceurs (Cp independant de T !! )
	1121	IF(callnlte.or.callthermos) THEN
	1122	call compo_hedin83_mod(pplay,rmu0,
	1123	& co2vmr_gcm,covmr_gcm,ovmr_gcm,n2vmr_gcm,nvmr_gcm)
	1124	ENDIF
	1125
	1126	if(callthermos) then
	1127	call concentrations2(pplay,t_seri,d_t,co2vmr_gcm, n2vmr_gcm,
	1128	& covmr_gcm, ovmr_gcm,nvmr_gcm,pdtphys)
	1129	endif
	1130
	1131
	1132	c NLTE cooling from CO2 emission
	1133	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1134
	1135	IF(callnlte) THEN
	1136	if(nltemodel.eq.0.or.nltemodel.eq.1) then
	1137	CALL nltecool(klon, klev, pplay*9.869e-6, t_seri,
	1138	$ co2vmr_gcm,n2vmr_gcm, covmr_gcm, ovmr_gcm,
	1139	$ d_t_nlte)
	1140	else if(nltemodel.eq.2) then
	1141	CALL nlte_tcool(klon,klev,pplay*9.869e-6,
	1142	$ t_seri,zzlay,co2vmr_gcm, n2vmr_gcm, covmr_gcm,
	1143	$ ovmr_gcm,d_t_nlte,ierr_nlte,varerr )
	1144	if(ierr_nlte.gt.0) then
	1145	write(,)
	1146	$ 'WARNING: nlte_tcool output with error message',
	1147	$ 'ierr_nlte=',ierr_nlte,'varerr=',varerr
	1148	write(,)'I will continue anyway'
	1149	endif
	1150
	1151	endif
	1152
	1153	ELSE
	1154
	1155	d_t_nlte(:,:)=0.
	1156
	1157	ENDIF
	1158
	1159	c Find number of layers for LTE radiation calculations
	1160
	1161	IF(callnlte .or. callnirco2)
	1162	$ CALL nlthermeq(klon, klev, paprs, pplay)
	1163
	1164
	1165	c LTE radiative transfert / solar / IR matrix
	1166	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1167
[3]	1168	CALL radlwsw
[1301]	1169	e (dist, rmu0, fract, zzlev,
[1310]	1170	e paprs, pplay,ftsol, t_seri)
[808]	1171
[1301]	1172
[1310]	1173	c CO2 near infrared absorption
	1174	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1175
	1176	d_t_nirco2(:,:)=0.
	1177	if (callnirco2) then
	1178	call nirco2abs (klon, klev, pplay, dist, nqmax, qx,
	1179	. rmu0, fract, d_t_nirco2,
	1180	. co2vmr_gcm, ovmr_gcm)
	1181	endif
	1182
	1183
	1184	c Net atmospheric radiative heating rate (K.s-1)
	1185	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1186
	1187	IF(callnlte.or.callnirco2) THEN
	1188	CALL blendrad(klon, klev, pplay,heat,
	1189	& cool, d_t_nirco2,d_t_nlte, dtsw, dtlw)
	1190	ELSE
	1191	dtsw(:,:)=heat(:,:)
	1192	dtlw(:,:)=-1*cool(:,:)
	1193	ENDIF
	1194
	1195	DO k=1,klev
	1196	DO i=1,klon
	1197	d_t_rad(i,k) = dtsw(i,k) + dtlw(i,k) ! K/s
	1198	ENDDO
	1199	ENDDO
	1200
	1201
	1202	cc---------------------------------------------
	1203
	1204	c EUV heating rate (K.s-1)
	1205	c ~~~~~~~~~~~~~~~~~~~~~~~~
	1206
	1207	d_t_euv(:,:)=0.
	1208
	1209	IF (callthermos) THEN
	1210
	1211	c call thermosphere(zplev,zplay,dist_sol,
	1212	c $ mu0,ptimestep,ptime,zday,tsurf,zzlev,zzlay,
	1213	c & pt,pq,pu,pv,pdt,pdq,
	1214	c $ zdteuv,zdtconduc,zdumolvis,zdvmolvis,zdqmoldiff)
	1215
	1216	call euvheat(klon, klev, t_seri,paprs,pplay,zzlay,
	1217	$ rmu0,pdtphys,gmtime,rjourvrai,
	1218	C $ pq,pdq,zdteuv)
	1219	$ co2vmr_gcm, n2vmr_gcm, covmr_gcm,
	1220	$ ovmr_gcm,nvmr_gcm,d_t_euv )
	1221
	1222	DO k=1,klev
	1223	DO ig=1,klon
	1224	d_t_rad(ig,k)=d_t_rad(ig,k)+d_t_euv(ig,k)
	1225
	1226	ENDDO
	1227	ENDDO
	1228
	1229	ENDIF ! callthermos
	1230
[3]	1231	c====================================================================
[1310]	1232	itaprad = 0
	1233	ENDIF ! radpas
	1234	c====================================================================
[3]	1235	c
	1236	c Ajouter la tendance des rayonnements (tous les pas)
	1237	c
	1238	DO k = 1, klev
	1239	DO i = 1, klon
[1310]	1240	t_seri(i,k) = t_seri(i,k) + d_t_rad(i,k) * dtime
[3]	1241	ENDDO
	1242	ENDDO
[973]	1243
[1310]	1244	! CONDUCTION and MOLECULAR VISCOSITY
	1245	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1246
	1247	d_t_conduc(:,:)=0.
	1248	d_u_molvis(:,:)=0.
	1249	d_v_molvis(:,:)=0.
	1250
	1251	IF (callthermos) THEN
	1252
	1253	tsurf(:)=t_seri(:,1)
	1254	call conduction(klon, klev,pdtphys,
	1255	$ pplay,paprs,t_seri,
	1256	$ tsurf,zzlev,zzlay,d_t_conduc)
	1257
	1258	call molvis(klon, klev,pdtphys,
	1259	$ pplay,paprs,t_seri,
	1260	$ u,tsurf,zzlev,zzlay,d_u_molvis)
	1261
	1262	call molvis(klon, klev, pdtphys,
	1263	$ pplay,paprs,t_seri,
	1264	$ v,tsurf,zzlev,zzlay,d_u_molvis)
	1265
	1266	DO k=1,klev
	1267	DO ig=1,klon
	1268	t_seri(ig,k)= t_seri(ig,k)+ d_t_conduc(ig,k)*dtime ! [K]
	1269	u_seri(ig,k)= u_seri(ig,k)+ d_u_molvis(ig,k)*dtime ! m/s
	1270	v_seri(ig,k)= v_seri(ig,k)+ d_v_molvis(ig,k)*dtime ! m/s
	1271	ENDDO
	1272	ENDDO
	1273
	1274	ENDIF ! callthermos
	1275
	1276	c====================================================================
[973]	1277	! tests: output tendencies
	1278	! call writefield_phy('physiq_dtrad',dtrad,klev)
[3]	1279
	1280	IF (if_ebil.ge.2) THEN
	1281	ztit='after rad'
	1282	CALL diagetpq(airephy,ztit,ip_ebil,2,2,dtime
	1283	e , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
	1284	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
	1285	call diagphy(airephy,ztit,ip_ebil
	1286	e , topsw, toplw, solsw, sollw, zero_v
	1287	e , zero_v, zero_v, zero_v, ztsol
	1288	e , d_h_vcol, d_qt, d_ec
	1289	s , fs_bound, fq_bound )
	1290	END IF
	1291	c
	1292
	1293	c====================================================================
	1294	c Calcul des gravity waves FLOTT
	1295	c====================================================================
	1296	c
	1297	if (ok_orodr.or.ok_gw_nonoro) then
	1298	c CALCUL DE N2
	1299	do i=1,klon
	1300	do k=2,klev
	1301	ztlev(i,k) = (t_seri(i,k)+t_seri(i,k-1))/2.
	1302	zpklev(i,k) = sqrt(ppk(i,k)*ppk(i,k-1))
	1303	enddo
	1304	enddo
	1305	call t2tpot(klon*klev,ztlev, ztetalev,zpklev)
	1306	call t2tpot(klon*klev,t_seri,ztetalay,ppk)
	1307	do i=1,klon
	1308	do k=2,klev
	1309	zdtetalev(i,k) = ztetalay(i,k)-ztetalay(i,k-1)
	1310	zdzlev(i,k) = (zphi(i,k)-zphi(i,k-1))/RG
	1311	zn2(i,k) = RGzdtetalev(i,k)/(ztetalev(i,k)zdzlev(i,k))
	1312	zn2(i,k) = max(zn2(i,k),1.e-12) ! securite
	1313	enddo
[808]	1314	zn2(i,1) = 1.e-12 ! securite
[3]	1315	enddo
	1316
	1317	endif
	1318
	1319	c ----------------------------ORODRAG
	1320	IF (ok_orodr) THEN
	1321	c
	1322	c selection des points pour lesquels le shema est actif:
	1323	igwd=0
	1324	DO i=1,klon
	1325	itest(i)=0
	1326	c IF ((zstd(i).gt.10.0)) THEN
	1327	IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
	1328	itest(i)=1
	1329	igwd=igwd+1
	1330	idx(igwd)=i
	1331	ENDIF
	1332	ENDDO
	1333	c igwdim=MAX(1,igwd)
	1334	c
	1335	c A ADAPTER POUR VENUS!!!
[1306]	1336	CALL drag_noro(klon,klev,dtime,paprs,pplay,pphi,zn2,
[3]	1337	e zmea,zstd, zsig, zgam, zthe,zpic,zval,
	1338	e igwd,idx,itest,
	1339	e t_seri, u_seri, v_seri,
	1340	s zulow, zvlow, zustrdr, zvstrdr,
	1341	s d_t_oro, d_u_oro, d_v_oro)
	1342
[808]	1343	c print*,"d_u_oro=",d_u_oro(klon/2,:)
[3]	1344	c ajout des tendances
	1345	t_seri(:,:) = t_seri(:,:) + d_t_oro(:,:)
	1346	d_t_oro(:,:)= d_t_oro(:,:)/dtime ! K/s
	1347	u_seri(:,:) = u_seri(:,:) + d_u_oro(:,:)
	1348	d_u_oro(:,:)= d_u_oro(:,:)/dtime ! (m/s)/s
	1349	v_seri(:,:) = v_seri(:,:) + d_v_oro(:,:)
	1350	d_v_oro(:,:)= d_v_oro(:,:)/dtime ! (m/s)/s
	1351	c
	1352	ELSE
	1353	d_t_oro = 0.
	1354	d_u_oro = 0.
	1355	d_v_oro = 0.
	1356	zustrdr = 0.
	1357	zvstrdr = 0.
	1358	c
	1359	ENDIF ! fin de test sur ok_orodr
	1360	c
[973]	1361	! tests: output tendencies
	1362	! call writefield_phy('physiq_d_t_oro',d_t_oro,klev)
	1363	! call writefield_phy('physiq_d_u_oro',d_u_oro,klev)
	1364	! call writefield_phy('physiq_d_v_oro',d_v_oro,klev)
	1365
[3]	1366	c ----------------------------OROLIFT
	1367	IF (ok_orolf) THEN
[808]	1368	print*,"ok_orolf NOT IMPLEMENTED !"
	1369	stop
[3]	1370	c
	1371	c selection des points pour lesquels le shema est actif:
	1372	igwd=0
	1373	DO i=1,klon
	1374	itest(i)=0
	1375	IF ((zpic(i)-zmea(i)).GT.100.) THEN
	1376	itest(i)=1
	1377	igwd=igwd+1
	1378	idx(igwd)=i
	1379	ENDIF
	1380	ENDDO
	1381	c igwdim=MAX(1,igwd)
	1382	c
	1383	c A ADAPTER POUR VENUS!!!
	1384	c CALL lift_noro(klon,klev,dtime,paprs,pplay,
	1385	c e rlatd,zmea,zstd,zpic,zgam,zthe,zpic,zval,
	1386	c e igwd,idx,itest,
	1387	c e t_seri, u_seri, v_seri,
	1388	c s zulow, zvlow, zustrli, zvstrli,
	1389	c s d_t_lif, d_u_lif, d_v_lif )
	1390
	1391	c
	1392	c ajout des tendances
	1393	t_seri(:,:) = t_seri(:,:) + d_t_lif(:,:)
	1394	d_t_lif(:,:)= d_t_lif(:,:)/dtime ! K/s
	1395	u_seri(:,:) = u_seri(:,:) + d_u_lif(:,:)
	1396	d_u_lif(:,:)= d_u_lif(:,:)/dtime ! (m/s)/s
	1397	v_seri(:,:) = v_seri(:,:) + d_v_lif(:,:)
	1398	d_v_lif(:,:)= d_v_lif(:,:)/dtime ! (m/s)/s
	1399	c
	1400	ELSE
	1401	d_t_lif = 0.
	1402	d_u_lif = 0.
	1403	d_v_lif = 0.
	1404	zustrli = 0.
	1405	zvstrli = 0.
	1406	c
	1407	ENDIF ! fin de test sur ok_orolf
	1408
	1409	c ---------------------------- NON-ORO GRAVITY WAVES
	1410	IF(ok_gw_nonoro) then
	1411
	1412	call flott_gwd_ran(klon,klev,dtime,pplay,zn2,
	1413	e t_seri, u_seri, v_seri,
	1414	o zustrhi,zvstrhi,
	1415	o d_t_hin, d_u_hin, d_v_hin)
	1416
	1417	c ajout des tendances
	1418
	1419	t_seri(:,:) = t_seri(:,:) + d_t_hin(:,:)
	1420	d_t_hin(:,:)= d_t_hin(:,:)/dtime ! K/s
	1421	u_seri(:,:) = u_seri(:,:) + d_u_hin(:,:)
	1422	d_u_hin(:,:)= d_u_hin(:,:)/dtime ! (m/s)/s
	1423	v_seri(:,:) = v_seri(:,:) + d_v_hin(:,:)
	1424	d_v_hin(:,:)= d_v_hin(:,:)/dtime ! (m/s)/s
	1425
	1426	ELSE
	1427	d_t_hin = 0.
	1428	d_u_hin = 0.
	1429	d_v_hin = 0.
	1430	zustrhi = 0.
	1431	zvstrhi = 0.
	1432
	1433	ENDIF ! fin de test sur ok_gw_nonoro
	1434
[973]	1435	! tests: output tendencies
	1436	! call writefield_phy('physiq_d_t_hin',d_t_hin,klev)
	1437	! call writefield_phy('physiq_d_u_hin',d_u_hin,klev)
	1438	! call writefield_phy('physiq_d_v_hin',d_v_hin,klev)
	1439
[3]	1440	c====================================================================
	1441	c Transport de ballons
	1442	c====================================================================
	1443	if (ballons.eq.1) then
[953]	1444	CALL ballon(30,pdtphys,rjourvrai,gmtime*RDAY,rlatd,rlond,
[3]	1445	c C t,pplay,u,v,pphi) ! alt above surface (smoothed for GCM)
	1446	C t,pplay,u,v,zphi) ! alt above planet average radius
	1447	endif !ballons
	1448
	1449	c====================================================================
	1450	c Bilan de mmt angulaire
	1451	c====================================================================
	1452	if (bilansmc.eq.1) then
	1453	CMODDEB FLOTT
	1454	C CALCULER LE BILAN DE MOMENT ANGULAIRE (DIAGNOSTIQUE)
	1455	C STRESS NECESSAIRES: COUCHE LIMITE ET TOUTE LA PHYSIQUE
	1456
	1457	DO i = 1, klon
	1458	zustrph(i)=0.
	1459	zvstrph(i)=0.
	1460	zustrcl(i)=0.
	1461	zvstrcl(i)=0.
	1462	ENDDO
	1463	DO k = 1, klev
	1464	DO i = 1, klon
	1465	zustrph(i)=zustrph(i)+(u_seri(i,k)-u(i,k))/dtime*
	1466	c (paprs(i,k)-paprs(i,k+1))/rg
	1467	zvstrph(i)=zvstrph(i)+(v_seri(i,k)-v(i,k))/dtime*
	1468	c (paprs(i,k)-paprs(i,k+1))/rg
	1469	zustrcl(i)=zustrcl(i)+d_u_vdf(i,k)*
	1470	c (paprs(i,k)-paprs(i,k+1))/rg
	1471	zvstrcl(i)=zvstrcl(i)+d_v_vdf(i,k)*
	1472	c (paprs(i,k)-paprs(i,k+1))/rg
	1473	ENDDO
	1474	ENDDO
	1475
[953]	1476	CALL aaam_bud (27,klon,klev,rjourvrai,gmtime*RDAY,
[3]	1477	C ra,rg,romega,
	1478	C rlatd,rlond,pphis,
	1479	C zustrdr,zustrli,zustrcl,
	1480	C zvstrdr,zvstrli,zvstrcl,
	1481	C paprs,u,v)
	1482
	1483	CCMODFIN FLOTT
	1484	endif !bilansmc
	1485
	1486	c====================================================================
	1487	c====================================================================
	1488	c Calculer le transport de l'eau et de l'energie (diagnostique)
	1489	c
	1490	c A REVOIR POUR VENUS...
	1491	c
	1492	c CALL transp (paprs,ftsol,
	1493	c e t_seri, q_seri, u_seri, v_seri, zphi,
	1494	c s ve, vq, ue, uq)
	1495	c
	1496	c====================================================================
	1497	c+jld ec_conser
	1498	DO k = 1, klev
	1499	DO i = 1, klon
	1500	d_t_ec(i,k)=0.5/cpdet(t_seri(i,k))
	1501	$ (u(i,k)2+v(i,k)2-u_seri(i,k)2-v_seri(i,k)*2)
	1502	t_seri(i,k)=t_seri(i,k)+d_t_ec(i,k)
	1503	d_t_ec(i,k) = d_t_ec(i,k)/dtime
	1504	END DO
	1505	END DO
	1506	do i=1,klon
	1507	flux_ec(i,1) = 0.0
	1508	do j=2,klev
	1509	flux_ec(i,j) = flux_ec(i,j-1)
	1510	. +cpdet(t_seri(i,j-1))/RGd_t_ec(i,j-1)(paprs(i,j-1)-paprs(i,j))
	1511	enddo
	1512	enddo
	1513
	1514	c-jld ec_conser
	1515	c====================================================================
	1516	IF (if_ebil.ge.1) THEN
	1517	ztit='after physic'
	1518	CALL diagetpq(airephy,ztit,ip_ebil,1,1,dtime
	1519	e , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
	1520	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
	1521	C Comme les tendances de la physique sont ajoute dans la dynamique,
	1522	C on devrait avoir que la variation d'entalpie par la dynamique
	1523	C est egale a la variation de la physique au pas de temps precedent.
	1524	C Donc la somme de ces 2 variations devrait etre nulle.
	1525	call diagphy(airephy,ztit,ip_ebil
	1526	e , topsw, toplw, solsw, sollw, sens
	1527	e , zero_v, zero_v, zero_v, ztsol
	1528	e , d_h_vcol, d_qt, d_ec
	1529	s , fs_bound, fq_bound )
	1530	C
	1531	d_h_vcol_phy=d_h_vcol
	1532	C
	1533	END IF
	1534	C
	1535	c=======================================================================
	1536	c SORTIES
	1537	c=======================================================================
	1538
	1539	c Convertir les incrementations en tendances
	1540	c
	1541	DO k = 1, klev
	1542	DO i = 1, klon
	1543	d_u(i,k) = ( u_seri(i,k) - u(i,k) ) / dtime
	1544	d_v(i,k) = ( v_seri(i,k) - v(i,k) ) / dtime
	1545	d_t(i,k) = ( t_seri(i,k) - t(i,k) ) / dtime
	1546	ENDDO
	1547	ENDDO
	1548	c
	1549	DO iq = 1, nqmax
	1550	DO k = 1, klev
	1551	DO i = 1, klon
	1552	d_qx(i,k,iq) = ( tr_seri(i,k,iq) - qx(i,k,iq) ) / dtime
	1553	ENDDO
	1554	ENDDO
	1555	ENDDO
	1556
	1557	c------------------------
	1558	c Calcul moment cinetique
	1559	c------------------------
	1560	c TEST VENUS...
	1561	c mangtot = 0.0
	1562	c DO k = 1, klev
	1563	c DO i = 1, klon
	1564	c mang(i,k) = RAcos(rlatd(i)RPI/180.)
	1565	c . (u_seri(i,k)+RAcos(rlatd(i)RPI/180.)ROMEGA)
	1566	c . airephy(i)(paprs(i,k)-paprs(i,k+1))/RG
	1567	c mangtot=mangtot+mang(i,k)
	1568	c ENDDO
	1569	c ENDDO
	1570	c print*,"Moment cinetique total = ",mangtot
	1571	c
	1572	c------------------------
	1573	c
	1574	c Sauvegarder les valeurs de t et u a la fin de la physique:
	1575	c
	1576	DO k = 1, klev
	1577	DO i = 1, klon
	1578	u_ancien(i,k) = u_seri(i,k)
	1579	t_ancien(i,k) = t_seri(i,k)
	1580	ENDDO
	1581	ENDDO
	1582	c
	1583	c=============================================================
	1584	c Ecriture des sorties
	1585	c=============================================================
	1586
	1587	#ifdef CPP_IOIPSL
	1588
	1589	#ifdef histhf
	1590	#include "write_histhf.h"
	1591	#endif
	1592
	1593	#ifdef histday
	1594	#include "write_histday.h"
	1595	#endif
	1596
	1597	#ifdef histmth
	1598	#include "write_histmth.h"
	1599	#endif
	1600
	1601	#ifdef histins
	1602	#include "write_histins.h"
	1603	#endif
	1604
	1605	#endif
	1606
	1607	c====================================================================
	1608	c Si c'est la fin, il faut conserver l'etat de redemarrage
	1609	c====================================================================
	1610	c
	1611	IF (lafin) THEN
	1612	itau_phy = itau_phy + itap
[892]	1613	CALL phyredem ("restartphy.nc")
[3]	1614
	1615	c--------------FLOTT
	1616	CMODEB LOTT
	1617	C FERMETURE DU FICHIER FORMATTE CONTENANT LES COMPOSANTES
	1618	C DU BILAN DE MOMENT ANGULAIRE.
	1619	if (bilansmc.eq.1) then
	1620	write(,)'Fermeture de aaam_bud.out (FL Vous parle)'
	1621	close(27)
	1622	close(28)
	1623	endif !bilansmc
	1624	CMODFIN
	1625	c-------------
	1626	c--------------SLEBONNOIS
	1627	C FERMETURE DES FICHIERS FORMATTES CONTENANT LES POSITIONS ET VITESSES
	1628	C DES BALLONS
	1629	if (ballons.eq.1) then
	1630	write(,)'Fermeture des ballons*.out'
	1631	close(30)
	1632	close(31)
	1633	close(32)
	1634	close(33)
	1635	close(34)
	1636	endif !ballons
	1637	c-------------
	1638	ENDIF
	1639
	1640	RETURN
	1641	END
	1642
	1643
	1644
	1645	***********************************************************************
	1646	***********************************************************************
	1647	***********************************************************************
	1648	***********************************************************************
	1649	***********************************************************************
	1650	***********************************************************************
	1651	***********************************************************************
	1652	***********************************************************************
	1653
	1654	SUBROUTINE gr_fi_ecrit(nfield,nlon,iim,jjmp1,fi,ecrit)
	1655	IMPLICIT none
	1656	c
	1657	c Tranformer une variable de la grille physique a
	1658	c la grille d'ecriture
	1659	c
	1660	INTEGER nfield,nlon,iim,jjmp1, jjm
	1661	REAL fi(nlon,nfield), ecrit(iim*jjmp1,nfield)
	1662	c
	1663	INTEGER i, n, ig
	1664	c
	1665	jjm = jjmp1 - 1
	1666	DO n = 1, nfield
	1667	DO i=1,iim
	1668	ecrit(i,n) = fi(1,n)
	1669	ecrit(i+jjm*iim,n) = fi(nlon,n)
	1670	ENDDO
	1671	DO ig = 1, nlon - 2
	1672	ecrit(iim+ig,n) = fi(1+ig,n)
	1673	ENDDO
	1674	ENDDO
	1675	RETURN
	1676	END

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

Download in other formats:

Original Format