Context Navigation

physiq.F @ 814

Last change on this file since 814 was 808, checked in by slebonnois, 13 years ago
SL: Many changes for VENUS (related to newstart) and TITAN (related to clouds). Please read DOC/chantiers/commit_importants.log (cf v808).
File size: 51.7 KB

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

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

Download in other formats:

Original Format