Context Navigation

gcm.F @ 5426

Last change on this file since 5426 was 1520, checked in by Ehouarn Millour, 14 years ago

Implementation of a different vertical discretization (from/for planets, but
can in principle also be used for Earth).
Choice of vertical discretization is set by flag 'disvert_type';
'disvert_type=1' is Earth standard (default; ie set to 1 if
planet_type=="earth") case.
With 'disvert_type=2', approximate altitude of layers and reference atmospheric
scale height must be given using an input file ("z2sig.def", first line
should give scale height, in km, following lines must specify the altitude,
in km above surface, of mid-layers, one per line; see disvert_noterre.F).

Checked that these changes do not impact on 'bench' results, on Vargas.

EM.

Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 16.7 KB

Line
1	!
2	! $Id: gcm.F 1520 2011-05-23 11:37:09Z fhourdin $
3	!
4	c
5	c
6	PROGRAM gcm
7
8	#ifdef CPP_IOIPSL
9	USE IOIPSL
10	#endif
11
12	USE mod_const_mpi, ONLY: init_const_mpi
13	USE parallel
14	USE infotrac
15	USE mod_interface_dyn_phys
16	USE mod_hallo
17	USE Bands
18	USE getparam
19	USE filtreg_mod
20	USE control_mod
21
22	! Ehouarn: for now these only apply to Earth:
23	#ifdef CPP_EARTH
24	USE mod_grid_phy_lmdz
25	USE mod_phys_lmdz_para, ONLY : klon_mpi_para_nb
26	USE mod_phys_lmdz_omp_data, ONLY: klon_omp
27	USE dimphy
28	USE comgeomphy
29	#endif
30	IMPLICIT NONE
31
32	c ...... Version du 10/01/98 ..........
33
34	c avec coordonnees verticales hybrides
35	c avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 )
36
37	c=======================================================================
38	c
39	c Auteur: P. Le Van /L. Fairhead/F.Hourdin
40	c -------
41	c
42	c Objet:
43	c ------
44	c
45	c GCM LMD nouvelle grille
46	c
47	c=======================================================================
48	c
49	c ... Dans inigeom , nouveaux calculs pour les elongations cu , cv
50	c et possibilite d'appeler une fonction f(y) a derivee tangente
51	c hyperbolique a la place de la fonction a derivee sinusoidale.
52	c ... Possibilite de choisir le schema pour l'advection de
53	c q , en modifiant iadv dans traceur.def (MAF,10/02) .
54	c
55	c Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron,10/99)
56	c Pour Van-Leer iadv=10
57	c
58	c-----------------------------------------------------------------------
59	c Declarations:
60	c -------------
61	#include "dimensions.h"
62	#include "paramet.h"
63	#include "comconst.h"
64	#include "comdissnew.h"
65	#include "comvert.h"
66	#include "comgeom.h"
67	#include "logic.h"
68	#include "temps.h"
69	#include "ener.h"
70	#include "description.h"
71	#include "serre.h"
72	!#include "com_io_dyn.h"
73	#include "iniprint.h"
74	#include "tracstoke.h"
75
76	#ifdef INCA
77	! Only INCA needs these informations (from the Earth's physics)
78	#include "indicesol.h"
79	#endif
80
81	INTEGER longcles
82	PARAMETER ( longcles = 20 )
83	REAL clesphy0( longcles )
84	SAVE clesphy0
85
86
87
88	REAL zdtvr
89	c INTEGER nbetatmoy, nbetatdem,nbetat
90	INTEGER nbetatmoy, nbetatdem
91
92	c variables dynamiques
93	REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
94	REAL teta(ip1jmp1,llm) ! temperature potentielle
95	REAL, ALLOCATABLE, DIMENSION(:,:,:) :: q ! champs advectes
96	REAL ps(ip1jmp1) ! pression au sol
97	c REAL p (ip1jmp1,llmp1 ) ! pression aux interfac.des couches
98	c REAL pks(ip1jmp1) ! exner au sol
99	c REAL pk(ip1jmp1,llm) ! exner au milieu des couches
100	c REAL pkf(ip1jmp1,llm) ! exner filt.au milieu des couches
101	REAL masse(ip1jmp1,llm) ! masse d'air
102	REAL phis(ip1jmp1) ! geopotentiel au sol
103	c REAL phi(ip1jmp1,llm) ! geopotentiel
104	c REAL w(ip1jmp1,llm) ! vitesse verticale
105
106	c variables dynamiques intermediaire pour le transport
107
108	c variables pour le fichier histoire
109	REAL dtav ! intervalle de temps elementaire
110
111	REAL time_0
112
113	LOGICAL lafin
114	c INTEGER ij,iq,l,i,j
115	INTEGER i,j
116
117
118	real time_step, t_wrt, t_ops
119
120
121	LOGICAL call_iniphys
122	data call_iniphys/.true./
123
124	c REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm)
125	c+jld variables test conservation energie
126	c REAL ecin(ip1jmp1,llm),ecin0(ip1jmp1,llm)
127	C Tendance de la temp. potentiel d (theta)/ d t due a la
128	C tansformation d'energie cinetique en energie thermique
129	C cree par la dissipation
130	c REAL dhecdt(ip1jmp1,llm)
131	c REAL vcont(ip1jm,llm),ucont(ip1jmp1,llm)
132	c REAL d_h_vcol, d_qt, d_qw, d_ql, d_ec
133	c CHARACTER (len=15) :: ztit
134	c-jld
135
136
137	character (len=80) :: dynhist_file, dynhistave_file
138	character (len=20) :: modname
139	character (len=80) :: abort_message
140	! locales pour gestion du temps
141	INTEGER :: an, mois, jour
142	REAL :: heure
143
144
145	c-----------------------------------------------------------------------
146	c variables pour l'initialisation de la physique :
147	c ------------------------------------------------
148	INTEGER ngridmx
149	PARAMETER( ngridmx = 2+(jjm-1)*iim - 1/jjm )
150	REAL zcufi(ngridmx),zcvfi(ngridmx)
151	REAL latfi(ngridmx),lonfi(ngridmx)
152	REAL airefi(ngridmx)
153	SAVE latfi, lonfi, airefi
154
155	INTEGER :: ierr
156
157
158	c-----------------------------------------------------------------------
159	c Initialisations:
160	c ----------------
161
162	abort_message = 'last timestep reached'
163	modname = 'gcm'
164	descript = 'Run GCM LMDZ'
165	lafin = .FALSE.
166	dynhist_file = 'dyn_hist'
167	dynhistave_file = 'dyn_hist_ave'
168
169
170
171	c----------------------------------------------------------------------
172	c lecture des fichiers gcm.def ou run.def
173	c ---------------------------------------
174	c
175	! Ehouarn: dump possibility of using defrun
176	!#ifdef CPP_IOIPSL
177	CALL conf_gcm( 99, .TRUE. , clesphy0 )
178	!#else
179	! CALL defrun( 99, .TRUE. , clesphy0 )
180	!#endif
181	c
182	c
183	c------------------------------------
184	c Initialisation partie parallele
185	c------------------------------------
186	CALL init_const_mpi
187
188	call init_parallel
189	call ini_getparam("out.def")
190	call Read_Distrib
191	! Ehouarn : temporarily (?) keep this only for Earth
192	if (planet_type.eq."earth") then
193	#ifdef CPP_EARTH
194	CALL Init_Phys_lmdz(iim,jjp1,llm,mpi_size,distrib_phys)
195	#endif
196	endif ! of if (planet_type.eq."earth")
197	CALL set_bands
198	#ifdef CPP_EARTH
199	! Ehouarn: For now only Earth physics is parallel
200	CALL Init_interface_dyn_phys
201	#endif
202	CALL barrier
203
204	if (mpi_rank==0) call WriteBands
205	call SetDistrib(jj_Nb_Caldyn)
206
207	c$OMP PARALLEL
208	call Init_Mod_hallo
209	c$OMP END PARALLEL
210
211	! Ehouarn : temporarily (?) keep this only for Earth
212	if (planet_type.eq."earth") then
213	#ifdef CPP_EARTH
214	c$OMP PARALLEL
215	call InitComgeomphy
216	c$OMP END PARALLEL
217	#endif
218	endif ! of if (planet_type.eq."earth")
219
220	c-----------------------------------------------------------------------
221	c Choix du calendrier
222	c -------------------
223
224	c calend = 'earth_365d'
225
226	#ifdef CPP_IOIPSL
227	if (calend == 'earth_360d') then
228	call ioconf_calendar('360d')
229	write(lunout,*)'CALENDRIER CHOISI: Terrestre a 360 jours/an'
230	else if (calend == 'earth_365d') then
231	call ioconf_calendar('noleap')
232	write(lunout,*)'CALENDRIER CHOISI: Terrestre a 365 jours/an'
233	else if (calend == 'earth_366d') then
234	call ioconf_calendar('gregorian')
235	write(lunout,*)'CALENDRIER CHOISI: Terrestre bissextile'
236	else
237	abort_message = 'Mauvais choix de calendrier'
238	call abort_gcm(modname,abort_message,1)
239	endif
240	#endif
241
242	IF (config_inca /= 'none') THEN
243	#ifdef INCA
244	call init_const_lmdz(
245	$ nbtr,anneeref,dayref,
246	$ iphysiq,day_step,nday,
247	$ nbsrf, is_oce,is_sic,
248	$ is_ter,is_lic)
249
250	call init_inca_para(
251	$ iim,jjm+1,llm,klon_glo,mpi_size,
252	$ distrib_phys,COMM_LMDZ)
253	#endif
254	END IF
255
256	c-----------------------------------------------------------------------
257	c Initialisation des traceurs
258	c ---------------------------
259	c Choix du nombre de traceurs et du schema pour l'advection
260	c dans fichier traceur.def, par default ou via INCA
261	call infotrac_init
262
263	c Allocation de la tableau q : champs advectes
264	ALLOCATE(q(ip1jmp1,llm,nqtot))
265
266	c-----------------------------------------------------------------------
267	c Lecture de l'etat initial :
268	c ---------------------------
269
270	c lecture du fichier start.nc
271	if (read_start) then
272	! we still need to run iniacademic to initialize some
273	! constants & fields, if we run the 'newtonian' or 'SW' cases:
274	if (iflag_phys.ne.1) then
275	CALL iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
276	endif
277
278	! if (planet_type.eq."earth") then
279	! Load an Earth-format start file
280	CALL dynetat0("start.nc",vcov,ucov,
281	& teta,q,masse,ps,phis, time_0)
282	! endif ! of if (planet_type.eq."earth")
283
284	c write(73,*) 'ucov',ucov
285	c write(74,*) 'vcov',vcov
286	c write(75,*) 'teta',teta
287	c write(76,*) 'ps',ps
288	c write(77,*) 'q',q
289
290	endif ! of if (read_start)
291
292	c le cas echeant, creation d un etat initial
293	IF (prt_level > 9) WRITE(lunout,*)
294	. 'GCM: AVANT iniacademic AVANT AVANT AVANT AVANT'
295	if (.not.read_start) then
296	CALL iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
297	endif
298
299	c-----------------------------------------------------------------------
300	c Lecture des parametres de controle pour la simulation :
301	c -------------------------------------------------------
302	c on recalcule eventuellement le pas de temps
303
304	IF(MOD(day_step,iperiod).NE.0) THEN
305	abort_message =
306	. 'Il faut choisir un nb de pas par jour multiple de iperiod'
307	call abort_gcm(modname,abort_message,1)
308	ENDIF
309
310	IF(MOD(day_step,iphysiq).NE.0) THEN
311	abort_message =
312	* 'Il faut choisir un nb de pas par jour multiple de iphysiq'
313	call abort_gcm(modname,abort_message,1)
314	ENDIF
315
316	zdtvr = daysec/REAL(day_step)
317	IF(dtvr.NE.zdtvr) THEN
318	WRITE(lunout,*)
319	. 'WARNING!!! changement de pas de temps',dtvr,'>',zdtvr
320	ENDIF
321
322	C
323	C on remet le calendrier à zero si demande
324	c
325	IF (raz_date == 1) THEN
326	annee_ref = anneeref
327	day_ref = dayref
328	day_ini = dayref
329	itau_dyn = 0
330	itau_phy = 0
331	time_0 = 0.
332	write(lunout,*)
333	. 'GCM: On reinitialise a la date lue dans gcm.def'
334	ELSE IF (annee_ref .ne. anneeref .or. day_ref .ne. dayref) THEN
335	write(lunout,*)
336	. 'GCM: Attention les dates initiales lues dans le fichier'
337	write(lunout,*)
338	. ' restart ne correspondent pas a celles lues dans '
339	write(lunout,*)' gcm.def'
340	write(lunout,*)' annee_ref=',annee_ref," anneeref=",anneeref
341	write(lunout,*)' day_ref=',day_ref," dayref=",dayref
342	write(lunout,*)' Pas de remise a zero'
343	ENDIF
344	c if (annee_ref .ne. anneeref .or. day_ref .ne. dayref) then
345	c write(lunout,*)
346	c . 'GCM: Attention les dates initiales lues dans le fichier'
347	c write(lunout,*)
348	c . ' restart ne correspondent pas a celles lues dans '
349	c write(lunout,*)' gcm.def'
350	c write(lunout,*)' annee_ref=',annee_ref," anneeref=",anneeref
351	c write(lunout,*)' day_ref=',day_ref," dayref=",dayref
352	c if (raz_date .ne. 1) then
353	c write(lunout,*)
354	c . 'GCM: On garde les dates du fichier restart'
355	c else
356	c annee_ref = anneeref
357	c day_ref = dayref
358	c day_ini = dayref
359	c itau_dyn = 0
360	c itau_phy = 0
361	c time_0 = 0.
362	c write(lunout,*)
363	c . 'GCM: On reinitialise a la date lue dans gcm.def'
364	c endif
365	c ELSE
366	c raz_date = 0
367	c endif
368
369	#ifdef CPP_IOIPSL
370	mois = 1
371	heure = 0.
372	call ymds2ju(annee_ref, mois, day_ref, heure, jD_ref)
373	jH_ref = jD_ref - int(jD_ref)
374	jD_ref = int(jD_ref)
375
376	call ioconf_startdate(INT(jD_ref), jH_ref)
377
378	write(lunout,*)'DEBUG'
379	write(lunout,*)'annee_ref, mois, day_ref, heure, jD_ref'
380	write(lunout,*)annee_ref, mois, day_ref, heure, jD_ref
381	call ju2ymds(jD_ref+jH_ref,an, mois, jour, heure)
382	write(lunout,*)'jD_ref+jH_ref,an, mois, jour, heure'
383	write(lunout,*)jD_ref+jH_ref,an, mois, jour, heure
384	#else
385	! Ehouarn: we still need to define JD_ref and JH_ref
386	! and since we don't know how many days there are in a year
387	! we set JD_ref to 0 (this should be improved ...)
388	jD_ref=0
389	jH_ref=0
390	#endif
391
392	c nombre d'etats dans les fichiers demarrage et histoire
393	nbetatdem = nday / iecri
394	nbetatmoy = nday / periodav + 1
395
396	if (iflag_phys.eq.1) then
397	! these initialisations have already been done (via iniacademic)
398	! if running in SW or Newtonian mode
399	c-----------------------------------------------------------------------
400	c Initialisation des constantes dynamiques :
401	c ------------------------------------------
402	dtvr = zdtvr
403	CALL iniconst
404
405	c-----------------------------------------------------------------------
406	c Initialisation de la geometrie :
407	c --------------------------------
408	CALL inigeom
409
410	c-----------------------------------------------------------------------
411	c Initialisation du filtre :
412	c --------------------------
413	CALL inifilr
414	endif ! of if (iflag_phys.eq.1)
415	c
416	c-----------------------------------------------------------------------
417	c Initialisation de la dissipation :
418	c ----------------------------------
419
420	CALL inidissip( lstardis, nitergdiv, nitergrot, niterh ,
421	* tetagdiv, tetagrot , tetatemp )
422
423	c-----------------------------------------------------------------------
424	c Initialisation de la physique :
425	c -------------------------------
426	IF (call_iniphys.and.iflag_phys.eq.1) THEN
427	latfi(1)=rlatu(1)
428	lonfi(1)=0.
429	zcufi(1) = cu(1)
430	zcvfi(1) = cv(1)
431	DO j=2,jjm
432	DO i=1,iim
433	latfi((j-2)*iim+1+i)= rlatu(j)
434	lonfi((j-2)*iim+1+i)= rlonv(i)
435	zcufi((j-2)iim+1+i) = cu((j-1)iip1+i)
436	zcvfi((j-2)iim+1+i) = cv((j-1)iip1+i)
437	ENDDO
438	ENDDO
439	latfi(ngridmx)= rlatu(jjp1)
440	lonfi(ngridmx)= 0.
441	zcufi(ngridmx) = cu(ip1jm+1)
442	zcvfi(ngridmx) = cv(ip1jm-iim)
443	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,aire,airefi)
444
445	WRITE(lunout,*)
446	. 'GCM: WARNING!!! vitesse verticale nulle dans la physique'
447	! Earth:
448	if (planet_type.eq."earth") then
449	#ifdef CPP_EARTH
450	CALL iniphysiq(ngridmx,llm,daysec,day_ini,dtphys/nsplit_phys ,
451	, latfi,lonfi,airefi,zcufi,zcvfi,rad,g,r,cpp )
452	#endif
453	endif ! of if (planet_type.eq."earth")
454	call_iniphys=.false.
455	ENDIF ! of IF (call_iniphys.and.(iflag_phys.eq.1))
456
457
458	c-----------------------------------------------------------------------
459	c Initialisation des dimensions d'INCA :
460	c --------------------------------------
461	IF (config_inca /= 'none') THEN
462	!$OMP PARALLEL
463	#ifdef INCA
464	CALL init_inca_dim(klon_omp,llm,iim,jjm,
465	$ rlonu,rlatu,rlonv,rlatv)
466	#endif
467	!$OMP END PARALLEL
468	END IF
469
470	c-----------------------------------------------------------------------
471	c Initialisation des I/O :
472	c ------------------------
473
474
475	day_end = day_ini + nday
476	WRITE(lunout,300)day_ini,day_end
477	300 FORMAT('1'/,15x,'run du jour',i7,2x,'au jour',i7//)
478
479	#ifdef CPP_IOIPSL
480	call ju2ymds(jD_ref + day_ini - day_ref, an, mois, jour, heure)
481	write (lunout,301)jour, mois, an
482	call ju2ymds(jD_ref + day_end - day_ref, an, mois, jour, heure)
483	write (lunout,302)jour, mois, an
484	301 FORMAT('1'/,15x,'run du ', i2,'/',i2,'/',i4)
485	302 FORMAT('1'/,15x,' au ', i2,'/',i2,'/',i4)
486	#endif
487
488	! if (planet_type.eq."earth") then
489	! Write an Earth-format restart file
490	CALL dynredem0_p("restart.nc", day_end, phis)
491	! endif
492
493	ecripar = .TRUE.
494
495	#ifdef CPP_IOIPSL
496	time_step = zdtvr
497	IF (mpi_rank==0) then
498	if (ok_dyn_ins) then
499	! initialize output file for instantaneous outputs
500	! t_ops = iecri * daysec ! do operations every t_ops
501	t_ops =((1.0iecri)/day_step) daysec
502	t_wrt = daysec ! iecri * daysec ! write output every t_wrt
503	t_wrt = daysec ! iecri * daysec ! write output every t_wrt
504	CALL inithist(day_ref,annee_ref,time_step,
505	& t_ops,t_wrt)
506	endif
507
508	IF (ok_dyn_ave) THEN
509	! initialize output file for averaged outputs
510	t_ops = iperiod * time_step ! do operations every t_ops
511	t_wrt = periodav * daysec ! write output every t_wrt
512	CALL initdynav(day_ref,annee_ref,time_step,
513	& t_ops,t_wrt)
514	! CALL initdynav_p(dynhistave_file,day_ref,annee_ref,time_step,
515	! . t_ops, t_wrt, histaveid)
516	END IF
517	ENDIF
518	dtav = iperiod*dtvr/daysec
519	#endif
520	! #endif of #ifdef CPP_IOIPSL
521
522	c Choix des frequences de stokage pour le offline
523	c istdyn=day_step/4 ! stockage toutes les 6h=1jour/4
524	c istdyn=day_step/12 ! stockage toutes les 2h=1jour/12
525	istdyn=day_step/4 ! stockage toutes les 6h=1jour/12
526	istphy=istdyn/iphysiq
527
528
529	c
530	c-----------------------------------------------------------------------
531	c Integration temporelle du modele :
532	c ----------------------------------
533
534	c write(78,*) 'ucov',ucov
535	c write(78,*) 'vcov',vcov
536	c write(78,*) 'teta',teta
537	c write(78,*) 'ps',ps
538	c write(78,*) 'q',q
539
540	c$OMP PARALLEL DEFAULT(SHARED) COPYIN(/temps/,/logici/,/logicl/)
541	CALL leapfrog_p(ucov,vcov,teta,ps,masse,phis,q,clesphy0,
542	. time_0)
543	c$OMP END PARALLEL
544
545
546	END
547

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

Download in other formats:

Original Format