- Timestamp:
- Dec 22, 2009, 12:07:26 PM (14 years ago)
- File:
-
- 1 copied
Legend:
- Unmodified
- Added
- Removed
-
LMDZ4/branches/LMDZ4V5.0-dev/libf/dyn3d/inter_barxy_m.F90
r1281 r1293 1 ! 2 ! $Header$ 3 ! 4 SUBROUTINE inter_barxy ( interfd,jnterfd,dlonid,dlatid , 5 , champ,imod,jmod,rlonimod,rlatimod, jsort,champint ) 6 7 c Auteur : P. Le Van 8 c 9 INTEGER interfd,jnterfd,imod,jmod 10 REAL champ(interfd,jnterfd +1 ),dlonid(interfd),dlatid(jnterfd), 11 , champint(imod,jsort) 12 REAL rlonimod(imod),rlatimod(jmod) 13 14 #include "dimensions.h" 15 #include "paramet.h" 16 #include "comgeom2.h" 17 18 REAL champx(imod),champy(jnterfd +1,imod),chpn(imod),chps(imod) 19 REAL chhpn,chhps 20 REAL fmody(jjp1) 21 c 22 23 DO j = 1, jnterfd + 1 24 CALL inter_barx( interfd, dlonid, champ( 1,j ), 25 , imod, rlonimod , champx ) 26 DO i = 1,imod 27 champy(j,i) = champx(i) 28 ENDDO 29 ENDDO 30 31 DO i = 1, imod 32 CALL inter_bary( jjm,jnterfd,dlatid,champy(1,i), 33 , jmod ,rlatimod, fmody ) 34 DO j = 1, jsort 35 champint(i,j) = fmody(j) 36 ENDDO 37 ENDDO 38 39 IF( jsort.EQ.jjp1) THEN 40 41 c .... Valeurs uniques aux poles .... 42 c 43 DO i = 1,imod 44 chpn(i) = aire( i, 1 ) * champint( i, 1 ) 45 chps(i) = aire( i, jjp1 ) * champint( i,jjp1 ) 46 ENDDO 47 chhpn = SSUM(imod,chpn,1)/apoln 48 chhps = SSUM(imod,chps,1)/apols 49 50 DO i = 1, imod 51 champint( i, 1 ) = chhpn 52 champint( i, jjp1) = chhps 53 ENDDO 54 c 55 ENDIF 56 57 RETURN 58 END 59 1 module inter_barxy_m 2 3 ! Authors: Robert SADOURNY, Phu LE VAN, Lionel GUEZ 4 5 implicit none 6 7 private 8 public inter_barxy 9 10 contains 11 12 SUBROUTINE inter_barxy(dlonid, dlatid, champ, rlonimod, rlatimod, champint) 13 14 use assert_eq_m, only: assert_eq 15 use assert_m, only: assert 16 17 include "dimensions.h" 18 ! (for "iim", "jjm") 19 20 include "paramet.h" 21 ! (for other included files) 22 23 include "comgeom2.h" 24 ! (for "aire", "apoln", "apols") 25 26 REAL, intent(in):: dlonid(:) 27 ! (longitude from input file, in rad, from -pi to pi) 28 29 REAL, intent(in):: dlatid(:), champ(:, :), rlonimod(:) 30 31 REAL, intent(in):: rlatimod(:) 32 ! (latitude angle, in degrees or rad, in strictly decreasing order) 33 34 real, intent(out):: champint(:, :) 35 ! Si taille de la seconde dim = jjm + 1, on veut interpoler sur les 36 ! jjm+1 latitudes rlatu du modele (latitudes des scalaires et de U) 37 ! Si taille de la seconde dim = jjm, on veut interpoler sur les 38 ! jjm latitudes rlatv du modèle (latitudes de V) 39 40 ! Variables local to the procedure: 41 42 REAL champy(iim, size(champ, 2)) 43 integer j, i, jnterfd, jmods 44 45 REAL yjmod(size(champint, 2)) 46 ! (angle, in degrees, in strictly increasing order) 47 48 REAL yjdat(size(dlatid) + 1) ! angle, in degrees, in increasing order 49 LOGICAL decrois ! "dlatid" is in decreasing order 50 51 !----------------------------------- 52 53 jnterfd = assert_eq(size(champ, 2) - 1, size(dlatid), & 54 "inter_barxy jnterfd") 55 jmods = size(champint, 2) 56 call assert(size(champ, 1) == size(dlonid), "inter_barxy size(champ, 1)") 57 call assert((/size(rlonimod), size(champint, 1)/) == iim, & 58 "inter_barxy iim") 59 call assert(any(jmods == (/jjm, jjm + 1/)), 'inter_barxy jmods') 60 call assert(size(rlatimod) == jjm, "inter_barxy size(rlatimod)") 61 62 ! Check decreasing order for "rlatimod": 63 DO i = 2, jjm 64 IF (rlatimod(i) >= rlatimod(i-1)) stop & 65 '"inter_barxy": "rlatimod" should be strictly decreasing' 66 ENDDO 67 68 yjmod(:jjm) = ord_coordm(rlatimod) 69 IF (jmods == jjm + 1) THEN 70 IF (90. - yjmod(jjm) < 0.01) stop & 71 '"inter_barxy": with jmods = jjm + 1, yjmod(jjm) should be < 90.' 72 ELSE 73 ! jmods = jjm 74 IF (ABS(yjmod(jjm) - 90.) > 0.01) stop & 75 '"inter_barxy": with jmods = jjm, yjmod(jjm) should be 90.' 76 ENDIF 77 78 if (jmods == jjm + 1) yjmod(jjm + 1) = 90. 79 80 DO j = 1, jnterfd + 1 81 champy(:, j) = inter_barx(dlonid, champ(:, j), rlonimod) 82 ENDDO 83 84 CALL ord_coord(dlatid, yjdat, decrois) 85 IF (decrois) champy(:, :) = champy(:, jnterfd + 1:1:-1) 86 DO i = 1, iim 87 champint(i, :) = inter_bary(yjdat, champy(i, :), yjmod) 88 ENDDO 89 champint(:, :) = champint(:, jmods:1:-1) 90 91 IF (jmods == jjm + 1) THEN 92 ! Valeurs uniques aux poles 93 champint(:, 1) = SUM(aire(:iim, 1) * champint(:, 1)) / apoln 94 champint(:, jjm + 1) = SUM(aire(:iim, jjm + 1) & 95 * champint(:, jjm + 1)) / apols 96 ENDIF 97 98 END SUBROUTINE inter_barxy 99 100 !****************************** 101 102 function inter_barx(dlonid, fdat, rlonimod) 103 104 ! INTERPOLATION BARYCENTRIQUE BASEE SUR LES AIRES 105 ! VERSION UNIDIMENSIONNELLE , EN LONGITUDE . 106 107 ! idat : indice du champ de donnees, de 1 a idatmax 108 ! imod : indice du champ du modele, de 1 a imodmax 109 ! fdat(idat) : champ de donnees (entrees) 110 ! inter_barx(imod) : champ du modele (sorties) 111 ! dlonid(idat): abscisses des interfaces des mailles donnees 112 ! rlonimod(imod): abscisses des interfaces des mailles modele 113 ! ( L'indice 1 correspond a l'interface mailLE 1 / maille 2) 114 ! ( Les abscisses sont exprimées en degres) 115 116 use assert_eq_m, only: assert_eq 117 118 IMPLICIT NONE 119 120 REAL, intent(in):: dlonid(:) 121 real, intent(in):: fdat(:) 122 real, intent(in):: rlonimod(:) 123 124 real inter_barx(size(rlonimod)) 125 126 ! ... Variables locales ... 127 128 INTEGER idatmax, imodmax 129 REAL xxid(size(dlonid)+1), xxd(size(dlonid)+1), fdd(size(dlonid)+1) 130 REAL fxd(size(dlonid)+1), xchan(size(dlonid)+1), fdchan(size(dlonid)+1) 131 REAL xxim(size(rlonimod)) 132 133 REAL x0, xim0, dx, dxm 134 REAL chmin, chmax, pi 135 136 INTEGER imod, idat, i, ichang, id0, id1, nid, idatmax1 137 138 !----------------------------------------------------- 139 140 idatmax = assert_eq(size(dlonid), size(fdat), "inter_barx idatmax") 141 imodmax = size(rlonimod) 142 143 pi = 2. * ASIN(1.) 144 145 ! REDEFINITION DE L'ORIGINE DES ABSCISSES 146 ! A L'INTERFACE OUEST DE LA PREMIERE MAILLE DU MODELE 147 DO imod = 1, imodmax 148 xxim(imod) = rlonimod(imod) 149 ENDDO 150 151 CALL minmax( imodmax, xxim, chmin, chmax) 152 IF( chmax.LT.6.50 ) THEN 153 DO imod = 1, imodmax 154 xxim(imod) = xxim(imod) * 180./pi 155 ENDDO 156 ENDIF 157 158 xim0 = xxim(imodmax) - 360. 159 160 DO imod = 1, imodmax 161 xxim(imod) = xxim(imod) - xim0 162 ENDDO 163 164 idatmax1 = idatmax +1 165 166 DO idat = 1, idatmax 167 xxd(idat) = dlonid(idat) 168 ENDDO 169 170 CALL minmax( idatmax, xxd, chmin, chmax) 171 IF( chmax.LT.6.50 ) THEN 172 DO idat = 1, idatmax 173 xxd(idat) = xxd(idat) * 180./pi 174 ENDDO 175 ENDIF 176 177 DO idat = 1, idatmax 178 xxd(idat) = AMOD( xxd(idat) - xim0, 360. ) 179 fdd(idat) = fdat (idat) 180 ENDDO 181 182 i = 2 183 DO while (xxd(i) >= xxd(i-1) .and. i < idatmax) 184 i = i + 1 185 ENDDO 186 IF (xxd(i) < xxd(i-1)) THEN 187 ichang = i 188 ! *** reorganisation des longitudes entre 0. et 360. degres **** 189 nid = idatmax - ichang +1 190 DO i = 1, nid 191 xchan (i) = xxd(i+ichang -1 ) 192 fdchan(i) = fdd(i+ichang -1 ) 193 ENDDO 194 DO i=1, ichang -1 195 xchan (i+ nid) = xxd(i) 196 fdchan(i+nid) = fdd(i) 197 ENDDO 198 DO i =1, idatmax 199 xxd(i) = xchan(i) 200 fdd(i) = fdchan(i) 201 ENDDO 202 end IF 203 204 ! translation des champs de donnees par rapport 205 ! a la nouvelle origine, avec redondance de la 206 ! maille a cheval sur les bords 207 208 id0 = 0 209 id1 = 0 210 211 DO idat = 1, idatmax 212 IF ( xxd( idatmax1- idat ).LT.360.) exit 213 id1 = id1 + 1 214 ENDDO 215 216 DO idat = 1, idatmax 217 IF (xxd(idat).GT.0.) exit 218 id0 = id0 + 1 219 END DO 220 221 IF( id1 /= 0 ) then 222 DO idat = 1, id1 223 xxid(idat) = xxd(idatmax - id1 + idat) - 360. 224 fxd (idat) = fdd(idatmax - id1 + idat) 225 END DO 226 DO idat = 1, idatmax - id1 227 xxid(idat + id1) = xxd(idat) 228 fxd (idat + id1) = fdd(idat) 229 END DO 230 end IF 231 232 IF(id0 /= 0) then 233 DO idat = 1, idatmax - id0 234 xxid(idat) = xxd(idat + id0) 235 fxd (idat) = fdd(idat + id0) 236 END DO 237 238 DO idat = 1, id0 239 xxid (idatmax - id0 + idat) = xxd(idat) + 360. 240 fxd (idatmax - id0 + idat) = fdd(idat) 241 END DO 242 else 243 DO idat = 1, idatmax 244 xxid(idat) = xxd(idat) 245 fxd (idat) = fdd(idat) 246 ENDDO 247 end IF 248 xxid(idatmax1) = xxid(1) + 360. 249 fxd (idatmax1) = fxd(1) 250 251 ! initialisation du champ du modele 252 253 inter_barx(:) = 0. 254 255 ! iteration 256 257 x0 = xim0 258 dxm = 0. 259 imod = 1 260 idat = 1 261 262 do while (imod <= imodmax) 263 do while (xxim(imod).GT.xxid(idat)) 264 dx = xxid(idat) - x0 265 dxm = dxm + dx 266 inter_barx(imod) = inter_barx(imod) + dx * fxd(idat) 267 x0 = xxid(idat) 268 idat = idat + 1 269 end do 270 IF (xxim(imod).LT.xxid(idat)) THEN 271 dx = xxim(imod) - x0 272 dxm = dxm + dx 273 inter_barx(imod) = (inter_barx(imod) + dx * fxd(idat)) / dxm 274 x0 = xxim(imod) 275 dxm = 0. 276 imod = imod + 1 277 ELSE 278 dx = xxim(imod) - x0 279 dxm = dxm + dx 280 inter_barx(imod) = (inter_barx(imod) + dx * fxd(idat)) / dxm 281 x0 = xxim(imod) 282 dxm = 0. 283 imod = imod + 1 284 idat = idat + 1 285 END IF 286 end do 287 288 END function inter_barx 289 290 !****************************** 291 292 function inter_bary(yjdat, fdat, yjmod) 293 294 ! Interpolation barycentrique basée sur les aires. 295 ! Version unidimensionnelle, en latitude. 296 ! L'indice 1 correspond à l'interface maille 1 -- maille 2. 297 298 use assert_m, only: assert 299 300 IMPLICIT NONE 301 302 REAL, intent(in):: yjdat(:) 303 ! (angles, ordonnées des interfaces des mailles des données, in 304 ! degrees, in increasing order) 305 306 REAL, intent(in):: fdat(:) ! champ de données 307 308 REAL, intent(in):: yjmod(:) 309 ! (ordonnées des interfaces des mailles du modèle) 310 ! (in degrees, in strictly increasing order) 311 312 REAL inter_bary(size(yjmod)) ! champ du modèle 313 314 ! Variables local to the procedure: 315 316 REAL y0, dy, dym 317 INTEGER jdat ! indice du champ de données 318 integer jmod ! indice du champ du modèle 319 320 !------------------------------------ 321 322 call assert(size(yjdat) == size(fdat), "inter_bary") 323 324 ! Initialisation des variables 325 inter_bary(:) = 0. 326 y0 = -90. 327 dym = 0. 328 jmod = 1 329 jdat = 1 330 331 do while (jmod <= size(yjmod)) 332 do while (yjmod(jmod) > yjdat(jdat)) 333 dy = yjdat(jdat) - y0 334 dym = dym + dy 335 inter_bary(jmod) = inter_bary(jmod) + dy * fdat(jdat) 336 y0 = yjdat(jdat) 337 jdat = jdat + 1 338 end do 339 IF (yjmod(jmod) < yjdat(jdat)) THEN 340 dy = yjmod(jmod) - y0 341 dym = dym + dy 342 inter_bary(jmod) = (inter_bary(jmod) + dy * fdat(jdat)) / dym 343 y0 = yjmod(jmod) 344 dym = 0. 345 jmod = jmod + 1 346 ELSE 347 ! {yjmod(jmod) == yjdat(jdat)} 348 dy = yjmod(jmod) - y0 349 dym = dym + dy 350 inter_bary(jmod) = (inter_bary(jmod) + dy * fdat(jdat)) / dym 351 y0 = yjmod(jmod) 352 dym = 0. 353 jmod = jmod + 1 354 jdat = jdat + 1 355 END IF 356 end do 357 ! Le test de fin suppose que l'interface 0 est commune aux deux 358 ! grilles "yjdat" et "yjmod". 359 360 END function inter_bary 361 362 !****************************** 363 364 SUBROUTINE ord_coord(xi, xo, decrois) 365 366 ! This procedure receives an array of latitudes. 367 ! It converts them to degrees if they are in radians. 368 ! If the input latitudes are in decreasing order, the procedure 369 ! reverses their order. 370 ! Finally, the procedure adds 90° as the last value of the array. 371 372 use assert_eq_m, only: assert_eq 373 374 IMPLICIT NONE 375 376 include "comconst.h" 377 ! (for "pi") 378 379 REAL, intent(in):: xi(:) 380 ! (latitude, in degrees or radians, in increasing or decreasing order) 381 ! ("xi" should contain latitudes from pole to pole. 382 ! "xi" should contain the latitudes of the boundaries of grid 383 ! cells, not the centers of grid cells. 384 ! So the extreme values should not be 90° and -90°.) 385 386 REAL, intent(out):: xo(:) ! angles in degrees 387 LOGICAL, intent(out):: decrois 388 389 ! Variables local to the procedure: 390 INTEGER nmax, i 391 392 !-------------------- 393 394 nmax = assert_eq(size(xi), size(xo) - 1, "ord_coord") 395 396 ! Check monotonicity: 397 decrois = xi(2) < xi(1) 398 DO i = 3, nmax 399 IF (decrois .neqv. xi(i) < xi(i-1)) stop & 400 '"ord_coord": latitudes are not monotonic' 401 ENDDO 402 403 IF (abs(xi(1)) < pi) then 404 ! "xi" contains latitudes in radians 405 xo(:nmax) = xi(:) * 180. / pi 406 else 407 ! "xi" contains latitudes in degrees 408 xo(:nmax) = xi(:) 409 end IF 410 411 IF (ABS(abs(xo(1)) - 90) < 0.001 .or. ABS(abs(xo(nmax)) - 90) < 0.001) THEN 412 print *, "ord_coord" 413 PRINT *, '"xi" should contain the latitudes of the boundaries of ' & 414 // 'grid cells, not the centers of grid cells.' 415 STOP 416 ENDIF 417 418 IF (decrois) xo(:nmax) = xo(nmax:1:- 1) 419 xo(nmax + 1) = 90. 420 421 END SUBROUTINE ord_coord 422 423 !*********************************** 424 425 function ord_coordm(xi) 426 427 ! This procedure converts to degrees, if necessary, and inverts the 428 ! order. 429 430 IMPLICIT NONE 431 432 include "comconst.h" 433 ! (for "pi") 434 435 REAL, intent(in):: xi(:) ! angle, in rad or degrees 436 REAL ord_coordm(size(xi)) ! angle, in degrees 437 438 !----------------------------- 439 440 IF (xi(1) < 6.5) THEN 441 ! "xi" is in rad 442 ord_coordm(:) = xi(size(xi):1:-1) * 180. / pi 443 else 444 ! "xi" is in degrees 445 ord_coordm(:) = xi(size(xi):1:-1) 446 ENDIF 447 448 END function ord_coordm 449 450 end module inter_barxy_m
Note: See TracChangeset
for help on using the changeset viewer.
