Context Navigation

physiq.F @ 104

Last change on this file since 104 was 104, checked in by slebonnois, 14 years ago

SLebonnois: modification de makelmdz et create_make_gcm pour pouvoir
compiler la chimie titan. Pas de raison que ca gene les autres.
Dans cette version, les compilations de Venus et Titan fonctionnent.

Phytitan: modifications pour pouvoir compiler correctement.
Il ne manque plus que physiq.F a faire.

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

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

Download in other formats:

Original Format