Context Navigation

startvar.F90 @ 1

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

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

Download in other formats:

Original Format