Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

startvar.F90 @ 1243

Last change on this file since 1243 was 1056, checked in by slebonnois, 11 years ago
SL: Titan runs ! see DOC/chantiers/commit_importants.log
File size: 30.6 KB

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

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

Download in other formats:

Original Format