Context Navigation

physiq.F @ 1242

Last change on this file since 1242 was 1126, checked in by slebonnois, 11 years ago
SL: update of Titan photochemical module to include computation of chemistry up to 1300 km
File size: 50.4 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,
	9	. omega,
	10	. d_u, d_v, d_t, d_qx, d_ps)
	11
	12	c======================================================================
	13	c
	14	c Modifications pour la physique de Titan
[1056]	15	c S. Lebonnois (LMD/CNRS) Juin 2013: Parallelisation
[3]	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
[1056]	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
	35	c gmtime--input-R-temps universel dans la journee (0 a RDAY s)
	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)
	48	c omega---input-R-vitesse verticale en Pa/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======================================================================
[102]	56	USE ioipsl
[849]	57	! USE histcom ! not needed; histcom is included in ioipsl
[102]	58	USE infotrac
	59	USE control_mod
	60	use dimphy
	61	USE comgeomphy
[1048]	62	use cpdet_mod, only: cpdet, t2tpot
[1056]	63	USE mod_phys_lmdz_para, only : is_parallel,jj_nb
	64	USE phys_state_var_mod ! Variables sauvegardees de la physique
	65	USE iophy
	66	USE common_mod, only: rmcbar,xfbar,ncount,
	67	& flxesp_i,tau_drop,tau_aer,solesp,precip,
	68	& evapch4,occcld_m,occcld,satch4,satc2h6,satc2h2,rmcloud,
	69	& TauHID,TauHVD,TauGID,TauGVD,TauCID,TauCVD,NSPECV,NSPECI,
	70	& common_init
	71
	72	USE moyzon_mod
	73	USE write_field_phy
[102]	74	IMPLICIT none
	75	c======================================================================
	76	c CLEFS CPP POUR LES IO
	77	c =====================
[1056]	78	#define histday
[102]	79	#define histmth
	80	#define histins
	81	c======================================================================
[3]	82	#include "dimensions.h"
	83	integer jjmp1
	84	parameter (jjmp1=jjm+1-1/jjm)
	85	#include "dimsoil.h"
	86	#include "clesphys.h"
	87	#include "temps.h"
	88	#include "iniprint.h"
	89	#include "logic.h"
[1056]	90	#include "tabcontrol.h"
[3]	91	#include "comorbit.h"
[175]	92	#include "microtab.h"
	93	#include "itemps.h"
[3]	94	c======================================================================
[1056]	95	LOGICAL ok_journe ! sortir le fichier journalier
	96	save ok_journe
	97	c PARAMETER (ok_journe=.true.)
	98	c
[3]	99	LOGICAL ok_mensuel ! sortir le fichier mensuel
	100	save ok_mensuel
	101	c PARAMETER (ok_mensuel=.true.)
	102	c
	103	LOGICAL ok_instan ! sortir le fichier instantane
	104	save ok_instan
	105	c PARAMETER (ok_instan=.true.)
	106	c
	107	c======================================================================
	108	c
	109	c Variables argument:
	110	c
	111	INTEGER nlon
	112	INTEGER nlev
	113	INTEGER nqmax
	114	REAL rjourvrai
	115	REAL gmtime
	116	REAL pdtphys
	117	LOGICAL debut, lafin
	118	REAL paprs(klon,klev+1)
	119	REAL pplay(klon,klev)
	120	REAL pphi(klon,klev)
	121	REAL pphis(klon)
	122	REAL presnivs(klev)
	123
	124	! ADAPTATION GCM POUR CP(T)
	125	REAL ppk(klon,klev)
	126
	127	REAL u(klon,klev)
	128	REAL v(klon,klev)
	129	REAL t(klon,klev)
	130	REAL qx(klon,klev,nqmax)
	131
	132	REAL d_u_dyn(klon,klev)
	133	REAL d_t_dyn(klon,klev)
	134
	135	REAL omega(klon,klev)
	136
	137	REAL d_u(klon,klev)
	138	REAL d_v(klon,klev)
	139	REAL d_t(klon,klev)
	140	REAL d_qx(klon,klev,nqmax)
	141	REAL d_ps(klon)
	142
	143	c Variables propres a la physique
	144	c
[105]	145	REAL,save,allocatable :: rlev(:,:) ! altitude a chaque niveau (interface inferieure de la couche)
	146	INTEGER,save :: itap ! compteur pour la physique
[3]	147	REAL delp(klon,klev) ! epaisseur d'une couche
	148
	149	INTEGER igwd,idx(klon),itest(klon)
	150	c
	151	c Diagnostiques 2D de drag_noro, lift_noro et gw_nonoro
	152
	153	REAL zulow(klon),zvlow(klon)
	154	REAL zustrdr(klon), zvstrdr(klon)
	155	REAL zustrli(klon), zvstrli(klon)
	156	REAL zustrhi(klon), zvstrhi(klon)
	157
	158	c Pour calcul GW drag oro et nonoro: CALCUL de N2:
	159	real zdzlev(klon,klev)
	160	real ztlev(klon,klev),zpklev(klon,klev)
	161	real ztetalay(klon,klev),ztetalev(klon,klev)
	162	real zdtetalev(klon,klev)
	163	real zn2(klon,klev) ! BV^2 at plev
	164
	165	c Pour les bilans de moment angulaire,
	166	integer bilansmc
	167	c Pour le transport de ballons
	168	integer ballons
	169	c j'ai aussi besoin
	170	c du stress de couche limite a la surface:
	171
	172	REAL zustrcl(klon),zvstrcl(klon)
	173
	174	c et du stress total c de la physique:
	175
	176	REAL zustrph(klon),zvstrph(klon)
	177
	178	c Variables locales:
	179	c
	180	REAL cdragh(klon) ! drag coefficient pour T and Q
	181	REAL cdragm(klon) ! drag coefficient pour vent
	182	c
[1056]	183	cAA Pour TRACEURS
[3]	184	cAA
[105]	185	REAL,save,allocatable :: source(:,:)
[3]	186	integer nmicro
[105]	187	save nmicro
[1056]	188	character*8 nom
	189	REAL qaer(klon,klev,nqmax)
[3]	190
	191	REAL ycoefh(klon,klev) ! coef d'echange pour phytrac
	192	REAL yu1(klon) ! vents dans la premiere couche U
	193	REAL yv1(klon) ! vents dans la premiere couche V
	194
	195	REAL sens(klon), dsens(klon) ! chaleur sensible et sa derivee
	196	REAL ve(klon) ! integr. verticale du transport meri. de l'energie
	197	REAL vq(klon) ! integr. verticale du transport meri. de l'eau
	198	REAL ue(klon) ! integr. verticale du transport zonal de l'energie
	199	REAL uq(klon) ! integr. verticale du transport zonal de l'eau
	200	c
	201
	202	c======================================================================
	203	c
	204	c Declaration des procedures appelees
	205	c
	206	EXTERNAL ajsec ! ajustement sec
	207	EXTERNAL clmain ! couche limite
	208	EXTERNAL hgardfou ! verifier les temperatures
	209	EXTERNAL orbite ! calculer l'orbite
	210	EXTERNAL phyetat0 ! lire l'etat initial de la physique
	211	EXTERNAL phyredem ! ecrire l'etat de redemarrage de la physique
	212	EXTERNAL radlwsw ! rayonnements solaire et infrarouge
	213	EXTERNAL suphec ! initialiser certaines constantes
[1056]	214	c EXTERNAL transp ! transport total de l'eau et de l'energie
[3]	215	EXTERNAL abort_gcm
	216	EXTERNAL printflag
	217	EXTERNAL zenang
	218	EXTERNAL diagetpq
	219	EXTERNAL conf_phys
	220	EXTERNAL diagphy
	221	EXTERNAL mucorr
	222	EXTERNAL phytrac
	223	c
	224	c Variables locales
	225	c
	226	CXXX PB
	227	REAL fluxt(klon,klev) ! flux turbulent de chaleur
	228	REAL fluxu(klon,klev) ! flux turbulent de vitesse u
	229	REAL fluxv(klon,klev) ! flux turbulent de vitesse v
	230	c
	231	REAL flux_dyn(klon,klev) ! flux de chaleur produit par la dynamique
	232	REAL flux_ajs(klon,klev) ! flux de chaleur ajustement sec
	233	REAL flux_ec(klon,klev) ! flux de chaleur Ec
	234	c
[815]	235	REAL dtimerad
[3]	236	INTEGER itaprad
[815]	237	SAVE itaprad,dtimerad
[3]	238	REAL zdtime
	239	c
	240	c CHIMIE
	241
	242	REAL dtimechim
	243	INTEGER itapchim,appel_chim
	244	SAVE itapchim,dtimechim
	245
	246	c ORBITE
	247
	248	REAL dist, rmu0(klon), fract(klon), pdecli
