Context Navigation

physiq_mod.F @ 1572

Last change on this file since 1572 was 1572, checked in by emillour, 8 years ago

All GCMs:
Further adaptations to keep up with changes in LMDZ5 concerning
physics/dynamics separation (up to rev r2500 of LMDZ5)

arch:
remove ifort debug option '-check all', replace it with '-check bounds,format,output_conversion,pointers,uninit' (i.e. get it to stop complaining about copying into temporary arrays)

dyn3d_common:
comconst_mod.F90 : add ngroup

dyn3d:
gcm.F90 : minor bug fix (arguments to a call_abort())
leapfrog.F90 : recompute geopotential for bilan_dyn outputs
conf_gcm.F90 : read "ngroup" from run.def
groupe.F , groupeun.F : ngroup no longer a local parameter

dyn3d_par:
conf_gcm.F90 : read "ngroup" from run.def
groupe_p.F , groupeun_p.F : ngroup no longer a local parameter

misc:
regr1_step_av_m.F90 : removed (not used)

phy_common:
mod_phys_lmdz_mpi_transfert.F90 , mod_phys_lmdz_mpi_data.F90 : change is_north_pole and is_south_pole to is_north_pole_dyn and is_south_pole_dyn
mod_phys_lmdz_omp_data.F90 : introduce is_nort_pole_phy and is_south_pole_phy

dynphy_lonlat:
mod_interface_dyn_phys.F90 : use is_north_pole_dyn and is_south_pole_dyn
calfis_p.F : use is_north_pole_dyn and is_south_pole_dyn

phyvenus:
physiq_mod , write_hist*.h : use is_north_pole_phy and is_south_pole_phy to correctly compute mesh area at poles to send to hist*nc files.

phytitan:
physiq_mod , write_hist*.h : use is_north_pole_phy and is_south_pole_phy to correctly compute mesh area at poles to send to hist*nc files.

File size: 50.1 KB

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

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

Download in other formats:

Original Format