| 1 | MODULE grid_noro_m |
|---|
| 2 | ! |
|---|
| 3 | !******************************************************************************* |
|---|
| 4 | |
|---|
| 5 | USE print_control_mod, ONLY: lunout |
|---|
| 6 | USE assert_eq_m, ONLY: assert_eq |
|---|
| 7 | PRIVATE |
|---|
| 8 | PUBLIC :: grid_noro, grid_noro0 |
|---|
| 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 | !=============================================================================== |
|---|
| 49 | IMPLICIT NONE |
|---|
| 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") |
|---|
| 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) |
|---|
| 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+2*iext)=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/zdeltay**2 |
|---|
| 158 | zxtzxusn(i,j)=(zusn(i+1,j)-zusn(i-1,j))**2/zdeltax**2 |
|---|
| 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)*rad*COS(yusn(j)) |
|---|
| 180 | zbordoue=(b(ii)+xincr-xusn(i))*rad*COS(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)*weighx*weighy |
|---|
| 187 | zytzy(ii,jj)=zytzy(ii,jj)+zytzyusn(i,j)*weighx*weighy |
|---|
| 188 | zxtzy(ii,jj)=zxtzy(ii,jj)+zxtzyusn(i,j)*weighx*weighy |
|---|
| 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)*weighx*weighy !--- 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(xl**2+xm**2) |
|---|
| 251 | xq=xk+SQRT(xl**2+xm**2) |
|---|
| 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.29577951*ATAN2(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 | ! |
|---|
| 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+2*iext)=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)*rad*COS(yusn(j)) |
|---|
| 412 | zbordoue=(b(ii)+xincr-xusn(i))*rad*COS(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)*weighx*weighy !--- 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 | ! |
|---|
| 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 | |
|---|
