Context Navigation

startvar.F90 @ 1638

Last change on this file since 1638 was 1638, checked in by slebonnois, 8 years ago
SL: corrections for Venus newstart
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	USE comconst_mod
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 "paramet.h"
275	#include "comgeom2.h"
276	INTEGER :: iml, jml
277	INTEGER :: lml
278	INTEGER :: jml2
279	REAL, DIMENSION(:,:,:), POINTER :: v3d=>NULL()
280	CHARACTER(LEN=10) :: vname
281	INTEGER :: il
282	REAL :: xppn, xpps
283	!-------------------------------------------------------------------------------
284	NULLIFY(v3d)
285	IF(MINVAL(champ)==MAXVAL(champ).AND.MINVAL(champ)==val_exp) THEN
286
287	iml = assert_eq((/size(lon_in), size(pls, 1), size(workvar, 1), &
288	& size(champ, 1), size(lon_in2)/), "startget_dyn iml")
289	jml = assert_eq(size(lat_in), size(pls, 2), size(workvar, 2), &
290	& size(champ, 2), "startget_dyn jml")
291	lml = assert_eq(size(pls, 3), size(workvar, 3), size(champ, 3), &
292	& "startget_dyn lml")
293	jml2 = size(lat_in2)
294
295	!--- READING UNALLOCATED FILES
296	IF(.NOT.ALLOCATED(psol_dyn)) &
297	CALL start_init_dyn(iml,jml,lon_in,lat_in,jml2,lon_in2,lat_in2,ibar)
298
299	!--- CHECKING IF THE FIELD IS KNOWN AND INTERPOLATING 3D FIELDS
300	SELECT CASE(varname)
301	CASE('u'); vname='U'
302	CASE('v'); vname='V'
303	CASE('t','tpot'); vname='TEMP'
304	CASE('q'); vname='R'
305	CASE DEFAULT
306	WRITE(lunout,*)'startget_dyn'
307	WRITE(lunout,*)'No rule is present to extract variable '//TRIM(varname)&
308	//' from any data set'; STOP
309	END SELECT
310	CALL start_inter_3d(TRIM(vname), iml, jml, lml, lon_in, lat_in, jml2, &
311	lon_in2, lat_in2, pls, champ,ibar )
312
313	!--- COMPUTING THE REQUIRED FILED
314	SELECT CASE(varname)
315	CASE('u') !--- Eastward wind
316	DO il=1,lml; champ(:,:,il)=champ(:,:,il)*cu(:,1:jml); END DO
317	champ(iml,:,:)=champ(1,:,:)
318
319	CASE('v') !--- Northward wind
320	DO il=1,lml; champ(:,:,il)=champ(:,:,il)*cv(:,1:jml); END DO
321	champ(iml,:,:)=champ(1,:,:)
322
323	CASE('tpot') !--- Temperature
324	IF(MINVAL(workvar)/=MAXVAL(workvar)) THEN
325	champ=champ*cpp/workvar
326	DO il=1,lml
327	xppn = SUM(aire(:,1 )*champ(:,1 ,il))/apoln
328	xpps = SUM(aire(:,jml)*champ(:,jml,il))/apols
329	champ(:,1 ,il) = xppn
330	champ(:,jml,il) = xpps
331	END DO
332	ELSE
333	WRITE(lunout,*)'Could not compute potential temperature as the'
334	WRITE(lunout,*)'Exner function is missing or constant.'; STOP
335	END IF
336
337	CASE('q') !--- Relat. humidity
338	IF(MINVAL(workvar)/=MAXVAL(workvar)) THEN
339	champ=0.01champworkvar
340	WHERE(champ<0.) champ=1.0E-10
341	DO il=1,lml
342	xppn = SUM(aire(:,1 )*champ(:,1 ,il))/apoln
343	xpps = SUM(aire(:,jml)*champ(:,jml,il))/apols
344	champ(:,1 ,il) = xppn
345	champ(:,jml,il) = xpps
346	END DO
347	ELSE
348	WRITE(lunout,*)'Could not compute specific humidity as the'
349	WRITE(lunout,*)'saturated humidity is missing or constant.'; STOP
350	END IF
351
352	END SELECT
353
354	END IF
355
356	END SUBROUTINE startget_dyn
357	!
358	!-------------------------------------------------------------------------------
359
360
361	!-------------------------------------------------------------------------------
362	!
363	SUBROUTINE start_init_orog(iml,jml,lon_in,lat_in)
364
365	USE comconst_mod
366	USE grid_atob_m, ONLY: rugsoro
367
368	!-------------------------------------------------------------------------------
369	! Arguments:
370	INTEGER, INTENT(IN) :: iml, jml
371	REAL, DIMENSION(iml), INTENT(IN) :: lon_in
372	REAL, DIMENSION(jml), INTENT(IN) :: lat_in
373	!-------------------------------------------------------------------------------
374	! Local variables:
375	#include "iniprint.h"
376	CHARACTER(LEN=25) :: title
377	CHARACTER(LEN=120) :: orofname
378	LOGICAL :: check=.TRUE.
379	REAL, DIMENSION(1) :: lev
380	REAL :: date, dt
381	INTEGER, DIMENSION(1) :: itau
382	INTEGER :: fid, llm_tmp, ttm_tmp
383	REAL, DIMENSION(:,:), ALLOCATABLE :: relief_hi, tmp_var
384	REAL, DIMENSION(:), ALLOCATABLE :: lon_rad, lat_rad, lon_ini, lat_ini
385	!-------------------------------------------------------------------------------
386	pi=2.0*ASIN(1.0); deg2rad=pi/180.0
387
388	orofname = 'Relief.nc'; title='RELIEF'
389	IF(check) WRITE(lunout,*)'Reading the high resolution orography'
390	CALL flininfo(orofname, iml_rel, jml_rel, llm_tmp, ttm_tmp, fid)
391
392	ALLOCATE(lat_rel(iml_rel,jml_rel),lon_rel(iml_rel,jml_rel))
393	CALL flinopen(orofname, .FALSE., iml_rel, jml_rel, llm_tmp, lon_rel, lat_rel,&
394	lev, ttm_tmp, itau, date, dt, fid)
395	ALLOCATE(relief_hi(iml_rel,jml_rel))
396	CALL flinget(fid, title, iml_rel, jml_rel, llm_tmp, ttm_tmp, 1, 1, relief_hi)
397	CALL flinclo(fid)
398
399	!--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS
400	ALLOCATE(lon_ini(iml_rel),lat_ini(jml_rel))
401	lon_ini(:)=lon_rel(:,1); IF(MAXVAL(lon_rel)>pi) lon_ini=lon_ini*deg2rad
402	lat_ini(:)=lat_rel(1,:); IF(MAXVAL(lat_rel)>pi) lat_ini=lat_ini*deg2rad
403
404	!--- FIELDS ARE PROCESSED TO BE ON STANDARD ANGULAR DOMAINS
405	ALLOCATE(lon_rad(iml_rel),lat_rad(jml_rel))
406	CALL conf_dat2d(title, iml_rel, jml_rel, lon_ini, lat_ini, lon_rad, lat_rad, &
407	relief_hi, .FALSE.)
408	DEALLOCATE(lon_ini,lat_ini)
409
410	!--- COMPUTING THE REQUIRED FIELDS USING ROUTINE grid_noro
411	IF(check) WRITE(lunout,*)'Computes all parameters needed for gravity wave dra&
412	&g code'
413
414	ALLOCATE(phis(iml,jml)) ! Geopotentiel au sol
415	ALLOCATE(zstd(iml,jml)) ! Deviation standard de l'orographie sous-maille
416	ALLOCATE(zsig(iml,jml)) ! Pente de l'orographie sous-maille
417	ALLOCATE(zgam(iml,jml)) ! Anisotropie de l'orographie sous maille
418	ALLOCATE(zthe(iml,jml)) ! Orientation axe +grande pente d'oro sous maille
419	ALLOCATE(zpic(iml,jml)) ! Hauteur pics de la SSO
420	ALLOCATE(zval(iml,jml)) ! Hauteur vallees de la SSO
421	ALLOCATE(relief(iml,jml)) ! Orographie moyenne
422
423	CALL grid_noro(iml_rel, jml_rel, lon_rad, lat_rad, relief_hi, iml-1, jml, &
424	lon_in, lat_in, phis, relief, zstd, zsig, zgam, zthe, zpic, zval)
425	phis = phis * g
426
427	!--- SURFACE ROUGHNESS COMPUTATION (UNUSED FOR THE MOMENT !!! )
428	IF(check) WRITE(lunout,*)'Compute surface roughness induced by the orography'
429	ALLOCATE(rugo (iml ,jml))
430	CALL rugsoro(lon_rad, lat_rad, relief_hi, &
431	lon_in, lat_in, rugo)
432	DEALLOCATE(relief_hi,lon_rad,lat_rad)
433	RETURN
434
435	END SUBROUTINE start_init_orog
436	!
437	!-------------------------------------------------------------------------------
438
439
440	!-------------------------------------------------------------------------------
441	!
442	SUBROUTINE start_init_phys(iml,jml,lon_in,lat_in,jml2,lon_in2,lat_in2,ibar)
443	!
444	!-------------------------------------------------------------------------------
445	! Arguments:
446	INTEGER, INTENT(IN) :: iml, jml
447	REAL, DIMENSION(iml), INTENT(IN) :: lon_in
448	REAL, DIMENSION(jml), INTENT(IN) :: lat_in
449	INTEGER, INTENT(IN) :: jml2
450	REAL, DIMENSION(iml), INTENT(IN) :: lon_in2
451	REAL, DIMENSION(jml2), INTENT(IN) :: lat_in2
452	LOGICAL, INTENT(IN) :: ibar
453	!-------------------------------------------------------------------------------
454	! Local variables:
455	#include "iniprint.h"
456	CHARACTER(LEN=25) :: title
457	CHARACTER(LEN=120) :: physfname
458	LOGICAL :: check=.TRUE.
459	REAL :: date, dt
460	INTEGER, DIMENSION(1) :: itau
461	INTEGER :: llm_tmp, ttm_tmp
462	REAL, DIMENSION(:,:), ALLOCATABLE :: var_ana
463	REAL, DIMENSION(:), ALLOCATABLE :: lon_rad, lat_rad, lon_ini, lat_ini
464	REAL, DIMENSION(:), ALLOCATABLE :: levphys_ini
465	!-------------------------------------------------------------------------------
466	physfname = 'ECPHY.nc'; pi=2.0*ASIN(1.0); deg2rad=pi/180.0
467	IF(check) WRITE(lunout,*)'Opening the surface analysis'
468	CALL flininfo(physfname, iml_phys, jml_phys, llm_tmp, ttm_tmp, fid_phys)
469
470	ALLOCATE(lat_phys(iml_phys,jml_phys))
471	ALLOCATE(lon_phys(iml_phys,jml_phys))
472	ALLOCATE(levphys_ini(llm_tmp))
473	CALL flinopen(physfname, .FALSE., iml_phys, jml_phys, llm_tmp, lon_phys, &
474	lat_phys, levphys_ini, ttm_tmp, itau, date, dt, fid_phys)
475	DEALLOCATE(levphys_ini)
476
477	!--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS
478	ALLOCATE(lon_ini(iml_phys),lat_ini(jml_phys))
479	lon_ini(:)=lon_phys(:,1); IF(MAXVAL(lon_phys)>pi) lon_ini=lon_ini*deg2rad
480	lat_ini(:)=lat_phys(1,:); IF(MAXVAL(lat_phys)>pi) lat_ini=lat_ini*deg2rad
481
482	ALLOCATE(var_ana(iml_phys,jml_phys),lon_rad(iml_phys),lat_rad(jml_phys))
483
484	!--- SURFACE TEMPERATURE
485	title='ST'
486	ALLOCATE(tsol(iml,jml))
487	CALL flinget(fid_phys,title,iml_phys,jml_phys,llm_tmp,ttm_tmp,1,1,var_ana)
488	CALL conf_dat2d(title,iml_phys, jml_phys, lon_ini, lat_ini, lon_rad, lat_rad,&
489	var_ana , ibar )
490	CALL interp_startvar(title, ibar, .TRUE., &
491	iml_phys, jml_phys, lon_rad, lat_rad, var_ana, iml, jml, jml-1, &
492	lon_in, lat_in, lon_in2, lat_in2, tsol)
493
494	!--- SOIL MOISTURE
495	title='CDSW'
496	ALLOCATE(qsol(iml,jml))
497	CALL flinget(fid_phys,title,iml_phys,jml_phys,llm_tmp,ttm_tmp,1,1,var_ana)
498	CALL conf_dat2d(title,iml_phys, jml_phys, lon_ini, lat_ini, lon_rad, lat_rad,&
499	var_ana, ibar )
500	CALL interp_startvar(title, ibar, .TRUE., &
501	iml_phys, jml_phys, lon_rad, lat_rad, var_ana, iml, jml, jml-1, &
502	lon_in, lat_in, lon_in2, lat_in2, qsol)
503
504	CALL flinclo(fid_phys)
505
506	DEALLOCATE(var_ana,lon_rad,lat_rad,lon_ini,lat_ini)
507
508	END SUBROUTINE start_init_phys
509	!
510	!-------------------------------------------------------------------------------
511
512
513	!-------------------------------------------------------------------------------
514	!
515	SUBROUTINE start_init_dyn(iml,jml,lon_in,lat_in,jml2,lon_in2,lat_in2,ibar)
516	!
517	!-------------------------------------------------------------------------------
518	! Arguments:
519	INTEGER, INTENT(IN) :: iml, jml
520	REAL, DIMENSION(iml), INTENT(IN) :: lon_in
521	REAL, DIMENSION(jml), INTENT(IN) :: lat_in
522	INTEGER, INTENT(IN) :: jml2
523	REAL, DIMENSION(iml), INTENT(IN) :: lon_in2
524	REAL, DIMENSION(jml2), INTENT(IN) :: lat_in2
525	LOGICAL, INTENT(IN) :: ibar
526	!-------------------------------------------------------------------------------
527	! Local variables:
528	#include "iniprint.h"
529	#include "dimensions.h"
530	#include "paramet.h"
531	#include "comgeom2.h"
532	CHARACTER(LEN=25) :: title
533	CHARACTER(LEN=120) :: physfname
534	LOGICAL :: check=.TRUE.
535	REAL :: date, dt
536	INTEGER, DIMENSION(1) :: itau
537	INTEGER :: i, j
538	REAL, DIMENSION(:,:), ALLOCATABLE :: var_ana, z
539	REAL, DIMENSION(:), ALLOCATABLE :: lon_rad, lat_rad, lon_ini, lat_ini
540	REAL, DIMENSION(:), ALLOCATABLE :: xppn, xpps
541	!-------------------------------------------------------------------------------
542
543	!--- KINETIC ENERGY
544	physfname = 'ECDYN.nc'; pi=2.0*ASIN(1.0); deg2rad=pi/180.0
545	IF(check) WRITE(lunout,*) 'Opening the surface analysis'
546	CALL flininfo(physfname, iml_dyn, jml_dyn, llm_dyn, ttm_dyn, fid_dyn)
547	IF(check) WRITE(lunout,*) 'Values read: ', iml_dyn, jml_dyn, llm_dyn, ttm_dyn
548
549	ALLOCATE(lat_dyn(iml_dyn,jml_dyn))
550	ALLOCATE(lon_dyn(iml_dyn,jml_dyn))
551	ALLOCATE(levdyn_ini(llm_dyn))
552	CALL flinopen(physfname, .FALSE., iml_dyn, jml_dyn, llm_dyn, lon_dyn,lat_dyn,&
553	levdyn_ini, ttm_dyn, itau, date, dt, fid_dyn)
554
555	!--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS
556	ALLOCATE(lon_ini(iml_dyn),lat_ini(jml_dyn))
557	lon_ini(:)=lon_dyn(:,1); IF(MAXVAL(lon_dyn)>pi) lon_ini=lon_ini*deg2rad
558	lat_ini(:)=lat_dyn(1,:); IF(MAXVAL(lat_dyn)>pi) lat_ini=lat_ini*deg2rad
559
560	ALLOCATE(var_ana(iml_dyn,jml_dyn),lon_rad(iml_dyn),lat_rad(jml_dyn))
561
562	!--- SURFACE GEOPOTENTIAL
563	title='Z'
564	ALLOCATE(z(iml,jml))
565	CALL flinget(fid_dyn, title, iml_dyn, jml_dyn, 0, ttm_dyn, 1, 1, var_ana)
566	CALL conf_dat2d(title, iml_dyn, jml_dyn, lon_ini, lat_ini, lon_rad, lat_rad, &
567	var_ana, ibar)
568	CALL interp_startvar(title, ibar, .TRUE., &
569	iml_dyn, jml_dyn, lon_rad, lat_rad, var_ana, iml, jml, jml-1, &
570	lon_in, lat_in, lon_in2, lat_in2, z)
571
572	!--- SURFACE PRESSURE
573	title='SP'
574	ALLOCATE(psol_dyn(iml,jml))
575	CALL flinget(fid_dyn, title, iml_dyn, jml_dyn, 0, ttm_dyn, 1, 1, var_ana)
576	CALL conf_dat2d(title, iml_dyn, jml_dyn, lon_ini, lat_ini, lon_rad, lat_rad, &
577	var_ana, ibar)
578	CALL interp_startvar(title, ibar, .TRUE., &
579	iml_dyn, jml_dyn, lon_rad, lat_rad, var_ana, iml, jml, jml-1, &
580	lon_in, lat_in, lon_in2, lat_in2, psol_dyn)
581
582	DEALLOCATE(var_ana,lon_rad,lat_rad,lon_ini,lat_ini)
583
584	!--- ALLOCATION OF VARIABLES CREATED IN OR COMING FROM RESTART FILE
585	IF(.NOT.ALLOCATED(tsol)) THEN
586	CALL start_init_phys(iml,jml,lon_in,lat_in,jml2,lon_in2,lat_in2,ibar)
587	ELSE
588	IF(SIZE(tsol)/=SIZE(psol_dyn)) THEN
589	WRITE(lunout,*)'start_init_dyn :'
590	WRITE(lunout,*)'The temperature field we have does not have the right size'
591	STOP
592	END IF
593	END IF
594
595	IF(.NOT.ALLOCATED(phis)) THEN
596	CALL start_init_orog(iml,jml,lon_in,lat_in)
597	ELSE
598	IF(SIZE(phis)/=SIZE(psol_dyn)) THEN
599	WRITE(lunout,*)'start_init_dyn :'
600	WRITE(lunout,*)'The orography field we have does not have the right size'
601	STOP
602	END IF
603	END IF
604
605	!--- PSOL IS COMPUTED IN PASCALS
606	DO j = 1, jml
607	DO i = 1, iml-1
608	psol_dyn(i,j) = psol_dyn(i,j)*(1.0+(z(i,j)-phis(i,j))/287.0/tsol(i,j))
609	END DO
610	psol_dyn(iml,j) = psol_dyn(1,j)
611	END DO
612	DEALLOCATE(z)
613
614	ALLOCATE(xppn(iml-1),xpps(iml-1))
615	DO i = 1, iml-1
616	xppn(i) = aire( i,1) * psol_dyn( i,1)
617	xpps(i) = aire( i,jml) * psol_dyn( i,jml)
618	END DO
619	DO i = 1, iml
620	psol_dyn(i,1 ) = SUM(xppn)/apoln
621	psol_dyn(i,jml) = SUM(xpps)/apols
622	END DO
623	DEALLOCATE(xppn,xpps)
624
625	RETURN
626
627	END SUBROUTINE start_init_dyn
628	!
629	!-------------------------------------------------------------------------------
630
631
632	!-------------------------------------------------------------------------------
633	!
634	SUBROUTINE start_inter_3d(varname, iml, jml, lml, lon_in, lat_in, jml2, &
635	lon_in2, lat_in2, pls_in, var3d, ibar)
636
637	use pchsp_95_m, only: pchsp_95
638	use pchfe_95_m, only: pchfe_95
639
640	! Arguments:
641	CHARACTER(LEN=*), INTENT(IN) :: varname
642	INTEGER, INTENT(IN) :: iml, jml, lml
643	REAL, DIMENSION(iml), INTENT(IN) :: lon_in
644	REAL, DIMENSION(jml), INTENT(IN) :: lat_in
645	INTEGER, INTENT(IN) :: jml2
646	REAL, DIMENSION(iml), INTENT(IN) :: lon_in2
647	REAL, DIMENSION(jml2), INTENT(IN) :: lat_in2
648	REAL, DIMENSION(iml, jml, lml), INTENT(IN) :: pls_in
649	REAL, DIMENSION(iml, jml, lml), INTENT(OUT) :: var3d
650	LOGICAL, INTENT(IN) :: ibar
651	!----------------------------------------------------------------------------
652	! Local variables:
653	#include "iniprint.h"
654	LOGICAL:: check=.TRUE., skip
655	REAL chmin, chmax
656	INTEGER ii, ij, il, ierr
657	integer n_extrap ! number of extrapolated points
658	REAL, DIMENSION(iml, jml, llm_dyn):: var_tmp3d
659	REAL, DIMENSION(:), ALLOCATABLE :: lon_rad, lat_rad, lon_ini, lat_ini
660	REAL, DIMENSION(llm_dyn):: lev_dyn, ax, ay, yder
661
662	!---------------------------------------------------------------------------
663	IF(check) WRITE(lunout, *)'Going into flinget to extract the 3D field.'
664	IF(check) WRITE(lunout, *) fid_dyn, varname, iml_dyn, jml_dyn, llm_dyn, &
665	ttm_dyn
666	IF(check) WRITE(lunout, *) 'Allocating space for interpolation', iml, jml, &
667	llm_dyn
668
669	IF(.NOT.ALLOCATED(var_ana3d)) ALLOCATE(var_ana3d(iml_dyn, jml_dyn, llm_dyn))
670	CALL flinget(fid_dyn, varname, iml_dyn, jml_dyn, llm_dyn, ttm_dyn, 1, 1, &
671	var_ana3d)
672
673	!--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS
674	ALLOCATE(lon_ini(iml_dyn), lat_ini(jml_dyn))
675	lon_ini(:)=lon_dyn(:, 1); IF(MAXVAL(lon_dyn)>pi) lon_ini=lon_ini*deg2rad
676	lat_ini(:)=lat_dyn(1, :); IF(MAXVAL(lat_dyn)>pi) lat_ini=lat_ini*deg2rad
677
678	!--- FIELDS ARE PROCESSED TO BE ON STANDARD ANGULAR DOMAINS
679	ALLOCATE(lon_rad(iml_dyn), lat_rad(jml_dyn))
680	CALL conf_dat3d (varname, iml_dyn, jml_dyn, llm_dyn, lon_ini, lat_ini, &
681	levdyn_ini, lon_rad, lat_rad, lev_dyn, var_ana3d, ibar)
682	DEALLOCATE(lon_ini, lat_ini)
683
684	!--- COMPUTING THE REQUIRED FIELDS USING ROUTINE grid_noro
685	DO il=1, llm_dyn
686	CALL interp_startvar(varname, ibar, il==1, iml_dyn, jml_dyn, lon_rad, &
687	lat_rad, var_ana3d(:, :, il), iml, jml, jml2, lon_in, lat_in, &
688	lon_in2, lat_in2, var_tmp3d(:, :, il))
689	END DO
690	DEALLOCATE(lon_rad, lat_rad)
691
692	!--- VERTICAL INTERPOLATION IS PERFORMED FROM TOP OF ATMOSPHERE TO GROUND
693	ax = lev_dyn(llm_dyn:1:-1)
694	skip = .false.
695	n_extrap = 0
696	DO ij=1, jml
697	DO ii=1, iml-1
698	ay = var_tmp3d(ii, ij, llm_dyn:1:-1)
699	yder = pchsp_95(ax, ay, ibeg=2, iend=2, vc_beg=0., vc_end=0.)
700	CALL pchfe_95(ax, ay, yder, skip, pls_in(ii, ij, lml:1:-1), &
701	var3d(ii, ij, lml:1:-1), ierr)
702	if (ierr < 0) stop 1
703	n_extrap = n_extrap + ierr
704	END DO
705	END DO
706	if (n_extrap /= 0) then
707	print *, "start_inter_3d pchfe_95: n_extrap = ", n_extrap
708	end if
709	var3d(iml, :, :) = var3d(1, :, :)
710
711	DO il=1, lml
712	CALL minmax(iml*jml, var3d(1, 1, il), chmin, chmax)
713	WRITE(lunout, *)' '//TRIM(varname)//' min max l ', il, chmin, chmax
714	END DO
715
716	END SUBROUTINE start_inter_3d
717	!
718	!-------------------------------------------------------------------------------
719
720
721	!-------------------------------------------------------------------------------
722	!
723	SUBROUTINE interp_startvar(vname, ibar, ibeg, ii, jj, lon, lat, vari, &
724	i1, j1, j2, lon1, lat1, lon2, lat2, varo)
725	!
726	!-------------------------------------------------------------------------------
727
728	USE inter_barxy_m, only: inter_barxy
729	USE grid_atob_m, only: grille_m
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 (lon, lat, vari, 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