[1056]	249	REAL rmu0bar(klon), fractbar(klon)
[3]	250	REAL zday
[854]	251	REAL zls,zlsdeg,zlsm1
[1048]	252	save zlsm1
[3]	253	c
	254	INTEGER i, k, iq, ig, j, ll, l
	255	c
	256	REAL zphi(klon,klev)
	257	REAL zzlev(klon,klev+1),zzlay(klon,klev),z1,z2
	258	c
	259	c Variables du changement
	260	c
	261	c ajs: ajustement sec
	262	c vdf: couche limite (Vertical DiFfusion)
	263	c mph: microphysique
	264	c kim: chimie
	265	REAL d_t_ajs(klon,klev), d_tr_ajs(klon,klev,nqmax)
	266	REAL d_u_ajs(klon,klev), d_v_ajs(klon,klev)
	267	c
	268	REAL d_ts(klon)
	269	c
	270	REAL d_u_vdf(klon,klev), d_v_vdf(klon,klev)
	271	REAL d_t_vdf(klon,klev), d_tr_vdf(klon,klev,nqmax)
	272	c
	273	REAL d_tr_mph(klon,klev,nqmax),d_tr_kim(klon,klev,nqmax)
	274
	275	CMOD LOTT: Tendances Orography Sous-maille
	276	REAL d_u_oro(klon,klev), d_v_oro(klon,klev)
	277	REAL d_t_oro(klon,klev)
	278	REAL d_u_lif(klon,klev), d_v_lif(klon,klev)
	279	REAL d_t_lif(klon,klev)
	280	C Tendances Ondes de G non oro (runs strato).
	281	REAL d_u_hin(klon,klev), d_v_hin(klon,klev)
	282	REAL d_t_hin(klon,klev)
	283
	284	c
	285	c Variables liees a l'ecriture de la bande histoire physique
	286	c
	287	INTEGER ecrit_mth
	288	SAVE ecrit_mth ! frequence d'ecriture (fichier mensuel)
	289	c
	290	INTEGER ecrit_day
	291	SAVE ecrit_day ! frequence d'ecriture (fichier journalier)
	292	c
	293	INTEGER ecrit_ins
	294	SAVE ecrit_ins ! frequence d'ecriture (fichier instantane)
	295	c
	296	integer itau_w ! pas de temps ecriture = itap + itau_phy
	297
	298	c Variables locales pour effectuer les appels en serie
	299	c
	300	REAL t_seri(klon,klev)
	301	REAL u_seri(klon,klev), v_seri(klon,klev)
	302	c
	303	REAL tr_seri(klon,klev,nqmax)
[1056]	304	REAL d_tr(klon,klev,nqmax)
[3]	305	c
	306	c pour ioipsl
	307	INTEGER nid_day, nid_mth, nid_ins
	308	SAVE nid_day, nid_mth, nid_ins
	309	INTEGER nhori, nvert, idayref
[97]	310	REAL zsto, zout, zsto1, zsto2, zero
[3]	311	parameter (zero=0.0e0)
	312	real zjulian
	313	save zjulian
[1056]	314	REAL tmpout(klon,klev) ! pour sorties
[3]	315
[1056]	316	CHARACTER*1 str1
[3]	317	CHARACTER*2 str2
	318	character*20 modname
	319	character*80 abort_message
	320	logical ok_sync
	321
	322	character*30 nom_fichier
	323	character*10 varname
	324	character*40 vartitle
	325	character*20 varunits
	326	C Variables liees au bilan d'energie et d'enthalpi
	327	REAL ztsol(klon)
	328	REAL h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
	329	$ , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
	330	SAVE h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
	331	$ , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
	332	REAL d_h_vcol, d_h_dair, d_qt, d_qw, d_ql, d_qs, d_ec
	333	REAL d_h_vcol_phy
	334	REAL fs_bound, fq_bound
	335	SAVE d_h_vcol_phy
	336	REAL zero_v(klon),zero_v2(klon,klev)
	337	CHARACTER*15 ztit
	338	INTEGER ip_ebil ! PRINT level for energy conserv. diag.
	339	SAVE ip_ebil
	340	DATA ip_ebil/2/
	341	INTEGER if_ebil ! level for energy conserv. dignostics
	342	SAVE if_ebil
	343	c+jld ec_conser
	344	REAL d_t_ec(klon,klev) ! tendance du a la conversion Ec -> E thermique
	345	c-jld ec_conser
	346
	347	c TEST VENUS...
	348	REAL mang(klon,klev) ! moment cinetique
	349	REAL mangtot ! moment cinetique total
	350
[97]	351	c Temporaire avant de trouver mieux :
	352	c Recuperation des TAU du TR
[105]	353	REAL t_tauhvd(klon,klev),t_khvd(klon,klev)
[175]	354	REAL t_tcld(klon,klev),t_kcld(klon,klev)
	355	REAL t_kcvd(klon,klev)
[97]	356
[175]	357	REAL ch4(klon,jjm+1),dch4(jjm+1)
	358	INTEGER ig0
	359	integer ich4
	360	common/ch4ind/ich4
	361
	362	c flux de chaleur latente d'evaporation CH4
	363	REAL fclat(klon)
	364	c reservoir de surface
	365	REAL,save,allocatable :: reservoir(:)
	366
[3]	367	c Declaration des constantes et des fonctions thermodynamiques
	368	c
	369	#include "YOMCST.h"
	370
	371	c======================================================================
	372	c INITIALISATIONS
[1056]	373	c================
[3]	374
	375	modname = 'physiq'
	376	ok_sync=.TRUE.
	377
[106]	378	bilansmc = 0
[3]	379	ballons = 0
[1056]	380	! NE FONCTIONNENT PAS ENCORE EN PARALLELE !!!
	381	if (is_parallel) then
	382	bilansmc = 0
	383	ballons = 0
	384	endif
[3]	385
	386	IF (if_ebil.ge.1) THEN
	387	DO i=1,klon
	388	zero_v(i)=0.
	389	END DO
	390	DO i=1,klon
	391	DO j=1,klev
	392	zero_v2(i,j)=0.
	393	END DO
	394	END DO
	395	END IF
	396
	397	c PREMIER APPEL SEULEMENT
	398	c========================
	399	IF (debut) THEN
[1056]	400	allocate(rlev(klon,klevp1))
[105]	401	allocate(source(klon,nqmax))
[175]	402	allocate(reservoir(klon))
[3]	403
	404	CALL suphec ! initialiser constantes et parametres phys.
	405
	406	IF (if_ebil.ge.1) d_h_vcol_phy=0.
	407	c
	408	c appel a la lecture du physiq.def
	409	c
[1071]	410	call conf_phys(ok_mensuel,ok_journe,
[1056]	411	. ok_instan,
	412	. if_ebil)
[3]	413
[1056]	414	call phys_state_var_init
	415	call common_init
[3]	416	c
	417	c Initialiser les compteurs:
	418	c
[1056]	419	itap = 0
	420	itaprad = 0
[1126]	421	itapchim = 1
[1056]	422
	423	c init rnuabar
[175]	424	ncount(:,:) = 0
[1056]	425	rmcbar = 0.
	426	xfbar = 0.
	427
[3]	428	c
	429	c Lecture startphy.nc :
	430	c
[1056]	431	CALL phyetat0 ("startphy.nc")
[3]	432
[1056]	433	c dtime est defini dans tabcontrol.h et lu dans startphy
[149]	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----------------
[149]	446	dtime=pdtphys
	447	ENDIF
	448
[3]	449	radpas = NINT( RDAY/pdtphys/nbapp_rad)
	450	chimpas = radpas*nbapp_rad/nbapp_chim
	451
[1056]	452	CALL printflag( ok_mensuel,ok_journe,ok_instan )
[3]	453	c
	454	c Initialiser les pas de temps:
	455	c
[1056]	456	dtimerad = dtime*REAL(radpas) ! pas de temps du rayonnement (s)
[3]	457	c PRINT*,'dtimerad,dtime,radpas',dtimerad,dtime,radpas
	458
[1056]	459	dtimechim = dtime*REAL(chimpas) ! pas de temps de la chimie (s)
[3]	460	c PRINT*,'dtimechim,dtime,chimpas',dtimechim,dtime,chimpas
	461
	462
	463	c INITIALISATION ORBITE
	464
	465	CALL iniorbit(aphelie,periheli,year_day,peri_day,obliquit)
	466
	467	c---------
	468	c FLOTT
[1056]	469	IF (ok_orodr) THEN
	470	DO i=1,klon
	471	rugoro(i) = MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
	472	ENDDO
	473	CALL SUGWD(klon,klev,paprs,pplay)
	474	DO i=1,klon
	475	zuthe(i)=0.
	476	zvthe(i)=0.
	477	if(zstd(i).gt.10.)then
	478	zuthe(i)=(1.-zgam(i))*cos(zthe(i))
	479	zvthe(i)=(1.-zgam(i))*sin(zthe(i))
	480	endif
	481	ENDDO
	482	ENDIF
