Context Navigation

grid_noro_m.F90 @ 2645

Last change on this file since 2645 was 2576, checked in by dcugnet, 8 years ago

Use "traditional" sst/sic files without additional records from contiguous years,
but check for previous and next files. Example for Amip Sea Surface Temperature:

amipbc_sst_1x1.nc must be present.
if amipbc_sst_1x1_m.nc is present (previous year), its last record is used

for interpolation ; otherwise, the last record of amipbc_sst_1x1.nc is used.

if amipbc_sst_1x1_p.nc is present (next year), its first record is used

for interpolation ; otherwise, the first record of amipbc_sst_1x1.nc is used.
If both contiguous fiels are missing, the interpolation then implies periodicity.
Note: grid_noro0 has been moved from limit_netcd to grid_noro_m where it belongs.

File size: 20.6 KB

Rev	Line
[2293]	1	MODULE grid_noro_m
	2	!
	3	!*******************************************************************************
	4
[2576]	5	USE print_control_mod, ONLY: lunout
	6	USE assert_eq_m, ONLY: assert_eq
[2293]	7	PRIVATE
[2576]	8	PUBLIC :: grid_noro, grid_noro0
[2293]	9
	10
	11	CONTAINS
	12
	13
	14	!-------------------------------------------------------------------------------
	15	!
	16	SUBROUTINE grid_noro(xd,yd,zd,x,y,zphi,zmea,zstd,zsig,zgam,zthe,zpic,zval,mask)
	17	!
	18	!-------------------------------------------------------------------------------
	19	! Author: F. Lott (see also Z.X. Li, A. Harzallah et L. Fairhead)
	20	!-------------------------------------------------------------------------------
	21	! Purpose: Compute the Parameters of the SSO scheme as described in LOTT &MILLER
	22	! (1997) and LOTT(1999).
	23	!-------------------------------------------------------------------------------
	24	! Comments:
	25	! * Target points are on a rectangular grid:
	26	! iim+1 latitudes including North and South Poles;
	27	! jjm+1 longitudes, with periodicity jjm+1=1.
	28	! * At the poles, the fields value is repeated jjm+1 time.
	29	! * The parameters a,b,c,d represent the limits of the target gridpoint region.
	30	! The means over this region are calculated from USN data, ponderated by a
	31	! weight proportional to the surface occupated by the data inside the model
	32	! gridpoint area. In most circumstances, this weight is the ratio between the
	33	! surfaces of the USN gridpoint area and the model gridpoint area.
	34	!
	35	! (c)
	36	! ----d-----
	37	! \| . . . .\|
	38	! \| \|
	39	! (b)a . * . .b(a)
	40	! \| \|
	41	! \| . . . .\|
	42	! ----c-----
	43	! (d)
	44	! * Hard-coded US Navy dataset dimensions (imdp=2160 ; jmdp=1080) have been
	45	! removed (ALLOCATABLE used).
	46	! * iext (currently 10% of imdp) represents the margin to ensure output cells
	47	! on the edge are contained in input cells.
	48	!===============================================================================
