[2293] | 1 | !******************************************************************************* |
---|
| 2 | ! |
---|
| 3 | MODULE grid_atob_m |
---|
| 4 | ! |
---|
| 5 | !******************************************************************************* |
---|
| 6 | |
---|
| 7 | USE assert_eq_m, ONLY: assert_eq |
---|
| 8 | |
---|
| 9 | PRIVATE |
---|
| 10 | PUBLIC :: grille_m, rugosite, sea_ice, rugsoro |
---|
| 11 | |
---|
| 12 | CONTAINS |
---|
| 13 | |
---|
| 14 | !------------------------------------------------------------------------------- |
---|
| 15 | ! |
---|
| 16 | SUBROUTINE fine2coarse(x_i, y_i, x_o, y_o, d_o1, d_i, msk, d_o2) |
---|
| 17 | ! |
---|
| 18 | !------------------------------------------------------------------------------- |
---|
| 19 | IMPLICIT NONE |
---|
| 20 | !------------------------------------------------------------------------------- |
---|
| 21 | ! Arguments: |
---|
[2299] | 22 | REAL, INTENT(IN) :: x_i(:), y_i(:) !-- INPUT X&Y COOR. (mi)(ni) |
---|
[2361] | 23 | REAL, INTENT(IN) :: x_o(:), y_o(:) !-- OUTPUT X&Y COOR. (mo)(no) |
---|
[2299] | 24 | REAL, INTENT(OUT) :: d_o1(:,:) !-- OUTPUT FIELD (mo,no) |
---|
| 25 | REAL, OPTIONAL, INTENT(IN) :: d_i (:,:) !-- INPUT FIELD (mi,ni) |
---|
[2361] | 26 | LOGICAL, OPTIONAL, INTENT(IN) :: msk (:,:) !-- MASK (mo,no) |
---|
[2299] | 27 | REAL, OPTIONAL, INTENT(OUT) :: d_o2(:,:) !-- OUTPUT FOR d_i^2 (mo,no) |
---|
[2293] | 28 | !------------------------------------------------------------------------------- |
---|
| 29 | ! Local variables: |
---|
| 30 | CHARACTER(LEN=256) :: modname="fine2coarse" |
---|
[2361] | 31 | INTEGER :: mi, ni, ii, ji, mo, no, io, jo, nr(2), m1, n1, m2, n2, mx, my, nn |
---|
| 32 | LOGICAL :: li, lo |
---|
| 33 | REAL :: inc |
---|
| 34 | INTEGER, DIMENSION(SIZE(x_o),SIZE(y_o)) :: num_tot |
---|
| 35 | LOGICAL, DIMENSION(SIZE(x_o),SIZE(y_o)) :: found, mask |
---|
| 36 | REAL, DIMENSION(SIZE(x_o),SIZE(y_o)) :: dist |
---|
| 37 | REAL, DIMENSION(SIZE(x_o)) :: a, b |
---|
| 38 | REAL, DIMENSION(SIZE(y_o)) :: c, d |
---|
| 39 | REAL, PARAMETER :: thresh=1.E-5 |
---|
[2293] | 40 | !------------------------------------------------------------------------------- |
---|
[2361] | 41 | mi=SIZE(x_i); ni=SIZE(y_i); mo=SIZE(x_o); no=SIZE(y_o) |
---|
| 42 | m1=m1; m2=mo; mx=mo; IF(PRESENT(msk)) mx=SIZE(msk,1) |
---|
| 43 | n1=ni; n2=no; my=no; IF(PRESENT(msk)) my=SIZE(msk,2) |
---|
[2293] | 44 | li=PRESENT(d_i ); IF(li) THEN; m1=SIZE(d_i ,1); n1=SIZE(d_i ,2); END IF |
---|
| 45 | lo=PRESENT(d_o2); IF(lo) THEN; m2=SIZE(d_o2,1); n2=SIZE(d_o2,2); END IF |
---|
[2361] | 46 | mi=assert_eq(mi,m1,TRIM(modname)//" mi") |
---|
| 47 | ni=assert_eq(ni,n1,TRIM(modname)//" ni") |
---|
| 48 | mo=assert_eq(mo,m2,mx,SIZE(d_o1,1),TRIM(modname)//" mo") |
---|
| 49 | no=assert_eq(no,n2,my,SIZE(d_o1,2),TRIM(modname)//" no") |
---|
| 50 | mask(:,:)=.TRUE.; IF(PRESENT(msk)) mask(:,:)=msk(:,:) |
---|
[2293] | 51 | |
---|
| 52 | !--- COMPUTE CELLS INTERFACES COORDINATES OF OUTPUT GRID |
---|
| 53 | b(mo)=x_o(mo)+(x_o(mo)-x_o(mo-1))/2.0; b(1:mo-1)=(x_o(1:mo-1)+x_o(2:mo))/2.0 |
---|
| 54 | d(no)=y_o(no)+(y_o(no)-y_o(no-1))/2.0; d(1:no-1)=(y_o(1:no-1)+y_o(2:no))/2.0 |
---|
| 55 | a(1 )=x_o(1 )-(x_o(2 )-x_o(1 ))/2.0; a(2:mo )= b(1:mo-1) |
---|
| 56 | c(1 )=y_o(1 )-(y_o(2 )-y_o(1 ))/2.0; c(2:no )= d(1:no-1) |
---|
| 57 | |
---|
| 58 | !--- ACCUMULATE INPUT POINTS ON OUTPUT GRID |
---|
[2299] | 59 | d_o1(:,:)=0.; num_tot(:,:)=0; inc=1.0 |
---|
[2293] | 60 | IF(lo) d_o2(:,:)=0. |
---|
| 61 | DO ji = 1, ni |
---|
| 62 | DO ii = 1, mi |
---|
[2299] | 63 | IF(li) inc=d_i(ii,ji) |
---|
[2293] | 64 | DO jo = 1, no |
---|
| 65 | IF((y_i(ji)-c(jo)<thresh.OR.y_i(ji)-d(jo)>thresh).AND. & |
---|
| 66 | (y_i(ji)-c(jo)>thresh.OR.y_i(ji)-d(jo)<thresh)) CYCLE |
---|
| 67 | DO io = 1, mo |
---|
| 68 | IF((x_i(ii)-a(io)<thresh.OR.x_i(ii)-b(io)>thresh).AND. & |
---|
| 69 | (x_i(ii)-a(io)>thresh.OR.x_i(ii)-b(io)<thresh)) CYCLE |
---|
| 70 | num_tot(io,jo)=num_tot(io,jo)+1 |
---|
[2361] | 71 | IF(mask(io,jo)) d_o1(io,jo)=d_o1(io,jo)+inc |
---|
[2293] | 72 | IF(.NOT.lo) CYCLE |
---|
[2361] | 73 | IF(mask(io,jo)) d_o2(io,jo)=d_o2(io,jo)+inc*inc |
---|
[2293] | 74 | END DO |
---|
| 75 | END DO |
---|
| 76 | END DO |
---|
| 77 | END DO |
---|
| 78 | |
---|
| 79 | !--- CHECK INPUT POINTS HAVE BEEN FOUND IN EACH OUTPUT CELL |
---|
| 80 | found(:,:)=num_tot(:,:)/=0 |
---|
[2299] | 81 | WHERE(found.AND.mask) d_o1(:,:)=d_o1(:,:)/REAL(num_tot(:,:)) |
---|
[2293] | 82 | IF(PRESENT(d_o2)) THEN |
---|
[2299] | 83 | WHERE(found.AND.mask) d_o2(:,:)=d_o2(:,:)/REAL(num_tot(:,:)) |
---|
[2293] | 84 | RETURN |
---|
| 85 | END IF |
---|
[2294] | 86 | nn=COUNT(found); IF(nn==0) RETURN |
---|
[2293] | 87 | |
---|
| 88 | !--- MISSING POINTS ; USE DISTANCE ON THE SPHERE TO FIND NEAREST POINT nr(2) |
---|
| 89 | DO io = 1, mo |
---|
| 90 | DO jo = 1, no |
---|
| 91 | IF(found(io,jo)) CYCLE |
---|
| 92 | ! IF(prt_level>=1) PRINT*, "Problem: point out of domain (i,j)=", io,jo |
---|
| 93 | CALL dist_sphe(x_o(io),y_o(jo),x_i,y_i,dist(:,:)) |
---|
| 94 | nr=MINLOC(dist(:,:))!; IF(prt_level>=1) PRINT*, "Solution: ", nr |
---|
[2299] | 95 | inc=1.0; IF(li) inc=d_i(nr(1),nr(2)) |
---|
[2361] | 96 | IF(mask(io,jo)) d_o1(io,jo)=inc |
---|
[2293] | 97 | END DO |
---|
| 98 | END DO |
---|
| 99 | |
---|
| 100 | END SUBROUTINE fine2coarse |
---|
| 101 | ! |
---|
| 102 | !------------------------------------------------------------------------------- |
---|
| 103 | |
---|
| 104 | |
---|
| 105 | !------------------------------------------------------------------------------- |
---|
| 106 | ! |
---|
| 107 | SUBROUTINE grille_m(xdata, ydata, entree, x, y, sortie) |
---|
| 108 | ! |
---|
| 109 | !------------------------------------------------------------------------------- |
---|
| 110 | ! Author: Z.X. Li (april 1st 1994) (see also A. Harzallah and L. Fairhead) |
---|
| 111 | !------------------------------------------------------------------------------- |
---|
| 112 | ! Purpose: Naive method to regrid on a coarser grid. |
---|
| 113 | ! Value at a new point is an average of the old points lcoated in a zone |
---|
| 114 | ! surrounding the considered new point. |
---|
| 115 | ! No ponderation is used (see grille_p) |
---|
| 116 | ! |
---|
| 117 | ! (c) |
---|
| 118 | ! ----d----- |
---|
| 119 | ! | . . . .| |
---|
| 120 | ! | | |
---|
| 121 | ! (b)a . * . .b(a) |
---|
| 122 | ! | | |
---|
| 123 | ! | . . . .| |
---|
| 124 | ! ----c----- |
---|
| 125 | ! (d) |
---|
| 126 | ! |
---|
| 127 | !------------------------------------------------------------------------------- |
---|
| 128 | IMPLICIT NONE |
---|
| 129 | !------------------------------------------------------------------------------- |
---|
| 130 | ! Arguments: |
---|
| 131 | REAL, INTENT(IN) :: xdata(:), ydata(:) !--- INPUT FIELD X AND Y COORD. |
---|
| 132 | REAL, INTENT(IN) :: entree(SIZE(xdata),SIZE(ydata)) !--- INPUT FIELD |
---|
| 133 | REAL, INTENT(IN) :: x(:), y(:) !--- OUTPUT FIELD X AND Y COORD. |
---|
| 134 | REAL, INTENT(OUT) :: sortie(SIZE(x),SIZE(y)) !--- OUTPUT FIELD |
---|
| 135 | !------------------------------------------------------------------------------- |
---|
[2299] | 136 | CALL fine2coarse(xdata,ydata,x,y,sortie,entree) |
---|
[2293] | 137 | |
---|
| 138 | END SUBROUTINE grille_m |
---|
| 139 | ! |
---|
| 140 | !------------------------------------------------------------------------------- |
---|
| 141 | |
---|
| 142 | |
---|
| 143 | !------------------------------------------------------------------------------- |
---|
| 144 | ! |
---|
| 145 | SUBROUTINE rugosite(xdata, ydata, entree, x, y, sortie, mask) |
---|
| 146 | ! |
---|
| 147 | !------------------------------------------------------------------------------- |
---|
| 148 | ! Author: Z.X. Li (april 1st 1994) |
---|
| 149 | !------------------------------------------------------------------------------- |
---|
[2299] | 150 | ! Purpose: Remap rugosity length ; constant value (0.001) on oceans. |
---|
[2293] | 151 | ! Naive method (see grille_m) |
---|
| 152 | !------------------------------------------------------------------------------- |
---|
| 153 | IMPLICIT NONE |
---|
| 154 | !------------------------------------------------------------------------------- |
---|
| 155 | ! Arguments: |
---|
| 156 | REAL, INTENT(IN) :: xdata(:), ydata(:) !--- INPUT FIELD X AND Y COORD. |
---|
| 157 | REAL, INTENT(IN) :: entree(SIZE(xdata),SIZE(ydata)) !--- INPUT FIELD |
---|
| 158 | REAL, INTENT(IN) :: x(:), y(:) !--- OUTPUT FIELD X AND Y COORD. |
---|
| 159 | REAL, INTENT(OUT) :: sortie(SIZE(x),SIZE(y)) !--- OUTPUT FIELD |
---|
| 160 | REAL, INTENT(IN) :: mask (SIZE(x),SIZE(y)) !--- MASK |
---|
| 161 | !------------------------------------------------------------------------------- |
---|
[2299] | 162 | CALL fine2coarse(xdata,ydata,x,y,sortie,LOG(entree)) |
---|
[2293] | 163 | WHERE(NINT(mask)==1) |
---|
[2299] | 164 | sortie(:,:)=EXP(sortie(:,:)) |
---|
[2293] | 165 | ELSE WHERE |
---|
[2299] | 166 | sortie(:,:)=0.001 |
---|
[2293] | 167 | END WHERE |
---|
| 168 | |
---|
| 169 | END SUBROUTINE rugosite |
---|
| 170 | ! |
---|
| 171 | !------------------------------------------------------------------------------- |
---|
| 172 | |
---|
| 173 | |
---|
| 174 | !------------------------------------------------------------------------------- |
---|
| 175 | ! |
---|
| 176 | SUBROUTINE sea_ice(xdata, ydata, glace01, x, y, frac_ice) |
---|
| 177 | ! |
---|
| 178 | !------------------------------------------------------------------------------- |
---|
| 179 | ! Author: Z.X. Li (april 1st 1994) |
---|
| 180 | ! Purpose: Regrid ice indicator (0 or 1) on a coarser grid to get an ice fract. |
---|
| 181 | ! field (between 0 and 1). |
---|
| 182 | ! Naive method (see grille_m) |
---|
| 183 | !------------------------------------------------------------------------------- |
---|
| 184 | IMPLICIT NONE |
---|
| 185 | !------------------------------------------------------------------------------- |
---|
| 186 | ! Arguments: |
---|
| 187 | REAL, INTENT(IN) :: xdata(:), ydata(:) !--- INPUT FIELD X AND Y COORD. |
---|
| 188 | REAL, INTENT(IN) :: glace01(:,:) !--- INPUT FIELD |
---|
| 189 | REAL, INTENT(IN) :: x(:), y(:) !--- OUTPUT FIELD X AND Y COORD. |
---|
| 190 | REAL, INTENT(OUT) :: frac_ice(SIZE(x),SIZE(y)) !--- OUTPUT FIELD |
---|
| 191 | !------------------------------------------------------------------------------- |
---|
[2299] | 192 | CALL fine2coarse(xdata,ydata,x,y,frac_ice,msk=NINT(glace01)==1) |
---|
[2293] | 193 | |
---|
| 194 | END SUBROUTINE sea_ice |
---|
| 195 | ! |
---|
| 196 | !------------------------------------------------------------------------------- |
---|
| 197 | |
---|
| 198 | |
---|
| 199 | !------------------------------------------------------------------------------- |
---|
| 200 | ! |
---|
| 201 | SUBROUTINE rugsoro(xrel, yrel, relief, xmod, ymod, rugs) |
---|
| 202 | ! |
---|
| 203 | !------------------------------------------------------------------------------- |
---|
| 204 | ! Author: Z.X. Li (april 1st 1994) ; Modifications: august 23rd 1995. |
---|
| 205 | ! Purpose: Compute rugosity due to orography by using standard dev in a 1x1 cell |
---|
| 206 | !------------------------------------------------------------------------------- |
---|
| 207 | IMPLICIT NONE |
---|
| 208 | !------------------------------------------------------------------------------- |
---|
| 209 | ! Arguments: |
---|
| 210 | REAL, INTENT(IN) :: xrel(:), yrel(:) !--- INPUT FIELD X AND Y COORD. |
---|
| 211 | REAL, INTENT(IN) :: relief(:,:) !--- INPUT FIELD |
---|
| 212 | REAL, INTENT(IN) :: xmod(:), ymod(:) !--- OUTPUT FIELD X AND Y COORD. |
---|
| 213 | REAL, INTENT(OUT) :: rugs(SIZE(xmod),SIZE(ymod)) !--- OUTPUT FIELD |
---|
| 214 | !------------------------------------------------------------------------------- |
---|
| 215 | ! Local variable: |
---|
| 216 | INTEGER :: k, nn |
---|
| 217 | INTEGER, PARAMETER:: itmp=360, jtmp=180 |
---|
[2299] | 218 | REAL :: out(SIZE(xmod),SIZE(xmod)), amin, amax |
---|
| 219 | REAL :: cham1tmp(itmp,jtmp), cham2tmp(itmp,jtmp), xtmp(itmp), ytmp(jtmp) |
---|
[2293] | 220 | !------------------------------------------------------------------------------- |
---|
| 221 | |
---|
| 222 | !--- TEST INPUT FILE FITS FOR ONE DEGREE RESOLUTION |
---|
| 223 | xtmp(:)=4.0*ATAN(1.0)*[(-1.0+REAL(2*k-1)/REAL(itmp),k=1,itmp)] |
---|
| 224 | ytmp(:)=2.0*ATAN(1.0)*[(-1.0+REAL(2*k-1)/REAL(jtmp),k=1,jtmp)] |
---|
| 225 | CALL fine2coarse(xrel,yrel,xtmp,ytmp,cham1tmp,relief,d_o2=cham2tmp) |
---|
| 226 | cham2tmp(:,:)=cham2tmp(:,:)-cham1tmp(:,:)**2 |
---|
| 227 | nn=COUNT(cham2tmp<=-7.5) |
---|
| 228 | IF(nn/=0) THEN |
---|
| 229 | PRINT*,'Problem for rugsoro ; std**2 < -7.5 for several points: ',nn |
---|
| 230 | CALL ABORT_GCM("", "", 1) |
---|
| 231 | END IF |
---|
| 232 | nn=COUNT(cham2tmp<0) |
---|
| 233 | IF(nn/=0) PRINT*,'Problem for rugsoro ; std**2 < 0. for several points: ',nn |
---|
[2299] | 234 | WHERE(cham2tmp<0.0) cham2tmp=0.0 |
---|
[2293] | 235 | cham2tmp(:,:)=SQRT(cham2tmp(:,:)) |
---|
| 236 | amin=MINVAL(cham2tmp); amax=MAXVAL(cham2tmp) |
---|
| 237 | PRINT*, 'Ecart-type 1x1:', amin, amax |
---|
| 238 | |
---|
| 239 | !--- COMPUTE RUGOSITY AT REQUIRED SCALE |
---|
[2299] | 240 | WHERE(cham2tmp<0.001) cham2tmp=0.001 |
---|
[2293] | 241 | CALL fine2coarse(xtmp,ytmp,xmod,ymod,out,REAL(LOG(cham2tmp))) |
---|
| 242 | out=EXP(out) |
---|
| 243 | amin=MINVAL(out); amax=MAXVAL(out) |
---|
| 244 | PRINT*, 'Ecart-type du modele:', amin, amax |
---|
[2299] | 245 | out=out/amax*20.0 |
---|
[2293] | 246 | amin=MINVAL(out); amax=MAXVAL(out) |
---|
| 247 | PRINT*, 'Longueur de rugosite du modele:', amin, amax |
---|
| 248 | rugs=REAL(out) |
---|
| 249 | |
---|
| 250 | END SUBROUTINE rugsoro |
---|
| 251 | ! |
---|
| 252 | !------------------------------------------------------------------------------- |
---|
| 253 | |
---|
| 254 | |
---|
| 255 | !------------------------------------------------------------------------------- |
---|
| 256 | ! |
---|
| 257 | SUBROUTINE dist_sphe(rf_lon,rf_lat,rlon,rlat,distance) |
---|
| 258 | ! |
---|
| 259 | !------------------------------------------------------------------------------- |
---|
| 260 | ! Author: Laurent Li (december 30th 1996). |
---|
[2361] | 261 | ! Purpose: Compute min. distance (along big circle) between 2 points in radians. |
---|
[2293] | 262 | !------------------------------------------------------------------------------- |
---|
| 263 | IMPLICIT NONE |
---|
| 264 | !------------------------------------------------------------------------------- |
---|
| 265 | ! Arguments: |
---|
[2361] | 266 | REAL, INTENT(IN) :: rf_lon, rf_lat !--- Reference point coordinates (radians) |
---|
| 267 | REAL, INTENT(IN) :: rlon(:), rlat(:)!--- Points longitudes/latitudes (radians) |
---|
| 268 | REAL, INTENT(OUT):: distance(SIZE(rlon),SIZE(rlat)) !--- Distance (radians) |
---|
[2293] | 269 | !------------------------------------------------------------------------------- |
---|
| 270 | ! Local variables: |
---|
| 271 | LOGICAL, SAVE :: first=.TRUE. |
---|
[2361] | 272 | REAL, SAVE :: pi, hpi |
---|
| 273 | REAL :: pa, pb, cpa, spa, crlo(SIZE(rlon)) |
---|
[2293] | 274 | INTEGER :: i, j |
---|
| 275 | !------------------------------------------------------------------------------- |
---|
[2361] | 276 | IF(first) THEN; pi=4.0*ATAN(1.0); hpi=pi/2.0; first=.FALSE.; END IF |
---|
| 277 | crlo(:)=COS(rf_lon-rlon(:)) !--- COS of points 1 and 2 angle |
---|
| 278 | pa=hpi-rf_lat !--- North Pole - Point 1 distance |
---|
[2293] | 279 | cpa=COS(pa); spa=SIN(pa) |
---|
| 280 | DO j=1,SIZE(rlat) |
---|
[2361] | 281 | pb=hpi-rlat(j) !--- North Pole - Point 2 distance |
---|
| 282 | distance(:,j)=ACOS(cpa*COS(pb)+spa*SIN(pb)*crlo(:)) |
---|
[2293] | 283 | END DO |
---|
| 284 | |
---|
| 285 | END SUBROUTINE dist_sphe |
---|
| 286 | ! |
---|
| 287 | !------------------------------------------------------------------------------- |
---|
| 288 | |
---|
| 289 | END MODULE grid_atob_m |
---|
| 290 | ! |
---|
| 291 | !******************************************************************************* |
---|
| 292 | |
---|