Context Navigation

startvar.F90 @ 1319

Last change on this file since 1319 was 1319, checked in by Laurent Fairhead, 14 years ago

Modifications to the start and limit creation routines to account for different

calendars

Modification to phyetat0 to force the mask read in the start files to match the

surface fractions read in the limit file

Force readaerosol.F90 to read in aerosols file with 12 timesteps

Modifications aux routines de créations des fichiers start et limit pour prendre

en compte différents calendriers

Modification à phyetat0 pour forcer le masque lu dans le fichier start à être

compatible avec les fractions de surface lu dans le fichier limit

Forcer readaerosol à ne lire que des fichiers à 12 pas de temps

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

File size: 30.9 KB

Line
1	!
2	! $Id: startvar.F90 1319 2010-02-23 21:29:54Z fairhead $
3	!
4	!*******************************************************************************
5	!
6	MODULE startvar
7	!
8	!*******************************************************************************
9	!
10	!-------------------------------------------------------------------------------
11	! Purpose: Access data from the database of atmospheric to initialize the model.
12	!-------------------------------------------------------------------------------
13	! Comments:
14	!
15	! * This module is designed to work for Earth (and with ioipsl)
16	!
17	! * There are three ways to acces data, depending on the type of field
18	! which needs to be extracted. In any case the call should come after a restget
19	! and should be of the type : CALL startget(...)
20	!
21	! - A 1D variable on the physical grid :
22	! CALL startget(varname, iml, jml, lon_in, lat_in, nbindex, &
23	! champ, val_exp, jml2, lon_in2, lat_in2, ibar )
24	!
25	! - A 2D variable on the dynamical grid :
26	! CALL startget(varname, iml, jml, lon_in, lat_in, &
27	! champ, val_exp, jml2, lon_in2, lat_in2, ibar )
28	!
29	! - A 3D variable on the dynamical grid :
30	! CALL startget(varname, iml, jml, lon_in, lat_in, lml, pls, workvar, &
31	! champ, val_exp, jml2, lon_in2, lat_in2, ibar )
32	!
33	! * Data needs to be in NetCDF format
34	!
35	! * Variables should have the following names in the files:
36	! 'RELIEF' : High resolution orography
37	! 'ST' : Surface temperature
38	! 'CDSW' : Soil moisture
39	! 'Z' : Surface geopotential
40	! 'SP' : Surface pressure
41	! 'U' : East ward wind
42	! 'V' : Northward wind
43	! 'TEMP' : Temperature
44	! 'R' : Relative humidity
45	!
46	! * There is a big mess with the longitude size. Should it be iml or iml+1 ?
47	! I have chosen to use the iml+1 as an argument to this routine and we declare
48	! internaly smaller fields when needed. This needs to be cleared once and for
49	! all in LMDZ. A convention is required.
50	!-------------------------------------------------------------------------------
51	#ifdef CPP_EARTH
52	USE ioipsl
53	IMPLICIT NONE
54
55	PRIVATE
56	PUBLIC startget
57	INTERFACE startget
58	MODULE PROCEDURE startget_phys1d, startget_phys2d, startget_dyn
59	END INTERFACE
60
61	REAL, SAVE :: deg2rad, pi
62	INTEGER, SAVE :: iml_rel, jml_rel
63	INTEGER, SAVE :: fid_phys, iml_phys, jml_phys
64	INTEGER, SAVE :: fid_dyn, iml_dyn, jml_dyn, llm_dyn, ttm_dyn
65	REAL, DIMENSION(:,:), ALLOCATABLE, TARGET, SAVE :: lon_phys, lon_dyn
66	REAL, DIMENSION(:,:), ALLOCATABLE, TARGET, SAVE :: lat_phys, lat_dyn
67	REAL, DIMENSION(:,:), ALLOCATABLE, TARGET, SAVE :: lon_rug, lon_alb, lon_rel
68	REAL, DIMENSION(:,:), ALLOCATABLE, TARGET, SAVE :: lat_rug, lat_alb, lat_rel
69	REAL, DIMENSION(:), ALLOCATABLE, TARGET, SAVE :: levdyn_ini
70	REAL, DIMENSION(:,:), ALLOCATABLE, TARGET, SAVE :: relief, zstd, zsig, zgam
71	REAL, DIMENSION(:,:), ALLOCATABLE, TARGET, SAVE :: masque, zthe, zpic, zval
72	REAL, DIMENSION(:,:), ALLOCATABLE, TARGET, SAVE :: rugo, phis, tsol, qsol
73	REAL, DIMENSION(:,:), ALLOCATABLE, TARGET, SAVE :: psol_dyn
74	REAL, DIMENSION(:,:,:), ALLOCATABLE, TARGET, SAVE :: var_ana3d
75
76	CONTAINS
77
78	!-------------------------------------------------------------------------------
79	!
80	SUBROUTINE startget_phys1d(varname, iml, jml, lon_in, lat_in, nbindex, champ, &
81	val_exp ,jml2, lon_in2, lat_in2, ibar)
82	!
83	!-------------------------------------------------------------------------------
84	! Comment:
85	! This routine only works if the variable does not exist or is constant.
86	!-------------------------------------------------------------------------------
87	! Arguments:
88	CHARACTER(LEN=*), INTENT(IN) :: varname
89	INTEGER, INTENT(IN) :: iml, jml
90	REAL, DIMENSION(iml), INTENT(IN) :: lon_in
91	REAL, DIMENSION(jml), INTENT(IN) :: lat_in
92	INTEGER, INTENT(IN) :: nbindex
93	REAL, DIMENSION(nbindex), INTENT(INOUT) :: champ
94	REAL, INTENT(IN) :: val_exp
95	INTEGER, INTENT(IN) :: jml2
96	REAL, DIMENSION(iml), INTENT(IN) :: lon_in2
97	REAL, DIMENSION(jml2), INTENT(IN) :: lat_in2
98	LOGICAL, INTENT(IN) :: ibar
99	!-------------------------------------------------------------------------------
100	! Local variables:
101	#include "iniprint.h"
102	REAL, DIMENSION(:,:), POINTER :: v2d
103	!-------------------------------------------------------------------------------
104	v2d=>NULL()
105	IF(MINVAL(champ)==MAXVAL(champ).AND.MINVAL(champ)==val_exp) THEN
106
107	!--- CHECKING IF THE FIELD IS KNOWN ; READING UNALLOCATED FILES
108	SELECT CASE(varname)
109	CASE('tsol')
110	IF(.NOT.ALLOCATED(tsol)) &
111	CALL start_init_phys(iml,jml,lon_in,lat_in,jml2,lon_in2,lat_in2,ibar)
112	CASE('qsol')
113	IF(.NOT.ALLOCATED(qsol)) &
114	CALL start_init_phys(iml,jml,lon_in,lat_in,jml2,lon_in2,lat_in2,ibar)
115	CASE('psol')
116	IF(.NOT.ALLOCATED(psol_dyn)) &
117	CALL start_init_dyn (iml,jml,lon_in,lat_in,jml2,lon_in2,lat_in2,ibar)
118	CASE('zmea','zstd','zsig','zgam','zthe','zpic','zval')
119	IF(.NOT.ALLOCATED(relief)) &
120	CALL start_init_orog(iml,jml,lon_in,lat_in)
121	CASE('rads','snow','tslab','seaice','rugmer','agsno')
122	CASE DEFAULT
123	WRITE(lunout,*)'startget_phys1d'
124	WRITE(lunout,*)'No rule is present to extract variable '//TRIM(varname)&
125	//' from any data set'; STOP
126	END SELECT
127
128	!--- SELECTING v2d FOR WANTED VARIABLE AND CHEKING ITS SIZE
129	SELECT CASE(varname)
130	CASE('rads','snow','tslab','seaice'); champ=0.0
131	CASE('rugmer'); champ(:)=0.001
132	CASE('agsno'); champ(:)=50.0
133	CASE DEFAULT
134	SELECT CASE(varname)
135	CASE('tsol'); v2d=>tsol
136	CASE('qsol'); v2d=>qsol
137	CASE('psol'); v2d=>psol_dyn
138	CASE('zmea'); v2d=>relief
139	CASE('zstd'); v2d=>zstd
140	CASE('zsig'); v2d=>zsig
141	CASE('zgam'); v2d=>zgam
142	CASE('zthe'); v2d=>zthe
143	CASE('zpic'); v2d=>zpic
144	CASE('zval'); v2d=>zval
145	END SELECT
146	IF(SIZE(v2d)/=SIZE(lon_in)*SIZE(lat_in)) THEN
147	WRITE(lunout,*)'STARTVAR module has been initialized to the wrong size'
148	STOP
149	END IF
150	CALL gr_dyn_fi(1,iml,jml,nbindex,v2d,champ)
151	END SELECT
152
153	ELSE
154
155	!--- SOME FIELDS ARE CAUGHT: MAY BE NEEDED FOR A 3D INTEPROLATION
156	SELECT CASE(varname)
157	CASE('tsol')
158	IF(.NOT.ALLOCATED(tsol)) ALLOCATE(tsol(iml,jml))
159	CALL gr_fi_dyn(1,iml,jml,nbindex,champ,tsol)
160	END SELECT
161
162	END IF
163
164	END SUBROUTINE startget_phys1d
165	!
166	!-------------------------------------------------------------------------------
167
168
169	!-------------------------------------------------------------------------------
170	!
171	SUBROUTINE startget_phys2d(varname, iml, jml, lon_in, lat_in, champ, val_exp, &
172	jml2, lon_in2, lat_in2 , ibar, msk)
173	!
174	!-------------------------------------------------------------------------------
175	! Comment:
176	! This routine only works if the variable does not exist or is constant.
177	!-------------------------------------------------------------------------------
178	! Arguments:
179	CHARACTER(LEN=*), INTENT(IN) :: varname
180	INTEGER, INTENT(IN) :: iml, jml
181	REAL, DIMENSION(iml), INTENT(IN) :: lon_in
182	REAL, DIMENSION(jml), INTENT(IN) :: lat_in
183	REAL, DIMENSION(iml,jml), INTENT(INOUT) :: champ
184	REAL, INTENT(IN) :: val_exp
185	INTEGER, INTENT(IN) :: jml2
186	REAL, DIMENSION(iml), INTENT(IN) :: lon_in2
187	REAL, DIMENSION(jml2), INTENT(IN) :: lat_in2
188	LOGICAL, INTENT(IN) :: ibar
189	REAL, DIMENSION(iml,jml), INTENT(IN), OPTIONAL :: msk
190	!-------------------------------------------------------------------------------
191	! Local variables:
192	#include "iniprint.h"
193	REAL, DIMENSION(:,:), POINTER :: v2d=>NULL()
194	LOGICAL :: lrelief1, lrelief2
195	!-------------------------------------------------------------------------------
196	v2d=>NULL()
197	lrelief1=(.NOT.ALLOCATED(relief).AND. PRESENT(msk))
198	lrelief2=(.NOT.ALLOCATED(relief).AND..NOT.PRESENT(msk))
199	IF(MINVAL(champ)==MAXVAL(champ).AND.MINVAL(champ)==val_exp) THEN
200
201	!--- CHECKING IF THE FIELD IS KNOWN ; READING UNALLOCATED FILES
202	SELECT CASE(varname)
203	CASE('psol')
204	IF(.NOT.ALLOCATED(psol_dyn)) &
205	CALL start_init_dyn (iml,jml,lon_in,lat_in,jml2,lon_in2,lat_in2,ibar)
206	CASE('relief')
207	IF(lrelief1) CALL start_init_orog(iml,jml,lon_in,lat_in,msk)
208	IF(lrelief2) CALL start_init_orog(iml,jml,lon_in,lat_in)
209	CASE('surfgeo')
210	IF(.NOT.ALLOCATED(phis)) CALL start_init_orog(iml,jml,lon_in,lat_in)
211	CASE('rugosite','masque')
212	IF(.NOT.ALLOCATED(rugo)) CALL start_init_orog(iml,jml,lon_in,lat_in)
213	CASE DEFAULT
214	WRITE(lunout,*)'startget_phys2d'
215	WRITE(lunout,*)'No rule is present to extract variable '//TRIM(varname)&
216	//' from any data set'; STOP
217	END SELECT
218
219	!--- SELECTING v2d FOR WANTED VARIABLE AND CHEKING ITS SIZE
220	SELECT CASE(varname)
221	CASE('psol'); v2d=>psol_dyn
222	CASE('relief'); v2d=>relief
223	CASE('rugosite'); v2d=>rugo
224	CASE('masque'); v2d=>masque
225	CASE('surfgeo'); v2d=>phis
226	END SELECT
227	IF(SIZE(champ)/=SIZE(v2d)) THEN
228	WRITE(lunout,*) 'STARTVAR module has been initialized to the wrong size'
229	STOP
230	END IF
231
232	champ(:,:)=v2d(:,:)
233
234	ELSE
235
236	!--- SOME FIELDS ARE CAUGHT: MAY BE NEEDED FOR A 3D INTEPROLATION
237	SELECT CASE(varname)
238	CASE ('surfgeo')
239	IF(.NOT.ALLOCATED(phis)) ALLOCATE(phis(iml,jml))
240	IF(SIZE(phis)/=SIZE(champ)) THEN
241	WRITE(lunout,*)'STARTVAR module has been initialized to the wrong size'
242	STOP
243	END IF
244	phis(:,:)=champ(:,:)
245	END SELECT
246
247	END IF
248
249	END SUBROUTINE startget_phys2d
250	!
251	!-------------------------------------------------------------------------------
252
253
254	!-------------------------------------------------------------------------------
255	!
256	SUBROUTINE startget_dyn(varname, iml, jml, lon_in, lat_in, lml, pls,workvar,&
257	champ, val_exp, jml2, lon_in2, lat_in2, ibar)
258	!-------------------------------------------------------------------------------
259	! Comment:
260	! This routine only works if the variable does not exist or is constant.
261	!-------------------------------------------------------------------------------
262	! Arguments:
263	CHARACTER(LEN=*), INTENT(IN) :: varname
264	INTEGER, INTENT(IN) :: iml, jml
265	REAL, DIMENSION(iml), INTENT(IN) :: lon_in
266	REAL, DIMENSION(jml), INTENT(IN) :: lat_in
267	INTEGER, INTENT(IN) :: lml
268	REAL, DIMENSION(iml,jml,lml), INTENT(IN) :: pls, workvar
269	REAL, DIMENSION(iml,jml,lml), INTENT(INOUT) :: champ
270	REAL, INTENT(IN) :: val_exp
271	INTEGER, INTENT(IN) :: jml2
272	REAL, DIMENSION(iml), INTENT(IN) :: lon_in2
273	REAL, DIMENSION(jml2), INTENT(IN) :: lat_in2
274	LOGICAL, INTENT(IN) :: ibar
275	!-------------------------------------------------------------------------------
276	! Local variables:
277	#include "iniprint.h"
278	#include "dimensions.h"
279	#include "comconst.h"
280	#include "paramet.h"
281	#include "comgeom2.h"
282	REAL, DIMENSION(:,:,:), POINTER :: v3d=>NULL()
283	CHARACTER(LEN=10) :: vname
284	INTEGER :: il
285	REAL :: xppn, xpps
286	!-------------------------------------------------------------------------------
287	NULLIFY(v3d)
288	IF(MINVAL(champ)==MAXVAL(champ).AND.MINVAL(champ)==val_exp) THEN
289
290	!--- READING UNALLOCATED FILES
291	IF(.NOT.ALLOCATED(psol_dyn)) &
292	CALL start_init_dyn(iml,jml,lon_in,lat_in,jml2,lon_in2,lat_in2,ibar)
293
294	!--- CHECKING IF THE FIELD IS KNOWN AND INTERPOLATING 3D FIELDS
295	SELECT CASE(varname)
296	CASE('u'); vname='U'
297	CASE('v'); vname='V'
298	CASE('t','tpot'); vname='TEMP'
299	CASE('q'); vname='R'
300	CASE DEFAULT
301	WRITE(lunout,*)'startget_dyn'
302	WRITE(lunout,*)'No rule is present to extract variable '//TRIM(varname)&
303	//' from any data set'; STOP
304	END SELECT
305	CALL start_inter_3d(TRIM(vname), iml, jml, lml, lon_in, lat_in, jml2, &
306	lon_in2, lat_in2, pls, champ,ibar )
307
308	!--- COMPUTING THE REQUIRED FILED
309	SELECT CASE(varname)
310	CASE('u') !--- Eastward wind
311	DO il=1,lml; champ(:,:,il)=champ(:,:,il)*cu(:,1:jml); END DO
312	champ(iml,:,:)=champ(1,:,:)
313
314	CASE('v') !--- Northward wind
315	DO il=1,lml; champ(:,:,il)=champ(:,:,il)*cv(:,1:jml); END DO
316	champ(iml,:,:)=champ(1,:,:)
317
318	CASE('tpot') !--- Temperature
319	IF(MINVAL(workvar)/=MAXVAL(workvar)) THEN
320	champ=champ*cpp/workvar
321	DO il=1,lml
322	xppn = SUM(aire(:,1 )*champ(:,1 ,il))/apoln
323	xpps = SUM(aire(:,jml)*champ(:,jml,il))/apols
324	champ(:,1 ,il) = xppn
325	champ(:,jml,il) = xpps
326	END DO
327	ELSE
328	WRITE(lunout,*)'Could not compute potential temperature as the'
329	WRITE(lunout,*)'Exner function is missing or constant.'; STOP
330	END IF
331
332	CASE('q') !--- Relat. humidity
333	IF(MINVAL(workvar)/=MAXVAL(workvar)) THEN
334	champ=0.01champworkvar
335	WHERE(champ<0.) champ=1.0E-10
336	DO il=1,lml
337	xppn = SUM(aire(:,1 )*champ(:,1 ,il))/apoln
338	xpps = SUM(aire(:,jml)*champ(:,jml,il))/apols
339	champ(:,1 ,il) = xppn
340	champ(:,jml,il) = xpps
341	END DO
342	ELSE
343	WRITE(lunout,*)'Could not compute specific humidity as the'
344	WRITE(lunout,*)'saturated humidity is missing or constant.'; STOP
345	END IF
346
347	END SELECT
348
349	END IF
350
351	END SUBROUTINE startget_dyn
352	!
353	!-------------------------------------------------------------------------------
354
355
356	!-------------------------------------------------------------------------------
357	!
358	SUBROUTINE start_init_orog(iml,jml,lon_in,lat_in,masque_lu)
359	!
360	!-------------------------------------------------------------------------------
361	! Arguments:
362	INTEGER, INTENT(IN) :: iml, jml
363	REAL, DIMENSION(iml), INTENT(IN) :: lon_in
364	REAL, DIMENSION(jml), INTENT(IN) :: lat_in
365	REAL, DIMENSION(iml,jml), INTENT(IN), OPTIONAL :: masque_lu
366	!-------------------------------------------------------------------------------
367	! Local variables:
368	#include "iniprint.h"
369	CHARACTER(LEN=25) :: title
370	CHARACTER(LEN=120) :: orofname
371	LOGICAL :: check=.TRUE.
372	REAL, DIMENSION(1) :: lev
373	REAL :: date, dt
374	INTEGER, DIMENSION(1) :: itau
375	INTEGER :: fid, llm_tmp, ttm_tmp
376	REAL, DIMENSION(:,:), ALLOCATABLE :: relief_hi, tmp_var
377	REAL, DIMENSION(:), ALLOCATABLE :: lon_rad, lat_rad, lon_ini, lat_ini
378	!-------------------------------------------------------------------------------
379	pi=2.0*ASIN(1.0); deg2rad=pi/180.0
380
381	orofname = 'Relief.nc'; title='RELIEF'
382	IF(check) WRITE(lunout,*)'Reading the high resolution orography'
383	CALL flininfo(orofname, iml_rel, jml_rel, llm_tmp, ttm_tmp, fid)
384
385	ALLOCATE(lat_rel(iml_rel,jml_rel),lon_rel(iml_rel,jml_rel))
386	CALL flinopen(orofname, .FALSE., iml_rel, jml_rel, llm_tmp, lon_rel, lat_rel,&
387	lev, ttm_tmp, itau, date, dt, fid)
388	ALLOCATE(relief_hi(iml_rel,jml_rel))
389	CALL flinget(fid, title, iml_rel, jml_rel, llm_tmp, ttm_tmp, 1, 1, relief_hi)
390	CALL flinclo(fid)
391
392	!--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS
393	ALLOCATE(lon_ini(iml_rel),lat_ini(jml_rel))
394	lon_ini(:)=lon_rel(:,1); IF(MAXVAL(lon_rel)>pi) lon_ini=lon_ini*deg2rad
395	lat_ini(:)=lat_rel(1,:); IF(MAXVAL(lat_rel)>pi) lat_ini=lat_ini*deg2rad
396
397	!--- FIELDS ARE PROCESSED TO BE ON STANDARD ANGULAR DOMAINS
398	ALLOCATE(lon_rad(iml_rel),lat_rad(jml_rel))
399	CALL conf_dat2d(title, iml_rel, jml_rel, lon_ini, lat_ini, lon_rad, lat_rad, &
400	relief_hi, .FALSE.)
401	DEALLOCATE(lon_ini,lat_ini)
402
403	!--- COMPUTING THE REQUIRED FIELDS USING ROUTINE grid_noro
404	IF(check) WRITE(lunout,*)'Computes all parameters needed for gravity wave dra&
405	&g code'
406
407	ALLOCATE(phis(iml,jml)) ! Geopotentiel au sol
408	ALLOCATE(zstd(iml,jml)) ! Deviation standard de l'orographie sous-maille
409	ALLOCATE(zsig(iml,jml)) ! Pente de l'orographie sous-maille
410	ALLOCATE(zgam(iml,jml)) ! Anisotropie de l'orographie sous maille
411	ALLOCATE(zthe(iml,jml)) ! Orientation axe +grande pente d'oro sous maille
412	ALLOCATE(zpic(iml,jml)) ! Hauteur pics de la SSO
413	ALLOCATE(zval(iml,jml)) ! Hauteur vallees de la SSO
414	ALLOCATE(relief(iml,jml)) ! Orographie moyenne
415	ALLOCATE(masque(iml,jml)) ! Masque terre ocean
416	masque = -99999.
417	IF(PRESENT(masque_lu)) masque=masque_lu
418
419	CALL grid_noro(iml_rel, jml_rel, lon_rad, lat_rad, relief_hi, iml-1, jml, &
420	lon_in, lat_in, phis, relief, zstd, zsig, zgam, zthe, zpic, zval, masque)
421	phis = phis * 9.81
422
423	!--- SURFACE ROUGHNESS COMPUTATION (UNUSED FOR THE MOMENT !!! )
424	IF(check) WRITE(lunout,*)'Compute surface roughness induced by the orography'
425	ALLOCATE(rugo (iml ,jml))
426	ALLOCATE(tmp_var(iml-1,jml))
427	CALL rugsoro(iml_rel, jml_rel, lon_rad, lat_rad, relief_hi, iml-1, jml, &
428	lon_in, lat_in, tmp_var)
429	rugo(1:iml-1,:)=tmp_var; rugo(iml,:)=tmp_var(1,:)
430	DEALLOCATE(relief_hi,tmp_var,lon_rad,lat_rad)
431	RETURN
432
433	END SUBROUTINE start_init_orog
434	!
435	!-------------------------------------------------------------------------------
436
437
438	!-------------------------------------------------------------------------------
439	!
440	SUBROUTINE start_init_phys(iml,jml,lon_in,lat_in,jml2,lon_in2,lat_in2,ibar)
441	!
442	!-------------------------------------------------------------------------------
443	! Arguments:
444	INTEGER, INTENT(IN) :: iml, jml
445	REAL, DIMENSION(iml), INTENT(IN) :: lon_in
446	REAL, DIMENSION(jml), INTENT(IN) :: lat_in
447	INTEGER, INTENT(IN) :: jml2
448	REAL, DIMENSION(iml), INTENT(IN) :: lon_in2
449	REAL, DIMENSION(jml2), INTENT(IN) :: lat_in2
450	LOGICAL, INTENT(IN) :: ibar
451	!-------------------------------------------------------------------------------
452	! Local variables:
453	#include "iniprint.h"
454	CHARACTER(LEN=25) :: title
455	CHARACTER(LEN=120) :: physfname
456	LOGICAL :: check=.TRUE.
457	REAL :: date, dt
458	INTEGER, DIMENSION(1) :: itau
459	INTEGER :: llm_tmp, ttm_tmp
460	REAL, DIMENSION(:,:), ALLOCATABLE :: var_ana
461	REAL, DIMENSION(:), ALLOCATABLE :: lon_rad, lat_rad, lon_ini, lat_ini
462	REAL, DIMENSION(:), ALLOCATABLE :: levphys_ini
463	!-------------------------------------------------------------------------------
464	physfname = 'ECPHY.nc'; pi=2.0*ASIN(1.0); deg2rad=pi/180.0
465	IF(check) WRITE(lunout,*)'Opening the surface analysis'
466	CALL flininfo(physfname, iml_phys, jml_phys, llm_tmp, ttm_tmp, fid_phys)
467
468	ALLOCATE(lat_phys(iml_phys,jml_phys))
469	ALLOCATE(lon_phys(iml_phys,jml_phys))
470	ALLOCATE(levphys_ini(llm_tmp))
471	CALL flinopen(physfname, .FALSE., iml_phys, jml_phys, llm_tmp, lon_phys, &
472	lat_phys, levphys_ini, ttm_tmp, itau, date, dt, fid_phys)
473	DEALLOCATE(levphys_ini)
474
475	!--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS
476	ALLOCATE(lon_ini(iml_phys),lat_ini(jml_phys))
477	lon_ini(:)=lon_phys(:,1); IF(MAXVAL(lon_phys)>pi) lon_ini=lon_ini*deg2rad
478	lat_ini(:)=lat_phys(1,:); IF(MAXVAL(lat_phys)>pi) lat_ini=lat_ini*deg2rad
479
480	ALLOCATE(var_ana(iml_phys,jml_phys),lon_rad(iml_phys),lat_rad(jml_phys))
481
482	!--- SURFACE TEMPERATURE
483	title='ST'
484	ALLOCATE(tsol(iml,jml))
485	CALL flinget(fid_phys,title,iml_phys,jml_phys,llm_tmp,ttm_tmp,1,1,var_ana)
486	CALL conf_dat2d(title,iml_phys, jml_phys, lon_ini, lat_ini, lon_rad, lat_rad,&
487	var_ana , ibar )
488	CALL interp_startvar(title, ibar, .TRUE., &
489	iml_phys, jml_phys, lon_rad, lat_rad, var_ana, iml, jml, jml-1, &
490	lon_in, lat_in, lon_in2, lat_in2, tsol)
491
492	!--- SOIL MOISTURE
493	title='CDSW'
494	ALLOCATE(qsol(iml,jml))
495	CALL flinget(fid_phys,title,iml_phys,jml_phys,llm_tmp,ttm_tmp,1,1,var_ana)
496	CALL conf_dat2d(title,iml_phys, jml_phys, lon_ini, lat_ini, lon_rad, lat_rad,&
497	var_ana, ibar )
498	CALL interp_startvar(title, ibar, .TRUE., &
499	iml_phys, jml_phys, lon_rad, lat_rad, var_ana, iml, jml, jml-1, &
500	lon_in, lat_in, lon_in2, lat_in2, qsol)
501
502	CALL flinclo(fid_phys)
503
504	DEALLOCATE(var_ana,lon_rad,lat_rad,lon_ini,lat_ini)
505
506	END SUBROUTINE start_init_phys
507	!
508	!-------------------------------------------------------------------------------
509
510
511	!-------------------------------------------------------------------------------
512	!
513	SUBROUTINE start_init_dyn(iml,jml,lon_in,lat_in,jml2,lon_in2,lat_in2,ibar)
514	!
515	!-------------------------------------------------------------------------------
516	! Arguments:
517	INTEGER, INTENT(IN) :: iml, jml
518	REAL, DIMENSION(iml), INTENT(IN) :: lon_in
519	REAL, DIMENSION(jml), INTENT(IN) :: lat_in
520	INTEGER, INTENT(IN) :: jml2
521	REAL, DIMENSION(iml), INTENT(IN) :: lon_in2
522	REAL, DIMENSION(jml2), INTENT(IN) :: lat_in2
523	LOGICAL, INTENT(IN) :: ibar
524	!-------------------------------------------------------------------------------
525	! Local variables:
526	#include "iniprint.h"
527	#include "dimensions.h"
528	#include "paramet.h"
529	#include "comgeom2.h"
530	CHARACTER(LEN=25) :: title
531	CHARACTER(LEN=120) :: physfname
532	LOGICAL :: check=.TRUE.
533	REAL :: date, dt
534	INTEGER, DIMENSION(1) :: itau
535	INTEGER :: i, j
536	REAL, DIMENSION(:,:), ALLOCATABLE :: var_ana, z
537	REAL, DIMENSION(:), ALLOCATABLE :: lon_rad, lat_rad, lon_ini, lat_ini
538	REAL, DIMENSION(:), ALLOCATABLE :: xppn, xpps
539	!-------------------------------------------------------------------------------
540
541	!--- KINETIC ENERGY
542	physfname = 'ECDYN.nc'; pi=2.0*ASIN(1.0); deg2rad=pi/180.0
543	IF(check) WRITE(lunout,*) 'Opening the surface analysis'
544	CALL flininfo(physfname, iml_dyn, jml_dyn, llm_dyn, ttm_dyn, fid_dyn)
545	IF(check) WRITE(lunout,*) 'Values read: ', iml_dyn, jml_dyn, llm_dyn, ttm_dyn
546
547	ALLOCATE(lat_dyn(iml_dyn,jml_dyn))
548	ALLOCATE(lon_dyn(iml_dyn,jml_dyn))
549	ALLOCATE(levdyn_ini(llm_dyn))
550	CALL flinopen(physfname, .FALSE., iml_dyn, jml_dyn, llm_dyn, lon_dyn,lat_dyn,&
551	levdyn_ini, ttm_dyn, itau, date, dt, fid_dyn)
552
553	!--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS
554	ALLOCATE(lon_ini(iml_dyn),lat_ini(jml_dyn))
555	lon_ini(:)=lon_dyn(:,1); IF(MAXVAL(lon_dyn)>pi) lon_ini=lon_ini*deg2rad
556	lat_ini(:)=lat_dyn(1,:); IF(MAXVAL(lat_dyn)>pi) lat_ini=lat_ini*deg2rad
557
558	ALLOCATE(var_ana(iml_dyn,jml_dyn),lon_rad(iml_dyn),lat_rad(jml_dyn))
559
560	!--- SURFACE GEOPOTENTIAL
561	title='Z'
562	ALLOCATE(z(iml,jml))
563	CALL flinget(fid_dyn, title, iml_dyn, jml_dyn, 0, ttm_dyn, 1, 1, var_ana)
564	CALL conf_dat2d(title, iml_dyn, jml_dyn, lon_ini, lat_ini, lon_rad, lat_rad, &
565	var_ana, ibar)
566	CALL interp_startvar(title, ibar, .TRUE., &
567	iml_dyn, jml_dyn, lon_rad, lat_rad, var_ana, iml, jml, jml-1, &
568	lon_in, lat_in, lon_in2, lat_in2, z)
569
570	!--- SURFACE PRESSURE
571	title='SP'
572	ALLOCATE(psol_dyn(iml,jml))
573	CALL flinget(fid_dyn, title, iml_dyn, jml_dyn, 0, ttm_dyn, 1, 1, var_ana)
574	CALL conf_dat2d(title, iml_dyn, jml_dyn, lon_ini, lat_ini, lon_rad, lat_rad, &
575	var_ana, ibar)
576	CALL interp_startvar(title, ibar, .TRUE., &
577	iml_dyn, jml_dyn, lon_rad, lat_rad, var_ana, iml, jml, jml-1, &
578	lon_in, lat_in, lon_in2, lat_in2, psol_dyn)
579
580	DEALLOCATE(var_ana,lon_rad,lat_rad,lon_ini,lat_ini)
581
582	!--- ALLOCATION OF VARIABLES CREATED IN OR COMING FROM RESTART FILE
583	IF(.NOT.ALLOCATED(tsol)) THEN
584	CALL start_init_phys(iml,jml,lon_in,lat_in,jml2,lon_in2,lat_in2,ibar)
585	ELSE
586	IF(SIZE(tsol)/=SIZE(psol_dyn)) THEN
587	WRITE(lunout,*)'start_init_dyn :'
588	WRITE(lunout,*)'The temperature field we have does not have the right size'
589	STOP
590	END IF
591	END IF
592
593	IF(.NOT.ALLOCATED(phis)) THEN
594	CALL start_init_orog(iml,jml,lon_in,lat_in)
595	ELSE
596	IF(SIZE(phis)/=SIZE(psol_dyn)) THEN
597	WRITE(lunout,*)'start_init_dyn :'
598	WRITE(lunout,*)'The orography field we have does not have the right size'
599	STOP
600	END IF
601	END IF
602
603	!--- PSOL IS COMPUTED IN PASCALS
604	DO j = 1, jml
605	DO i = 1, iml-1
606	psol_dyn(i,j) = psol_dyn(i,j)*(1.0+(z(i,j)-phis(i,j))/287.0/tsol(i,j))
607	END DO
608	psol_dyn(iml,j) = psol_dyn(1,j)
609	END DO
610	DEALLOCATE(z)
611
612	ALLOCATE(xppn(iml-1),xpps(iml-1))
613	DO i = 1, iml-1
614	xppn(i) = aire( i,1) * psol_dyn( i,1)
615	xpps(i) = aire( i,jml) * psol_dyn( i,jml)
616	END DO
617	DO i = 1, iml
618	psol_dyn(i,1 ) = SUM(xppn)/apoln
619	psol_dyn(i,jml) = SUM(xpps)/apols
620	END DO
621	DEALLOCATE(xppn,xpps)
622
623	RETURN
624
625	END SUBROUTINE start_init_dyn
626	!
627	!-------------------------------------------------------------------------------
628
629
630	!-------------------------------------------------------------------------------
631	!
632	SUBROUTINE start_inter_3d(varname, iml, jml, lml, lon_in, lat_in, jml2,lon_in2,&
633	lat_in2, pls_in, var3d, ibar)
634	!
635	!-------------------------------------------------------------------------------
636	! Arguments:
637	CHARACTER(LEN=*), INTENT(IN) :: varname
638	INTEGER, INTENT(IN) :: iml, jml, lml
639	REAL, DIMENSION(iml), INTENT(IN) :: lon_in
640	REAL, DIMENSION(jml), INTENT(IN) :: lat_in
641	INTEGER, INTENT(IN) :: jml2
642	REAL, DIMENSION(iml), INTENT(IN) :: lon_in2
643	REAL, DIMENSION(jml2), INTENT(IN) :: lat_in2
644	REAL, DIMENSION(iml,jml,lml), INTENT(IN) :: pls_in
645	REAL, DIMENSION(iml,jml,lml), INTENT(OUT) :: var3d
646	LOGICAL, INTENT(IN) :: ibar
647	!-------------------------------------------------------------------------------
648	! Local variables:
649	#include "iniprint.h"
650	LOGICAL :: check=.TRUE.
651	REAL :: bx, by, chmin, chmax
652	INTEGER :: ii, ij, il
653	REAL, DIMENSION(:,:,:), ALLOCATABLE :: var_tmp3d
654	REAL, DIMENSION(:), ALLOCATABLE :: lon_rad, lat_rad, lon_ini, lat_ini
655	REAL, DIMENSION(:), ALLOCATABLE :: lev_dyn, ax, ay, yder
656	INTEGER, DIMENSION(:), ALLOCATABLE :: lind
657	!-------------------------------------------------------------------------------
658	IF(check) WRITE(lunout,*)'Going into flinget to extract the 3D field.'
659	IF(check) WRITE(lunout,*)fid_dyn, varname, iml_dyn, jml_dyn, llm_dyn, ttm_dyn
660	IF(check) WRITE(lunout,*)'Allocating space for interpolation',iml,jml,llm_dyn
661
662	IF(.NOT.ALLOCATED(var_ana3d)) ALLOCATE(var_ana3d(iml_dyn,jml_dyn,llm_dyn))
663	CALL flinget(fid_dyn,varname,iml_dyn,jml_dyn,llm_dyn,ttm_dyn,1,1,var_ana3d)
664
665	!--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS
666	ALLOCATE(lon_ini(iml_dyn),lat_ini(jml_dyn))
667	lon_ini(:)=lon_dyn(:,1); IF(MAXVAL(lon_dyn)>pi) lon_ini=lon_ini*deg2rad
668	lat_ini(:)=lat_dyn(1,:); IF(MAXVAL(lat_dyn)>pi) lat_ini=lat_ini*deg2rad
669
670	!--- FIELDS ARE PROCESSED TO BE ON STANDARD ANGULAR DOMAINS
671	ALLOCATE(lon_rad(iml_dyn),lat_rad(jml_dyn),lev_dyn(llm_dyn))
672	CALL conf_dat3d (varname, iml_dyn, jml_dyn, llm_dyn, lon_ini, lat_ini, &
673	levdyn_ini, lon_rad, lat_rad, lev_dyn, var_ana3d, ibar)
674	DEALLOCATE(lon_ini,lat_ini)
675
676	!--- COMPUTING THE REQUIRED FIELDS USING ROUTINE grid_noro
677	ALLOCATE(var_tmp3d(iml,jml,llm_dyn))
678	DO il=1,llm_dyn
679	CALL interp_startvar(varname, ibar, il==1, &
680	iml_dyn, jml_dyn, lon_rad, lat_rad, var_ana3d(:,:,il), iml, jml, jml2, &
681	lon_in, lat_in, lon_in2, lat_in2, var_tmp3d(:,:,il))
682	END DO
683	DEALLOCATE(lon_rad,lat_rad)
684
685	ALLOCATE(lind(llm_dyn))
686	DO il=1,llm_dyn
687	lind(il) = llm_dyn-il+1
688	END DO
689
690	!--- VERTICAL INTERPOLATION IS PERFORMED FROM TOP OF ATMOSPHERE TO GROUND
691	ALLOCATE(ax(llm_dyn),ay(llm_dyn),yder(llm_dyn))
692	DO ij=1,jml
693	DO ii=1,iml-1
694	ax(:)=lev_dyn(lind(:))
695	ay(:)=var_tmp3d(ii,ij,lind(:))
696	CALL SPLINE(ax, ay, llm_dyn, 1.e30, 1.e30, yder)
697	DO il=1,lml
698	bx=pls_in(ii,ij,il)
699	CALL SPLINT(ax, ay, yder, llm_dyn, bx, by)
700	var3d(ii,ij,il) = by
701	END DO
702	END DO
703	var3d(iml,ij,:) = var3d(1,ij,:)
704	END DO
705	DEALLOCATE(var_tmp3d,lev_dyn,ax,ay,yder,lind)
706
707	DO il=1,lml
708	CALL minmax(iml*jml,var3d(1,1,il),chmin,chmax)
709	WRITE(lunout,*)' '//TRIM(varname)//' min max l ',il,chmin,chmax
710	END DO
711
712	RETURN
713
714	END SUBROUTINE start_inter_3d
715	!
716	!-------------------------------------------------------------------------------
717
718
719	!-------------------------------------------------------------------------------
720	!
721	SUBROUTINE interp_startvar(vname, ibar, ibeg, ii, jj, lon, lat, vari, &
722	i1, j1, j2, lon1, lat1, lon2, lat2, varo)
723	!
724	!-------------------------------------------------------------------------------
725	! Arguments:
726	CHARACTER(LEN=*), INTENT(IN) :: vname
727	LOGICAL, INTENT(IN) :: ibar, ibeg
728	INTEGER, INTENT(IN) :: ii, jj
729	REAL, DIMENSION(ii), INTENT(IN) :: lon
730	REAL, DIMENSION(jj), INTENT(IN) :: lat
731	REAL, DIMENSION(ii,jj), INTENT(IN) :: vari
732	INTEGER, INTENT(IN) :: i1, j1, j2
733	REAL, DIMENSION(i1), INTENT(IN) :: lon1
734	REAL, DIMENSION(j1), INTENT(IN) :: lat1
735	REAL, DIMENSION(i1), INTENT(IN) :: lon2
736	REAL, DIMENSION(j2), INTENT(IN) :: lat2
737	REAL, DIMENSION(i1,j1), INTENT(OUT) :: varo
738	!-------------------------------------------------------------------------------
739	! Local variables:
740	#include "iniprint.h"
741	REAL, DIMENSION(i1-1,j1) :: vtmp
742	!-------------------------------------------------------------------------------
743	IF(ibar) THEN
744	IF(ibeg) THEN
745	WRITE(lunout,*) &
746	'---------------------------------------------------------------'
747	WRITE(lunout,*) &
748	'$$$ Utilisation de l interpolation barycentrique pour '//TRIM(vname)//' $$$'
749	WRITE(lunout,*) &
750	'---------------------------------------------------------------'
751	END IF
752	CALL inter_barxy(ii, jj-1, lon, lat, vari, i1-1, j2, lon2, lat2, j1, vtmp)
753	ELSE
754	CALL grille_m (ii, jj, lon, lat, vari, i1-1, j1, lon1, lat1, vtmp)
755	END IF
756	CALL gr_int_dyn(vtmp, varo, i1-1, j1)
757
758	END SUBROUTINE interp_startvar
759	!
760	!-------------------------------------------------------------------------------
761
762	#endif
763	! of #ifdef CPP_EARTH
764
765	END MODULE startvar
766	!
767	!*******************************************************************************

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

Download in other formats:

Original Format