[2576]	49	IMPLICIT NONE
[2293]	50	!-------------------------------------------------------------------------------
	51	! Arguments:
	52	REAL, INTENT(IN) :: xd(:), yd(:) !--- INPUT COORDINATES (imdp) (jmdp)
	53	REAL, INTENT(IN) :: zd(:,:) !--- INPUT FIELD (imdp,jmdp)
	54	REAL, INTENT(IN) :: x(:), y(:) !--- OUTPUT COORDINATES (imar+1) (jmar)
	55	REAL, INTENT(OUT) :: zphi(:,:) !--- GEOPOTENTIAL (imar+1,jmar)
	56	REAL, INTENT(OUT) :: zmea(:,:) !--- MEAN OROGRAPHY (imar+1,jmar)
	57	REAL, INTENT(OUT) :: zstd(:,:) !--- STANDARD DEVIATION (imar+1,jmar)
	58	REAL, INTENT(OUT) :: zsig(:,:) !--- SLOPE (imar+1,jmar)
	59	REAL, INTENT(OUT) :: zgam(:,:) !--- ANISOTROPY (imar+1,jmar)
	60	REAL, INTENT(OUT) :: zthe(:,:) !--- SMALL AXIS ORIENTATION (imar+1,jmar)
	61	REAL, INTENT(OUT) :: zpic(:,:) !--- MAXIMUM ALTITITUDE (imar+1,jmar)
	62	REAL, INTENT(OUT) :: zval(:,:) !--- MINIMUM ALTITITUDE (imar+1,jmar)
	63	REAL, INTENT(OUT) :: mask(:,:) !--- MASK (imar+1,jmar)
	64	!-------------------------------------------------------------------------------
	65	! Local variables:
	66	CHARACTER(LEN=256) :: modname="grid_noro"
	67	REAL, ALLOCATABLE :: xusn(:), yusn(:) ! dim (imdp+2*iext) (jmdp+2)
	68	REAL, ALLOCATABLE :: zusn(:,:) ! dim (imdp+2*iext,jmdp+2)
	69	! CORRELATIONS OF OROGRAPHY GRADIENT ! dim (imar+1,jmar)
	70	REAL, ALLOCATABLE :: ztz(:,:), zxtzx(:,:), zytzy(:,:), zxtzy(:,:), weight(:,:)
	71	! CORRELATIONS OF USN OROGRAPHY GRADIENTS ! dim (imar+2*iext,jmdp+2)
	72	REAL, ALLOCATABLE :: zxtzxusn(:,:), zytzyusn(:,:), zxtzyusn(:,:)
	73	REAL, ALLOCATABLE :: mask_tmp(:,:), zmea0(:,:) ! dim (imar+1,jmar)
	74	REAL, ALLOCATABLE :: num_tot(:,:), num_lan(:,:) ! dim (imax,jmax)
	75	REAL, ALLOCATABLE :: a(:), b(:) ! dim (imax)
	76	REAL, ALLOCATABLE :: c(:), d(:) ! dim (jmax)
	77	LOGICAL :: masque_lu
	78	INTEGER :: i, ii, imdp, imar, iext
	79	INTEGER :: j, jj, jmdp, jmar, nn
	80	REAL :: xpi, zdeltax, zlenx, weighx, xincr, zmeanor0
	81	REAL :: rad, zdeltay, zleny, weighy, masque, zmeasud0
	82	REAL :: zbordnor, zmeanor, zstdnor, zsignor, zweinor, zpicnor, zvalnor
	83	REAL :: zbordsud, zmeasud, zstdsud, zsigsud, zweisud, zpicsud, zvalsud
	84	REAL :: zbordest, zbordoue, xk, xl, xm, xp, xq, xw
	85	!-------------------------------------------------------------------------------
	86	imdp=assert_eq(SIZE(xd),SIZE(zd,1),TRIM(modname)//" imdp")
	87	jmdp=assert_eq(SIZE(yd),SIZE(zd,2),TRIM(modname)//" jmdp")
	88	imar=assert_eq([SIZE(x),SIZE(zphi,1),SIZE(zmea,1),SIZE(zstd,1),SIZE(zsig,1), &
	89	SIZE(zgam,1),SIZE(zthe,1),SIZE(zpic,1),SIZE(zval,1), &
	90	SIZE(mask,1)],TRIM(modname)//" imar")-1
	91	jmar=assert_eq([SIZE(y),SIZE(zphi,2),SIZE(zmea,2),SIZE(zstd,2),SIZE(zsig,2), &
	92	SIZE(zgam,2),SIZE(zthe,2),SIZE(zpic,2),SIZE(zval,2), &
	93	SIZE(mask,2)],TRIM(modname)//" jmar")
[2311]	94	! IF(imar/=iim) CALL abort_physic(TRIM(modname),'imar/=iim' ,1)
	95	! IF(jmar/=jjm+1) CALL abort_physic(TRIM(modname),'jmar/=jjm+1',1)
[2293]	96	iext=imdp/10 !--- OK up to 36 degrees cell
	97	xpi = ACOS(-1.)
	98	rad = 6371229.
	99	zdeltay=2.xpi/REAL(jmdp)rad
	100	WRITE(lunout,)" Orography parameters at sub-cell scale *"
	101
	102	!--- ARE WE USING A READ MASK ?
	103	masque_lu=ANY(mask/=-99999.); IF(.NOT.masque_lu) mask=0.0
	104	WRITE(lunout,*)'Masque lu: ',masque_lu
	105
	106	!--- EXTENSION OF THE INPUT DATABASE TO PROCEED COMPUTATIONS AT BOUNDARIES:
	107	ALLOCATE(xusn(imdp+2*iext))
	108	xusn(1 +iext:imdp +iext)=xd(:)
	109	xusn(1 : iext)=xd(1+imdp-iext:imdp)-2.*xpi
	110	xusn(1+imdp+iext:imdp+2iext)=xd(1 :iext)+2.xpi
	111
	112	ALLOCATE(yusn(jmdp+2))
	113	yusn(1 )=yd(1) +(yd(1) -yd(2))
	114	yusn(2:jmdp+1)=yd(:)
	115	yusn( jmdp+2)=yd(jmdp)+(yd(jmdp)-yd(jmdp-1))
	116
	117	ALLOCATE(zusn(imdp+2*iext,jmdp+2))
	118	zusn(1 +iext:imdp +iext,2:jmdp+1)=zd (: , :)
	119	zusn(1 : iext,2:jmdp+1)=zd (imdp-iext+1:imdp , :)
	120	zusn(1+imdp +iext:imdp+2*iext,2:jmdp+1)=zd (1:iext , :)
	121	zusn(1 :imdp/2+iext, 1)=zusn(1+imdp/2:imdp +iext, 2)
	122	zusn(1+imdp/2+iext:imdp+2*iext, 1)=zusn(1 :imdp/2+iext, 2)
	123	zusn(1 :imdp/2+iext, jmdp+2)=zusn(1+imdp/2:imdp +iext,jmdp+1)
	124	zusn(1+imdp/2+iext:imdp+2*iext, jmdp+2)=zusn(1 :imdp/2+iext,jmdp+1)
	125
	126	!--- COMPUTE LIMITS OF MODEL GRIDPOINT AREA (REGULAR GRID)
	127	ALLOCATE(a(imar+1),b(imar+1))
	128	b(1:imar)=(x(1:imar )+ x(2:imar+1))/2.0
	129	b(imar+1)= x( imar+1)+(x( imar+1)-x(imar))/2.0
	130	a(1)=x(1)-(x(2)-x(1))/2.0
	131	a(2:imar+1)= b(1:imar)
	132
	133	ALLOCATE(c(jmar),d(jmar))
	134	d(1:jmar-1)=(y(1:jmar-1)+ y(2:jmar))/2.0
	135	d( jmar )= y( jmar )+(y( jmar)-y(jmar-1))/2.0
	136	c(1)=y(1)-(y(2)-y(1))/2.0
	137	c(2:jmar)=d(1:jmar-1)
	138
	139	!--- INITIALIZATIONS:
	140	ALLOCATE(weight(imar+1,jmar)); weight(:,:)= 0.0
	141	ALLOCATE(zxtzx (imar+1,jmar)); zxtzx (:,:)= 0.0
	142	ALLOCATE(zytzy (imar+1,jmar)); zytzy (:,:)= 0.0
	143	ALLOCATE(zxtzy (imar+1,jmar)); zxtzy (:,:)= 0.0
	144	ALLOCATE(ztz (imar+1,jmar)); ztz (:,:)= 0.0
	145	zmea(:,:)= 0.0
	146	zpic(:,:)=-1.E+10
	147	zval(:,:)= 1.E+10
	148
	149	!--- COMPUTE SLOPES CORRELATIONS ON USN GRID
	150	! CORRELATIONS OF USN OROGRAPHY GRADIENTS ! dim (imdp+2*iext,jmdp+2)
	151	ALLOCATE(zytzyusn(imdp+2*iext,jmdp+2)); zytzyusn(:,:)=0.0
	152	ALLOCATE(zxtzxusn(imdp+2*iext,jmdp+2)); zxtzxusn(:,:)=0.0
	153	ALLOCATE(zxtzyusn(imdp+2*iext,jmdp+2)); zxtzyusn(:,:)=0.0
	154	DO j = 2, jmdp+1
	155	zdeltax=zdeltay*cos(yusn(j))
	156	DO i = 2, imdp+2*iext-1
	157	zytzyusn(i,j)=(zusn(i,j+1)-zusn(i,j-1))2/zdeltay2
	158	zxtzxusn(i,j)=(zusn(i+1,j)-zusn(i-1,j))2/zdeltax2
	159	zxtzyusn(i,j)=(zusn(i,j+1)-zusn(i,j-1)) /zdeltay &
	160	& *(zusn(i+1,j)-zusn(i-1,j)) /zdeltax
	161	END DO
	162	END DO
	163
	164	!--- SUMMATION OVER GRIDPOINT AREA
	165	zleny=xpi/REAL(jmdp)*rad
	166	xincr=xpi/REAL(jmdp)/2.
	167	ALLOCATE(num_tot(imar+1,jmar)); num_tot(:,:)=0.
	168	ALLOCATE(num_lan(imar+1,jmar)); num_lan(:,:)=0.
	169	DO ii = 1, imar+1
	170	DO jj = 1, jmar
	171	DO j = 2,jmdp+1
	172	zlenx =zleny *COS(yusn(j))
	173	zdeltax=zdeltay*COS(yusn(j))
	174	zbordnor=(xincr+c(jj)-yusn(j))*rad
	175	zbordsud=(xincr-d(jj)+yusn(j))*rad
	176	weighy=AMAX1(0.,AMIN1(zbordnor,zbordsud,zleny))
	177	IF(weighy==0.) CYCLE
	178	DO i = 2, imdp+2*iext-1
	179	zbordest=(xusn(i)-a(ii)+xincr)radCOS(yusn(j))
	180	zbordoue=(b(ii)+xincr-xusn(i))radCOS(yusn(j))
	181	weighx=AMAX1(0.,AMIN1(zbordest,zbordoue,zlenx))
	182	IF(weighx==0.) CYCLE
	183	num_tot(ii,jj)=num_tot(ii,jj)+1.0
	184	IF(zusn(i,j)>=1.)num_lan(ii,jj)=num_lan(ii,jj)+1.0
	185	weight(ii,jj)=weight(ii,jj)+weighx*weighy
	186	zxtzx(ii,jj)=zxtzx(ii,jj)+zxtzxusn(i,j)weighxweighy
	187	zytzy(ii,jj)=zytzy(ii,jj)+zytzyusn(i,j)weighxweighy
	188	zxtzy(ii,jj)=zxtzy(ii,jj)+zxtzyusn(i,j)weighxweighy
	189	ztz (ii,jj)= ztz(ii,jj)+zusn(i,j)zusn(i,j)weighx*weighy
	190	zmea (ii,jj)= zmea(ii,jj)+zusn(i,j)weighxweighy !--- MEAN
	191	zpic (ii,jj)=AMAX1(zpic(ii,jj),zusn(i,j)) !--- PEAKS
	192	zval (ii,jj)=AMIN1(zval(ii,jj),zusn(i,j)) !--- VALLEYS
	193	END DO
	194	END DO
	195	END DO
	196	END DO
	197
	198	!--- COMPUTE PARAMETERS NEEDED BY LOTT & MILLER (1997) AND LOTT (1999) SSO SCHEME
	199	IF(.NOT.masque_lu) THEN
	200	WHERE(weight(:,1:jmar-1)/=0.0) mask=num_lan(:,:)/num_tot(:,:)
	201	END IF
	202	nn=COUNT(weight(:,1:jmar-1)==0.0)
	203	IF(nn/=0) WRITE(lunout,*)'Problem with weight ; vanishing occurrences: ',nn
	204	WHERE(weight(:,:)/=0.0)
	205	zmea (:,:)=zmea (:,:)/weight(:,:)
	206	zxtzx(:,:)=zxtzx(:,:)/weight(:,:)
	207	zytzy(:,:)=zytzy(:,:)/weight(:,:)
	208	zxtzy(:,:)=zxtzy(:,:)/weight(:,:)
	209	ztz (:,:)=ztz (:,:)/weight(:,:)
	210	zstd (:,:)=ztz (:,:)-zmea(:,:)**2
	211	END WHERE
	212	WHERE(zstd(:,:)<0) zstd(:,:)=0.
	213	zstd (:,:)=SQRT(zstd(:,:))
	214
	215	!--- CORRECT VALUES OF HORIZONTAL SLOPE NEAR THE POLES:
	216	zxtzx(:, 1)=zxtzx(:, 2)
	217	zxtzx(:,jmar)=zxtzx(:,jmar-1)
	218	zxtzy(:, 1)=zxtzy(:, 2)
	219	zxtzy(:,jmar)=zxtzy(:,jmar-1)
	220	zytzy(:, 1)=zytzy(:, 2)
	221	zytzy(:,jmar)=zytzy(:,jmar-1)
	222
	223	!=== FILTERS TO SMOOTH OUT FIELDS FOR INPUT INTO SSO SCHEME.
	224	!--- FIRST FILTER, MOVING AVERAGE OVER 9 POINTS.
	225	!-------------------------------------------------------------------------------
	226	ALLOCATE(zmea0(imar+1,jmar))
	227	zmea0(:,:)=zmea(:,:) ! GK211005 (CG) UNSMOOTHED TOPO
	228	CALL MVA9(zmea); CALL MVA9(zstd); CALL MVA9(zpic); CALL MVA9(zval)
	229	CALL MVA9(zxtzx); CALL MVA9(zxtzy); CALL MVA9(zytzy)
	230
	231	!--- MASK BASED ON GROUND MAXIMUM, 10% THRESHOLD. (SURFACE PARAMS MEANINGLESS)
	232	ALLOCATE(mask_tmp(imar+1,jmar)); mask_tmp(:,:)=0.0
	233	WHERE(mask>=0.1) mask_tmp = 1.
	234	WHERE(weight(:,:)/=0.0)
	235	! zphi (:,:)= mask_tmp(:,:)*zmea (:,:) ! GK211005 (CG) not necessarly smoothed
	236	zphi (:,:)= mask_tmp(:,:)*zmea0(:,:)
	237	zmea0(:,:)= mask_tmp(:,:)*zmea0(:,:)
	238	zmea (:,:)= mask_tmp(:,:)*zmea (:,:)
	239	zpic (:,:)= mask_tmp(:,:)*zpic (:,:)
	240	zval (:,:)= mask_tmp(:,:)*zval (:,:)
	241	zstd (:,:)= mask_tmp(:,:)*zstd (:,:)
	242	END WHERE
	243	DO ii = 1, imar
	244	DO jj = 1, jmar
	245	IF (weight(ii,jj)/=0.0) THEN
	246	!--- Coefficients K, L et M:
	247	xk=(zxtzx(ii,jj)+zytzy(ii,jj))/2.
	248	xl=(zxtzx(ii,jj)-zytzy(ii,jj))/2.
	249	xm=zxtzy(ii,jj)
	250	xp=xk-SQRT(xl2+xm2)
	251	xq=xk+SQRT(xl2+xm2)
	252	xw=1.e-8
	253	IF(xp<=xw) xp=0.
	254	IF(xq<=xw) xq=xw
	255	IF(ABS(xm)<=xw) xm=xw*SIGN(1.,xm)
	256	!--- SLOPE
	257	zsig(ii,jj)=SQRT(xq)*mask_tmp(ii,jj)
	258	!---ISOTROPY
	259	zgam(ii,jj)=xp/xq*mask_tmp(ii,jj)
	260	!--- THETA ANGLE
	261	zthe(ii,jj)=57.29577951ATAN2(xm,xl)/2.mask_tmp(ii,jj)
	262	END IF
	263	END DO
	264	END DO
	265	WRITE(lunout,*)' MEAN ORO:' ,MAXVAL(zmea)
	266	WRITE(lunout,*)' ST. DEV.:' ,MAXVAL(zstd)
	267	WRITE(lunout,*)' PENTE:' ,MAXVAL(zsig)
	268	WRITE(lunout,*)' ANISOTROP:',MAXVAL(zgam)
	269	WRITE(lunout,*)' ANGLE:' ,MINVAL(zthe),MAXVAL(zthe)
	270	WRITE(lunout,*)' pic:' ,MAXVAL(zpic)
	271	WRITE(lunout,*)' val:' ,MAXVAL(zval)
	272
	273	!--- Values at poles
	274	zmea0(imar+1,:)=zmea0(1,:)
	275	zmea (imar+1,:)=zmea (1,:)
	276	zphi (imar+1,:)=zphi (1,:)
	277	zpic (imar+1,:)=zpic (1,:)
	278	zval (imar+1,:)=zval (1,:)
	279	zstd (imar+1,:)=zstd (1,:)
	280	zsig (imar+1,:)=zsig (1,:)
	281	zgam (imar+1,:)=zgam (1,:)
	282	zthe (imar+1,:)=zthe (1,:)
	283
	284	zweinor =SUM(weight(1:imar, 1),DIM=1)
	285	zweisud =SUM(weight(1:imar,jmar),DIM=1)
	286	zmeanor0=SUM(weight(1:imar, 1)*zmea0(1:imar, 1),DIM=1)
	287	zmeasud0=SUM(weight(1:imar,jmar)*zmea0(1:imar,jmar),DIM=1)
	288	zmeanor =SUM(weight(1:imar, 1)*zmea (1:imar, 1),DIM=1)
	289	zmeasud =SUM(weight(1:imar,jmar)*zmea (1:imar,jmar),DIM=1)
	290	zstdnor =SUM(weight(1:imar, 1)*zstd (1:imar, 1),DIM=1)
	291	zstdsud =SUM(weight(1:imar,jmar)*zstd (1:imar,jmar),DIM=1)
	292	zsignor =SUM(weight(1:imar, 1)*zsig (1:imar, 1),DIM=1)
	293	zsigsud =SUM(weight(1:imar,jmar)*zsig (1:imar,jmar),DIM=1)
	294	zpicnor =SUM(weight(1:imar, 1)*zpic (1:imar, 1),DIM=1)
	295	zpicsud =SUM(weight(1:imar,jmar)*zpic (1:imar,jmar),DIM=1)
	296	zvalnor =SUM(weight(1:imar, 1)*zval (1:imar, 1),DIM=1)
	297	zvalsud =SUM(weight(1:imar,jmar)*zval (1:imar,jmar),DIM=1)
	298
	299	zmea(:,1)=zmeanor /zweinor; zmea(:,jmar)=zmeasud /zweisud
	300	! zphi(:,1)=zmeanor0/zweinor; zphi(:,jmar)=zmeasud0/zweisud TO COMMIT
	301	zphi(:,1)=zmeanor /zweinor; zphi(:,jmar)=zmeasud /zweisud
	302	zpic(:,1)=zpicnor /zweinor; zpic(:,jmar)=zpicsud /zweisud
	303	zval(:,1)=zvalnor /zweinor; zval(:,jmar)=zvalsud /zweisud
	304	zstd(:,1)=zstdnor /zweinor; zstd(:,jmar)=zstdsud /zweisud
	305	zsig(:,1)=zsignor /zweinor; zsig(:,jmar)=zsigsud /zweisud
	306	zgam(:,1)=1.; zgam(:,jmar)=1.
	307	zthe(:,1)=0.; zthe(:,jmar)=0.
	308
	309	END SUBROUTINE grid_noro
	310	!
	311	!-------------------------------------------------------------------------------
	312
	313
	314	!-------------------------------------------------------------------------------
	315	!
[2576]	316	SUBROUTINE grid_noro0(xd,yd,zd,x,y,zphi,mask)
	317	!
	318	!===============================================================================
	319	! Purpose: Extracted from grid_noro to provide geopotential height for dynamics
	320	! without any call to physics subroutines.
	321	!===============================================================================
	322	IMPLICIT NONE
	323	!-------------------------------------------------------------------------------
	324	! Arguments:
	325	REAL, INTENT(IN) :: xd(:), yd(:) !--- INPUT COORDINATES (imdp) (jmdp)
	326	REAL, INTENT(IN) :: zd(:,:) !--- INPUT FIELD (imdp,jmdp)
	327	REAL, INTENT(IN) :: x(:), y(:) !--- OUTPUT COORDINATES (imar+1) (jmar)
	328	REAL, INTENT(OUT) :: zphi(:,:) !--- GEOPOTENTIAL (imar+1,jmar)
	329	REAL, INTENT(INOUT):: mask(:,:) !--- MASK (imar+1,jmar)
	330	!-------------------------------------------------------------------------------
	331	! Local variables:
	332	CHARACTER(LEN=256) :: modname="grid_noro0"
	333	REAL, ALLOCATABLE :: xusn(:), yusn(:) ! dim (imdp+2*iext) (jmdp+2)
	334	REAL, ALLOCATABLE :: zusn(:,:) ! dim (imdp+2*iext,jmdp+2)
	335	REAL, ALLOCATABLE :: weight(:,:) ! dim (imar+1,jmar)
	336	REAL, ALLOCATABLE :: mask_tmp(:,:), zmea(:,:) ! dim (imar+1,jmar)
	337	REAL, ALLOCATABLE :: num_tot(:,:), num_lan(:,:) ! dim (imax,jmax)
	338	REAL, ALLOCATABLE :: a(:), b(:) ! dim (imax)
	339	REAL, ALLOCATABLE :: c(:), d(:) ! dim (jmax)
	340	LOGICAL :: masque_lu
	341	INTEGER :: i, ii, imdp, imar, iext
	342	INTEGER :: j, jj, jmdp, jmar, nn
	343	REAL :: xpi, zlenx, weighx, xincr, zbordnor, zmeanor, zweinor, zbordest
	344	REAL :: rad, zleny, weighy, masque, zbordsud, zmeasud, zweisud, zbordoue
	345	!-------------------------------------------------------------------------------
	346	imdp=assert_eq(SIZE(xd),SIZE(zd,1),TRIM(modname)//" imdp")
	347	jmdp=assert_eq(SIZE(yd),SIZE(zd,2),TRIM(modname)//" jmdp")
	348	imar=assert_eq(SIZE(x),SIZE(zphi,1),SIZE(mask,1),TRIM(modname)//" imar")-1
	349	jmar=assert_eq(SIZE(y),SIZE(zphi,2),SIZE(mask,2),TRIM(modname)//" jmar")
	350	! IF(imar/=iim) CALL abort_gcm(TRIM(modname),'imar/=iim' ,1)
	351	! IF(jmar/=jjm+1) CALL abort_gcm(TRIM(modname),'jmar/=jjm+1',1)
	352	iext=imdp/10
	353	xpi = ACOS(-1.)
	354	rad = 6371229.
	355
	356	!--- ARE WE USING A READ MASK ?
	357	masque_lu=ANY(mask/=-99999.); IF(.NOT.masque_lu) mask=0.0
	358	WRITE(lunout,*)'Masque lu: ',masque_lu
	359
	360	!--- EXTENSION OF THE INPUT DATABASE TO PROCEED COMPUTATIONS AT BOUNDARIES:
	361	ALLOCATE(xusn(imdp+2*iext))
	362	xusn(1 +iext:imdp +iext)=xd(:)
	363	xusn(1 : iext)=xd(1+imdp-iext:imdp)-2.*xpi
	364	xusn(1+imdp+iext:imdp+2iext)=xd(1 :iext)+2.xpi
	365
	366	ALLOCATE(yusn(jmdp+2))
	367	yusn(1 )=yd(1) +(yd(1) -yd(2))
	368	yusn(2:jmdp+1)=yd(:)
	369	yusn( jmdp+2)=yd(jmdp)+(yd(jmdp)-yd(jmdp-1))
	370
	371	ALLOCATE(zusn(imdp+2*iext,jmdp+2))
	372	zusn(1 +iext:imdp +iext,2:jmdp+1)=zd (: , :)
	373	zusn(1 : iext,2:jmdp+1)=zd (imdp-iext+1:imdp , :)
	374	zusn(1+imdp +iext:imdp+2*iext,2:jmdp+1)=zd (1:iext , :)
	375	zusn(1 :imdp/2+iext, 1)=zusn(1+imdp/2:imdp +iext, 2)
	376	zusn(1+imdp/2+iext:imdp+2*iext, 1)=zusn(1 :imdp/2+iext, 2)
	377	zusn(1 :imdp/2+iext, jmdp+2)=zusn(1+imdp/2:imdp +iext,jmdp+1)
	378	zusn(1+imdp/2+iext:imdp+2*iext, jmdp+2)=zusn(1 :imdp/2+iext,jmdp+1)
	379
	380	!--- COMPUTE LIMITS OF MODEL GRIDPOINT AREA (REGULAR GRID)
	381	ALLOCATE(a(imar+1),b(imar+1))
	382	b(1:imar)=(x(1:imar )+ x(2:imar+1))/2.0
	383	b(imar+1)= x( imar+1)+(x( imar+1)-x(imar))/2.0
	384	a(1)=x(1)-(x(2)-x(1))/2.0
	385	a(2:imar+1)= b(1:imar)
	386
	387	ALLOCATE(c(jmar),d(jmar))
	388	d(1:jmar-1)=(y(1:jmar-1)+ y(2:jmar))/2.0
	389	d( jmar )= y( jmar )+(y( jmar)-y(jmar-1))/2.0
	390	c(1)=y(1)-(y(2)-y(1))/2.0
	391	c(2:jmar)=d(1:jmar-1)
	392
	393	!--- INITIALIZATIONS:
	394	ALLOCATE(weight(imar+1,jmar)); weight(:,:)= 0.0
	395	ALLOCATE(zmea (imar+1,jmar)); zmea (:,:)= 0.0
	396
	397	!--- SUMMATION OVER GRIDPOINT AREA
	398	zleny=xpi/REAL(jmdp)*rad
	399	xincr=xpi/REAL(jmdp)/2.
	400	ALLOCATE(num_tot(imar+1,jmar)); num_tot(:,:)=0.
	401	ALLOCATE(num_lan(imar+1,jmar)); num_lan(:,:)=0.
	402	DO ii = 1, imar+1
	403	DO jj = 1, jmar
	404	DO j = 2,jmdp+1
	405	zlenx =zleny *COS(yusn(j))
	406	zbordnor=(xincr+c(jj)-yusn(j))*rad
	407	zbordsud=(xincr-d(jj)+yusn(j))*rad
	408	weighy=AMAX1(0.,AMIN1(zbordnor,zbordsud,zleny))
	409	IF(weighy/=0) THEN
	410	DO i = 2, imdp+2*iext-1
	411	zbordest=(xusn(i)-a(ii)+xincr)radCOS(yusn(j))
	412	zbordoue=(b(ii)+xincr-xusn(i))radCOS(yusn(j))
	413	weighx=AMAX1(0.,AMIN1(zbordest,zbordoue,zlenx))
	414	IF(weighx/=0)THEN
	415	num_tot(ii,jj)=num_tot(ii,jj)+1.0
	416	IF(zusn(i,j)>=1.)num_lan(ii,jj)=num_lan(ii,jj)+1.0
	417	weight(ii,jj)=weight(ii,jj)+weighx*weighy
	418	zmea (ii,jj)=zmea (ii,jj)+zusn(i,j)weighxweighy !--- MEAN
	419	END IF
	420	END DO
	421	END IF
	422	END DO
	423	END DO
	424	END DO
	425
	426	!--- COMPUTE PARAMETERS NEEDED BY LOTT & MILLER (1997) AND LOTT (1999) SSO SCHEME
	427	IF(.NOT.masque_lu) THEN
	428	WHERE(weight(:,1:jmar-1)/=0.0) mask=num_lan(:,:)/num_tot(:,:)
	429	END IF
	430	nn=COUNT(weight(:,1:jmar-1)==0.0)
	431	IF(nn/=0) WRITE(lunout,*)'Problem with weight ; vanishing occurrences: ',nn
	432	WHERE(weight/=0.0) zmea(:,:)=zmea(:,:)/weight(:,:)
	433
	434	!--- MASK BASED ON GROUND MAXIMUM, 10% THRESHOLD (<10%: SURF PARAMS MEANINGLESS)
	435	ALLOCATE(mask_tmp(imar+1,jmar)); mask_tmp(:,:)=0.0
	436	WHERE(mask>=0.1) mask_tmp = 1.
	437	WHERE(weight(:,:)/=0.0)
	438	zphi(:,:)=mask_tmp(:,:)*zmea(:,:)
	439	zmea(:,:)=mask_tmp(:,:)*zmea(:,:)
	440	END WHERE
	441	WRITE(lunout,*)' MEAN ORO:' ,MAXVAL(zmea)
	442
	443	!--- Values at poles
	444	zphi(imar+1,:)=zphi(1,:)
	445
	446	zweinor=SUM(weight(1:imar, 1),DIM=1)
	447	zweisud=SUM(weight(1:imar,jmar),DIM=1)
	448	zmeanor=SUM(weight(1:imar, 1)*zmea(1:imar, 1),DIM=1)
	449	zmeasud=SUM(weight(1:imar,jmar)*zmea(1:imar,jmar),DIM=1)
	450	zphi(:,1)=zmeanor/zweinor; zphi(:,jmar)=zmeasud/zweisud
	451
	452	END SUBROUTINE grid_noro0
	453	!
	454	!-------------------------------------------------------------------------------
	455
	456
	457	!-------------------------------------------------------------------------------
	458	!
[2293]	459	SUBROUTINE MVA9(x)
	460	!
	461	!-------------------------------------------------------------------------------
	462	IMPLICIT NONE
	463	! MAKE A MOVING AVERAGE OVER 9 GRIDPOINTS OF THE X FIELDS
	464	!-------------------------------------------------------------------------------
	465	! Arguments:
	466	REAL, INTENT(INOUT) :: x(:,:)
	467	!-------------------------------------------------------------------------------
	468	! Local variables:
	469	REAL :: xf(SIZE(x,DIM=1),SIZE(x,DIM=2)), WEIGHTpb(-1:1,-1:1)
	470	INTEGER :: i, j, imar, jmar
	471	!-------------------------------------------------------------------------------
	472	WEIGHTpb=RESHAPE([((1./REAL((1+i*2)(1+j**2)),i=-1,1),j=-1,1)],SHAPE=[3,3])
	473	WEIGHTpb=WEIGHTpb/SUM(WEIGHTpb)
	474	imar=SIZE(X,DIM=1); jmar=SIZE(X,DIM=2)
	475	DO j=2,jmar-1
	476	DO i=2,imar-1
	477	xf(i,j)=SUM(x(i-1:i+1,j-1:j+1)*WEIGHTpb(:,:))
	478	END DO
	479	END DO
	480	DO j=2,jmar-1
	481	xf(1,j)=SUM(x(imar-1,j-1:j+1)*WEIGHTpb(-1,:))
	482	xf(1,j)=xf(1,j)+SUM(x(1:2,j-1:j+1)*WEIGHTpb(0:1,-1:1))
	483	xf(imar,j)=xf(1,j)
	484	END DO
	485	xf(:, 1)=xf(:, 2)
	486	xf(:,jmar)=xf(:,jmar-1)
	487	x(:,:)=xf(:,:)
	488
	489	END SUBROUTINE MVA9
	490	!
	491	!-------------------------------------------------------------------------------
	492
	493
	494	END MODULE grid_noro_m
	495

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

Download in other formats:

Original Format