Context Navigation

physiq.F @ 160

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

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

Download in other formats:

Original Format