[3]	483
	484	if (bilansmc.eq.1) then
	485	C OUVERTURE D'UN FICHIER FORMATTE POUR STOCKER LES COMPOSANTES
	486	C DU BILAN DE MOMENT ANGULAIRE.
	487	open(27,file='aaam_bud.out',form='formatted')
	488	open(28,file='fields_2d.out',form='formatted')
	489	write(,)'Ouverture de aaam_bud.out (FL Vous parle)'
	490	write(,)'Ouverture de fields_2d.out (FL Vous parle)'
	491	endif !bilansmc
	492
	493	c--------------SLEBONNOIS
	494	C OUVERTURE DES FICHIERS FORMATTES CONTENANT LES POSITIONS ET VITESSES
	495	C DES BALLONS
	496	if (ballons.eq.1) then
	497	open(30,file='ballons-lat.out',form='formatted')
	498	open(31,file='ballons-lon.out',form='formatted')
	499	open(32,file='ballons-u.out',form='formatted')
	500	open(33,file='ballons-v.out',form='formatted')
	501	open(34,file='ballons-alt.out',form='formatted')
	502	write(,)'Ouverture des ballons*.out'
	503	endif !ballons
	504	c-------------
	505
	506	c---------
	507	C TRACEURS
	508	C source dans couche limite
[1056]	509	source = 0.0 ! pas de source, pour l'instant
[3]	510	C
	511	c Si microphysique offline, pas besoin d'avoir de traceurs microphysiques
	512	c car on lit les profils verticaux des qaer dans une look-up table pour
	513	c le rayonnement.
	514
	515	c calcul de nmicro
	516	c !!!! Les traceurs microphysiques doivent etre toujours en premiers!!
	517
	518	nmicro = 0
	519	do iq=1,nqmax
	520	nom = tname(iq)
[102]	521	c print*,iq,"nom=",nom,"tname=",tname(iq)
[3]	522	print*,iq,"nom=",nom
	523	if (nom(1:1).eq."q") then
	524	nmicro = nmicro+1
	525	endif
	526	enddo
	527	print*,"nmicro=",nmicro
	528
[1056]	529	c --------------
[3]	530	c Verifications:
[1056]	531	c --------------
[175]	532	IF ((nmicro.eq.0).and.(microfi.eq.1)) THEN
	533	abort_message="MICROPHYSIQUE ONLINE, MAIS NMICRO=0..."
	534	call abort_gcm(modname,abort_message,1)
	535	ENDIF
	536	IF (microfi.lt.1.and.clouds.eq.1) THEN
	537	write(lunout,*)"microfi.lt.1.and.clouds.eq.1"
	538	abort_message =
	539	& "Impossible de faire des nuages sans microphysique..."
	540	call abort_gcm(modname,abort_message,1)
	541	ENDIF
[3]	542	IF (nlon .NE. klon) THEN
	543	WRITE(lunout,*)'nlon et klon ne sont pas coherents', nlon,
	544	. klon
	545	abort_message='nlon et klon ne sont pas coherents'
	546	call abort_gcm(modname,abort_message,1)
	547	ENDIF
	548	IF (nlev .NE. klev) THEN
	549	WRITE(lunout,*)'nlev et klev ne sont pas coherents', nlev,
	550	. klev
	551	abort_message='nlev et klev ne sont pas coherents'
	552	call abort_gcm(modname,abort_message,1)
	553	ENDIF
[1056]	554
	555	IF (((moyzon_mu).and.(microfi.ne.1)).or.
	556	. ((.not.moyzon_mu).and.(microfi.eq.1))) THEN
	557	abort_message="Microphysic 2D and moyzon_mu not compatible"
	558	write(lunout,*) "moyzon_mu=",moyzon_mu
	559	write(lunout,*) "microfi=",microfi
	560	call abort_gcm(modname,abort_message,1)
	561	ENDIF
	562	IF (((moyzon_ch).and.(.not.chimi)).or.
	563	. ((.not.moyzon_ch).and.(chimi))) THEN
	564	abort_message="Chemistry and moyzon_ch not compatible"
	565	write(lunout,*) "moyzon_ch=",moyzon_ch
	566	write(lunout,*) "chimi=",chimi
	567	call abort_gcm(modname,abort_message,1)
	568	ENDIF
	569
	570	IF (dtimeREAL(radpas).GT.(RDAY0.25).AND.cycle_diurne)
[3]	571	$ THEN
	572	WRITE(lunout,*)'Nbre d appels au rayonnement insuffisant'
	573	WRITE(lunout,*)"Au minimum 4 appels par jour si cycle diurne"
	574	abort_message='Nbre d appels au rayonnement insuffisant'
	575	call abort_gcm(modname,abort_message,1)
	576	ENDIF
	577	c
	578	WRITE(lunout,*)"Clef pour la convection seche, iflag_ajs=",
	579	. iflag_ajs
	580	c
[97]	581	ecrit_mth = NINT(RDAY/dtime) *nday ! tous les nday jours
[3]	582	IF (ok_mensuel) THEN
[1056]	583	WRITE(lunout,*)'La frequence de sortie mensuelle est de ',
[3]	584	. ecrit_mth
	585	ENDIF
	586
	587	ecrit_day = NINT(RDAY/dtime *1.0) ! tous les jours
	588	IF (ok_journe) THEN
	589	WRITE(lunout,*)'La frequence de sortie journaliere est de ',
	590	. ecrit_day
	591	ENDIF
	592
[1120]	593	ecrit_ins = NINT(RDAY/dtime*ecriphy) ! Fraction de jour reglable
[3]	594	IF (ok_instan) THEN
	595	WRITE(lunout,*)'La frequence de sortie instant. est de ',
	596	. ecrit_ins
	597	ENDIF
	598
	599	c Initialisation des sorties
	600	c===========================
	601
	602	#ifdef CPP_IOIPSL
	603
[1056]	604	#ifdef histday
	605	#include "ini_histday.h"
	606	#endif
	607
