Context Navigation

etat0dyn_netcdf.F90 @ 2301

Last change on this file since 2301 was 2301, checked in by Ehouarn Millour, 9 years ago
Moved ioipsl_getin_p_mod.F90 to "phylmd" beacause files in "misc" should rely on neither physics nor dynamics modules (or include) files. Moved etat0dyn_netcdf.F90 to "dynlonlat_phylonlat/phylmd" as it is only used by ce0l and moreover requires "grid_atob_m" which is in "dynlonlat_phylonlat". With these adjustements, compiling with '-nophys' works again. EM
File size: 32.3 KB

Line
1	MODULE etat0dyn
2	!
3	!*******************************************************************************
4	! Purpose: Create dynamical initial state using atmospheric fields from a
5	! database of atmospheric to initialize the model.
6	!-------------------------------------------------------------------------------
7	! Comments:
8	!
9	! * This module is designed to work for Earth (and with ioipsl)
10	!
11	! * etat0dyn_netcdf routine can access to NetCDF data through the following
12	! routine (to be called after restget):
13	! CALL startget_dyn3d(varname, lon_in, lat_in, pls, workvar,&
14	! champ, lon_in2, lat_in2, ibar)
15	!
16	! * Variables should have the following names in the NetCDF files:
17	! 'U' : East ward wind (in "ECDYN.nc")
18	! 'V' : Northward wind (in "ECDYN.nc")
19	! 'TEMP' : Temperature (in "ECDYN.nc")
20	! 'R' : Relative humidity (in "ECDYN.nc")
21	! 'RELIEF' : High resolution orography (in "Relief.nc")
22	!
23	! * The land mask and corresponding weights can be:
24	! 1) already known (in particular if etat0dyn has been called before) ;
25	! in this case, ANY(masque(:,:)/=-99999.) = .TRUE.
26	! 2) computed using the ocean mask from the ocean model (to ensure ocean
27	! fractions are the same for atmosphere and ocean) for coupled runs.
28	! File name: "o2a.nc" ; variable name: "OceMask"
29	! 3) computed from topography file "Relief.nc" for forced runs.
30	!
31	! * There is a big mess with the longitude size. Should it be iml or iml+1 ?
32	! I have chosen to use the iml+1 as an argument to this routine and we declare
33	! internaly smaller fields when needed. This needs to be cleared once and for
34	! all in LMDZ. A convention is required.
35	!-------------------------------------------------------------------------------
36	USE ioipsl, ONLY: flininfo, flinopen, flinget, flinclo, histclo
37	USE assert_eq_m, ONLY: assert_eq
38	#ifdef CPP_PHYS
39	USE indice_sol_mod, ONLY: epsfra
40	#endif
41	IMPLICIT NONE
42
43	PRIVATE
44	PUBLIC :: etat0dyn_netcdf
45
46	include "iniprint.h"
47	include "dimensions.h"
48	include "paramet.h"
49	include "comgeom2.h"
50	include "comvert.h"
51	include "comconst.h"
52	include "temps.h"
53	include "comdissnew.h"
54	include "serre.h"
55	REAL, SAVE :: deg2rad
56	#ifndef CPP_PHYS
57	REAL, SAVE :: epsfra= 1.E-5
58	#endif
59	INTEGER, SAVE :: iml_dyn, jml_dyn, llm_dyn, ttm_dyn, fid_dyn
60	REAL, ALLOCATABLE, SAVE :: lon_dyn(:,:), lat_dyn(:,:), levdyn_ini(:)
61	CHARACTER(LEN=120), PARAMETER :: dynfname='ECDYN.nc'
62	CHARACTER(LEN=120), PARAMETER :: orofname='Relief.nc'
63
64	CONTAINS
65
66	!-------------------------------------------------------------------------------
67	!
68	SUBROUTINE etat0dyn_netcdf(ib, masque, phis)
69	!
70	!-------------------------------------------------------------------------------
71	! Purpose: Create dynamical initial states.
72	!-------------------------------------------------------------------------------
73	! Notes: 1) This routine is designed to work for Earth
74	! 2) If masque(:,:)/=-99999., masque and phis are already known.
75	! Otherwise: read it from ocean model file (coupled run) or compute it.
76	!-------------------------------------------------------------------------------
77	USE control_mod
78	#ifdef CPP_PHYS
79	#ifdef CPP_EARTH
80	USE regr_lat_time_coefoz_m, ONLY: regr_lat_time_coefoz
81	USE regr_pr_o3_m, ONLY: regr_pr_o3
82	USE press_coefoz_m, ONLY: press_coefoz
83	#endif
84	#endif
85	USE exner_hyb_m, ONLY: exner_hyb
86	USE exner_milieu_m, ONLY: exner_milieu
87	USE infotrac
88	USE filtreg_mod
89	IMPLICIT NONE
90	!-------------------------------------------------------------------------------
91	! Arguments:
92	LOGICAL, INTENT(IN) :: ib !--- Barycentric interpolation
93	REAL, INTENT(INOUT) :: masque(iip1,jjp1) !--- Land-ocean mask
94	REAL, INTENT(INOUT) :: phis (iip1,jjp1) !--- Ground geopotential
95	!#ifndef CPP_EARTH
96	! WRITE(lunout,*)'etat0dyn_netcdf: Earth-specific routine, needs Earth physics'
97	!#else
98	!-------------------------------------------------------------------------------
99	! Local variables:
100	CHARACTER(LEN=256) :: modname, fmt
101	INTEGER :: i, j, l, ji, itau, iday
102	REAL :: xpn, xps, time, phystep
103	REAL, DIMENSION(iip1,jjp1) :: psol, masque_tmp
104	REAL, DIMENSION(iip1,jjp1,llm+1) :: p3d
105	REAL, DIMENSION(iip1,jjp1,llm) :: uvent, t3d, tpot, qsat, qd
106	REAL, DIMENSION(iip1,jjp1,llm) :: pk, pls, y, masse
107	REAL, DIMENSION(iip1,jjm ,llm) :: vvent
108	REAL, DIMENSION(ip1jm ,llm) :: pbarv
109	REAL, DIMENSION(ip1jmp1 ,llm) :: pbaru, phi, w
110	REAL, DIMENSION(ip1jmp1) :: pks
111	REAL, DIMENSION(iim) :: xppn, xpps
112	REAL, ALLOCATABLE :: q3d(:,:,:,:)
113	!-------------------------------------------------------------------------------
114	modname='etat0dyn_netcdf'
115
116	deg2rad = pi/180.0
117
118	! Initializations for tracers and filter
119	!*******************************************************************************
120	CALL infotrac_init
121	CALL inifilr()
122
123	! Compute ground geopotential and possibly the mask.
124	!*******************************************************************************
125	masque_tmp(:,:)=masque(:,:)
126	CALL start_init_orog0(rlonv, rlatu, phis, masque_tmp)
127	WRITE(fmt,"(i4,'i1)')")iip1 ; fmt='('//ADJUSTL(fmt)
128	IF(ALL(masque==-99999.)) THEN !--- KEEP NEW MASK
129	masque=masque_tmp
130	IF(prt_level>=1) THEN
131	WRITE(lunout,*)'BUILT MASK :'
132	WRITE(lunout,fmt) NINT(masque)
133	END IF
134	WHERE( masque(:,:)<EPSFRA) masque(:,:)=0.
135	WHERE(1.-masque(:,:)<EPSFRA) masque(:,:)=1.
136	END IF
137
138	! Compute psol AND tsol, knowing phis.
139	!*******************************************************************************
140	CALL start_init_dyn(rlonv, rlatu, rlonu, rlatv, ib, phis, psol)
141
142	! Mid-levels pressure computation
143	!*******************************************************************************
144	CALL pression(ip1jmp1, ap, bp, psol, p3d) !--- Update p3d
145	IF(pressure_exner) THEN !--- Update pk, pks
146	CALL exner_hyb (ip1jmp1,psol,p3d,pks,pk)
147	ELSE
148	CALL exner_milieu(ip1jmp1,psol,p3d,pks,pk)
149	END IF
150	pls(:,:,:)=preff(pk(:,:,:)/cpp)*(1./kappa) !--- Update pls
151
152	! Update uvent, vvent, t3d and tpot
153	!*******************************************************************************
154	uvent(:,:,:) = 0.0 ; vvent(:,:,:) = 0.0 ; t3d (:,:,:) = 0.0
155	CALL startget_dyn3d('u' ,rlonu,rlatu,pls,y ,uvent,rlonv,rlatv,ib)
156	CALL startget_dyn3d('v' ,rlonv,rlatv,pls(:,:jjm,:),y(:,:jjm,:),vvent, &
157	& rlonu,rlatu(:jjm),ib)
158	CALL startget_dyn3d('t' ,rlonv,rlatu,pls,y ,t3d ,rlonu,rlatv,ib)
159	tpot(:,:,:)=t3d(:,:,:)
160	CALL startget_dyn3d('tpot',rlonv,rlatu,pls,pk,tpot,rlonu,rlatv,ib)
161
162	WRITE(lunout,*) 'T3D min,max:',MINVAL(t3d(:,:,:)),MAXVAL(t3d(:,:,:))
163	WRITE(lunout,*) 'PLS min,max:',MINVAL(pls(:,:,:)),MAXVAL(pls(:,:,:))
164
165	! Humidity at saturation computation
166	!*******************************************************************************
167	WRITE(lunout,*) 'avant q_sat'
168	CALL q_sat(llmjjp1iip1, t3d, pls, qsat)
169	WRITE(lunout,*) 'apres q_sat'
170	WRITE(lunout,*) 'QSAT min,max:',MINVAL(qsat(:,:,:)),MAXVAL(qsat(:,:,:))
171	! WRITE(lunout,*) 'QSAT :',qsat(10,20,:)
172	qd (:,:,:) = 0.0
173	CALL startget_dyn3d('q',rlonv,rlatu,pls,qsat,qd,rlonu,rlatv,ib)
174	ALLOCATE(q3d(iip1,jjp1,llm,nqtot)); q3d(:,:,:,:)=0.0 ; q3d(:,:,:,1)=qd(:,:,:)
175	CALL flinclo(fid_dyn)
176
177	#ifdef CPP_PHYS
178	#ifdef CPP_EARTH
179	! Parameterization of ozone chemistry:
180	!*******************************************************************************
181	! Look for ozone tracer:
182	DO i=1,nqtot; IF(ANY(["O3","o3"]==tname(i))) EXIT; END DO
183	IF(i/=nqtot+1) THEN
184	CALL regr_lat_time_coefoz
185	CALL press_coefoz
186	CALL regr_pr_o3(p3d, q3d(:,:,:,i))
187	q3d(:,:,:,i)=q3d(:,:,:,i)*48./ 29. !--- Mole->mass fraction
188	END IF
189
190	#endif
191	#endif
192	q3d(iip1,:,:,:)=q3d(1,:,:,:)
193
194	! Intermediate computation
195	!*******************************************************************************
196	CALL massdair(p3d,masse)
197	WRITE(lunout,*)' ALPHAX ',alphax
198	DO l=1,llm
199	xppn(:)=aire(1:iim,1 )*masse(1:iim,1 ,l)
200	xpps(:)=aire(1:iim,jjp1)*masse(1:iim,jjp1,l)
201	xpn=SUM(xppn)/apoln
202	xps=SUM(xpps)/apols
203	masse(:,1 ,l)=xpn
204	masse(:,jjp1,l)=xps
205	END DO
206
207	! Writing
208	!*******************************************************************************
209	CALL inidissip(lstardis, nitergdiv, nitergrot, niterh, tetagdiv, tetagrot, &
210	tetatemp, vert_prof_dissip)
211	WRITE(lunout,*)'sortie inidissip'
212	itau=0
213	itau_dyn=0
214	itau_phy=0
215	iday=dayref+itau/day_step
216	time=FLOAT(itau-(iday-dayref)*day_step)/day_step
217	IF(time>1.) THEN
218	time=time-1
219	iday=iday+1
220	END IF
221	day_ref=dayref
222	annee_ref=anneeref
223	CALL geopot( ip1jmp1, tpot, pk, pks, phis, phi )
224	WRITE(lunout,*)'sortie geopot'
225	CALL caldyn0( itau, uvent, vvent, tpot, psol, masse, pk, phis, &
226	phi, w, pbaru, pbarv, time+iday-dayref)
227	WRITE(lunout,*)'sortie caldyn0'
228	CALL dynredem0( "start.nc", dayref, phis)
229	WRITE(lunout,*)'sortie dynredem0'
230	CALL dynredem1( "start.nc", 0.0, vvent, uvent, tpot, q3d, masse, psol)
231	WRITE(lunout,*)'sortie dynredem1'
232	CALL histclo()
233
234	!#endif
235	! of ifdef CPP_EARTH
236
237	END SUBROUTINE etat0dyn_netcdf
238	!
239	!-------------------------------------------------------------------------------
240
241
242	!#ifdef CPP_EARTH
243	!-------------------------------------------------------------------------------
244	!
245	SUBROUTINE startget_dyn3d(var, lon_in, lat_in, pls, workvar,&
246	champ, lon_in2, lat_in2, ibar)
247	!-------------------------------------------------------------------------------
248	IMPLICIT NONE
249	!===============================================================================
250	! Purpose: Compute some quantities (u,v,t,q,tpot) using variables U,V,TEMP and R
251	! (3D fields) of file dynfname.
252	!-------------------------------------------------------------------------------
253	! Note: An input auxilliary field "workvar" has to be specified in two cases:
254	! * for "q": the saturated humidity.
255	! * for "tpot": the Exner function.
256	!===============================================================================
257	! Arguments:
258	CHARACTER(LEN=*), INTENT(IN) :: var
259	REAL, INTENT(IN) :: lon_in(:) ! dim (iml)
260	REAL, INTENT(IN) :: lat_in(:) ! dim (jml)
261	REAL, INTENT(IN) :: pls (:, :, :) ! dim (iml, jml, lml)
262	REAL, INTENT(IN) :: workvar(:, :, :) ! dim (iml, jml, lml)
263	REAL, INTENT(INOUT) :: champ (:, :, :) ! dim (iml, jml, lml)
264	REAL, INTENT(IN) :: lon_in2(:) ! dim (iml)
265	REAL, INTENT(IN) :: lat_in2(:) ! dim (jml2)
266	LOGICAL, INTENT(IN) :: ibar
267	!-------------------------------------------------------------------------------
268	! Local variables:
269	CHARACTER(LEN=10) :: vname
270	CHARACTER(LEN=256) :: msg, modname="startget_dyn3d"
271	INTEGER :: iml, jml, jml2, lml, il
272	REAL :: xppn, xpps
273	!-------------------------------------------------------------------------------
274	iml=assert_eq([SIZE(lon_in),SIZE(pls,1),SIZE(workvar,1),SIZE(champ,1), &
275	& SIZE(lon_in2)], TRIM(modname)//" iml")
276	jml=assert_eq( SIZE(lat_in),SIZE(pls,2),SIZE(workvar,2),SIZE(champ,2), &
277	& TRIM(modname)//" jml")
278	lml=assert_eq( SIZE(pls,3),SIZE(workvar,3),SIZE(champ,3), &
279	& TRIM(modname)//" lml")
280	jml2=SIZE(lat_in2)
281
282	!--- CHECK IF THE FIELD IS KNOWN
283	SELECT CASE(var)
284	CASE('u'); vname='U'
285	CASE('v'); vname='V'
286	CASE('t'); vname='TEMP'
287	CASE('q'); vname='R'; msg='humidity as the saturated humidity'
288	CASE('tpot'); msg='potential temperature as the Exner function'
289	CASE DEFAULT; msg='No rule to extract variable '//TRIM(var)
290	CALL abort_gcm(modname,TRIM(msg)//' from any data set',1)
291	END SELECT
292
293	!--- CHECK IF SOMETHING IS MISSING
294	IF((var=='tpot'.OR.var=='q').AND.MINVAL(workvar)==MAXVAL(workvar)) THEN
295	msg='Could not compute '//TRIM(msg)//' is missing or constant.'
296	CALL abort_gcm(modname,TRIM(msg),1)
297	END IF
298
299	!--- INTERPOLATE 3D FIELD IF NEEDED
300	IF(var/='tpot') CALL start_inter_3d(TRIM(vname),lon_in,lat_in,lon_in2, &
301	lat_in2,pls,champ,ibar)
302
303	!--- COMPUTE THE REQUIRED FILED
304	SELECT CASE(var)
305	CASE('u'); DO il=1,lml; champ(:,:,il)=champ(:,:,il)*cu(:,1:jml); END DO
306	champ(iml,:,:)=champ(1,:,:) !--- Eastward wind
307
308	CASE('v'); DO il=1,lml; champ(:,:,il)=champ(:,:,il)*cv(:,1:jml); END DO
309	champ(iml,:,:)=champ(1,:,:) !--- Northward wind
310
311	CASE('tpot','q')
312	IF(var=='tpot') THEN; champ=champ*cpp/workvar !--- Potential temperature
313	ELSE; champ=champ.01workvar !--- Relative humidity
314	WHERE(champ<0.) champ=1.0E-10
315	END IF
316	DO il=1,lml
317	xppn = SUM(aire(:,1 )*champ(:,1 ,il))/apoln
318	xpps = SUM(aire(:,jml)*champ(:,jml,il))/apols
319	champ(:,1 ,il) = xppn
320	champ(:,jml,il) = xpps
321	END DO
322	END SELECT
323
324	END SUBROUTINE startget_dyn3d
325	!
326	!-------------------------------------------------------------------------------
327
328
329	!-------------------------------------------------------------------------------
330	!
331	SUBROUTINE start_init_orog0(lon_in,lat_in,phis,masque)
332	!
333	!-------------------------------------------------------------------------------
334	USE conf_dat_m, ONLY: conf_dat2d
335	IMPLICIT NONE
336	!===============================================================================
337	! Purpose: Compute "phis" just like it would be in start_init_orog.
338	!===============================================================================
339	! Arguments:
340	REAL, INTENT(IN) :: lon_in(:), lat_in(:) ! dim (iml) (jml)
341	REAL, INTENT(INOUT) :: phis(:,:), masque(:,:) ! dim (iml,jml)
342	!-------------------------------------------------------------------------------
343	! Local variables:
344	CHARACTER(LEN=256) :: modname="start_init_orog0"
345	CHARACTER(LEN=256) :: title="RELIEF"
346	INTEGER :: fid, llm_tmp,ttm_tmp, iml,jml, iml_rel,jml_rel, itau(1)
347	REAL :: lev(1), date, dt
348	REAL, ALLOCATABLE :: lon_rad(:), lon_ini(:), lon_rel(:,:), relief_hi(:,:)
349	REAL, ALLOCATABLE :: lat_rad(:), lat_ini(:), lat_rel(:,:)
350	!-------------------------------------------------------------------------------
351	iml=assert_eq(SIZE(lon_in),SIZE(phis,1),SIZE(masque,1),TRIM(modname)//" iml")
352	jml=assert_eq(SIZE(lat_in),SIZE(phis,2),SIZE(masque,2),TRIM(modname)//" jml")
353	IF(iml/=iip1) CALL abort_gcm(TRIM(modname),'iml/=iip1',1)
354	IF(jml/=jjp1) CALL abort_gcm(TRIM(modname),'jml/=jjp1',1)
355	pi=2.0*ASIN(1.0); deg2rad=pi/180.0
356	IF(ANY(phis/=-99999.)) RETURN !--- phis ALREADY KNOWN
357
358	!--- HIGH RESOLUTION OROGRAPHY
359	CALL flininfo(orofname, iml_rel, jml_rel, llm_tmp, ttm_tmp, fid)
360
361	ALLOCATE(lat_rel(iml_rel,jml_rel),lon_rel(iml_rel,jml_rel))
362	CALL flinopen(orofname, .FALSE., iml_rel, jml_rel, llm_tmp, lon_rel, lat_rel,&
363	lev, ttm_tmp, itau, date, dt, fid)
364	ALLOCATE(relief_hi(iml_rel,jml_rel))
365	CALL flinget(fid, title, iml_rel, jml_rel, llm_tmp, ttm_tmp, 1, 1, relief_hi)
366	CALL flinclo(fid)
367
368	!--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS
369	ALLOCATE(lon_ini(iml_rel),lat_ini(jml_rel))
370	lon_ini(:)=lon_rel(:,1); IF(MAXVAL(lon_rel)>pi) lon_ini=lon_ini*deg2rad
371	lat_ini(:)=lat_rel(1,:); IF(MAXVAL(lat_rel)>pi) lat_ini=lat_ini*deg2rad
372
373	!--- FIELDS ARE PROCESSED TO BE ON STANDARD ANGULAR DOMAINS
374	ALLOCATE(lon_rad(iml_rel),lat_rad(jml_rel))
375	CALL conf_dat2d(title, lon_ini, lat_ini, lon_rad, lat_rad, relief_hi, .FALSE.)
376	DEALLOCATE(lon_ini,lat_ini)
377
378	!--- COMPUTING SURFACE GEOPOTENTIAL USING ROUTINE grid_noro0
379	WRITE(lunout,*)
380	WRITE(lunout,)' Compute surface geopotential *'
381
382	!--- CALL OROGRAPHY MODULE (REDUCED VERSION) TO COMPUTE FIELDS
383	CALL grid_noro0(lon_rad, lat_rad, relief_hi, lon_in, lat_in, phis, masque)
384	phis = phis * 9.81
385	phis(iml,:) = phis(1,:)
386	DEALLOCATE(relief_hi,lon_rad,lat_rad)
387
388	END SUBROUTINE start_init_orog0
389	!
390	!-------------------------------------------------------------------------------
391
392
393	!-------------------------------------------------------------------------------
394	!
395	SUBROUTINE grid_noro0(xd,yd,zd,x,y,zphi,mask)
396	!
397	!===============================================================================
398	! Purpose: Extracted from grid_noro to provide geopotential height for dynamics
399	! without any call to physics subroutines.
400	!===============================================================================
401	IMPLICIT NONE
402	!-------------------------------------------------------------------------------
403	! Arguments:
404	REAL, INTENT(IN) :: xd(:), yd(:) !--- INPUT COORDINATES (imdp) (jmdp)
405	REAL, INTENT(IN) :: zd(:,:) !--- INPUT FIELD (imdp,jmdp)
406	REAL, INTENT(IN) :: x(:), y(:) !--- OUTPUT COORDINATES (imar+1) (jmar)
407	REAL, INTENT(OUT) :: zphi(:,:) !--- GEOPOTENTIAL (imar+1,jmar)
408	REAL, INTENT(INOUT):: mask(:,:) !--- MASK (imar+1,jmar)
409	!-------------------------------------------------------------------------------
410	! Local variables:
411	CHARACTER(LEN=256) :: modname="grid_noro0"
412	REAL, ALLOCATABLE :: xusn(:), yusn(:) ! dim (imdp+2*iext) (jmdp+2)
413	REAL, ALLOCATABLE :: zusn(:,:) ! dim (imdp+2*iext,jmdp+2)
414	REAL, ALLOCATABLE :: weight(:,:) ! dim (imar+1,jmar)
415	REAL, ALLOCATABLE :: mask_tmp(:,:), zmea(:,:) ! dim (imar+1,jmar)
416	REAL, ALLOCATABLE :: num_tot(:,:), num_lan(:,:) ! dim (imax,jmax)
417	REAL, ALLOCATABLE :: a(:), b(:) ! dim (imax)
418	REAL, ALLOCATABLE :: c(:), d(:) ! dim (jmax)
419	LOGICAL :: masque_lu
420	INTEGER :: i, ii, imdp, imar, iext
421	INTEGER :: j, jj, jmdp, jmar, nn
422	REAL :: xpi, zlenx, weighx, xincr, zbordnor, zmeanor, zweinor, zbordest
423	REAL :: rad, zleny, weighy, masque, zbordsud, zmeasud, zweisud, zbordoue
424	!-------------------------------------------------------------------------------
425	imdp=assert_eq(SIZE(xd),SIZE(zd,1),TRIM(modname)//" imdp")
426	jmdp=assert_eq(SIZE(yd),SIZE(zd,2),TRIM(modname)//" jmdp")
427	imar=assert_eq(SIZE(x),SIZE(zphi,1),SIZE(mask,1),TRIM(modname)//" imar")-1
428	jmar=assert_eq(SIZE(y),SIZE(zphi,2),SIZE(mask,2),TRIM(modname)//" jmar")
429	IF(imar/=iim) CALL abort_gcm(TRIM(modname),'imar/=iim' ,1)
430	IF(jmar/=jjm+1) CALL abort_gcm(TRIM(modname),'jmar/=jjm+1',1)
431	iext=imdp/10
432	xpi = ACOS(-1.)
433	rad = 6371229.
434
435	!--- ARE WE USING A READ MASK ?
436	masque_lu=ANY(mask/=-99999.); IF(.NOT.masque_lu) mask=0.0
437	WRITE(lunout,*)'Masque lu: ',masque_lu
438
439	!--- EXTENSION OF THE INPUT DATABASE TO PROCEED COMPUTATIONS AT BOUNDARIES:
440	ALLOCATE(xusn(imdp+2*iext))
441	xusn(1 +iext:imdp +iext)=xd(:)
442	xusn(1 : iext)=xd(1+imdp-iext:imdp)-2.*xpi
443	xusn(1+imdp+iext:imdp+2iext)=xd(1 :iext)+2.xpi
444
445	ALLOCATE(yusn(jmdp+2))
446	yusn(1 )=yd(1) +(yd(1) -yd(2))
447	yusn(2:jmdp+1)=yd(:)
448	yusn( jmdp+2)=yd(jmdp)+(yd(jmdp)-yd(jmdp-1))
449
450	ALLOCATE(zusn(imdp+2*iext,jmdp+2))
451	zusn(1 +iext:imdp +iext,2:jmdp+1)=zd (: , :)
452	zusn(1 : iext,2:jmdp+1)=zd (imdp-iext+1:imdp , :)
453	zusn(1+imdp +iext:imdp+2*iext,2:jmdp+1)=zd (1:iext , :)
454	zusn(1 :imdp/2+iext, 1)=zusn(1+imdp/2:imdp +iext, 2)
455	zusn(1+imdp/2+iext:imdp+2*iext, 1)=zusn(1 :imdp/2+iext, 2)
456	zusn(1 :imdp/2+iext, jmdp+2)=zusn(1+imdp/2:imdp +iext,jmdp+1)
457	zusn(1+imdp/2+iext:imdp+2*iext, jmdp+2)=zusn(1 :imdp/2+iext,jmdp+1)
458
459	!--- COMPUTE LIMITS OF MODEL GRIDPOINT AREA (REGULAR GRID)
460	ALLOCATE(a(imar+1),b(imar+1))
461	b(1:imar)=(x(1:imar )+ x(2:imar+1))/2.0
462	b(imar+1)= x( imar+1)+(x( imar+1)-x(imar))/2.0
463	a(1)=x(1)-(x(2)-x(1))/2.0
464	a(2:imar+1)= b(1:imar)
465
466	ALLOCATE(c(jmar),d(jmar))
467	d(1:jmar-1)=(y(1:jmar-1)+ y(2:jmar))/2.0
468	d( jmar )= y( jmar )+(y( jmar)-y(jmar-1))/2.0
469	c(1)=y(1)-(y(2)-y(1))/2.0
470	c(2:jmar)=d(1:jmar-1)
471
472	!--- INITIALIZATIONS:
473	ALLOCATE(weight(imar+1,jmar)); weight(:,:)= 0.0
474	ALLOCATE(zmea (imar+1,jmar)); zmea (:,:)= 0.0
475
476	!--- SUMMATION OVER GRIDPOINT AREA
477	zleny=xpi/REAL(jmdp)*rad
478	xincr=xpi/REAL(jmdp)/2.
479	ALLOCATE(num_tot(imar+1,jmar)); num_tot(:,:)=0.
480	ALLOCATE(num_lan(imar+1,jmar)); num_lan(:,:)=0.
481	DO ii = 1, imar+1
482	DO jj = 1, jmar
483	DO j = 2,jmdp+1
484	zlenx =zleny *COS(yusn(j))
485	zbordnor=(xincr+c(jj)-yusn(j))*rad
486	zbordsud=(xincr-d(jj)+yusn(j))*rad
487	weighy=AMAX1(0.,AMIN1(zbordnor,zbordsud,zleny))
488	IF(weighy/=0) THEN
489	DO i = 2, imdp+2*iext-1
490	zbordest=(xusn(i)-a(ii)+xincr)radCOS(yusn(j))
491	zbordoue=(b(ii)+xincr-xusn(i))radCOS(yusn(j))
492	weighx=AMAX1(0.,AMIN1(zbordest,zbordoue,zlenx))
493	IF(weighx/=0)THEN
494	num_tot(ii,jj)=num_tot(ii,jj)+1.0
495	IF(zusn(i,j)>=1.)num_lan(ii,jj)=num_lan(ii,jj)+1.0
496	weight(ii,jj)=weight(ii,jj)+weighx*weighy
497	zmea (ii,jj)=zmea (ii,jj)+zusn(i,j)weighxweighy !--- MEAN
498	END IF
499	END DO
500	END IF
501	END DO
502	END DO
503	END DO
504
505	!--- COMPUTE PARAMETERS NEEDED BY LOTT & MILLER (1997) AND LOTT (1999) SSO SCHEME
506	IF(.NOT.masque_lu) THEN
507	WHERE(weight(:,1:jmar-1)/=0.0) mask=num_lan(:,:)/num_tot(:,:)
508	END IF
509	nn=COUNT(weight(:,1:jmar-1)==0.0)
510	IF(nn/=0) WRITE(lunout,*)'Problem with weight ; vanishing occurrences: ',nn
511	WHERE(weight/=0.0) zmea(:,:)=zmea(:,:)/weight(:,:)
512
513	!--- MASK BASED ON GROUND MAXIMUM, 10% THRESHOLD (<10%: SURF PARAMS MEANINGLESS)
514	ALLOCATE(mask_tmp(imar+1,jmar)); mask_tmp(:,:)=0.0
515	WHERE(mask>=0.1) mask_tmp = 1.
516	WHERE(weight(:,:)/=0.0)
517	zphi(:,:)=mask_tmp(:,:)*zmea(:,:)
518	zmea(:,:)=mask_tmp(:,:)*zmea(:,:)
519	END WHERE
520	WRITE(lunout,*)' MEAN ORO:' ,MAXVAL(zmea)
521
522	!--- Values at poles
523	zphi(imar+1,:)=zphi(1,:)
524
525	zweinor=SUM(weight(1:imar, 1),DIM=1)
526	zweisud=SUM(weight(1:imar,jmar),DIM=1)
527	zmeanor=SUM(weight(1:imar, 1)*zmea(1:imar, 1),DIM=1)
528	zmeasud=SUM(weight(1:imar,jmar)*zmea(1:imar,jmar),DIM=1)
529	zphi(:,1)=zmeanor/zweinor; zphi(:,jmar)=zmeasud/zweisud
530
531	END SUBROUTINE grid_noro0
532	!
533	!-------------------------------------------------------------------------------
534
535
536	!-------------------------------------------------------------------------------
537	!
538	SUBROUTINE start_init_dyn(lon_in,lat_in,lon_in2,lat_in2,ibar,zs,psol)
539	!
540	!-------------------------------------------------------------------------------
541	IMPLICIT NONE
542	!===============================================================================
543	! Purpose: Compute psol, knowing phis.
544	!===============================================================================
545	! Arguments:
546	REAL, INTENT(IN) :: lon_in (:), lat_in (:) ! dim (iml) (jml)
547	REAL, INTENT(IN) :: lon_in2(:), lat_in2(:) ! dim (iml) (jml2)
548	LOGICAL, INTENT(IN) :: ibar
549	REAL, INTENT(IN) :: zs (:,:) ! dim (iml,jml)
550	REAL, INTENT(OUT) :: psol(:,:) ! dim (iml,jml)
551	!-------------------------------------------------------------------------------
552	! Local variables:
553	CHARACTER(LEN=256) :: modname='start_init_dyn'
554	REAL :: date, dt
555	INTEGER :: iml, jml, jml2, itau(1)
556	REAL, ALLOCATABLE :: lon_rad(:), lon_ini(:), var_ana(:,:)
557	REAL, ALLOCATABLE :: lat_rad(:), lat_ini(:)
558	REAL, ALLOCATABLE :: z(:,:), ps(:,:), ts(:,:)
559	!-------------------------------------------------------------------------------
560	iml=assert_eq(SIZE(lon_in),SIZE(zs,1),SIZE(psol,1),SIZE(lon_in2), &
561	& TRIM(modname)//" iml")
562	jml=assert_eq(SIZE(lat_in),SIZE(zs,2),SIZE(psol,2),TRIM(modname)//" jml")
563	jml2=SIZE(lat_in2)
564
565	WRITE(lunout,*) 'Opening the surface analysis'
566	CALL flininfo(dynfname, iml_dyn, jml_dyn, llm_dyn, ttm_dyn, fid_dyn)
567	WRITE(lunout,*) 'Values read: ', iml_dyn, jml_dyn, llm_dyn, ttm_dyn
568
569	ALLOCATE(lon_dyn(iml_dyn,jml_dyn), lat_dyn(iml_dyn,jml_dyn))
570	ALLOCATE(levdyn_ini(llm_dyn))
571	CALL flinopen(dynfname, .FALSE., iml_dyn, jml_dyn, llm_dyn, &
572	lon_dyn,lat_dyn,levdyn_ini,ttm_dyn,itau,date,dt,fid_dyn)
573
574	!--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS
575	ALLOCATE(lon_ini(iml_dyn),lat_ini(jml_dyn))
576	lon_ini(:)=lon_dyn(:,1); IF(MAXVAL(lon_dyn)>pi) lon_ini=lon_ini*deg2rad
577	lat_ini(:)=lat_dyn(1,:); IF(MAXVAL(lat_dyn)>pi) lat_ini=lat_ini*deg2rad
578
579	ALLOCATE(var_ana(iml_dyn,jml_dyn),lon_rad(iml_dyn),lat_rad(jml_dyn))
580	CALL get_var_dyn('Z',z) !--- SURFACE GEOPOTENTIAL
581	CALL get_var_dyn('SP',ps) !--- SURFACE PRESSURE
582	CALL get_var_dyn('ST',ts) !--- SURFACE TEMPERATURE
583	! CALL flinclo(fid_dyn)
584	DEALLOCATE(var_ana,lon_rad,lat_rad,lon_ini,lat_ini)
585
586	!--- PSOL IS COMPUTED IN PASCALS
587	psol(:iml-1,:) = ps(:iml-1,:)*(1.0+(z(:iml-1,:)-zs(:iml-1,:))/287.0 &
588	& /ts(:iml-1,:))
589	psol(iml,:)=psol(1,:)
590	DEALLOCATE(z,ps,ts)
591	psol(:,1 )=SUM(aire(1:iml-1,1 )*psol(1:iml-1,1 ))/apoln !--- NORTH POLE
592	psol(:,jml)=SUM(aire(1:iml-1,jml)*psol(1:iml-1,jml))/apols !--- SOUTH POLE
593
594	CONTAINS
595
596	!-------------------------------------------------------------------------------
597	!
598	SUBROUTINE get_var_dyn(title,field)
599	!
600	!-------------------------------------------------------------------------------
601	USE conf_dat_m, ONLY: conf_dat2d
602	IMPLICIT NONE
603	!-------------------------------------------------------------------------------
604	! Arguments:
605	CHARACTER(LEN=*), INTENT(IN) :: title
606	REAL, ALLOCATABLE, INTENT(INOUT) :: field(:,:)
607	!-------------------------------------------------------------------------------
608	! Local variables:
609	CHARACTER(LEN=256) :: msg
610	INTEGER :: tllm
611	!-------------------------------------------------------------------------------
612	SELECT CASE(title)
613	CASE('Z'); tllm=0; msg='geopotential'
614	CASE('SP'); tllm=0; msg='surface pressure'
615	CASE('ST'); tllm=llm_dyn; msg='temperature'
616	END SELECT
617	IF(.NOT.ALLOCATED(field)) THEN
618	ALLOCATE(field(iml,jml))
619	CALL flinget(fid_dyn, title, iml_dyn,jml_dyn, tllm, ttm_dyn, 1, 1, var_ana)
620	CALL conf_dat2d(title, lon_ini, lat_ini, lon_rad, lat_rad, var_ana, ibar)
621	CALL interp_startvar(title, ibar, .TRUE., lon_rad, lat_rad, var_ana, &
622	lon_in, lat_in, lon_in2, lat_in2, field)
623	ELSE IF(SIZE(field)/=SIZE(z)) THEN
624	msg='The '//TRIM(msg)//' field we have does not have the right size'
625	CALL abort_gcm(TRIM(modname),msg,1)
626	END IF
627
628	END SUBROUTINE get_var_dyn
629	!
630	!-------------------------------------------------------------------------------
631
632	END SUBROUTINE start_init_dyn
633	!
634	!-------------------------------------------------------------------------------
635
636
637	!-------------------------------------------------------------------------------
638	!
639	SUBROUTINE start_inter_3d(var,lon_in,lat_in,lon_in2,lat_in2,pls_in,var3d,ibar)
640	!
641	!-------------------------------------------------------------------------------
642	USE conf_dat_m, ONLY: conf_dat3d
643	USE pchsp_95_m, ONLY: pchsp_95
644	USE pchfe_95_m, ONLY: pchfe_95
645	IMPLICIT NONE
646	!-------------------------------------------------------------------------------
647	! Arguments:
648	CHARACTER(LEN=*), INTENT(IN) :: var
649	REAL, INTENT(IN) :: lon_in(:), lat_in(:) ! dim (iml) (jml)
650	REAL, INTENT(IN) :: lon_in2(:), lat_in2(:) ! dim (iml) (jml2)
651	REAL, INTENT(IN) :: pls_in(:,:,:) ! dim (iml,jml,lml)
652	REAL, INTENT(OUT) :: var3d (:,:,:) ! dim (iml,jml,lml)
653	LOGICAL, INTENT(IN) :: ibar
654	!-------------------------------------------------------------------------------
655	! Local variables:
656	CHARACTER(LEN=256) :: modname='start_inter_3d'
657	LOGICAL :: skip
658	REAL :: chmin, chmax
659	INTEGER :: iml, jml, lml, jml2, ii, ij, il, ierr
660	INTEGER :: n_extrap ! Extrapolated points number
661	REAL, ALLOCATABLE :: ax(:), lon_rad(:), lon_ini(:), lev_dyn(:), yder(:)
662	REAL, ALLOCATABLE :: ay(:), lat_rad(:), lat_ini(:), var_tmp3d(:,:,:)
663	REAL, ALLOCATABLE, SAVE :: var_ana3d(:,:,:)
664	!-------------------------------------------------------------------------------
665	iml=assert_eq(SIZE(lon_in),SIZE(lon_in2),SIZE(pls_in,1),SIZE(var3d,1),TRIM(modname)//" iml")
666	jml=assert_eq(SIZE(lat_in), SIZE(pls_in,2),SIZE(var3d,2),TRIM(modname)//" jml")
667	lml=assert_eq(SIZE(pls_in,3),SIZE(var3d,3),TRIM(modname)//" lml"); jml2=SIZE(lat_in2)
668
669	WRITE(lunout, *)'Going into flinget to extract the 3D field.'
670	IF(.NOT.ALLOCATED(var_ana3d)) ALLOCATE(var_ana3d(iml_dyn, jml_dyn, llm_dyn))
671	CALL flinget(fid_dyn,var,iml_dyn,jml_dyn,llm_dyn,ttm_dyn,1,1,var_ana3d)
672
673	!--- ANGLES IN DEGREES 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), lev_dyn(llm_dyn))
680	CALL conf_dat3d(var, lon_ini, lat_ini, levdyn_ini, &
681	lon_rad, lat_rad, lev_dyn, var_ana3d, ibar)
682	DEALLOCATE(lon_ini, lat_ini)
683
684	!--- COMPUTE THE REQUIRED FIELDS USING ROUTINE grid_noro
685	ALLOCATE(var_tmp3d(iml,jml,llm_dyn))
686	DO il = 1,llm_dyn
687	CALL interp_startvar(var,ibar,il==1,lon_rad,lat_rad,var_ana3d(:,:,il), &
688	lon_in,lat_in,lon_in2,lat_in2,var_tmp3d(:,:,il))
689	END DO
690	DEALLOCATE(lon_rad, lat_rad)
691
692	!--- VERTICAL INTERPOLATION FROM TOP OF ATMOSPHERE TO GROUND
693	ALLOCATE(ax(llm_dyn),ay(llm_dyn),yder(llm_dyn))
694	ax = lev_dyn(llm_dyn:1:-1)
695	skip = .FALSE.
696	n_extrap = 0
697	DO ij=1, jml
698	DO ii=1, iml-1
699	ay = var_tmp3d(ii, ij, llm_dyn:1:-1)
700	yder = pchsp_95(ax, ay, ibeg=2, iend=2, vc_beg=0., vc_end=0.)
701	CALL pchfe_95(ax, ay, yder, skip, pls_in(ii, ij, lml:1:-1), &
702	var3d(ii, ij, lml:1:-1), ierr)
703	IF(ierr<0) CALL abort_gcm(TRIM(modname),'error in pchfe_95',1)
704	n_extrap = n_extrap + ierr
705	END DO
706	END DO
707	IF(n_extrap/=0) WRITE(lunout,*)TRIM(modname)//" pchfe_95: n_extrap=", n_extrap
708	var3d(iml, :, :) = var3d(1, :, :)
709
710	DO il=1, lml
711	CALL minmax(iml*jml, var3d(1, 1, il), chmin, chmax)
712	WRITE(lunout, *)' '//TRIM(var)//' min max l ', il, chmin, chmax
713	END DO
714
715	END SUBROUTINE start_inter_3d
716	!
717	!-------------------------------------------------------------------------------
718
719
720	!-------------------------------------------------------------------------------
721	!
722	SUBROUTINE interp_startvar(nam,ibar,ibeg,lon,lat,vari,lon1,lat1,lon2,lat2,varo)
723	!
724	!-------------------------------------------------------------------------------
725	USE inter_barxy_m, ONLY: inter_barxy
726	USE grid_atob_m, ONLY: grille_m
727	IMPLICIT NONE
728	!-------------------------------------------------------------------------------
729	! Arguments:
730	CHARACTER(LEN=*), INTENT(IN) :: nam
731	LOGICAL, INTENT(IN) :: ibar, ibeg
732	REAL, INTENT(IN) :: lon(:), lat(:) ! dim (ii) (jj)
733	REAL, INTENT(IN) :: vari(:,:) ! dim (ii,jj)
734	REAL, INTENT(IN) :: lon1(:), lat1(:) ! dim (i1) (j1)
735	REAL, INTENT(IN) :: lon2(:), lat2(:) ! dim (i1) (j2)
736	REAL, INTENT(OUT) :: varo(:,:) ! dim (i1) (j1)
737	!-------------------------------------------------------------------------------
738	! Local variables:
739	CHARACTER(LEN=256) :: modname="interp_startvar"
740	INTEGER :: ii, jj, i1, j1, j2
741	REAL, ALLOCATABLE :: vtmp(:,:)
742	!-------------------------------------------------------------------------------
743	ii=assert_eq(SIZE(lon), SIZE(vari,1),TRIM(modname)//" ii")
744	jj=assert_eq(SIZE(lat), SIZE(vari,2),TRIM(modname)//" jj")
745	i1=assert_eq(SIZE(lon1),SIZE(lon2),SIZE(varo,1),TRIM(modname)//" i1")
746	j1=assert_eq(SIZE(lat1), SIZE(varo,2),TRIM(modname)//" j1")
747	j2=SIZE(lat2)
748	ALLOCATE(vtmp(i1-1,j1))
749	IF(ibar) THEN
750	IF(ibeg.AND.prt_level>1) THEN
751	WRITE(lunout,*)"---------------------------------------------------------"
752	WRITE(lunout,*)"$$$ Interpolation barycentrique pour "//TRIM(nam)//" $$$"
753	WRITE(lunout,*)"---------------------------------------------------------"
754	END IF
755	CALL inter_barxy(lon, lat(:jj-1), vari, lon2(:i1-1), lat2, vtmp)
756	ELSE
757	CALL grille_m (lon, lat, vari, lon1, lat1, vtmp)
758	END IF
759	CALL gr_int_dyn(vtmp, varo, i1-1, j1)
760
761	END SUBROUTINE interp_startvar
762	!
763	!-------------------------------------------------------------------------------
764
765	!#endif
766	! of ifdef CPP_EARTH
767
768	END MODULE etat0dyn
769	!
770	!*******************************************************************************

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

Download in other formats:

Original Format