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