Context Navigation

physiq.F @ 1442

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

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

Download in other formats:

Original Format