Context Navigation

startvar.F90 @ 4089

Last change on this file since 4089 was 1425, checked in by lguez, 14 years ago

Replaced Numerical Recipes procedures for spline interpolation (not in
the public domain) by procedures from the Pchip package in the Slatec
library. This only affects the program "ce0l", not the program
"gcm". Tested on Brodie SX8 with "-debug" and "-prod", "-parallel
none" and "-parallel mpi". "start.nc" and "limit.nc" are
changed. "startphy.nc" is not changed. The relative change is of order
1e-7 or less. The revision makes the program faster (tested on Brodie
with "-prod -d 144x142x39", CPU time is 38 s, instead of 54
s). Procedures from Slatec are untouched, except for
"i1mach.F". Created procedures "pchfe_95" and "pchsp_95" which are
wrappers for "pchfe" and "pchsp" from Slatec. "pchfe_95" and
"pchsp_95" have a safer and simpler interface.

Replaced "make" by "sxgmake" in "arch-SX8_BRODIE.fcm". Added files for
compilation by FCM with "g95".

In "arch-linux-32bit.fcm", replaced "pgf90" by "pgf95". There was no
difference between "dev" and "debug" so added "-O1" to "dev". Added
debugging options. Removed "-Wl,-Bstatic
-L/usr/lib/gcc-lib/i386-linux/2.95.2", which usually produces an error
at link-time.

Bash is now ubiquitous while KornShell? is not so use Bash instead of
KornShell? in FCM.

Replaced some statements "write(6,*)" by "write(lunout,*)". Replaced
"stop" by "stop 1" in the case where "abort_gcm" is called with "ierr
/= 0". Removed "stop" statements at the end of procedures
"limit_netcdf" and main program "ce0l" (why not let the program end
normally?).

Made some arrays automatic instead of allocatable in "start_inter_3d".

Zeroed "wake_pe", "fm_therm", "entr_therm" and "detr_therm" in
"dyn3dpar/etat0_netcdf.F90". The parallel and sequential results of
"ce0l" are thus identical.

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

File size: 31.2 KB

Line
1	!
2	! $Id: startvar.F90 1425 2010-09-02 13:45:23Z fhourdin $
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