Context Navigation

startvar.F @ 5326

Last change on this file since 5326 was 2, checked in by lmdz, 25 years ago
Initial revision
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 32.6 KB

Line
1	MODULE startvar
2	!
3	!
4	! There are three ways to access data from the database of atmospheric data which
5	! can be used to initialize the model. This depends on the type of field which needs
6	! to be extracted. In any case the call should come after a restget and should be of the type :
7	! CALL startget(...)
8	!
9	! We will details the possible arguments to startget here :
10	!
11	! - A 2D variable on the dynamical grid :
12	! CALL startget(varname, iml, jml, lon_in, lat_in, champ, val_ex)
13	!
14	!
15	! - A 1D variable on the physical grid :
16	! CALL startget(varname, iml, jml, lon_in, lat_in, nbindex, champ, val_exp)
17	!
18	!
19	! - A 3D variable on the dynamical grid :
20	! CALL startget(varname, iml, jml, lon_in, lat_in, lml, pls, workvar, champ, val_exp)
21	!
22	!
23	! There is special constraint on the atmospheric data base except that the
24	! the data needs to be in netCDF and the variables should have the the following
25	! names in the file :
26	!
27	! 'RELIEF' : High resolution orography
28	! 'ST' : Surface temperature
29	! 'CDSW' : Soil moisture
30	! 'Z' : Surface geopotential
31	! 'SP' : Surface pressure
32	! 'U' : East ward wind
33	! 'V' : Northward wind
34	! 'TEMP' : Temperature
35	! 'R' : Relative humidity
36	!
37	USE ioipsl
38	!
39	!
40	IMPLICIT NONE
41	!
42	!
43	PRIVATE
44	PUBLIC startget
45	!
46	!
47	INTERFACE startget
48	MODULE PROCEDURE startget_phys2d, startget_phys1d, startget_dyn
49	END INTERFACE
50	!
51	INTEGER, SAVE :: fid_phys, fid_dyn
52	INTEGER, SAVE :: iml_phys, iml_rel, iml_dyn
53	INTEGER, SAVE :: jml_phys, jml_rel, jml_dyn
54	INTEGER, SAVE :: llm_dyn, ttm_dyn
55	REAL, ALLOCATABLE, SAVE, DIMENSION (:,:) :: lon_phys, lon_rug,
56	. lon_alb, lon_rel, lon_dyn
57	REAL, ALLOCATABLE, SAVE, DIMENSION (:,:) :: lat_phys, lat_rug,
58	. lat_alb, lat_rel, lat_dyn
59	REAL, ALLOCATABLE, SAVE, DIMENSION (:) :: lev_dyn
60	REAL, ALLOCATABLE, SAVE, DIMENSION (:,:) :: relief, zstd, zsig,
61	. zgam, zthe, zpic, zval
62	REAL, ALLOCATABLE, SAVE, DIMENSION (:,:) :: rugo, masque, phis
63	!
64	REAL, ALLOCATABLE, SAVE, DIMENSION (:,:) :: tsol, qsol, psol_dyn
65	REAL, ALLOCATABLE, SAVE, DIMENSION (:,:,:) :: var_ana3d
66	!
67	CONTAINS
68	!
69	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
70	!
71	SUBROUTINE startget_phys2d(varname, iml, jml, lon_in, lat_in,
72	. champ, val_exp)
73	!
74	! There is a big mess with the size in logitude, should it be iml or iml+1.
75	! I have chosen to use the iml+1 as an argument to this routine and we declare
76	! internaly smaler fields when needed. This needs to be cleared once and for all in LMDZ.
77	! A convention is required.
78	!
79	!
80	CHARACTER(), INTENT(in) :: varname
81	INTEGER, INTENT(in) :: iml, jml
82	REAL, INTENT(in) :: lon_in(iml), lat_in(jml)
83	REAL, INTENT(inout) :: champ(iml,jml)
84	REAL, INTENT(in) :: val_exp
85	!
86	! This routine only works if the variable does not exist or is constant
87	!
88	IF ( MINVAL(champ(:,:)).EQ.MAXVAL(champ(:,:)) .AND.
89	.MINVAL(champ(:,:)).EQ.val_exp ) THEN
90	!
91	SELECTCASE(varname)
92	!
93	CASE ('relief')
94	!
95	! If we do not have the orography we need to get it
96	!
97	IF ( .NOT.ALLOCATED(relief)) THEN
98	!
99	CALL start_init_orog( iml, jml, lon_in, lat_in)
100	!
101	ENDIF
102	!
103	IF (SIZE(relief) .NE. SIZE(champ)) THEN
104	!
105	WRITE(,) 'STARTVAR module has been',
106	.' initialized to the wrong size'
107	STOP
108	!
109	ENDIF
110	!
111	champ(:,:) = relief(:,:)
112	!
113	CASE ('rugosite')
114	!
115	! If we do not have the orography we need to get it
116	!
117	IF ( .NOT.ALLOCATED(rugo)) THEN
118	!
119	CALL start_init_orog( iml, jml, lon_in, lat_in)
120	!
121	ENDIF
122	!
123	IF (SIZE(rugo) .NE. SIZE(champ)) THEN
124	!
125	WRITE(,)
126	. 'STARTVAR module has been initialized to the wrong size'
127	STOP
128	!
129	ENDIF
130	!
131	champ(:,:) = rugo(:,:)
132	!
133	CASE ('masque')
134	!
135	! If we do not have the orography we need to get it
136	!
137	IF ( .NOT.ALLOCATED(masque)) THEN
138	!
139	CALL start_init_orog( iml, jml, lon_in, lat_in)
140	!
141	ENDIF
142	!
143	IF (SIZE(masque) .NE. SIZE(champ)) THEN
144	!
145	WRITE(,)
146	. 'STARTVAR module has been initialized to the wrong size'
147	STOP
148	!
149	ENDIF
150	!
151	champ(:,:) = masque(:,:)
152	!
153	CASE ('surfgeo')
154	!
155	! If we do not have the orography we need to get it
156	!
157	IF ( .NOT.ALLOCATED(phis)) THEN
158	!
159	CALL start_init_orog( iml, jml, lon_in, lat_in)
160	!
161	ENDIF
162	!
163	IF (SIZE(phis) .NE. SIZE(champ)) THEN
164	!
165	WRITE(,)
166	. 'STARTVAR module has been initialized to the wrong size'
167	STOP
168	!
169	ENDIF
170	!
171	champ(:,:) = phis(:,:)
172	!
173	CASE ('psol')
174	!
175	! If we do not have the orography we need to get it
176	!
177	IF ( .NOT.ALLOCATED(psol_dyn)) THEN
178	!
179	CALL start_init_dyn( iml, jml, lon_in, lat_in)
180	!
181	ENDIF
182	!
183	IF (SIZE(psol_dyn) .NE. SIZE(champ)) THEN
184	!
185	WRITE(,)
186	. 'STARTVAR module has been initialized to the wrong size'
187	STOP
188	!
189	ENDIF
190	!
191	champ(:,:) = psol_dyn(:,:)
192	!
193	CASE DEFAULT
194	!
195	WRITE(,) 'startget_phys2d'
196	WRITE(,) 'No rule is present to extract variable',
197	. varname(:LEN_TRIM(varname)),' from any data set'
198	STOP
199	!
200	END SELECT
201	!
202	ELSE
203	!
204	! There are a few fields we might need if we need to interpolate 3D filed. Thus if they come through here we
205	! will catch them
206	!
207	SELECTCASE(varname)
208	!
209	CASE ('surfgeo')
210	!
211	IF ( .NOT.ALLOCATED(phis)) THEN
212	ALLOCATE(phis(iml,jml))
213	ENDIF
214	!
215	IF (SIZE(phis) .NE. SIZE(champ)) THEN
216	!
217	WRITE(,)
218	. 'STARTVAR module has been initialized to the wrong size'
219	STOP
220	!
221	ENDIF
222	!
223	phis(:,:) = champ(:,:)
224	!
225	END SELECT
226	!
227	ENDIF
228	!
229	END SUBROUTINE startget_phys2d
230	!
231	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
232	!
233	SUBROUTINE start_init_orog( iml, jml, lon_in, lat_in)
234	!
235	INTEGER, INTENT(in) :: iml, jml
236	REAL, INTENT(in) :: lon_in(iml), lat_in(jml)
237	!
238	! LOCAL
239	!
240	REAL :: lev(1), date, dt
241	INTEGER :: itau(1), fid
242	INTEGER :: llm_tmp, ttm_tmp
243	INTEGER :: i, j
244	INTEGER :: iret
245	REAL, ALLOCATABLE :: relief_hi(:,:)
246	REAL, ALLOCATABLE :: lon_rad(:), lat_rad(:)
247	REAL, ALLOCATABLE :: tmp_var(:,:)
248	INTEGER, ALLOCATABLE :: tmp_int(:,:)
249	!
250	CHARACTER*120 :: orogfname
251	LOGICAL :: check=.TRUE.
252	!
253	!
254	orogfname = 'Relief.nc'
255	!
256	IF ( check ) WRITE(,) 'Reading the high resolution orography'
257	!
258	CALL flininfo(orogfname,iml_rel, jml_rel, llm_tmp, ttm_tmp, fid)
259	!
260	ALLOCATE (lat_rel(iml_rel,jml_rel), stat=iret)
261	ALLOCATE (lon_rel(iml_rel,jml_rel), stat=iret)
262	ALLOCATE (relief_hi(iml_rel,jml_rel), stat=iret)
263	!
264	CALL flinopen(orogfname, .FALSE., iml_rel, jml_rel,
265	.llm_tmp, lon_rel, lat_rel, lev, ttm_tmp,
266	. itau, date, dt, fid)
267	!
268	CALL flinget(fid, 'RELIEF', iml_rel, jml_rel, llm_tmp,
269	. ttm_tmp, 1, 1, relief_hi)
270	!
271	CALL flinclo(fid)
272	!
273	! In case we have a file which is in degrees we do the transformation
274	!
275	ALLOCATE(lon_rad(iml_rel))
276	IF ( MAXVAL(lon_rel(:,:)) .GT. 2.0 * ASIN(1.0) ) THEN
277	lon_rad(:) = lon_rel(:,1) * 2.0 * ASIN(1.0) / 180.0
278	ELSE
279	lon_rad(:) = lon_rel(:,1)
280	ENDIF
281	ALLOCATE(lat_rad(jml_rel))
282	IF ( MAXVAL(lat_rel(:,:)) .GT. 2.0 * ASIN(1.0) ) THEN
283	lat_rad(:) = lat_rel(1,:) * 2.0 * ASIN(1.0) / 180.0
284	ELSE
285	lat_rad(:) = lat_rel(1,:)
286	ENDIF
287	!
288	!
289	IF ( check ) WRITE(,) 'Computes all the parameters needed',
290	.' for the gravity wave drag code'
291	!
292	! Allocate the data we need to put in the interpolated fields
293	!
294	! RELIEF: orographie moyenne
295	ALLOCATE(relief(iml,jml))
296	! zphi : orographie moyenne
297	ALLOCATE(phis(iml,jml))
298	! zstd: deviation standard de l'orographie sous-maille
299	ALLOCATE(zstd(iml,jml))
300	! zsig: pente de l'orographie sous-maille
301	ALLOCATE(zsig(iml,jml))
302	! zgam: anisotropy de l'orographie sous maille
303	ALLOCATE(zgam(iml,jml))
304	! zthe: orientation de l'axe oriente dans la direction
305	! de plus grande pente de l'orographie sous maille
306	ALLOCATE(zthe(iml,jml))
307	! zpic: hauteur pics de la SSO
308	ALLOCATE(zpic(iml,jml))
309	! zval: hauteur vallees de la SSO
310	ALLOCATE(zval(iml,jml))
311	! masque : Masque terre ocean
312	ALLOCATE(tmp_int(iml,jml))
313	ALLOCATE(masque(iml,jml))
314	!
315	CALL grid_noro(iml_rel, jml_rel, lon_rad, lat_rad, relief_hi,
316	. iml-1, jml, lon_in, lat_in,
317	. phis, relief, zstd, zsig, zgam, zthe, zpic, zval, tmp_int)
318	!
319	masque(:,:) = FLOAT(tmp_int(:,:))
320	!
321	! Compute surface roughness
322	!
323	IF ( check ) WRITE(,)
324	.'Compute surface roughness induced by the orography'
325	!
326	ALLOCATE(rugo(iml,jml))
327	ALLOCATE(tmp_var(iml-1,jml))
328	!
329	CALL rugsoro(iml_rel, jml_rel, lon_rad, lat_rad, relief_hi,
330	. iml-1, jml, lon_in, lat_in, tmp_var)
331	!
332	DO j = 1, jml
333	DO i = 1, iml-1
334	rugo(i,j) = tmp_var(i,j)
335	ENDDO
336	rugo(iml,j) = tmp_var(1,j)
337	ENDDO
338	!
339	! Build land-sea mask
340	!
341	!
342	RETURN
343	!
344	END SUBROUTINE start_init_orog
345	!
346	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
347	!
348	SUBROUTINE startget_phys1d(varname, iml, jml, lon_in,
349	.lat_in, nbindex, champ, val_exp)
350	!
351	CHARACTER(), INTENT(in) :: varname
352	INTEGER, INTENT(in) :: iml, jml, nbindex
353	REAL, INTENT(in) :: lon_in(iml), lat_in(jml)
354	REAL, INTENT(inout) :: champ(nbindex)
355	REAL, INTENT(in) :: val_exp
356	!
357	!
358	! This routine only works if the variable does not exist or is constant
359	!
360	IF ( MINVAL(champ(:)).EQ.MAXVAL(champ(:)) .AND.
361	.MINVAL(champ(:)).EQ.val_exp ) THEN
362	SELECTCASE(varname)
363	CASE ('tsol')
364	IF ( .NOT.ALLOCATED(tsol)) THEN
365	CALL start_init_phys( iml, jml, lon_in, lat_in)
366	ENDIF
367	IF ( SIZE(tsol) .NE. SIZE(lon_in)*SIZE(lat_in) ) THEN
368	WRITE(,)
369	. 'STARTVAR module has been initialized to the wrong size'
370	STOP
371	ENDIF
372	CALL gr_dyn_fi(1, iml, jml, nbindex, tsol, champ)
373	CASE ('qsol')
374	IF ( .NOT.ALLOCATED(qsol)) THEN
375	CALL start_init_phys( iml, jml, lon_in, lat_in)
376	ENDIF
377	IF ( SIZE(qsol) .NE. SIZE(lon_in)*SIZE(lat_in) ) THEN
378	WRITE(,)
379	. 'STARTVAR module has been initialized to the wrong size'
380	STOP
381	ENDIF
382	CALL gr_dyn_fi(1, iml, jml, nbindex, qsol, champ)
383	CASE ('psol')
384	IF ( .NOT.ALLOCATED(psol_dyn)) THEN
385	CALL start_init_dyn( iml, jml, lon_in, lat_in)
386	ENDIF
387	IF (SIZE(psol_dyn) .NE. SIZE(lon_in)*SIZE(lat_in)) THEN
388	WRITE(,)
389	. 'STARTVAR module has been initialized to the wrong size'
390	STOP
391	ENDIF
392	CALL gr_dyn_fi(1, iml, jml, nbindex, psol_dyn, champ)
393	CASE ('zmea')
394	IF ( .NOT.ALLOCATED(relief)) THEN
395	CALL start_init_orog( iml, jml, lon_in, lat_in)
396	ENDIF
397	IF ( SIZE(relief) .NE. SIZE(lon_in)*SIZE(lat_in) ) THEN
398	WRITE(,)
399	. 'STARTVAR module has been initialized to the wrong size'
400	STOP
401	ENDIF
402	CALL gr_dyn_fi(1, iml, jml, nbindex, relief, champ)
403	CASE ('zstd')
404	IF ( .NOT.ALLOCATED(zstd)) THEN
405	CALL start_init_orog( iml, jml, lon_in, lat_in)
406	ENDIF
407	IF ( SIZE(zstd) .NE. SIZE(lon_in)*SIZE(lat_in) ) THEN
408	WRITE(,)
409	. 'STARTVAR module has been initialized to the wrong size'
410	STOP
411	ENDIF
412	CALL gr_dyn_fi(1, iml, jml, nbindex,zstd, champ)
413	CASE ('zsig')
414	IF ( .NOT.ALLOCATED(zsig)) THEN
415	CALL start_init_orog( iml, jml, lon_in, lat_in)
416	ENDIF
417	IF ( SIZE(zsig) .NE. SIZE(lon_in)*SIZE(lat_in) ) THEN
418	WRITE(,)
419	. 'STARTVAR module has been initialized to the wrong size'
420	STOP
421	ENDIF
422	CALL gr_dyn_fi(1, iml, jml, nbindex,zsig, champ)
423	CASE ('zgam')
424	IF ( .NOT.ALLOCATED(zgam)) THEN
425	CALL start_init_orog( iml, jml, lon_in, lat_in)
426	ENDIF
427	IF ( SIZE(zgam) .NE. SIZE(lon_in)*SIZE(lat_in) ) THEN
428	WRITE(,)
429	. 'STARTVAR module has been initialized to the wrong size'
430	STOP
431	ENDIF
432	CALL gr_dyn_fi(1, iml, jml, nbindex,zgam, champ)
433	CASE ('zthe')
434	IF ( .NOT.ALLOCATED(zthe)) THEN
435	CALL start_init_orog( iml, jml, lon_in, lat_in)
436	ENDIF
437	IF ( SIZE(zthe) .NE. SIZE(lon_in)*SIZE(lat_in) ) THEN
438	WRITE(,)
439	. 'STARTVAR module has been initialized to the wrong size'
440	STOP
441	ENDIF
442	CALL gr_dyn_fi(1, iml, jml, nbindex,zthe, champ)
443	CASE ('zpic')
444	IF ( .NOT.ALLOCATED(zpic)) THEN
445	CALL start_init_orog( iml, jml, lon_in, lat_in)
446	ENDIF
447	IF ( SIZE(zpic) .NE. SIZE(lon_in)*SIZE(lat_in) ) THEN
448	WRITE(,)
449	. 'STARTVAR module has been initialized to the wrong size'
450	STOP
451	ENDIF
452	CALL gr_dyn_fi(1, iml, jml, nbindex,zpic, champ)
453	CASE ('zval')
454	IF ( .NOT.ALLOCATED(zval)) THEN
455	CALL start_init_orog( iml, jml, lon_in, lat_in)
456	ENDIF
457	IF ( SIZE(zval) .NE. SIZE(lon_in)*SIZE(lat_in) ) THEN
458	WRITE(,)
459	. 'STARTVAR module has been initialized to the wrong size'
460	STOP
461	ENDIF
462	CALL gr_dyn_fi(1, iml, jml, nbindex,zval, champ)
463	CASE ('rads')
464	champ(:) = 0.0
465	CASE ('snow')
466	champ(:) = 0.0
467	CASE ('deltat')
468	champ(:) = 0.0
469	CASE ('rugmer')
470	champ(:) = 0.001
471	CASE ('agsno')
472	champ(:) = 50.0
473	CASE DEFAULT
474	WRITE(,) 'startget_phys1d'
475	WRITE(,) 'No rule is present to extract variable ',
476	. varname(:LEN_TRIM(varname)),' from any data set'
477	STOP
478	END SELECT
479	ELSE
480	!
481	! If we see tsol we catch it as we may need it for a 3D interpolation
482	!
483	SELECTCASE(varname)
484	CASE ('tsol')
485	IF ( .NOT.ALLOCATED(tsol)) THEN
486	ALLOCATE(tsol(SIZE(lon_in),SIZE(lat_in) ))
487	ENDIF
488	CALL gr_fi_dyn(1, iml, jml, nbindex, champ, tsol)
489	END SELECT
490	ENDIF
491	END SUBROUTINE startget_phys1d
492	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
493	!
494	SUBROUTINE start_init_phys( iml, jml, lon_in, lat_in)
495	!
496	INTEGER, INTENT(in) :: iml, jml
497	REAL, INTENT(in) :: lon_in(iml), lat_in(jml)
498	!
499	! LOCAL
500	!
501	REAL :: lev(1), date, dt
502	INTEGER :: itau(1)
503	INTEGER :: llm_tmp, ttm_tmp
504	INTEGER :: i, j
505	!
506	CHARACTER*120 :: physfname
507	LOGICAL :: check=.TRUE.
508	!
509	REAL, ALLOCATABLE :: lon_rad(:), lat_rad(:)
510	REAL, ALLOCATABLE :: var_ana(:,:), tmp_var(:,:)
511	!
512	physfname = 'ECPHY.nc'
513	!
514	IF ( check ) WRITE(,) 'Opening the surface analysis'
515	!
516	CALL flininfo(physfname, iml_phys, jml_phys, llm_tmp,
517	. ttm_tmp, fid_phys)
518	!
519	!
520	ALLOCATE (lat_phys(iml_phys,jml_phys))
521	ALLOCATE (lon_phys(iml_phys,jml_phys))
522	!
523	CALL flinopen(physfname, .FALSE., iml_phys, jml_phys,
524	. llm_tmp, lon_phys, lat_phys, lev, ttm_tmp,
525	. itau, date, dt, fid_phys)
526	!
527	! Allocate the space we will need to get the data out of this file
528	!
529	ALLOCATE(var_ana(iml_phys, jml_phys))
530	!
531	! In case we have a file which is in degrees we do the transformation
532	!
533	ALLOCATE(lon_rad(iml_phys))
534	IF ( MAXVAL(lon_phys(:,:)) .GT. 2.0 * ASIN(1.0) ) THEN
535	lon_rad(:) = lon_phys(:,1) * 2.0 * ASIN(1.0) / 180.0
536	ELSE
537	lon_rad(:) = lon_phys(:,1)
538	ENDIF
539	ALLOCATE(lat_rad(jml_phys))
540	IF ( MAXVAL(lat_phys(:,:)) .GT. 2.0 * ASIN(1.0) ) THEN
541	lat_rad(:) = lat_phys(1,:) * 2.0 * ASIN(1.0) / 180.0
542	ELSE
543	lat_rad(:) = lat_phys(1,:)
544	ENDIF
545	!
546	! We get the two standard varibales
547	! Surface temperature
548	!
549	ALLOCATE(tsol(iml,jml))
550	ALLOCATE(tmp_var(iml-1,jml))
551	!
552	!
553	CALL flinget(fid_phys, 'ST', iml_phys, jml_phys,
554	.llm_tmp, ttm_tmp, 1, 1, var_ana)
555	CALL grille_m(iml_phys, jml_phys, lon_rad, lat_rad,
556	. var_ana, iml-1, jml, lon_in, lat_in, tmp_var)
557	CALL gr_int_dyn(tmp_var, tsol, iml-1, jml)
558	!
559	! Soil moisture
560	!
561	ALLOCATE(qsol(iml,jml))
562	CALL flinget(fid_phys, 'CDSW', iml_phys, jml_phys,
563	. llm_tmp, ttm_tmp, 1, 1, var_ana)
564	CALL grille_m(iml_phys, jml_phys, lon_rad, lat_rad,
565	. var_ana, iml-1, jml, lon_in, lat_in, tmp_var)
566	CALL gr_int_dyn(tmp_var, qsol, iml-1, jml)
567	!
568	CALL flinclo(fid_phys)
569	!
570	END SUBROUTINE start_init_phys
571	!
572	!
573	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
574	!
575	!
576	SUBROUTINE startget_dyn(varname, iml, jml, lon_in, lat_in,
577	. lml, pls, workvar, champ, val_exp)
578	!
579	! ARGUMENTS
580	!
581	CHARACTER(), INTENT(in) :: varname
582	INTEGER, INTENT(in) :: iml, jml, lml
583	REAL, INTENT(in) :: lon_in(iml), lat_in(jml)
584	REAL, INTENT(in) :: pls(iml, jml, lml)
585	REAL, INTENT(in) :: workvar(iml, jml, lml)
586	REAL, INTENT(inout) :: champ(iml, jml, lml)
587	REAL, INTENT(in) :: val_exp
588	!
589	! LOCAL
590	!
591	INTEGER :: il, ij, ii
592	REAL :: xppn, xpps
593	!
594	! C'est vraiment une galere de devoir rajouter tant de commons just pour avoir les aires.
595	! Il faudrait mettre une structure plus flexible et moins dangereuse.
596	!
597	#include "dimensions.h"
598	#include "paramet.h"
599	#include "comgeom2.h"
600	#include "comconst.h"
601	!
602	! This routine only works if the variable does not exist or is constant
603	!
604	IF ( MINVAL(champ(:,:,:)).EQ.MAXVAL(champ(:,:,:)) .AND.
605	. MINVAL(champ(:,:,:)).EQ.val_exp ) THEN
606	!
607	SELECTCASE(varname)
608	CASE ('u')
609	IF ( .NOT.ALLOCATED(psol_dyn)) THEN
610	CALL start_init_dyn( iml, jml, lon_in, lat_in)
611	ENDIF
612	CALL start_inter_3d('U', iml, jml, lml, lon_in,
613	. lat_in, pls, champ)
614	DO il=1,lml
615	DO ij=1,jml
616	DO ii=1,iml-1
617	champ(ii,ij,il) = champ(ii,ij,il) * cu(ii,ij)
618	ENDDO
619	champ(iml,ij, il) = champ(1,ij, il)
620	ENDDO
621	ENDDO
622	CASE ('v')
623	IF ( .NOT.ALLOCATED(psol_dyn)) THEN
624	CALL start_init_dyn( iml, jml, lon_in, lat_in)
625	ENDIF
626	CALL start_inter_3d('V', iml, jml, lml, lon_in,
627	. lat_in, pls, champ)
628	DO il=1,lml
629	DO ij=1,jml
630	DO ii=1,iml-1
631	champ(ii,ij,il) = champ(ii,ij,il) * cv(ii,ij)
632	ENDDO
633	champ(iml,ij, il) = champ(1,ij, il)
634	ENDDO
635	ENDDO
636	CASE ('t')
637	IF ( .NOT.ALLOCATED(psol_dyn)) THEN
638	CALL start_init_dyn( iml, jml, lon_in, lat_in)
639	ENDIF
640	CALL start_inter_3d('TEMP', iml, jml, lml, lon_in,
641	. lat_in, pls, champ)
642
643	CASE ('tpot')
644	IF ( .NOT.ALLOCATED(psol_dyn)) THEN
645	CALL start_init_dyn( iml, jml, lon_in, lat_in)
646	ENDIF
647	CALL start_inter_3d('TEMP', iml, jml, lml, lon_in,
648	. lat_in, pls, champ)
649	IF ( MINVAL(workvar(:,:,:)) .NE. MAXVAL(workvar(:,:,:)) )
650	. THEN
651	DO il=1,lml
652	DO ij=1,jml
653	DO ii=1,iml-1
654	champ(ii,ij,il) = champ(ii,ij,il) * cpp
655	. / workvar(ii,ij,il)
656	ENDDO
657	champ(iml,ij,il) = champ(1,ij,il)
658	ENDDO
659	ENDDO
660	DO il=1,lml
661	xppn = SUM(aire(:,1)*champ(:,1,il))/apoln
662	xpps = SUM(aire(:,jml)*champ(:,jml,il))/apoln
663	champ(:,1,il) = xppn
664	champ(:,jml,il) = xpps
665	ENDDO
666	ELSE
667	WRITE(,)'Could not compute potential temperature as the'
668	WRITE(,)'Exner function is missing or constant.'
669	STOP
670	ENDIF
671	CASE ('q')
672	IF ( .NOT.ALLOCATED(psol_dyn)) THEN
673	CALL start_init_dyn( iml, jml, lon_in, lat_in)
674	ENDIF
675	CALL start_inter_3d('R', iml, jml, lml, lon_in, lat_in,
676	. pls, champ)
677	IF ( MINVAL(workvar(:,:,:)) .NE. MAXVAL(workvar(:,:,:)) )
678	. THEN
679	DO il=1,lml
680	DO ij=1,jml
681	DO ii=1,iml-1
682	champ(ii,ij,il) = 0.01 * champ(ii,ij,il) *
683	. workvar(ii,ij,il)
684	ENDDO
685	champ(iml,ij,il) = champ(1,ij,il)
686	ENDDO
687	ENDDO
688	WHERE ( champ .LT. 0.) champ = 1.0E-10
689	DO il=1,lml
690	xppn = SUM(aire(:,1)*champ(:,1,il))/apoln
691	xpps = SUM(aire(:,jml)*champ(:,jml,il))/apoln
692	champ(:,1,il) = xppn
693	champ(:,jml,il) = xpps
694	ENDDO
695	ELSE
696	WRITE(,)'Could not compute specific humidity as the'
697	WRITE(,)'saturated humidity is missing or constant.'
698	STOP
699	ENDIF
700	CASE DEFAULT
701	WRITE(,) 'startget_dyn'
702	WRITE(,) 'No rule is present to extract variable ',
703	. varname(:LEN_TRIM(varname)),' from any data set'
704	STOP
705	END SELECT
706	ENDIF
707	END SUBROUTINE startget_dyn
708	!
709	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
710	!
711	SUBROUTINE start_init_dyn( iml, jml, lon_in, lat_in)
712	!
713	INTEGER, INTENT(in) :: iml, jml
714	REAL, INTENT(in) :: lon_in(iml), lat_in(jml)
715	!
716	! LOCAL
717	!
718	REAL :: lev(1), date, dt
719	INTEGER :: itau(1)
720	INTEGER :: i, j
721	integer :: iret
722	!
723	CHARACTER*120 :: physfname
724	LOGICAL :: check=.TRUE.
725	!
726	REAL, ALLOCATABLE :: lon_rad(:), lat_rad(:)
727	REAL, ALLOCATABLE :: var_ana(:,:), tmp_var(:,:), z(:,:)
728	REAL, ALLOCATABLE :: xppn(:), xpps(:)
729	LOGICAL :: allo
730	!
731	! Ce n'est pas tres pratique d'avoir a charger 3 include pour avoir la grille du modele
732	!
733	#include "dimensions.h"
734	#include "paramet.h"
735	#include "comgeom2.h"
736	!
737	physfname = 'ECDYN.nc'
738	!
739	IF ( check ) WRITE(,) 'Opening the surface analysis'
740	!
741	CALL flininfo(physfname, iml_dyn, jml_dyn, llm_dyn,
742	. ttm_dyn, fid_dyn)
743	IF ( check ) WRITE(,) 'Values read: ', iml_dyn, jml_dyn,
744	. llm_dyn, ttm_dyn
745	!
746	ALLOCATE (lat_dyn(iml_dyn,jml_dyn), stat=iret)
747	ALLOCATE (lon_dyn(iml_dyn,jml_dyn), stat=iret)
748	ALLOCATE (lev_dyn(llm_dyn), stat=iret)
749	!
750	CALL flinopen(physfname, .FALSE., iml_dyn, jml_dyn, llm_dyn,
751	. lon_dyn, lat_dyn, lev_dyn, ttm_dyn,
752	. itau, date, dt, fid_dyn)
753	!
754
755	allo = allocated (var_ana)
756	if (allo) then
757	DEALLOCATE(var_ana, stat=iret)
758	endif
759	ALLOCATE(var_ana(iml_dyn, jml_dyn), stat=iret)
760
761	allo = allocated (lon_rad)
762	if (allo) then
763	DEALLOCATE(lon_rad, stat=iret)
764	endif
765	ALLOCATE(lon_rad(iml_dyn), stat=iret)
766
767	IF ( MAXVAL(lon_dyn(:,:)) .GT. 2.0 * ASIN(1.0) ) THEN
768	lon_rad(:) = lon_dyn(:,1) * 2.0 * ASIN(1.0) / 180.0
769	ELSE
770	lon_rad(:) = lon_dyn(:,1)
771	ENDIF
772	ALLOCATE(lat_rad(jml_dyn))
773	IF ( MAXVAL(lat_dyn(:,:)) .GT. 2.0 * ASIN(1.0) ) THEN
774	lat_rad(:) = lat_dyn(1,:) * 2.0 * ASIN(1.0) / 180.0
775	ELSE
776	lat_rad(:) = lat_dyn(1,:)
777	ENDIF
778	!
779	ALLOCATE(z(iml, jml))
780	ALLOCATE(tmp_var(iml-1,jml))
781	!
782	CALL flinget(fid_dyn, 'Z', iml_dyn, jml_dyn, 0, ttm_dyn,
783	. 1, 1, var_ana)
784	CALL grille_m(iml_dyn, jml_dyn , lon_rad, lat_rad, var_ana,
785	. iml-1, jml, lon_in, lat_in, tmp_var)
786	CALL gr_int_dyn(tmp_var, z, iml-1, jml)
787	!
788	ALLOCATE(psol_dyn(iml, jml))
789	!
790	CALL flinget(fid_dyn, 'SP', iml_dyn, jml_dyn, 0, ttm_dyn,
791	. 1, 1, var_ana)
792	CALL grille_m(iml_dyn, jml_dyn , lon_rad, lat_rad, var_ana,
793	. iml-1, jml, lon_in, lat_in, tmp_var)
794	CALL gr_int_dyn(tmp_var, psol_dyn, iml-1, jml)
795	!
796	IF ( .NOT.ALLOCATED(tsol)) THEN
797	! These variables may have been allocated by the need to
798	! create a start field for them or by the varibale
799	! coming out of the restart file. In case we dor have it we will initialize it.
800	!
801	CALL start_init_phys( iml, jml, lon_in, lat_in)
802	ELSE
803	IF ( SIZE(tsol) .NE. SIZE(psol_dyn) ) THEN
804	WRITE(,) 'start_init_dyn :'
805	WRITE(,) 'The temperature field we have does not ',
806	. 'have the right size'
807	STOP
808	ENDIF
809	ENDIF
810	IF ( .NOT.ALLOCATED(phis)) THEN
811	!
812	! These variables may have been allocated by the need to create a start field for them or by the varibale
813	! coming out of the restart file. In case we dor have it we will initialize it.
814	!
815	CALL start_init_orog( iml, jml, lon_in, lat_in)
816	!
817	ELSE
818	!
819	IF (SIZE(phis) .NE. SIZE(psol_dyn)) THEN
820	!
821	WRITE(,) 'start_init_dyn :'
822	WRITE(,) 'The orography field we have does not ',
823	. ' have the right size'
824	STOP
825	ENDIF
826	!
827	ENDIF
828	!
829	! PSOL is computed in Pascals
830	!
831	!
832	DO j = 1, jml
833	DO i = 1, iml-1
834	psol_dyn(i,j) = psol_dyn(i,j)*(1.0+(z(i,j)-phis(i,j))
835	. /287.0/tsol(i,j))
836	ENDDO
837	psol_dyn(iml,j) = psol_dyn(1,j)
838	ENDDO
839	!
840	!
841	ALLOCATE(xppn(iml-1))
842	ALLOCATE(xpps(iml-1))
843	!
844	DO i = 1, iml-1
845	xppn(i) = aire( i,1) * psol_dyn( i,1)
846	xpps(i) = aire( i,jml) * psol_dyn( i,jml)
847	ENDDO
848	!
849	DO i = 1, iml
850	psol_dyn(i,1 ) = SUM(xppn)/apoln
851	psol_dyn(i,jml) = SUM(xpps)/apols
852	ENDDO
853	!
854	RETURN
855	!
856	END SUBROUTINE start_init_dyn
857	!
858	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
859	!
860	SUBROUTINE start_inter_3d(varname, iml, jml, lml, lon_in,
861	. lat_in, pls_in, var3d)
862	!
863	! This subroutine gets a variables from a 3D file and does the interpolations needed
864	!
865	!
866	! ARGUMENTS
867	!
868	CHARACTER() :: varname
869	INTEGER :: iml, jml, lml
870	REAL :: lon_in(iml), lat_in(jml), pls_in(iml, jml, lml)
871	REAL :: var3d(iml, jml, lml)
872	!
873	! LOCAL
874	!
875	INTEGER :: ii, ij, il
876	REAL :: bx, by
877	REAL, ALLOCATABLE :: lon_rad(:), lat_rad(:)
878	REAL, ALLOCATABLE :: var_tmp2d(:,:), var_tmp3d(:,:,:)
879	REAL, ALLOCATABLE :: ax(:), ay(:), yder(:)
880	INTEGER, ALLOCATABLE :: lind(:)
881	!
882	LOGICAL :: check = .TRUE.
883	!
884	IF ( .NOT. ALLOCATED(var_ana3d)) THEN
885	ALLOCATE(var_ana3d(iml_dyn, jml_dyn, llm_dyn))
886	ENDIF
887	!
888	!
889	IF ( check) WRITE(,) 'Going into flinget to extract the 3D ',
890	. ' field.', fid_dyn
891	IF ( check) WRITE(,) fid_dyn, varname, iml_dyn, jml_dyn,
892	. llm_dyn,ttm_dyn
893	!
894	CALL flinget(fid_dyn, varname, iml_dyn, jml_dyn, llm_dyn,
895	. ttm_dyn, 1, 1, var_ana3d)
896	!
897	IF ( check) WRITE(,) 'Allocating space for the interpolation',
898	. iml, jml, llm_dyn
899	!
900	ALLOCATE(lon_rad(iml_dyn))
901	IF ( MAXVAL(lon_dyn(:,:)) .GT. 2.0 * ASIN(1.0) ) THEN
902	lon_rad(:) = lon_dyn(:,1) * 2.0 * ASIN(1.0) / 180.0
903	ELSE
904	lon_rad(:) = lon_dyn(:,1)
905	ENDIF
906	ALLOCATE(lat_rad(jml_dyn))
907	IF ( MAXVAL(lat_dyn(:,:)) .GT. 2.0 * ASIN(1.0) ) THEN
908	lat_rad(:) = lat_dyn(1,:) * 2.0 * ASIN(1.0) / 180.0
909	ELSE
910	lat_rad(:) = lat_dyn(1,:)
911	ENDIF
912	!
913	ALLOCATE(var_tmp2d(iml-1, jml))
914	ALLOCATE(var_tmp3d(iml, jml, llm_dyn))
915	ALLOCATE(ax(llm_dyn))
916	ALLOCATE(ay(llm_dyn))
917	ALLOCATE(yder(llm_dyn))
918	ALLOCATE(lind(llm_dyn))
919	!
920	DO il=1,llm_dyn
921	!
922	CALL grille_m(iml_dyn, jml_dyn, lon_rad, lat_rad,
923	.var_ana3d(:,:,il), iml-1, jml, lon_in, lat_in, var_tmp2d)
924	!
925	CALL gr_int_dyn(var_tmp2d, var_tmp3d(:,:,il), iml-1, jml)
926	!
927	ENDDO
928	!
929	! IF needed we return the vertical axis. The spline interpolation
930	! Requires the coordinate to be in increasing order.
931	!
932	IF ( lev_dyn(1) .LT. lev_dyn(llm_dyn)) THEN
933	DO il=1,llm_dyn
934	lind(il) = il
935	ENDDO
936	ELSE
937	DO il=1,llm_dyn
938	lind(il) = llm_dyn-il+1
939	ENDDO
940	ENDIF
941	!
942	DO ij=1,jml
943	DO ii=1,iml-1
944	!
945	ax(:) = lev_dyn(lind(:)) * 100
946	ay(:) = var_tmp3d(ii, ij, lind(:))
947	!
948	CALL SPLINE(ax, ay, llm_dyn, 1.e30, 1.e30, yder)
949	!
950	DO il=1,lml
951	bx = pls_in(ii, ij, il)
952	CALL SPLINT(ax, ay, yder, llm_dyn, bx, by)
953	var3d(ii, ij, il) = by
954	ENDDO
955	!
956	ENDDO
957	var3d(iml, ij, :) = var3d(1, ij, :)
958	ENDDO
959
960	DEALLOCATE(lon_rad)
961	DEALLOCATE(lat_rad)
962	DEALLOCATE(var_tmp2d)
963	DEALLOCATE(var_tmp3d)
964	DEALLOCATE(ax)
965	DEALLOCATE(ay)
966	DEALLOCATE(yder)
967	DEALLOCATE(lind)
968
969	!
970	RETURN
971	!
972	END SUBROUTINE start_inter_3d
973	!
974	END MODULE startvar

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

Download in other formats:

Original Format