Context Navigation

dynetat0.F90 @ 1980

Last change on this file since 1980 was 1850, checked in by jvatant, 7 years ago
Follow-up of r1849 : understood the black magic ! --JVO
File size: 13.1 KB

Line
1	!
2	! $Id $
3	!
4	SUBROUTINE dynetat0(fichnom,vcov,ucov,teta,q,masse,ps,phis,time0)
5
6	USE infotrac, only: tname, nqtot, zone_num, iso_indnum,&
7	iso_num, phase_num, alpha_ideal, iqiso, &
8	ok_isotopes, iqpere, tnat
9	use netcdf, only: nf90_open,NF90_NOWRITE,nf90_noerr,nf90_strerror, &
10	nf90_get_var, nf90_inq_varid, nf90_inq_dimid, &
11	nf90_inquire_dimension,nf90_close
12
13	use control_mod, only : planet_type, timestart
14	USE comvert_mod, ONLY: pa,preff
15	USE comconst_mod, ONLY: im,jm,lllm,daysec,dtvr, &
16	rad,omeg,g,cpp,kappa,pi
17	USE logic_mod, ONLY: fxyhypb,ysinus
18	USE serre_mod, ONLY: clon,clat,grossismx,grossismy
19	USE temps_mod, ONLY: annee_ref,day_ref,itau_dyn, &
20	start_time,day_ini,hour_ini
21	USE ener_mod, ONLY: etot0,ptot0,ztot0,stot0,ang0
22
23	IMPLICIT NONE
24
25	!=======================================================================
26	!
27	! Read initial confitions file
28	!
29	!=======================================================================
30
31	include "dimensions.h"
32	include "paramet.h"
33	include "comgeom2.h"
34	include "iniprint.h"
35
36	!===============================================================================
37	! Arguments:
38	CHARACTER(LEN=*), INTENT(IN) :: fichnom !--- FILE NAME
39	REAL, INTENT(OUT) :: vcov(iip1,jjm, llm) !--- V COVARIANT WIND
40	REAL, INTENT(OUT) :: ucov(iip1,jjp1,llm) !--- U COVARIANT WIND
41	REAL, INTENT(OUT) :: teta(iip1,jjp1,llm) !--- POTENTIAL TEMP.
42	REAL, INTENT(OUT) :: q(iip1,jjp1,llm,nqtot) !--- TRACERS
43	REAL, INTENT(OUT) :: masse(iip1,jjp1,llm) !--- MASS PER CELL
44	REAL, INTENT(OUT) :: ps(iip1,jjp1) !--- GROUND PRESSURE
45	REAL, INTENT(OUT) :: phis(iip1,jjp1) !--- GEOPOTENTIAL
46	REAL,INTENT(OUT) :: time0
47	!===============================================================================
48	! Local Variables
49	CHARACTER(LEN=256) :: msg, var, modname
50	INTEGER,PARAMETER :: length=100
51	INTEGER :: iq, fID, vID, idecal
52	REAL :: tab_cntrl(length) ! array containing run parameters
53	INTEGER :: ierr
54	CHARACTER(len=12) :: start_file_type="earth" ! default start file type
55
56	REAL,ALLOCATABLE :: time(:) ! times stored in start
57	INTEGER :: timelen ! number of times stored in the file
58	INTEGER :: indextime ! index of selected time
59	!REAL hour_ini ! fraction of day of stored date. Equivalent of day_ini, but 0=<hour_ini<1
60
61	INTEGER :: edges(4),corner(4)
62	INTEGER :: i
63
64	!-----------------------------------------------------------------------
65	modname="dynetat0"
66
67	! Open initial state NetCDF file
68	var=fichnom
69	CALL err(NF90_OPEN(var,NF90_NOWRITE,fID),"open",var)
70	!
71	CALL get_var1("controle",tab_cntrl)
72
73	!!! AS: idecal is a hack to be able to read planeto starts...
74	!!! .... while keeping everything OK for LMDZ EARTH
75	if ((planet_type.eq."generic").or.(planet_type.eq."mars")) then
76	write(lunout,*)'dynetat0 : Planeto-like start file'
77	start_file_type="planeto"
78	idecal = 4
79	annee_ref = 2000
80	else
81	if (planet_type.eq."titan") then
82	! Titan inherited Earth-like start files with idecal=5
83	write(lunout,*)'dynetat0 : Titan start file'
84	else
85	write(lunout,*)'dynetat0 : Earth-like start file'
86	endif
87	idecal = 5
88	annee_ref = tab_cntrl(5)
89	endif
90
91
92	im = tab_cntrl(1)
93	jm = tab_cntrl(2)
94	lllm = tab_cntrl(3)
95	if (start_file_type.eq."earth") then
96	day_ref = tab_cntrl(4)
97	else
98	day_ini = tab_cntrl(4)
99	day_ref=0
100	endif
101	rad = tab_cntrl(idecal+1)
102	omeg = tab_cntrl(idecal+2)
103	g = tab_cntrl(idecal+3)
104	cpp = tab_cntrl(idecal+4)
105	kappa = tab_cntrl(idecal+5)
106	daysec = tab_cntrl(idecal+6)
107	dtvr = tab_cntrl(idecal+7)
108	etot0 = tab_cntrl(idecal+8)
109	ptot0 = tab_cntrl(idecal+9)
110	ztot0 = tab_cntrl(idecal+10)
111	stot0 = tab_cntrl(idecal+11)
112	ang0 = tab_cntrl(idecal+12)
113	pa = tab_cntrl(idecal+13)
114	preff = tab_cntrl(idecal+14)
115	!
116	clon = tab_cntrl(idecal+15)
117	clat = tab_cntrl(idecal+16)
118	grossismx = tab_cntrl(idecal+17)
119	grossismy = tab_cntrl(idecal+18)
120	!
121	IF ( tab_cntrl(idecal+19).EQ.1. ) THEN
122	fxyhypb = .TRUE.
123	! dzoomx = tab_cntrl(25)
124	! dzoomy = tab_cntrl(26)
125	! taux = tab_cntrl(28)
126	! tauy = tab_cntrl(29)
127	ELSE
128	fxyhypb = .FALSE.
129	ysinus = .FALSE.
130	IF( tab_cntrl(idecal+22).EQ.1. ) ysinus = .TRUE.
131	ENDIF
132
133	if (planet_type=="mars") then ! so far this is only for Mars
134	hour_ini = tab_cntrl(29)
135	else
136	hour_ini=0
137	endif
138
139	if (start_file_type.eq."earth") then
140	day_ini = tab_cntrl(30)
141	itau_dyn = tab_cntrl(31)
142	start_time = tab_cntrl(32)
143	else
144	day_ini=tab_cntrl(4)
145	itau_dyn=0
146	start_time=0
147	endif
148	! .................................................................
149	!
150	!
151	WRITE(lunout,*)trim(modname)//': rad,omeg,g,cpp,kappa ', &
152	rad,omeg,g,cpp,kappa
153
154	CALL check_dim(im,iim,'im','iim')
155	CALL check_dim(jm,jjm,'jm','jjm')
156	CALL check_dim(llm,lllm,'llm','lllm')
157
158	CALL get_var1("rlonu",rlonu)
159	CALL get_var1("rlatu",rlatu)
160	CALL get_var1("rlonv",rlonv)
161	CALL get_var1("rlatv",rlatv)
162
163	CALL get_var2("cu" ,cu)
164	CALL get_var2("cv" ,cv)
165
166	CALL get_var2("aire" ,aire)
167	CALL get_var2("phisinit",phis)
168
169	! read time axis
170	ierr = nf90_inq_varid (fID, "temps", vID)
171	IF (ierr .NE. nf90_noerr) THEN
172	write(lunout,*)"dynetat0: Le champ <temps> est absent"
173	write(lunout,*)"dynetat0: J essaie <Time>"
174	ierr = nf90_inq_varid (fID, "Time", vID)
175	IF (ierr .NE. nf90_noerr) THEN
176	write(lunout,*)"dynetat0: Le champ <Time> est absent"
177	write(lunout,*)trim(nf90_strerror(ierr))
178	CALL ABORT_gcm("dynetat0", "", 1)
179	ENDIF
180	! Get the length of the "Time" dimension
181	ierr = nf90_inq_dimid(fID,"Time",vID)
182	ierr = nf90_inquire_dimension(fID,vID,len=timelen)
183	allocate(time(timelen))
184	! Then look for the "Time" variable
185	ierr =nf90_inq_varid(fID,"Time",vID)
186	ierr = nf90_get_var(fID, vID, time)
187	IF (ierr .NE. nf90_noerr) THEN
188	write(lunout,*)"dynetat0: Lecture echouee <Time>"
189	write(lunout,*)trim(nf90_strerror(ierr))
190	CALL ABORT_gcm("dynetat0", "", 1)
191	ENDIF
192	ELSE
193	! Get the length of the "temps" dimension
194	ierr = nf90_inq_dimid(fID,"temps",vID)
195	ierr = nf90_inquire_dimension(fID,vID,len=timelen)
196	allocate(time(timelen))
197	! Then look for the "temps" variable
198	ierr = nf90_inq_varid (fID, "temps", vID)
199	ierr = nf90_get_var(fID, vID, time)
200	IF (ierr .NE. nf90_noerr) THEN
201	write(lunout,*)"dynetat0: Lecture echouee <temps>"
202	write(lunout,*)trim(nf90_strerror(ierr))
203	CALL ABORT_gcm("dynetat0", "", 1)
204	ENDIF
205	ENDIF
206
207	! select the desired time
208	IF (timestart .lt. 0) THEN ! default: we use the last time value
209	indextime = timelen
210	ELSE ! else we look for the desired value in the time axis
211	indextime = 0
212	DO i=1,timelen
213	IF (abs(time(i) - timestart) .lt. 0.01) THEN
214	indextime = i
215	EXIT
216	ENDIF
217	ENDDO
218	IF (indextime .eq. 0) THEN
219	write(lunout,*)"Time", timestart," is not in " &
220	//trim(fichnom)//"!!"
221	write(lunout,*)"Stored times are:"
222	DO i=1,timelen
223	PRINT*, time(i)
224	ENDDO
225	CALL ABORT_gcm("dynetat0", "", 1)
226	ENDIF
227	ENDIF
228
229	if (planet_type=="mars") then
230	! In start the absolute date is day_ini + hour_ini + time
231	! For now on, in the GCM dynamics, it is day_ini + time0
232	time0 = time(indextime) + hour_ini
233	day_ini = day_ini + INT(time0)
234	time0 = time0 - INT(time0) ! time0 devient le nouveau hour_ini
235	hour_ini = time0
236	else
237	time0 = time(indextime)
238	endif
239
240	PRINT*, "dynetat0: Selected time ",time(indextime), &
241	" at index ",indextime
242
243	DEALLOCATE(time)
244
245	! read vcov
246	CALL get_var3v_t("vcov",vcov,indextime)
247
248	! read ucov
249	CALL get_var3u_t("ucov",ucov,indextime)
250
251	! read teta (same corner/edges as ucov)
252	CALL get_var3u_t("teta",teta,indextime)
253
254	! read tracers (same corner/edges as ucov)
255	corner(1)=1
256	corner(2)=1
257	corner(3)=1
258	corner(4)=indextime
259	edges(1)=iip1
260	edges(2)=jjp1
261	edges(3)=llm
262	edges(4)=1
263	IF(nqtot.GE.1) THEN
264	DO iq=1,nqtot
265	ierr= nf90_inq_varid(fID,tname(iq),vID)
266	IF (ierr .NE. nf90_noerr) THEN
267	write(lunout,*)TRIM(modname)//": Tracer <"//TRIM(var)//"> is missing"
268	write(lunout,*)" It is hence initialized to zero"
269	q(:,:,:,iq)=0.
270	IF (planet_type=="earth") THEN
271	!--- CRisi: for isotops, theoretical initialization using very simplified
272	! Rayleigh distillation las.
273	IF(ok_isotopes.AND.iso_num(iq)>0) THEN
274	IF(zone_num(iq)==0) q(:,:,:,iq)=q(:,:,:,iqpere(iq))*tnat(iso_num(iq)) &
275	& (q(:,:,:,iqpere(iq))/30.e-3)*(alpha_ideal(iso_num(iq))-1)
276	IF(zone_num(iq)==1) q(:,:,:,iq)=q(:,:,:,iqiso(iso_indnum(iq),phase_num(iq)))
277	END IF
278	ENDIF
279	ELSE
280	ierr=nf90_get_var(fID,vID,q(:,:,:,iq),corner,edges)
281	IF (ierr .NE. nf90_noerr) THEN
282	write(lunout,*)"dynetat0: Lecture echouee pour " &
283	//trim(tname(iq))
284	write(lunout,*)trim(nf90_strerror(ierr))
285	CALL ABORT_gcm("dynetat0", "", 1)
286	ENDIF
287	ENDIF
288	ENDDO
289	ENDIF
290
291	!read masse (same corner/edges as ucov)
292	CALL get_var3u_t("masse",masse,indextime)
293
294	!read ps
295	CALL get_var2_t("ps",ps,indextime)
296
297	CALL err(NF90_CLOSE(fID),"close",fichnom)
298
299	if (planet_type/="mars") then
300	day_ini=day_ini+INT(time0) ! obsolete stuff ; 0<time<1 anyways
301	time0=time0-INT(time0)
302	endif
303
304
305	CONTAINS
306
307	SUBROUTINE check_dim(n1,n2,str1,str2)
308	INTEGER, INTENT(IN) :: n1, n2
309	CHARACTER(LEN=*), INTENT(IN) :: str1, str2
310	CHARACTER(LEN=256) :: s1, s2
311	IF(n1/=n2) THEN
312	s1='value of '//TRIM(str1)//' ='
313	s2=' read in starting file differs from parametrized '//TRIM(str2)//' ='
314	WRITE(msg,'(10x,a,i4,2x,a,i4)')TRIM(s1),n1,TRIM(s2),n2
315	CALL ABORT_gcm(TRIM(modname),TRIM(msg),1)
316	END IF
317	END SUBROUTINE check_dim
318
319
320	SUBROUTINE get_var1(var,v)
321	CHARACTER(LEN=*), INTENT(IN) :: var
322	REAL, INTENT(OUT) :: v(:)
323	CALL err(NF90_INQ_VARID(fID,var,vID),"inq",var)
324	CALL err(NF90_GET_VAR(fID,vID,v),"get",var)
325	END SUBROUTINE get_var1
326
327
328	SUBROUTINE get_var2(var,v)
329	CHARACTER(LEN=*), INTENT(IN) :: var
330	REAL, INTENT(OUT) :: v(:,:)
331	CALL err(NF90_INQ_VARID(fID,var,vID),"inq",var)
332	CALL err(NF90_GET_VAR(fID,vID,v),"get",var)
333	END SUBROUTINE get_var2
334
335	SUBROUTINE get_var2_t(var,v,indextime)
336	CHARACTER(LEN=*), INTENT(IN) :: var
337	REAL, INTENT(OUT) :: v(:,:)
338	INTEGER, INTENT(IN) :: indextime
339	corner(1)=1
340	corner(2)=1
341	corner(3)=indextime
342	edges(1)=iip1
343	edges(2)=jjp1
344	edges(3)=1
345	CALL err(NF90_INQ_VARID(fID,var,vID),"inq",var)
346	CALL err(NF90_GET_VAR(fID,vID,v,corner,edges),"get",var)
347	END SUBROUTINE get_var2_t
348
349
350	SUBROUTINE get_var3(var,v) ! on U grid
351	CHARACTER(LEN=*), INTENT(IN) :: var
352	REAL, INTENT(OUT) :: v(:,:,:)
353	CALL err(NF90_INQ_VARID(fID,var,vID),"inq",var)
354	CALL err(NF90_GET_VAR(fID,vID,v),"get",var)
355	END SUBROUTINE get_var3
356
357	SUBROUTINE get_var3u_t(var,v,indextime) ! on U grid
358	CHARACTER(LEN=*), INTENT(IN) :: var
359	REAL, INTENT(OUT) :: v(:,:,:)
360	INTEGER, INTENT(IN) :: indextime
361	corner(1)=1
362	corner(2)=1
363	corner(3)=1
364	corner(4)=indextime
365	edges(1)=iip1
366	edges(2)=jjp1
367	edges(3)=llm
368	edges(4)=1
369	CALL err(NF90_INQ_VARID(fID,var,vID),"inq",var)
370	CALL err(NF90_GET_VAR(fID,vID,v,corner,edges),"get",var)
371	END SUBROUTINE get_var3u_t
372
373	SUBROUTINE get_var3v_t(var,v,indextime) ! on V grid
374	CHARACTER(LEN=*), INTENT(IN) :: var
375	REAL, INTENT(OUT) :: v(:,:,:)
376	INTEGER, INTENT(IN) :: indextime
377	corner(1)=1
378	corner(2)=1
379	corner(3)=1
380	corner(4)=indextime
381	edges(1)=iip1
382	edges(2)=jjm
383	edges(3)=llm
384	edges(4)=1
385	CALL err(NF90_INQ_VARID(fID,var,vID),"inq",var)
386	CALL err(NF90_GET_VAR(fID,vID,v,corner,edges),"get",var)
387	END SUBROUTINE get_var3v_t
388
389	SUBROUTINE err(ierr,typ,nam)
390	INTEGER, INTENT(IN) :: ierr !--- NetCDF ERROR CODE
391	CHARACTER(LEN=*), INTENT(IN) :: typ !--- TYPE OF OPERATION
392	CHARACTER(LEN=*), INTENT(IN) :: nam !--- FIELD/FILE NAME
393	IF(ierr==NF90_NoERR) RETURN
394	SELECT CASE(typ)
395	CASE('inq'); msg="Field <"//TRIM(nam)//"> is missing"
396	CASE('get'); msg="Reading failed for <"//TRIM(nam)//">"
397	CASE('open'); msg="File opening failed for <"//TRIM(nam)//">"
398	CASE('close'); msg="File closing failed for <"//TRIM(nam)//">"
399	END SELECT
400	CALL ABORT_gcm(TRIM(modname),TRIM(msg),ierr)
401	END SUBROUTINE err
402
403	END SUBROUTINE dynetat0

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

Download in other formats:

Original Format