- Timestamp:
- Jul 24, 2024, 2:54:37 PM (4 months ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/branches/Amaury_dev/libf/dyn3d_common/inter_barxy_m.F90
r5113 r5116 1 2 1 ! $Id$ 3 2 … … 15 14 SUBROUTINE inter_barxy(dlonid, dlatid, champ, rlonimod, rlatimod, champint) 16 15 17 use lmdz_assert_eq, only: assert_eq18 use lmdz_assert, only: assert16 use lmdz_assert_eq, ONLY: assert_eq 17 use lmdz_assert, ONLY: assert 19 18 20 19 include "dimensions.h" … … 27 26 ! (for "aire", "apoln", "apols") 28 27 29 REAL, intent(in) :: dlonid(:)28 REAL, intent(in) :: dlonid(:) 30 29 ! (longitude from input file, in rad, from -pi to pi) 31 30 32 REAL, intent(in) :: dlatid(:), champ(:, :), rlonimod(:)33 34 REAL, intent(in) :: rlatimod(:)31 REAL, intent(in) :: dlatid(:), champ(:, :), rlonimod(:) 32 33 REAL, intent(in) :: rlatimod(:) 35 34 ! (latitude angle, in degrees or rad, in strictly decreasing order) 36 35 37 real, intent(out) :: champint(:, :)36 real, intent(out) :: champint(:, :) 38 37 ! Si taille de la seconde dim = jjm + 1, on veut interpoler sur les 39 38 ! jjm+1 latitudes rlatu du modele (latitudes des scalaires et de U) … … 55 54 56 55 jnterfd = assert_eq(size(champ, 2) - 1, size(dlatid), & 57 "inter_barxy jnterfd")56 "inter_barxy jnterfd") 58 57 jmods = size(champint, 2) 59 58 CALL assert(size(champ, 1) == size(dlonid), "inter_barxy size(champ, 1)") 60 59 CALL assert((/size(rlonimod), size(champint, 1)/) == iim, & 61 "inter_barxy iim")60 "inter_barxy iim") 62 61 CALL assert(any(jmods == (/jjm, jjm + 1/)), 'inter_barxy jmods') 63 62 CALL assert(size(rlatimod) == jjm, "inter_barxy size(rlatimod)") … … 65 64 ! Check decreasing order for "rlatimod": 66 65 DO i = 2, jjm 67 IF (rlatimod(i) >= rlatimod(i-1)) stop &68 '"inter_barxy": "rlatimod" should be strictly decreasing'66 IF (rlatimod(i) >= rlatimod(i - 1)) stop & 67 '"inter_barxy": "rlatimod" should be strictly decreasing' 69 68 ENDDO 70 69 71 70 yjmod(:jjm) = ord_coordm(rlatimod) 72 71 IF (jmods == jjm + 1) THEN 73 74 '"inter_barxy": with jmods = jjm + 1, yjmod(jjm) should be < 90.'72 IF (90. - yjmod(jjm) < 0.01) stop & 73 '"inter_barxy": with jmods = jjm + 1, yjmod(jjm) should be < 90.' 75 74 ELSE 76 77 78 '"inter_barxy": with jmods = jjm, yjmod(jjm) should be 90.'75 ! jmods = jjm 76 IF (ABS(yjmod(jjm) - 90.) > 0.01) stop & 77 '"inter_barxy": with jmods = jjm, yjmod(jjm) should be 90.' 79 78 ENDIF 80 79 … … 82 81 83 82 DO j = 1, jnterfd + 1 84 85 ENDDO 86 87 CALL ord_coord(dlatid, yjdat, decrois) 83 champy(:, j) = inter_barx(dlonid, champ(:, j), rlonimod) 84 ENDDO 85 86 CALL ord_coord(dlatid, yjdat, decrois) 88 87 IF (decrois) champy(:, :) = champy(:, jnterfd + 1:1:-1) 89 88 DO i = 1, iim 90 89 champint(i, :) = inter_bary(yjdat, champy(i, :), yjmod) 91 90 ENDDO 92 91 champint(:, :) = champint(:, jmods:1:-1) 93 92 94 93 IF (jmods == jjm + 1) THEN 95 96 champint(:, 1) = SUM(aire(:iim,1) * champint(:, 1)) / apoln97 98 * champint(:, jjm + 1)) / apols94 ! Valeurs uniques aux poles 95 champint(:, 1) = SUM(aire(:iim, 1) * champint(:, 1)) / apoln 96 champint(:, jjm + 1) = SUM(aire(:iim, jjm + 1) & 97 * champint(:, jjm + 1)) / apols 99 98 ENDIF 100 99 … … 103 102 !****************************** 104 103 105 function inter_barx(dlonid, fdat, rlonimod) 104 function inter_barx(dlonid, fdat, rlonimod) 106 105 107 106 ! INTERPOLATION BARYCENTRIQUE BASEE SUR LES AIRES … … 117 116 ! ( Les abscisses sont exprimees en degres) 118 117 119 use lmdz_assert_eq, only: assert_eq 118 USE lmdz_assert_eq, ONLY: assert_eq 119 USE lmdz_libmath, ONLY: minmax 120 120 121 121 IMPLICIT NONE 122 122 123 REAL, intent(in) :: dlonid(:)124 real, intent(in) :: fdat(:)125 real, intent(in) :: rlonimod(:)123 REAL, intent(in) :: dlonid(:) 124 real, intent(in) :: fdat(:) 125 real, intent(in) :: rlonimod(:) 126 126 127 127 real inter_barx(size(rlonimod)) … … 130 130 131 131 INTEGER idatmax, imodmax 132 REAL xxid(size(dlonid) +1), xxd(size(dlonid)+1), fdd(size(dlonid)+1)133 REAL fxd(size(dlonid) +1), xchan(size(dlonid)+1), fdchan(size(dlonid)+1)132 REAL xxid(size(dlonid) + 1), xxd(size(dlonid) + 1), fdd(size(dlonid) + 1) 133 REAL fxd(size(dlonid) + 1), xchan(size(dlonid) + 1), fdchan(size(dlonid) + 1) 134 134 REAL xxim(size(rlonimod)) 135 135 … … 149 149 ! A L'INTERFACE OUEST DE LA PREMIERE MAILLE DU MODELE 150 150 DO imod = 1, imodmax 151 152 ENDDO 153 154 CALL minmax( 155 IF( chmax<6.50) THEN156 157 xxim(imod) = xxim(imod) * 180./pi158 151 xxim(imod) = rlonimod(imod) 152 ENDDO 153 154 CALL minmax(imodmax, xxim, chmin, chmax) 155 IF(chmax<6.50) THEN 156 DO imod = 1, imodmax 157 xxim(imod) = xxim(imod) * 180. / pi 158 ENDDO 159 159 ENDIF 160 160 … … 162 162 163 163 DO imod = 1, imodmax 164 165 ENDDO 166 167 idatmax1 = idatmax + 1164 xxim(imod) = xxim(imod) - xim0 165 ENDDO 166 167 idatmax1 = idatmax + 1 168 168 169 169 DO idat = 1, idatmax 170 171 ENDDO 172 173 CALL minmax( 174 IF( chmax<6.50) THEN175 176 xxd(idat) = xxd(idat) * 180./pi177 170 xxd(idat) = dlonid(idat) 171 ENDDO 172 173 CALL minmax(idatmax, xxd, chmin, chmax) 174 IF(chmax<6.50) THEN 175 DO idat = 1, idatmax 176 xxd(idat) = xxd(idat) * 180. / pi 177 ENDDO 178 178 ENDIF 179 179 180 180 DO idat = 1, idatmax 181 xxd(idat) = MOD( xxd(idat) - xim0, 360.)182 181 xxd(idat) = MOD(xxd(idat) - xim0, 360.) 182 fdd(idat) = fdat (idat) 183 183 ENDDO 184 184 185 185 i = 2 186 DO while (xxd(i) >= xxd(i -1) .and. i < idatmax)187 188 ENDDO 189 IF (xxd(i) < xxd(i -1)) THEN190 191 192 nid = idatmax - ichang +1193 194 xchan (i) = xxd(i+ichang -1)195 fdchan(i) = fdd(i+ichang -1)196 197 DO i=1, ichang -1198 xchan (i+ nid) = xxd(i)199 fdchan(i+nid) = fdd(i)200 201 DO i =1, idatmax202 203 204 186 DO while (xxd(i) >= xxd(i - 1) .and. i < idatmax) 187 i = i + 1 188 ENDDO 189 IF (xxd(i) < xxd(i - 1)) THEN 190 ichang = i 191 ! *** reorganisation des longitudes entre 0. et 360. degres **** 192 nid = idatmax - ichang + 1 193 DO i = 1, nid 194 xchan (i) = xxd(i + ichang - 1) 195 fdchan(i) = fdd(i + ichang - 1) 196 ENDDO 197 DO i = 1, ichang - 1 198 xchan (i + nid) = xxd(i) 199 fdchan(i + nid) = fdd(i) 200 ENDDO 201 DO i = 1, idatmax 202 xxd(i) = xchan(i) 203 fdd(i) = fdchan(i) 204 ENDDO 205 205 end IF 206 206 … … 213 213 214 214 DO idat = 1, idatmax 215 IF ( xxd( idatmax1- idat)<360.) exit216 215 IF (xxd(idatmax1 - idat)<360.) exit 216 id1 = id1 + 1 217 217 ENDDO 218 218 219 219 DO idat = 1, idatmax 220 221 220 IF (xxd(idat)>0.) exit 221 id0 = id0 + 1 222 222 END DO 223 223 224 IF( id1 /= 0 ) then225 226 227 fxd (idat) = fdd(idatmax - id1 + idat)228 229 230 231 232 224 IF(id1 /= 0) THEN 225 DO idat = 1, id1 226 xxid(idat) = xxd(idatmax - id1 + idat) - 360. 227 fxd (idat) = fdd(idatmax - id1 + idat) 228 END DO 229 DO idat = 1, idatmax - id1 230 xxid(idat + id1) = xxd(idat) 231 fxd (idat + id1) = fdd(idat) 232 END DO 233 233 end IF 234 234 235 IF(id0 /= 0) then236 237 238 239 240 241 242 xxid (idatmax - id0 + idat) =xxd(idat) + 360.243 fxd (idatmax - id0 + idat) = fdd(idat)244 245 else 246 247 xxid(idat)= xxd(idat)248 fxd (idat)= fdd(idat)249 235 IF(id0 /= 0) THEN 236 DO idat = 1, idatmax - id0 237 xxid(idat) = xxd(idat + id0) 238 fxd (idat) = fdd(idat + id0) 239 END DO 240 241 DO idat = 1, id0 242 xxid (idatmax - id0 + idat) = xxd(idat) + 360. 243 fxd (idatmax - id0 + idat) = fdd(idat) 244 END DO 245 else 246 DO idat = 1, idatmax 247 xxid(idat) = xxd(idat) 248 fxd (idat) = fdd(idat) 249 ENDDO 250 250 end IF 251 251 xxid(idatmax1) = xxid(1) + 360. … … 258 258 ! iteration 259 259 260 x0 261 dxm 260 x0 = xim0 261 dxm = 0. 262 262 imod = 1 263 263 idat = 1 264 264 265 265 do while (imod <= imodmax) 266 267 dx= xxid(idat) - x0268 dxm= dxm + dx269 270 x0= xxid(idat)271 272 273 274 dx= xxim(imod) - x0275 dxm= dxm + dx276 277 x0= xxim(imod)278 dxm= 0.279 280 281 dx= xxim(imod) - x0282 dxm= dxm + dx283 284 x0= xxim(imod)285 dxm= 0.286 287 288 266 do while (xxim(imod)>xxid(idat)) 267 dx = xxid(idat) - x0 268 dxm = dxm + dx 269 inter_barx(imod) = inter_barx(imod) + dx * fxd(idat) 270 x0 = xxid(idat) 271 idat = idat + 1 272 END DO 273 IF (xxim(imod)<xxid(idat)) THEN 274 dx = xxim(imod) - x0 275 dxm = dxm + dx 276 inter_barx(imod) = (inter_barx(imod) + dx * fxd(idat)) / dxm 277 x0 = xxim(imod) 278 dxm = 0. 279 imod = imod + 1 280 ELSE 281 dx = xxim(imod) - x0 282 dxm = dxm + dx 283 inter_barx(imod) = (inter_barx(imod) + dx * fxd(idat)) / dxm 284 x0 = xxim(imod) 285 dxm = 0. 286 imod = imod + 1 287 idat = idat + 1 288 END IF 289 289 END DO 290 290 … … 299 299 ! L'indice 1 correspond à l'interface maille 1 -- maille 2. 300 300 301 use lmdz_assert, only: assert301 use lmdz_assert, ONLY: assert 302 302 303 303 IMPLICIT NONE 304 304 305 REAL, intent(in) :: yjdat(:)305 REAL, intent(in) :: yjdat(:) 306 306 ! (angles, ordonnées des interfaces des mailles des données, in 307 307 ! degrees, in increasing order) 308 308 309 REAL, intent(in) :: fdat(:) ! champ de données310 311 REAL, intent(in) :: yjmod(:)309 REAL, intent(in) :: fdat(:) ! champ de données 310 311 REAL, intent(in) :: yjmod(:) 312 312 ! (ordonnées des interfaces des mailles du modèle) 313 313 ! (in degrees, in strictly increasing order) … … 317 317 ! Variables local to the procedure: 318 318 319 REAL y0, dy, dym 319 REAL y0, dy, dym 320 320 INTEGER jdat ! indice du champ de données 321 321 integer jmod ! indice du champ du modèle … … 327 327 ! Initialisation des variables 328 328 inter_bary(:) = 0. 329 y0 330 dym 331 jmod 332 jdat 329 y0 = -90. 330 dym = 0. 331 jmod = 1 332 jdat = 1 333 333 334 334 do while (jmod <= size(yjmod)) 335 336 dy= yjdat(jdat) - y0337 dym= dym + dy338 339 y0= yjdat(jdat)340 jdat= jdat + 1341 342 343 dy= yjmod(jmod) - y0344 dym= dym + dy345 346 y0= yjmod(jmod)347 dym= 0.348 jmod= jmod + 1349 350 351 dy= yjmod(jmod) - y0352 dym= dym + dy353 354 y0= yjmod(jmod)355 dym= 0.356 jmod= jmod + 1357 jdat= jdat + 1358 335 do while (yjmod(jmod) > yjdat(jdat)) 336 dy = yjdat(jdat) - y0 337 dym = dym + dy 338 inter_bary(jmod) = inter_bary(jmod) + dy * fdat(jdat) 339 y0 = yjdat(jdat) 340 jdat = jdat + 1 341 END DO 342 IF (yjmod(jmod) < yjdat(jdat)) THEN 343 dy = yjmod(jmod) - y0 344 dym = dym + dy 345 inter_bary(jmod) = (inter_bary(jmod) + dy * fdat(jdat)) / dym 346 y0 = yjmod(jmod) 347 dym = 0. 348 jmod = jmod + 1 349 ELSE 350 ! {yjmod(jmod) == yjdat(jdat)} 351 dy = yjmod(jmod) - y0 352 dym = dym + dy 353 inter_bary(jmod) = (inter_bary(jmod) + dy * fdat(jdat)) / dym 354 y0 = yjmod(jmod) 355 dym = 0. 356 jmod = jmod + 1 357 jdat = jdat + 1 358 END IF 359 359 END DO 360 360 ! Le test de fin suppose que l'interface 0 est commune aux deux … … 373 373 ! Finally, the procedure adds 90° as the last value of the array. 374 374 375 use lmdz_assert_eq, only: assert_eq376 use comconst_mod, only: pi375 use lmdz_assert_eq, ONLY: assert_eq 376 use comconst_mod, ONLY: pi 377 377 378 378 IMPLICIT NONE 379 379 380 REAL, intent(in) :: xi(:)380 REAL, intent(in) :: xi(:) 381 381 ! (latitude, in degrees or radians, in increasing or decreasing order) 382 382 ! ("xi" should contain latitudes from pole to pole. … … 385 385 ! So the extreme values should not be 90° and -90°.) 386 386 387 REAL, intent(out) :: xo(:) ! angles in degrees388 LOGICAL, intent(out) :: decrois387 REAL, intent(out) :: xo(:) ! angles in degrees 388 LOGICAL, intent(out) :: decrois 389 389 390 390 ! Variables local to the procedure: … … 398 398 decrois = xi(2) < xi(1) 399 399 DO i = 3, nmax 400 IF (decrois .neqv. xi(i) < xi(i-1)) stop &401 '"ord_coord": latitudes are not monotonic'402 ENDDO 403 404 IF (abs(xi(1)) < pi) then405 406 400 IF (decrois .neqv. xi(i) < xi(i - 1)) stop & 401 '"ord_coord": latitudes are not monotonic' 402 ENDDO 403 404 IF (abs(xi(1)) < pi) THEN 405 ! "xi" contains latitudes in radians 406 xo(:nmax) = xi(:) * 180. / pi 407 407 else 408 409 408 ! "xi" contains latitudes in degrees 409 xo(:nmax) = xi(:) 410 410 end IF 411 411 412 412 IF (ABS(abs(xo(1)) - 90) < 0.001 .or. ABS(abs(xo(nmax)) - 90) < 0.001) THEN 413 414 415 // 'grid cells, not the centers of grid cells.'416 413 print *, "ord_coord" 414 PRINT *, '"xi" should contain the latitudes of the boundaries of ' & 415 // 'grid cells, not the centers of grid cells.' 416 STOP 417 417 ENDIF 418 418 … … 429 429 ! order. 430 430 431 use comconst_mod, only: pi431 use comconst_mod, ONLY: pi 432 432 433 433 IMPLICIT NONE 434 434 435 REAL, intent(in) :: xi(:) ! angle, in rad or degrees435 REAL, intent(in) :: xi(:) ! angle, in rad or degrees 436 436 REAL ord_coordm(size(xi)) ! angle, in degrees 437 437 … … 439 439 440 440 IF (xi(1) < 6.5) THEN 441 442 441 ! "xi" is in rad 442 ord_coordm(:) = xi(size(xi):1:-1) * 180. / pi 443 443 else 444 445 444 ! "xi" is in degrees 445 ord_coordm(:) = xi(size(xi):1:-1) 446 446 ENDIF 447 447
Note: See TracChangeset
for help on using the changeset viewer.