[3]	608	#ifdef histmth
	609	#include "ini_histmth.h"
	610	#endif
	611
	612	#ifdef histins
	613	#include "ini_histins.h"
	614	#endif
	615
	616	#endif
	617
	618	c
	619	c Initialiser les valeurs de u pour calculs tendances
	620	c (pour T, c'est fait dans phyetat0)
	621	c
	622	DO k = 1, klev
	623	DO i = 1, klon
	624	u_ancien(i,k) = u(i,k)
	625	ENDDO
	626	ENDDO
	627
	628	ENDIF ! debut
	629	c====================================================================
	630	c======================================================================
	631
[175]	632	c Creer un reservoir de surface infini
	633	c
	634	reservoir(:) = 2.
	635
[3]	636	c Mettre a zero des variables de sortie (pour securite)
	637	c
	638	DO i = 1, klon
	639	d_ps(i) = 0.0
	640	ENDDO
	641	DO k = 1, klev
	642	DO i = 1, klon
	643	d_t(i,k) = 0.0
	644	d_u(i,k) = 0.0
	645	d_v(i,k) = 0.0
	646	ENDDO
	647	ENDDO
	648	DO iq = 1, nqmax
	649	DO k = 1, klev
	650	DO i = 1, klon
	651	d_qx(i,k,iq) = 0.0
	652	ENDDO
	653	ENDDO
	654	ENDDO
	655	c
	656	c Ne pas affecter les valeurs entrees de u, v, h, et q
	657	c
	658	DO k = 1, klev
	659	DO i = 1, klon
	660	t_seri(i,k) = t(i,k)
	661	u_seri(i,k) = u(i,k)
	662	v_seri(i,k) = v(i,k)
	663	ENDDO
	664	ENDDO
	665	DO iq = 1, nqmax
	666	DO k = 1, klev
	667	DO i = 1, klon
	668	tr_seri(i,k,iq) = qx(i,k,iq)
	669	ENDDO
	670	ENDDO
	671	ENDDO
	672	C
	673	DO i = 1, klon
[1056]	674	ztsol(i) = ftsol(i)
[3]	675	ENDDO
	676	C
	677	IF (if_ebil.ge.1) THEN
	678	ztit='after dynamic'
	679	CALL diagetpq(airephy,ztit,ip_ebil,1,1,dtime
	680	e , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
	681	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
	682	C Comme les tendances de la physique sont ajoute dans la dynamique,
	683	C on devrait avoir que la variation d'entalpie par la dynamique
	684	C est egale a la variation de la physique au pas de temps precedent.
	685	C Donc la somme de ces 2 variations devrait etre nulle.
	686	call diagphy(airephy,ztit,ip_ebil
	687	e , zero_v, zero_v, zero_v, zero_v, zero_v
	688	e , zero_v, zero_v, zero_v, ztsol
	689	e , d_h_vcol+d_h_vcol_phy, d_qt, 0.
	690	s , fs_bound, fq_bound )
	691	END IF
	692
	693	c====================================================================
	694	c Diagnostiquer la tendance dynamique
	695	c
	696	IF (ancien_ok) THEN
	697	DO k = 1, klev
	698	DO i = 1, klon
	699	d_u_dyn(i,k) = (u_seri(i,k)-u_ancien(i,k))/dtime
	700	d_t_dyn(i,k) = (t_seri(i,k)-t_ancien(i,k))/dtime
	701	ENDDO
	702	ENDDO
	703
	704	! ADAPTATION GCM POUR CP(T)
	705	do i=1,klon
	706	flux_dyn(i,1) = 0.0
	707	do j=2,klev
	708	flux_dyn(i,j) = flux_dyn(i,j-1)
	709	. +cpdet(t_seri(i,j-1))/RGd_t_dyn(i,j-1)(paprs(i,j-1)-paprs(i,j))
	710	enddo
	711	enddo
	712
	713	ELSE
	714	DO k = 1, klev
	715	DO i = 1, klon
	716	d_u_dyn(i,k) = 0.0
	717	d_t_dyn(i,k) = 0.0
	718	ENDDO
	719	ENDDO
	720	ancien_ok = .TRUE.
	721	ENDIF
	722	c====================================================================
	723	c
	724	c Ajouter le geopotentiel du sol:
	725	c
	726	DO k = 1, klev
	727	DO i = 1, klon
	728	zphi(i,k) = pphi(i,k) + pphis(i)
	729	ENDDO
	730	ENDDO
	731
[1056]	732	c call WriteField_phy('physiq_pphi',pphi,klev)
	733	c call WriteField_phy('physiq_pphis',pphis,1)
	734
[3]	735	c calcul du geopotentiel aux niveaux intercouches
	736	c ponderation des altitudes au niveau des couches en dp/p
	737
	738	DO l=1,klev
	739	DO i=1,klon
[887]	740	c zzlay(i,l)=zphi(i,l)/RG
[106]	741	c SI ON TIENT COMPTE DE LA VARIATION DE G AVEC L'ALTITUDE:
[887]	742	zzlay(i,l)=RGRARA/(RG*RA-zphi(i,l))-RA
[3]	743	ENDDO
	744	ENDDO
	745	DO i=1,klon
[106]	746	c zzlev(i,1)=0.
	747	c CORRECTION 13/01/2011
	748	c (correspond a la position de la surface en ce point vs RA)
	749	zzlev(i,1)=pphis(i)/RG
[3]	750	ENDDO
	751	DO l=2,klev
	752	DO i=1,klon
	753	z1=(pplay(i,l-1)+paprs(i,l))/(pplay(i,l-1)-paprs(i,l))
[1056]	754	z2=(paprs(i,l) +pplay(i,l))/(paprs(i,l) -pplay(i,l))
[3]	755	zzlev(i,l)=(z1zzlay(i,l-1)+z2zzlay(i,l))/(z1+z2)
	756	ENDDO
	757	ENDDO
	758	DO i=1,klon
	759	zzlev(i,klev+1)=zzlay(i,klev)+(zzlay(i,klev)-zzlev(i,klev))
	760	ENDDO
	761
[1056]	762	! zonal averages needed
	763	if (moyzon_ch.or.moyzon_mu) then
	764
	765	c zzlaybar(1,:)=(zphibar(1,:)+zphisbar(1))/RG
	766	c SI ON TIENT COMPTE DE LA VARIATION DE G AVEC L'ALTITUDE:
	767	zzlaybar(1,:)=RGRARA/(RG*RA-(zphibar(1,:)+zphisbar(1)))-RA
	768	zzlevbar(1,1)=zphisbar(1)/RG
	769	DO l=2,klev
	770	z1=(zplaybar(1,l-1)+zplevbar(1,l))/
	771	. (zplevbar(1,l-1)-zplevbar(1,l))
	772	z2=(zplevbar(1,l) +zplaybar(1,l))/
	773	. (zplevbar(1,l) -zplaybar(1,l))
	774	zzlevbar(1,l)=(z1zzlaybar(1,l-1)+z2zzlaybar(1,l))/(z1+z2)
	775	ENDDO
	776	zzlevbar(1,klev+1)=zzlaybar(1,klev)+
	777	. (zzlaybar(1,klev)-zzlevbar(1,klev))
	778
	779	DO i=2,klon
	780	if (rlatd(i).ne.rlatd(i-1)) then
	781	DO l=1,klev
	782	c zzlaybar(i,l)=(zphibar(i,l)+zphisbar(i))/RG
	783	c SI ON TIENT COMPTE DE LA VARIATION DE G AVEC L'ALTITUDE:
	784	zzlaybar(i,l)=RGRARA/(RG*RA-(zphibar(i,l)+zphisbar(i)))-RA
	785	ENDDO
	786	zzlevbar(i,1)=zphisbar(i)/RG
	787	DO l=2,klev
	788	z1=(zplaybar(i,l-1)+zplevbar(i,l))/
	789	. (zplevbar(i,l-1)-zplevbar(i,l))
	790	z2=(zplevbar(i,l) +zplaybar(i,l))/
	791	. (zplevbar(i,l) -zplaybar(i,l))
	792	zzlevbar(i,l)=(z1zzlaybar(i,l-1)+z2zzlaybar(i,l))/(z1+z2)
	793	ENDDO
	794	zzlevbar(i,klev+1)=zzlaybar(i,klev)+
	795	. (zzlaybar(i,klev)-zzlevbar(i,klev))
	796	else
	797	zzlaybar(i,:)=zzlaybar(i-1,:)
	798	zzlevbar(i,:)=zzlevbar(i-1,:)
	799	endif
	800	ENDDO
	801
	802	endif ! moyzon
	803
	804	c call WriteField_phy('physiq_zphi',zphi,klev)
	805	c call WriteField_phy('physiq_zzlay',zzlay,klev)
	806	c call WriteField_phy('physiq_zzlev',zzlev,klev+1)
[106]	807	c- - - - - - - - - - - - - - - -
	808	c DIAGNOSTIQUE GRILLE VERTICALE
	809	c- - - - - - - - - - - - - - - -
	810	c print*,"DIAGNOSTIQUE GRILLE VERTICALE"
	811	c i=klon/2
	812	c print*,"Niveau Pression Altitude (lev puis lay)"
	813	c do l=1,klev
	814	c print*,l,paprs(i,l),zzlev(i,l)
	815	c print*,l,pplay(i,l),zzlay(i,l)
	816	c enddo
	817	c print*,klev+1,paprs(i,klev+1),zzlev(i,klev+1)
	818	c stop
	819
[3]	820	c====================================================================
	821	c
	822	c Verifier les temperatures
	823	c
	824	CALL hgardfou(t_seri,ftsol,'debutphy')
	825	c====================================================================
	826	c
	827	c Incrementer le compteur de la physique
	828	c
	829	itap = itap + 1
	830
	831	c====================================================================
	832	c
	833	c Epaisseurs couches
	834
	835	DO k = 1, klev
	836	DO i = 1, klon
	837	delp(i,k) = paprs(i,k)-paprs(i,k+1)
	838	ENDDO
	839	ENDDO
	840
	841	c====================================================================
[1056]	842	c Orbite et eclairement
[3]	843	c====================================================================
	844
	845	c Pour TITAN:
	846	c calcul de la longitude solaire
[123]	847	CALL solarlong(rjourvrai+gmtime,zls)
[815]	848	zlsdeg = zls*180./RPI ! zls est en radians !!
	849	print*,'Ls',zlsdeg
	850
[3]	851	CALL orbite(zls,dist,pdecli)
[854]	852	IF (debut) zlsm1=zls
[3]	853
	854	c dans zenang, Ls en degres ; dans mucorr, Ls en radians
[1056]	855	call mucorr(klon,zls,rlatd,rmu0bar,fractbar)
[3]	856	IF (cycle_diurne) THEN
[1056]	857	zdtime=dtime*REAL(radpas) ! pas de temps du rayonnement (s)
[815]	858	CALL zenang(zlsdeg,gmtime,zdtime,rlatd,rlond,rmu0,fract)
[3]	859	ELSE
[1056]	860	rmu0 = rmu0bar
	861	fract = fractbar
[3]	862	ENDIF
[1056]	863
[3]	864	c====================================================================
[1056]	865	c Appeler la diffusion verticale (programme de couche limite)
[3]	866	c====================================================================
	867
	868	c-------------------------------
	869	c TEST: on ne tient pas compte des calculs de clmain mais on force
	870	c l'equilibre radiatif du sol
	871	if (1.eq.0) then
	872	if (debut) then
	873	print*,"ATTENTION, CLMAIN SHUNTEE..."
	874	endif
	875
	876	DO i = 1, klon
	877	sens(i) = 0.0e0 ! flux de chaleur sensible au sol
	878	fder(i) = 0.0e0
	879	dlw(i) = 0.0e0
	880	ENDDO
	881
	882	c Incrementer la temperature du sol
	883	c
	884	DO i = 1, klon
	885	d_ts(i) = dtime * radsol(i)/22000. !valeur calculee par GCM pour I=200
	886	ftsol(i) = ftsol(i) + d_ts(i)
	887	do j=1,nsoilmx
	888	ftsoil(i,j)=ftsol(i)
	889	enddo
	890	ENDDO
	891
	892	c-------------------------------
	893	else
	894	c-------------------------------
	895
	896	fder = dlw
	897
	898	c print*,"radsol avant clmain=",radsol(klon/2)
	899	c print*,"solsw avant clmain=",solsw(klon/2)
	900	c print*,"sollw avant clmain=",sollw(klon/2)
	901
[1056]	902	! ADAPTATION GCM POUR CP(T)
[3]	903
	904	CALL clmain(dtime,itap,
	905	e t_seri,u_seri,v_seri,
	906	e rmu0,
	907	e ftsol,
	908	$ ftsoil,
	909	$ paprs,pplay,ppk,radsol,falbe,
	910	e solsw, sollw, sollwdown, fder,
[1056]	911	e rlond, rlatd, cuphy, cvphy,
[3]	912	e debut, lafin,
	913	s d_t_vdf,d_u_vdf,d_v_vdf,d_ts,
	914	s fluxt,fluxu,fluxv,cdragh,cdragm,
	915	s dsens,
	916	s ycoefh,yu1,yv1)
	917
	918	c print*,"radsol apres clmain=",radsol(klon/2)
	919	c print*,"solsw apres clmain=",solsw(klon/2)
	920	c print*,"sollw apres clmain=",sollw(klon/2)
	921
	922	CXXX Incrementation des flux
	923	DO i = 1, klon
	924	sens(i) = - fluxt(i,1) ! flux de chaleur sensible au sol
	925	fder(i) = dlw(i) + dsens(i)
	926	ENDDO
	927	CXXX
	928
	929	DO k = 1, klev
	930	DO i = 1, klon
	931	t_seri(i,k) = t_seri(i,k) + d_t_vdf(i,k)
	932	d_t_vdf(i,k)= d_t_vdf(i,k)/dtime ! K/s
	933	u_seri(i,k) = u_seri(i,k) + d_u_vdf(i,k)
	934	d_u_vdf(i,k)= d_u_vdf(i,k)/dtime ! (m/s)/s
	935	v_seri(i,k) = v_seri(i,k) + d_v_vdf(i,k)
	936	d_v_vdf(i,k)= d_v_vdf(i,k)/dtime ! (m/s)/s
	937	ENDDO
	938	ENDDO
[123]	939
[1056]	940	c call WriteField_phy('physiq_dtvdf',d_t_vdf,klev)
	941	c call WriteField_phy('physiq_duvdf',d_u_vdf,klev)
	942	c call WriteField_phy('physiq_dvvdf',d_v_vdf,klev)
[3]	943
	944	C TRACEURS
	945
	946	d_tr_vdf = 0.
	947	if (iflag_trac.eq.1) then
[1056]	948	DO iq=1, nqmax
	949	CALL cltrac(dtime,ycoefh,t_seri,
	950	s tr_seri(1,1,iq),source,
	951	e paprs, pplay,delp,
[3]	952	s d_tr_vdf(1,1,iq))
[1056]	953	tr_seri(:,:,iq) = tr_seri(:,:,iq) + d_tr_vdf(:,:,iq)
	954	d_tr_vdf(:,:,iq)= d_tr_vdf(:,:,iq)/dtime ! /s
	955	ENDDO
[3]	956	endif
[123]	957
[3]	958	IF (if_ebil.ge.2) THEN
	959	ztit='after clmain'
	960	CALL diagetpq(airephy,ztit,ip_ebil,2,1,dtime
	961	e , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
	962	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
	963	call diagphy(airephy,ztit,ip_ebil
	964	e , zero_v, zero_v, zero_v, zero_v, sens
	965	e , zero_v, zero_v, zero_v, ztsol
	966	e , d_h_vcol, d_qt, d_ec
	967	s , fs_bound, fq_bound )
	968	END IF
	969	C
	970	c
	971	c Incrementer la temperature du sol
	972	c
[1056]	973	c print*,'Tsol avant clmain:',ftsol(1)
[3]	974	DO i = 1, klon
	975	ftsol(i) = ftsol(i) + d_ts(i)
	976	ENDDO
	977	c print*,'DTsol apres clmain:',d_ts(klon/2)
[1056]	978	c print*,'Tsol apres clmain:',ftsol(1)
[3]	979
	980	c Calculer la derive du flux infrarouge
	981	c
	982	DO i = 1, klon
[495]	983	dlw(i) = - 4.0emisRSIGMAftsol(i)*3
[3]	984	ENDDO
	985
	986	c-------------------------------
	987	endif ! fin du TEST
	988
	989	c
	990	c Appeler l'ajustement sec
	991	c
	992	c===================================================================
	993	c Convection seche
	994	c===================================================================
	995	c
	996	d_t_ajs(:,:)=0.
	997	d_u_ajs(:,:)=0.
	998	d_v_ajs(:,:)=0.
	999	d_tr_ajs(:,:,:)=0.
	1000	c
	1001	IF(prt_level>9)WRITE(lunout,*)
	1002	. 'AVANT LA CONVECTION SECHE , iflag_ajs='
	1003	s ,iflag_ajs
	1004
	1005	if(iflag_ajs.eq.0) then
	1006	c Rien
	1007	c ====
	1008	IF(prt_level>9)WRITE(lunout,*)'pas de convection'
	1009
	1010	else if(iflag_ajs.eq.1) then
	1011
	1012	c Ajustement sec
	1013	c ==============
	1014	IF(prt_level>9)WRITE(lunout,*)'ajsec'
	1015
	1016	! ADAPTATION GCM POUR CP(T)
	1017	CALL ajsec(paprs, pplay, ppk, t_seri, u_seri, v_seri, nqmax,
	1018	. tr_seri, d_t_ajs, d_u_ajs, d_v_ajs, d_tr_ajs)
	1019
	1020	! ADAPTATION GCM POUR CP(T)
	1021	do i=1,klon
	1022	flux_ajs(i,1) = 0.0
	1023	do j=2,klev
	1024	flux_ajs(i,j) = flux_ajs(i,j-1)
	1025	. + cpdet(t_seri(i,j-1))/RG*d_t_ajs(i,j-1)/dtime
	1026	. *delp(i,j-1)
	1027	enddo
	1028	enddo
	1029
	1030	t_seri(:,:) = t_seri(:,:) + d_t_ajs(:,:)
	1031	d_t_ajs(:,:)= d_t_ajs(:,:)/dtime ! K/s
	1032	u_seri(:,:) = u_seri(:,:) + d_u_ajs(:,:)
	1033	d_u_ajs(:,:)= d_u_ajs(:,:)/dtime ! (m/s)/s
	1034	v_seri(:,:) = v_seri(:,:) + d_v_ajs(:,:)
	1035	d_v_ajs(:,:)= d_v_ajs(:,:)/dtime ! (m/s)/s
	1036	if (iflag_trac.eq.1) then
[1056]	1037	tr_seri(:,:,:) = tr_seri(:,:,:) + d_tr_ajs(:,:,:)
	1038	d_tr_ajs(:,:,:)= d_tr_ajs(:,:,:)/dtime ! /s
[3]	1039	endif
	1040
[1056]	1041	c call WriteField_phy('physiq_dtajs',d_t_ajs,klev)
	1042	c call WriteField_phy('physiq_duajs',d_u_ajs,klev)
	1043	c call WriteField_phy('physiq_dvajs',d_v_ajs,klev)
	1044
[3]	1045	endif
	1046	c
	1047	IF (if_ebil.ge.2) THEN
	1048	ztit='after dry_adjust'
	1049	CALL diagetpq(airephy,ztit,ip_ebil,2,2,dtime
	1050	e , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
	1051	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
	1052	call diagphy(airephy,ztit,ip_ebil
	1053	e , zero_v, zero_v, zero_v, zero_v, sens
	1054	e , zero_v, zero_v, zero_v, ztsol
	1055	e , d_h_vcol, d_qt, d_ec
	1056	s , fs_bound, fq_bound )
	1057	END IF
	1058
	1059	c====================================================================
	1060	c MICROPHYSIQUE ET CHIMIE
	1061	c====================================================================
	1062
	1063	d_tr_mph(:,:,:)=0.
	1064	d_tr_kim(:,:,:)=0.
	1065
	1066	c on recupere tr_seri inchange, d_tr_micro, d_tr_chim, tous les trois sur nqmax
[104]	1067	c on recupere aussi qaer pour le mettre dans les sorties
[3]	1068	c si microfi=1, sortie de qaer(1:nmicro)
	1069	c si nmicro != nqmax et si chimi, sortie de tr_seri(nmicro+1:nqmax)
	1070
	1071	c faire un test comme pour rayonnement, avec chimi en plus comme flag,
	1072	c pour voir si chimie appelee -> bouleen, qui passe dans phytrac.
	1073	c faut aussi le pas de temps chimique: dtimechim, a passer..
	1074
	1075	appel_chim = 0
	1076	IF (MOD(itapchim,chimpas).EQ.0) THEN
	1077	c print*,'CHIMIE ',
	1078	c $ ' (itapchim=',itapchim,'/chimpas=',chimpas,')'
	1079	appel_chim = 1
	1080	itapchim = 0
	1081	ENDIF
	1082	itapchim = itapchim + 1
	1083
	1084	if (iflag_trac.eq.1) then
[1056]	1085	c call WriteField_phy('physiq_qaer01',
	1086	c . qaer(:,:,1),klev)
	1087	c call WriteField_phy('physiq_qaer10',
	1088	c . qaer(:,:,10),klev)
	1089	c call WriteField_phy('physiq_tr_seri01',
	1090	c . tr_seri(:,:,1),klev)
	1091	c call WriteField_phy('physiq_tr_seri10',
	1092	c . tr_seri(:,:,10),klev)
	1093
[175]	1094	c call begintime(tt0)
[1056]	1095	c in phytrac call, mu0 and fract are only used in brume
	1096	c so we need to pass either rmu0 ou rmu0bar depending on
	1097	c moyzon_mu
	1098	if (moyzon_mu) then
[3]	1099	call phytrac (debut,lafin,
[1056]	1100	. nqmax,nmicro,dtime,appel_chim,dtimechim,
	1101	. paprs,pplay,delp,t,rmu0bar,fractbar,pdecli,zls,
	1102	. yu1,yv1,zzlev,zzlay,ftsol,
	1103	. tr_seri,qaer,d_tr_mph,d_tr_kim,
	1104	. fclat,reservoir)
	1105	else
	1106	call phytrac (debut,lafin,
	1107	. nqmax,nmicro,dtime,appel_chim,dtimechim,
	1108	. paprs,pplay,delp,t,rmu0,fract,pdecli,zls,
	1109	. yu1,yv1,zzlev,zzlay,ftsol,
	1110	. tr_seri,qaer,d_tr_mph,d_tr_kim,
	1111	. fclat,reservoir)
	1112	endif
[3]	1113
[175]	1114	c call endtime(tt0,tt1)
	1115	c ttphytra=ttphytra+tt1
	1116
	1117	c ----- ICI on ajuste radsol en tenant compte du flux de chaleur latente
	1118	c d'evaporation du reservoir.
	1119	c NOTE : c'est pas tres elegant mais ca permet d'eviter d'aller
	1120	c toucher a clmain.
	1121	if (clouds.eq.1) then
	1122	radsol(:) = radsol(:)+fclat(:) !test pas de flx de chaleur latente
	1123	endif
	1124
	1125	if (microfi.ge.1) then
[3]	1126	tr_seri(:,:,1:nmicro) = tr_seri(:,:,1:nmicro)
	1127	. + d_tr_mph(:,:,1:nmicro)*dtime
[1126]	1128	c call WriteField_phy('physiq_d_tr_mph01',
	1129	c . d_tr_mph(:,:,1),klev)
	1130	c call WriteField_phy('physiq_d_tr_mph10',
	1131	c . d_tr_mph(:,:,10),klev)
	1132	endif
[175]	1133	c PAS ELEGANT mais je n'ai pas trouve d'autres solutions :
	1134	c Il semblerait qu'il y ait un probleme lorsque les tendances de traceurs
	1135	c retourne des traceurs nuls et il y a parfois des valeurs negatives qui trainent.
	1136	c Pour ne diffuser le probleme, on force les valeurs negatives a ZERO.
	1137	DO iq=1,nmicro
	1138	DO i=1,klon
	1139	DO l=1,klev
	1140	if (tr_seri(i,l,iq).lt.0.) then
	1141	tr_seri(i,l,iq) = 0.
	1142	endif
	1143	ENDDO
	1144	ENDDO
	1145	ENDDO
	1146
	1147	c condensation:
	1148	c NE PAS OUBLIER LA CONDENSATION DES NUAGES !!!!
	1149	if ((clouds.eq.1.or.(chimi)).and.nqmax.gt.nmicro) then
	1150	tr_seri(:,:,nmicro+1:nqmax) = tr_seri(:,:,nmicro+1:nqmax)
	1151	. + d_tr_mph(:,:,nmicro+1:nqmax)*dtime
	1152	endif
[1126]	1153
	1154	c chimie:
[3]	1155	if ((chimi).and.(nqmax.gt.nmicro)) then
	1156	tr_seri(:,:,:) = tr_seri(:,:,:) + d_tr_kim(:,:,:)*dtime
	1157	endif
[123]	1158
[1126]	1159	endif !iflag_trac=1
[3]	1160
[175]	1161	c ch4=0.
	1162	c do l=1,llm
	1163	c ch4(1,l) = tr_seri(1,l,ich4)
	1164	c do j=2,jjm
	1165	c ig0=1+(j-2)*iim
	1166	c do i=1,iim
	1167	c ch4(j,l)= ch4(j,l) + tr_seri(ig0+i,l,ich4)/iim
	1168	c enddo
	1169	c enddo
	1170	c ch4(jjm+1,l) = tr_seri(klon,l,ich4)
	1171	c enddo
	1172	c do j=1,jjm+1
	1173	c write(501,*) j,ch4(j,1)
	1174	c enddo
	1175	c do l=1,llm
	1176	c write(502,'(I3,49(ES24.17,1X))') l,
	1177	c & (ch4(j,l),j=1,jjm+1)
	1178	c enddo
	1179	c write(501,*) ""
	1180	c write(502,*) ""
	1181
[3]	1182	c------------------
	1183	c test condensation
[97]	1184	c do i=1,nqmax
[3]	1185	c if(tname(i).eq."HCN") then
	1186	c print*,"HCN="
	1187	c do k=1,klev
	1188	c print,k,tr_seri(klon/2,k,i),d_tr_mph(klon/2,k,i)dtime
	1189	c v ,d_tr_kim(klon/2,k,i)*dtime
	1190	c enddo
	1191	c stop
	1192	c endif
	1193	c enddo
	1194	c------------------
	1195
	1196	c====================================================================
	1197	c RAYONNEMENT
	1198	c====================================================================
	1199
	1200	IF (MOD(itaprad,radpas).EQ.0) THEN
	1201	c print*,'RAYONNEMENT ',
	1202	c $ ' (itaprad=',itaprad,'/radpas=',radpas,')'
	1203
	1204	c ATTENTION, (klon/2) ne marche pas toujours............
	1205	c print*,"radsol avant radlwsw=",radsol(klon/2)
	1206	c print*,"solsw avant radlwsw=",solsw(klon/2)
	1207	c print*,"sollw avant radlwsw=",sollw(klon/2)
	1208	c print*,"avant radlwsw"
[175]	1209
	1210	c ----------------
	1211	c Calcul du rayon moyen des gouttes et des fractions volumique pour le TR
	1212	c ----------------
	1213	IF (clouds.eq.1) THEN
	1214	DO i=1,klon
	1215	DO j=1,klev
[1056]	1216	rmcbar(i,j)=rmcbar(i,j)/MAX(REAL(ncount(i,j)),1.)
	1217	xfbar(i,j,:)=xfbar(i,j,:)/MAX(REAL(ncount(i,j)),1.)
[175]	1218	ENDDO
	1219	ENDDO
	1220	ENDIF
	1221
	1222	c call begintime(tt0)
[3]	1223	CALL radlwsw
[1056]	1224	e (dist, rmu0, fract, zzlev,
[3]	1225	e paprs, pplay,ftsol, t_seri, nqmax, nmicro,
[1056]	1226	c tr_seri, qaer)
	1227	c print*,"apres radlwsw"
	1228
[175]	1229	c call endtime(tt0,tt1)
	1230	c ttrad=ttrad+tt1
	1231
[3]	1232	c print*,"apres radlwsw"
[175]	1233	c mise a zero du rayon moyen des gouttes et des fractions volumique
	1234	IF (clouds.eq.1) THEN
	1235	rmcbar(:,:) = 0.
	1236	xfbar(:,:,:) = 0.
	1237	ncount(:,:) = 0
	1238	ENDIF
[3]	1239
	1240	itaprad = 0
	1241	DO k = 1, klev
[175]	1242	DO i = 1, klon
[3]	1243	dtrad(i,k) = heat(i,k)-cool(i,k) !K/s
[175]	1244	ENDDO
[3]	1245	ENDDO
[1056]	1246
	1247	c call WriteField_phy('physiq_heat',heat,klev)
	1248	c call WriteField_phy('physiq_cool',cool,klev)
	1249
[3]	1250	ENDIF
	1251	itaprad = itaprad + 1
	1252	c====================================================================
	1253	c
	1254	c Ajouter la tendance des rayonnements (tous les pas)
	1255	c
	1256	DO k = 1, klev
	1257	DO i = 1, klon
	1258	t_seri(i,k) = t_seri(i,k) + dtrad(i,k) * dtime
	1259	ENDDO
	1260	ENDDO
	1261
	1262	IF (if_ebil.ge.2) THEN
	1263	ztit='after rad'
	1264	CALL diagetpq(airephy,ztit,ip_ebil,2,2,dtime
	1265	e , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
	1266	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
	1267	call diagphy(airephy,ztit,ip_ebil
	1268	e , topsw, toplw, solsw, sollw, zero_v
	1269	e , zero_v, zero_v, zero_v, ztsol
	1270	e , d_h_vcol, d_qt, d_ec
	1271	s , fs_bound, fq_bound )
	1272	END IF
	1273	c
	1274
	1275	c====================================================================
	1276	c Calcul des gravity waves FLOTT
	1277	c====================================================================
	1278	c
[1056]	1279	if (ok_orodr.or.ok_gw_nonoro) then
	1280	c CALCUL DE N2
	1281	do i=1,klon
	1282	do k=2,klev
	1283	ztlev(i,k) = (t_seri(i,k)+t_seri(i,k-1))/2.
	1284	zpklev(i,k) = sqrt(ppk(i,k)*ppk(i,k-1))
	1285	enddo
	1286	enddo
	1287	call t2tpot(klon*klev,ztlev, ztetalev,zpklev)
	1288	call t2tpot(klon*klev,t_seri,ztetalay,ppk)
	1289	do i=1,klon
	1290	do k=2,klev
	1291	zdtetalev(i,k) = ztetalay(i,k)-ztetalay(i,k-1)
	1292	zdzlev(i,k) = (zphi(i,k)-zphi(i,k-1))/RG
	1293	zn2(i,k) = RGzdtetalev(i,k)/(ztetalev(i,k)zdzlev(i,k))
	1294	zn2(i,k) = max(zn2(i,k),1.e-12) ! securite
	1295	enddo
	1296	zn2(i,1) = 1.e-12 ! securite
	1297	enddo
	1298
	1299	endif
	1300
	1301	c ----------------------------ORODRAG
	1302	IF (ok_orodr) THEN
[3]	1303	c
[1056]	1304	c selection des points pour lesquels le shema est actif:
	1305	igwd=0
	1306	DO i=1,klon
	1307	itest(i)=0
	1308	c IF ((zstd(i).gt.10.0)) THEN
	1309	IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
	1310	itest(i)=1
	1311	igwd=igwd+1
	1312	idx(igwd)=i
	1313	ENDIF
	1314	ENDDO
	1315	c igwdim=MAX(1,igwd)
	1316	c
	1317	c A ADAPTER POUR VENUS!!!
	1318	CALL drag_noro(klon,klev,dtime,paprs,pplay,zphi,zn2,
	1319	e zmea,zstd, zsig, zgam, zthe,zpic,zval,
	1320	e igwd,idx,itest,
	1321	e t_seri, u_seri, v_seri,
	1322	s zulow, zvlow, zustrdr, zvstrdr,
	1323	s d_t_oro, d_u_oro, d_v_oro)
	1324
	1325	c print*,"d_u_oro=",d_u_oro(klon/2,:)
	1326	c ajout des tendances
	1327	t_seri(:,:) = t_seri(:,:) + d_t_oro(:,:)
	1328	d_t_oro(:,:)= d_t_oro(:,:)/dtime ! K/s
	1329	u_seri(:,:) = u_seri(:,:) + d_u_oro(:,:)
	1330	d_u_oro(:,:)= d_u_oro(:,:)/dtime ! (m/s)/s
	1331	v_seri(:,:) = v_seri(:,:) + d_v_oro(:,:)
	1332	d_v_oro(:,:)= d_v_oro(:,:)/dtime ! (m/s)/s
	1333	c
	1334	ELSE
	1335	d_t_oro = 0.
	1336	d_u_oro = 0.
	1337	d_v_oro = 0.
	1338	zustrdr = 0.
	1339	zvstrdr = 0.
	1340	c
	1341	ENDIF ! fin de test sur ok_orodr
	1342	c
	1343	c ----------------------------OROLIFT
	1344	IF (ok_orolf) THEN
	1345	print*,"ok_orolf NOT IMPLEMENTED !"
	1346	stop
	1347	c
	1348	c selection des points pour lesquels le shema est actif:
	1349	igwd=0
	1350	DO i=1,klon
	1351	itest(i)=0
	1352	IF ((zpic(i)-zmea(i)).GT.100.) THEN
	1353	itest(i)=1
	1354	igwd=igwd+1
	1355	idx(igwd)=i
	1356	ENDIF
	1357	ENDDO
	1358	c igwdim=MAX(1,igwd)
	1359	c
	1360	c A ADAPTER POUR VENUS ET TITAN!!!
	1361	c CALL lift_noro(klon,klev,dtime,paprs,pplay,
	1362	c e rlatd,zmea,zstd,zpic,zgam,zthe,zpic,zval,
[3]	1363	c e igwd,idx,itest,
	1364	c e t_seri, u_seri, v_seri,
[1056]	1365	c s zulow, zvlow, zustrli, zvstrli,
	1366	c s d_t_lif, d_u_lif, d_v_lif )
	1367
[3]	1368	c
[1056]	1369	c ajout des tendances
	1370	t_seri(:,:) = t_seri(:,:) + d_t_lif(:,:)
	1371	d_t_lif(:,:)= d_t_lif(:,:)/dtime ! K/s
	1372	u_seri(:,:) = u_seri(:,:) + d_u_lif(:,:)
	1373	d_u_lif(:,:)= d_u_lif(:,:)/dtime ! (m/s)/s
	1374	v_seri(:,:) = v_seri(:,:) + d_v_lif(:,:)
	1375	d_v_lif(:,:)= d_v_lif(:,:)/dtime ! (m/s)/s
[3]	1376	c
[1056]	1377	ELSE
	1378	d_t_lif = 0.
	1379	d_u_lif = 0.
	1380	d_v_lif = 0.
	1381	zustrli = 0.
	1382	zvstrli = 0.
[3]	1383	c
[1056]	1384	ENDIF ! fin de test sur ok_orolf
	1385
	1386	c ---------------------------- NON-ORO GRAVITY WAVES
	1387	IF(ok_gw_nonoro) then
	1388
	1389	abort_message="Option non developpee pour Titan"
	1390	call abort_gcm(modname,abort_message,1)
	1391	c A FAIRE POUR TITAN
[3]	1392	c call flott_gwd_ran(klon,klev,dtime,pplay,zn2,
	1393	c e t_seri, u_seri, v_seri,
	1394	c o zustrhi,zvstrhi,
	1395	c o d_t_hin, d_u_hin, d_v_hin)
[1056]	1396
	1397	c ajout des tendances
	1398
[3]	1399	c t_seri(:,:) = t_seri(:,:) + d_t_hin(:,:)
	1400	c d_t_hin(:,:)= d_t_hin(:,:)/dtime ! K/s
	1401	c u_seri(:,:) = u_seri(:,:) + d_u_hin(:,:)
	1402	c d_u_hin(:,:)= d_u_hin(:,:)/dtime ! (m/s)/s
	1403	c v_seri(:,:) = v_seri(:,:) + d_v_hin(:,:)
	1404	c d_v_hin(:,:)= d_v_hin(:,:)/dtime ! (m/s)/s
[1056]	1405
	1406	ELSE
	1407	d_t_hin = 0.
	1408	d_u_hin = 0.
	1409	d_v_hin = 0.
	1410	zustrhi = 0.
	1411	zvstrhi = 0.
	1412
	1413	ENDIF ! fin de test sur ok_gw_nonoro
	1414
[3]	1415	c====================================================================
	1416	c Transport de ballons
	1417	c====================================================================
	1418	if (ballons.eq.1) then
[1056]	1419	CALL ballon(30,pdtphys,rjourvrai,gmtime,rlatd,rlond,
[3]	1420	c C t,pplay,u,v,pphi) ! alt above surface (smoothed for GCM)
	1421	C t,pplay,u,v,zphi) ! alt above planet average radius
	1422	endif !ballons
	1423
	1424	c====================================================================
	1425	c Bilan de mmt angulaire
	1426	c====================================================================
	1427	if (bilansmc.eq.1) then
	1428	CMODDEB FLOTT
	1429	C CALCULER LE BILAN DE MOMENT ANGULAIRE (DIAGNOSTIQUE)
	1430	C STRESS NECESSAIRES: COUCHE LIMITE ET TOUTE LA PHYSIQUE
	1431
	1432	DO i = 1, klon
	1433	zustrph(i)=0.
	1434	zvstrph(i)=0.
	1435	zustrcl(i)=0.
	1436	zvstrcl(i)=0.
	1437	ENDDO
	1438	DO k = 1, klev
	1439	DO i = 1, klon
	1440	zustrph(i)=zustrph(i)+(u_seri(i,k)-u(i,k))/dtime*
	1441	c (paprs(i,k)-paprs(i,k+1))/rg
	1442	zvstrph(i)=zvstrph(i)+(v_seri(i,k)-v(i,k))/dtime*
	1443	c (paprs(i,k)-paprs(i,k+1))/rg
	1444	zustrcl(i)=zustrcl(i)+d_u_vdf(i,k)*
	1445	c (paprs(i,k)-paprs(i,k+1))/rg
	1446	zvstrcl(i)=zvstrcl(i)+d_v_vdf(i,k)*
	1447	c (paprs(i,k)-paprs(i,k+1))/rg
	1448	ENDDO
	1449	ENDDO
	1450
[1056]	1451	CALL aaam_bud (27,klon,klev,rjourvrai,gmtime,
[3]	1452	C ra,rg,romega,
	1453	C rlatd,rlond,pphis,
	1454	C zustrdr,zustrli,zustrcl,
	1455	C zvstrdr,zvstrli,zvstrcl,
	1456	C paprs,u,v)
	1457
	1458	CCMODFIN FLOTT
	1459	endif !bilansmc
	1460
[1056]	1461	c====================================================================
	1462	c====================================================================
	1463	c Calculer le transport de l'eau et de l'energie (diagnostique)
	1464	c
	1465	c A REVOIR POUR VENUS ET TITAN...
	1466	c
	1467	c CALL transp (paprs,ftsol,
	1468	c e t_seri, q_seri, u_seri, v_seri, zphi,
	1469	c s ve, vq, ue, uq)
	1470	c
	1471	c====================================================================
	1472	c+jld ec_conser
	1473	DO k = 1, klev
	1474	DO i = 1, klon
	1475	d_t_ec(i,k)=0.5/cpdet(t_seri(i,k))
	1476	$ (u(i,k)2+v(i,k)2-u_seri(i,k)2-v_seri(i,k)*2)
	1477	t_seri(i,k)=t_seri(i,k)+d_t_ec(i,k)
	1478	d_t_ec(i,k) = d_t_ec(i,k)/dtime
	1479	END DO
	1480	END DO
	1481	do i=1,klon
	1482	flux_ec(i,1) = 0.0
	1483	do j=2,klev
	1484	flux_ec(i,j) = flux_ec(i,j-1)
	1485	. +cpdet(t_seri(i,j-1))/RGd_t_ec(i,j-1)delp(i,j-1)
	1486	enddo
	1487	enddo
	1488
	1489	c-jld ec_conser
	1490	c====================================================================
	1491	IF (if_ebil.ge.1) THEN
	1492	ztit='after physic'
	1493	CALL diagetpq(airephy,ztit,ip_ebil,1,1,dtime
	1494	e , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
	1495	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
	1496	C Comme les tendances de la physique sont ajoute dans la dynamique,
	1497	C on devrait avoir que la variation d'entalpie par la dynamique
	1498	C est egale a la variation de la physique au pas de temps precedent.
	1499	C Donc la somme de ces 2 variations devrait etre nulle.
	1500	call diagphy(airephy,ztit,ip_ebil
	1501	e , topsw, toplw, solsw, sollw, sens
	1502	e , zero_v, zero_v, zero_v, ztsol
	1503	e , d_h_vcol, d_qt, d_ec
	1504	s , fs_bound, fq_bound )
	1505	C
	1506	d_h_vcol_phy=d_h_vcol
	1507	C
	1508	END IF
	1509	C
[3]	1510	c=======================================================================
	1511	c SORTIES
	1512	c=======================================================================
	1513
	1514	c Convertir les incrementations en tendances
	1515	c
	1516	DO k = 1, klev
	1517	DO i = 1, klon
	1518	d_u(i,k) = ( u_seri(i,k) - u(i,k) ) / dtime
	1519	d_v(i,k) = ( v_seri(i,k) - v(i,k) ) / dtime
	1520	d_t(i,k) = ( t_seri(i,k) - t(i,k) ) / dtime
	1521	ENDDO
	1522	ENDDO
	1523	c
	1524	DO iq = 1, nqmax
	1525	DO k = 1, klev
	1526	DO i = 1, klon
	1527	d_qx(i,k,iq) = ( tr_seri(i,k,iq) - qx(i,k,iq) ) / dtime
	1528	ENDDO
	1529	ENDDO
	1530	ENDDO
	1531
	1532	c------------------------
	1533	c Calcul moment cinetique
	1534	c------------------------
[1056]	1535	c TEST...
[3]	1536	c mangtot = 0.0
	1537	c DO k = 1, klev
	1538	c DO i = 1, klon
	1539	c mang(i,k) = RAcos(rlatd(i)RPI/180.)
	1540	c . (u_seri(i,k)+RAcos(rlatd(i)RPI/180.)ROMEGA)
[1056]	1541	c . airephy(i)(paprs(i,k)-paprs(i,k+1))/RG
[3]	1542	c mangtot=mangtot+mang(i,k)
	1543	c ENDDO
	1544	c ENDDO
	1545	c print*,"Moment cinetique total = ",mangtot
	1546	c
	1547	c------------------------
	1548	c
	1549	c Sauvegarder les valeurs de t et u a la fin de la physique:
	1550	c
	1551	DO k = 1, klev
	1552	DO i = 1, klon
	1553	u_ancien(i,k) = u_seri(i,k)
	1554	t_ancien(i,k) = t_seri(i,k)
	1555	ENDDO
	1556	ENDDO
	1557	c
	1558	c=============================================================
	1559	c Ecriture des sorties
	1560	c=============================================================
	1561
	1562	#ifdef CPP_IOIPSL
	1563
[1056]	1564	#ifdef histday
	1565	#include "write_histday.h"
	1566	#endif
	1567
[3]	1568	#ifdef histmth
	1569	#include "write_histmth.h"
	1570	#endif
	1571
	1572	#ifdef histins
	1573	#include "write_histins.h"
	1574	#endif
[1056]	1575
[3]	1576	#endif
	1577
	1578	c====================================================================
	1579	c Si c'est la fin, il faut conserver l'etat de redemarrage
	1580	c====================================================================
	1581	c
	1582	IF (lafin) THEN
	1583	itau_phy = itau_phy + itap
[815]	1584	lsinit = zlsdeg
[1056]	1585	CALL phyredem ("restartphy.nc")
	1586
[3]	1587	c--------------FLOTT
	1588	CMODEB LOTT
	1589	C FERMETURE DU FICHIER FORMATTE CONTENANT LES COMPOSANTES
	1590	C DU BILAN DE MOMENT ANGULAIRE.
	1591	if (bilansmc.eq.1) then
	1592	write(,)'Fermeture de aaam_bud.out (FL Vous parle)'
	1593	close(27)
	1594	close(28)
	1595	endif !bilansmc
	1596	CMODFIN
	1597	c-------------
	1598	c--------------SLEBONNOIS
	1599	C FERMETURE DES FICHIERS FORMATTES CONTENANT LES POSITIONS ET VITESSES
	1600	C DES BALLONS
	1601	if (ballons.eq.1) then
	1602	write(,)'Fermeture des ballons*.out'
	1603	close(30)
	1604	close(31)
	1605	close(32)
	1606	close(33)
	1607	close(34)
	1608	endif !ballons
	1609	c-------------
	1610	ENDIF
	1611
	1612	RETURN
	1613	END
	1614
	1615
	1616
	1617	***********************************************************************
	1618	***********************************************************************
	1619	***********************************************************************
	1620	***********************************************************************
	1621	***********************************************************************
	1622	***********************************************************************
	1623	***********************************************************************
	1624	***********************************************************************
	1625
	1626	SUBROUTINE gr_fi_ecrit(nfield,nlon,iim,jjmp1,fi,ecrit)
	1627	IMPLICIT none
	1628	c
	1629	c Tranformer une variable de la grille physique a
	1630	c la grille d'ecriture
	1631	c
	1632	INTEGER nfield,nlon,iim,jjmp1, jjm
	1633	REAL fi(nlon,nfield), ecrit(iim*jjmp1,nfield)
	1634	c
	1635	INTEGER i, n, ig
	1636	c
	1637	jjm = jjmp1 - 1
	1638	DO n = 1, nfield
	1639	DO i=1,iim
	1640	ecrit(i,n) = fi(1,n)
	1641	ecrit(i+jjm*iim,n) = fi(nlon,n)
	1642	ENDDO
	1643	DO ig = 1, nlon - 2
	1644	ecrit(iim+ig,n) = fi(1+ig,n)
	1645	ENDDO
	1646	ENDDO
	1647	RETURN
	1648	END

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

Download in other formats:

Original Format