Context Navigation

etat0dyn_netcdf.F90 @ 2315

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

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

Download in other formats:

Original Format