4 |
|
|
5 |
implicit none |
implicit none |
6 |
|
|
|
private |
|
|
public inter_barxy |
|
|
|
|
7 |
contains |
contains |
8 |
|
|
9 |
SUBROUTINE inter_barxy(dlonid, dlatid, champ, rlonimod, rlatimod, champint) |
SUBROUTINE inter_barxy(dlonid, dlatid, champ, rlonimod, rlatimod, champint) |
10 |
|
|
11 |
! Author: P. Le Van |
! Author: P. Le Van |
12 |
|
|
|
use nr_util, only: assert_eq, assert |
|
|
use dimens_m, only: iim, jjm |
|
13 |
use comgeom, only: aire_2d, apoln, apols |
use comgeom, only: aire_2d, apoln, apols |
14 |
|
use dimens_m, only: iim, jjm |
15 |
|
use inter_barx_m, only: inter_barx |
16 |
|
use inter_bary_m, only: inter_bary |
17 |
|
use nr_util, only: assert_eq, assert |
18 |
|
use ord_coord_m, only: ord_coord |
19 |
|
use ord_coordm_m, only: ord_coordm |
20 |
|
|
21 |
REAL, intent(in):: dlonid(:) |
REAL, intent(in):: dlonid(:) |
22 |
! (longitude from input file, in rad, from -pi to pi) |
! (longitude from input file, in rad, from -pi to pi) |
30 |
! Si taille de la seconde dim = jjm + 1, on veut interpoler sur les |
! Si taille de la seconde dim = jjm + 1, on veut interpoler sur les |
31 |
! jjm+1 latitudes rlatu du modele (latitudes des scalaires et de U) |
! jjm+1 latitudes rlatu du modele (latitudes des scalaires et de U) |
32 |
! Si taille de la seconde dim = jjm, on veut interpoler sur les |
! Si taille de la seconde dim = jjm, on veut interpoler sur les |
33 |
! jjm latitudes rlatv du modèle (latitudes de V) |
! jjm latitudes rlatv du modèle (latitudes de V) |
34 |
|
|
35 |
! Local: |
! Local: |
36 |
|
|
83 |
champy(:, j) = inter_barx(dlonid, champ(:, j), rlonimod) |
champy(:, j) = inter_barx(dlonid, champ(:, j), rlonimod) |
84 |
ENDDO |
ENDDO |
85 |
|
|
86 |
CALL ord_coord(dlatid, yjdat, decrois) |
CALL ord_coord(dlatid, yjdat, decrois) |
87 |
IF (decrois) champy(:, :) = champy(:, jnterfd + 1:1:-1) |
IF (decrois) champy(:, :) = champy(:, jnterfd + 1:1:-1) |
88 |
DO i = 1, iim |
DO i = 1, iim |
89 |
champint(i, :) = inter_bary(yjdat, champy(i, :), yjmod) |
champint(i, :) = inter_bary(yjdat, champy(i, :), yjmod) |
99 |
|
|
100 |
END SUBROUTINE inter_barxy |
END SUBROUTINE inter_barxy |
101 |
|
|
|
!****************************** |
|
|
|
|
|
function inter_barx(dlonid, fdat, rlonimod) |
|
|
|
|
|
! From dyn3d/inter_barx.F, v 1.1.1.1 2004/05/19 12:53:06 |
|
|
|
|
|
! Auteurs : Robert Sadourny, P. Le Van |
|
|
|
|
|
! INTERPOLATION BARYCENTRIQUE BASEE SUR LES AIRES |
|
|
! VERSION UNIDIMENSIONNELLE , EN LONGITUDE . |
|
|
|
|
|
! idat : indice du champ de donnees, de 1 a idatmax |
|
|
! imod : indice du champ du modele, de 1 a imodmax |
|
|
! fdat(idat) : champ de donnees (entrees) |
|
|
! inter_barx(imod) : champ du modele (sorties) |
|
|
! dlonid(idat): abscisses des interfaces des mailles donnees |
|
|
! rlonimod(imod): abscisses des interfaces des mailles modele |
|
|
! ( L'indice 1 correspond a l'interface mailLE 1 / maille 2) |
|
|
! ( Les abscisses sont exprimées en degres) |
|
|
|
|
|
use nr_util, only: assert_eq |
|
|
|
|
|
IMPLICIT NONE |
|
|
|
|
|
REAL, intent(in):: dlonid(:) |
|
|
real, intent(in):: fdat(:) |
|
|
real, intent(in):: rlonimod(:) |
|
|
|
|
|
real inter_barx(size(rlonimod)) |
|
|
|
|
|
! ... Variables locales ... |
|
|
|
|
|
INTEGER idatmax, imodmax |
|
|
REAL xxid(size(dlonid)+1), xxd(size(dlonid)+1), fdd(size(dlonid)+1) |
|
|
REAL fxd(size(dlonid)+1), xchan(size(dlonid)+1), fdchan(size(dlonid)+1) |
|
|
REAL xxim(size(rlonimod)) |
|
|
|
|
|
REAL x0, xim0, dx, dxm |
|
|
REAL chmin, chmax, pi |
|
|
|
|
|
INTEGER imod, idat, i, ichang, id0, id1, nid, idatmax1 |
|
|
|
|
|
!----------------------------------------------------- |
|
|
|
|
|
idatmax = assert_eq(size(dlonid), size(fdat), "inter_barx idatmax") |
|
|
imodmax = size(rlonimod) |
|
|
|
|
|
pi = 2. * ASIN(1.) |
|
|
|
|
|
! REDEFINITION DE L'ORIGINE DES ABSCISSES |
|
|
! A L'INTERFACE OUEST DE LA PREMIERE MAILLE DU MODELE |
|
|
DO imod = 1, imodmax |
|
|
xxim(imod) = rlonimod(imod) |
|
|
ENDDO |
|
|
|
|
|
CALL minmax( imodmax, xxim, chmin, chmax) |
|
|
IF( chmax.LT.6.50 ) THEN |
|
|
DO imod = 1, imodmax |
|
|
xxim(imod) = xxim(imod) * 180./pi |
|
|
ENDDO |
|
|
ENDIF |
|
|
|
|
|
xim0 = xxim(imodmax) - 360. |
|
|
|
|
|
DO imod = 1, imodmax |
|
|
xxim(imod) = xxim(imod) - xim0 |
|
|
ENDDO |
|
|
|
|
|
idatmax1 = idatmax +1 |
|
|
|
|
|
DO idat = 1, idatmax |
|
|
xxd(idat) = dlonid(idat) |
|
|
ENDDO |
|
|
|
|
|
CALL minmax( idatmax, xxd, chmin, chmax) |
|
|
IF( chmax.LT.6.50 ) THEN |
|
|
DO idat = 1, idatmax |
|
|
xxd(idat) = xxd(idat) * 180./pi |
|
|
ENDDO |
|
|
ENDIF |
|
|
|
|
|
DO idat = 1, idatmax |
|
|
xxd(idat) = AMOD( xxd(idat) - xim0, 360. ) |
|
|
fdd(idat) = fdat (idat) |
|
|
ENDDO |
|
|
|
|
|
i = 2 |
|
|
DO while (xxd(i) >= xxd(i-1) .and. i < idatmax) |
|
|
i = i + 1 |
|
|
ENDDO |
|
|
IF (xxd(i) < xxd(i-1)) THEN |
|
|
ichang = i |
|
|
! *** reorganisation des longitudes entre 0. et 360. degres **** |
|
|
nid = idatmax - ichang +1 |
|
|
DO i = 1, nid |
|
|
xchan (i) = xxd(i+ichang -1 ) |
|
|
fdchan(i) = fdd(i+ichang -1 ) |
|
|
ENDDO |
|
|
DO i=1, ichang -1 |
|
|
xchan (i+ nid) = xxd(i) |
|
|
fdchan(i+nid) = fdd(i) |
|
|
ENDDO |
|
|
DO i =1, idatmax |
|
|
xxd(i) = xchan(i) |
|
|
fdd(i) = fdchan(i) |
|
|
ENDDO |
|
|
end IF |
|
|
|
|
|
! translation des champs de donnees par rapport |
|
|
! a la nouvelle origine, avec redondance de la |
|
|
! maille a cheval sur les bords |
|
|
|
|
|
id0 = 0 |
|
|
id1 = 0 |
|
|
|
|
|
DO idat = 1, idatmax |
|
|
IF ( xxd( idatmax1- idat ).LT.360.) exit |
|
|
id1 = id1 + 1 |
|
|
ENDDO |
|
|
|
|
|
DO idat = 1, idatmax |
|
|
IF (xxd(idat).GT.0.) exit |
|
|
id0 = id0 + 1 |
|
|
END DO |
|
|
|
|
|
IF( id1 /= 0 ) then |
|
|
DO idat = 1, id1 |
|
|
xxid(idat) = xxd(idatmax - id1 + idat) - 360. |
|
|
fxd (idat) = fdd(idatmax - id1 + idat) |
|
|
END DO |
|
|
DO idat = 1, idatmax - id1 |
|
|
xxid(idat + id1) = xxd(idat) |
|
|
fxd (idat + id1) = fdd(idat) |
|
|
END DO |
|
|
end IF |
|
|
|
|
|
IF(id0 /= 0) then |
|
|
DO idat = 1, idatmax - id0 |
|
|
xxid(idat) = xxd(idat + id0) |
|
|
fxd (idat) = fdd(idat + id0) |
|
|
END DO |
|
|
|
|
|
DO idat = 1, id0 |
|
|
xxid (idatmax - id0 + idat) = xxd(idat) + 360. |
|
|
fxd (idatmax - id0 + idat) = fdd(idat) |
|
|
END DO |
|
|
else |
|
|
DO idat = 1, idatmax |
|
|
xxid(idat) = xxd(idat) |
|
|
fxd (idat) = fdd(idat) |
|
|
ENDDO |
|
|
end IF |
|
|
xxid(idatmax1) = xxid(1) + 360. |
|
|
fxd (idatmax1) = fxd(1) |
|
|
|
|
|
! initialisation du champ du modele |
|
|
|
|
|
inter_barx(:) = 0. |
|
|
|
|
|
! iteration |
|
|
|
|
|
x0 = xim0 |
|
|
dxm = 0. |
|
|
imod = 1 |
|
|
idat = 1 |
|
|
|
|
|
do while (imod <= imodmax) |
|
|
do while (xxim(imod).GT.xxid(idat)) |
|
|
dx = xxid(idat) - x0 |
|
|
dxm = dxm + dx |
|
|
inter_barx(imod) = inter_barx(imod) + dx * fxd(idat) |
|
|
x0 = xxid(idat) |
|
|
idat = idat + 1 |
|
|
end do |
|
|
IF (xxim(imod).LT.xxid(idat)) THEN |
|
|
dx = xxim(imod) - x0 |
|
|
dxm = dxm + dx |
|
|
inter_barx(imod) = (inter_barx(imod) + dx * fxd(idat)) / dxm |
|
|
x0 = xxim(imod) |
|
|
dxm = 0. |
|
|
imod = imod + 1 |
|
|
ELSE |
|
|
dx = xxim(imod) - x0 |
|
|
dxm = dxm + dx |
|
|
inter_barx(imod) = (inter_barx(imod) + dx * fxd(idat)) / dxm |
|
|
x0 = xxim(imod) |
|
|
dxm = 0. |
|
|
imod = imod + 1 |
|
|
idat = idat + 1 |
|
|
END IF |
|
|
end do |
|
|
|
|
|
END function inter_barx |
|
|
|
|
|
!****************************** |
|
|
|
|
|
function inter_bary(yjdat, fdat, yjmod) |
|
|
|
|
|
! From dyn3d/inter_bary.F, version 1.1.1.1 2004/05/19 12:53:06 |
|
|
! Authors: R. Sadourny, P. Le Van |
|
|
|
|
|
! Interpolation barycentrique basée sur les aires. |
|
|
! Version unidimensionnelle, en latitude. |
|
|
! L'indice 1 correspond à l'interface maille 1 -- maille 2. |
|
|
|
|
|
use nr_util, only: assert |
|
|
|
|
|
IMPLICIT NONE |
|
|
|
|
|
REAL, intent(in):: yjdat(:) |
|
|
! (angles, ordonnées des interfaces des mailles des données, in |
|
|
! degrees, in increasing order) |
|
|
|
|
|
REAL, intent(in):: fdat(:) ! champ de données |
|
|
|
|
|
REAL, intent(in):: yjmod(:) |
|
|
! (ordonnées des interfaces des mailles du modèle) |
|
|
! (in degrees, in strictly increasing order) |
|
|
|
|
|
REAL inter_bary(size(yjmod)) ! champ du modèle |
|
|
|
|
|
! Variables local to the procedure: |
|
|
|
|
|
REAL y0, dy, dym |
|
|
INTEGER jdat ! indice du champ de données |
|
|
integer jmod ! indice du champ du modèle |
|
|
|
|
|
!------------------------------------ |
|
|
|
|
|
call assert(size(yjdat) == size(fdat), "inter_bary") |
|
|
|
|
|
! Initialisation des variables |
|
|
inter_bary(:) = 0. |
|
|
y0 = -90. |
|
|
dym = 0. |
|
|
jmod = 1 |
|
|
jdat = 1 |
|
|
|
|
|
do while (jmod <= size(yjmod)) |
|
|
do while (yjmod(jmod) > yjdat(jdat)) |
|
|
dy = yjdat(jdat) - y0 |
|
|
dym = dym + dy |
|
|
inter_bary(jmod) = inter_bary(jmod) + dy * fdat(jdat) |
|
|
y0 = yjdat(jdat) |
|
|
jdat = jdat + 1 |
|
|
end do |
|
|
IF (yjmod(jmod) < yjdat(jdat)) THEN |
|
|
dy = yjmod(jmod) - y0 |
|
|
dym = dym + dy |
|
|
inter_bary(jmod) = (inter_bary(jmod) + dy * fdat(jdat)) / dym |
|
|
y0 = yjmod(jmod) |
|
|
dym = 0. |
|
|
jmod = jmod + 1 |
|
|
ELSE |
|
|
! {yjmod(jmod) == yjdat(jdat)} |
|
|
dy = yjmod(jmod) - y0 |
|
|
dym = dym + dy |
|
|
inter_bary(jmod) = (inter_bary(jmod) + dy * fdat(jdat)) / dym |
|
|
y0 = yjmod(jmod) |
|
|
dym = 0. |
|
|
jmod = jmod + 1 |
|
|
jdat = jdat + 1 |
|
|
END IF |
|
|
end do |
|
|
! Le test de fin suppose que l'interface 0 est commune aux deux |
|
|
! grilles "yjdat" et "yjmod". |
|
|
|
|
|
END function inter_bary |
|
|
|
|
|
!****************************** |
|
|
|
|
|
SUBROUTINE ord_coord(xi, xo, decrois) |
|
|
|
|
|
! From dyn3d/ord_coord.F, version 1.1.1.1 2004/05/19 12:53:06 |
|
|
! Author : P. Le Van |
|
|
|
|
|
! This procedure receives an array of latitudes. |
|
|
! It converts them to degrees if they are in radians. |
|
|
! If the input latitudes are in decreasing order, the procedure |
|
|
! reverses their order. |
|
|
! Finally, the procedure adds 90° as the last value of the array. |
|
|
|
|
|
use nr_util, only: assert_eq, pi |
|
|
|
|
|
IMPLICIT NONE |
|
|
|
|
|
REAL, intent(in):: xi(:) |
|
|
! (latitude, in degrees or radians, in increasing or decreasing order) |
|
|
! ("xi" should contain latitudes from pole to pole. |
|
|
! "xi" should contain the latitudes of the boundaries of grid |
|
|
! cells, not the centers of grid cells. |
|
|
! So the extreme values should not be 90° and -90°.) |
|
|
|
|
|
REAL, intent(out):: xo(:) ! angles in degrees |
|
|
LOGICAL, intent(out):: decrois |
|
|
|
|
|
! Variables local to the procedure: |
|
|
INTEGER nmax, i |
|
|
|
|
|
!-------------------- |
|
|
|
|
|
nmax = assert_eq(size(xi), size(xo) - 1, "ord_coord") |
|
|
|
|
|
! Check monotonicity: |
|
|
decrois = xi(2) < xi(1) |
|
|
DO i = 3, nmax |
|
|
IF (decrois .neqv. xi(i) < xi(i-1)) then |
|
|
print *, '"ord_coord": latitudes are not monotonic' |
|
|
stop 1 |
|
|
end IF |
|
|
ENDDO |
|
|
|
|
|
IF (abs(xi(1)) < pi) then |
|
|
! "xi" contains latitudes in radians |
|
|
xo(:nmax) = xi(:) * 180. / pi |
|
|
else |
|
|
! "xi" contains latitudes in degrees |
|
|
xo(:nmax) = xi(:) |
|
|
end IF |
|
|
|
|
|
IF (ABS(abs(xo(1)) - 90) < 0.001 .or. ABS(abs(xo(nmax)) - 90) < 0.001) THEN |
|
|
print *, "ord_coord" |
|
|
PRINT *, '"xi" should contain the latitudes of the boundaries of ' & |
|
|
// 'grid cells, not the centers of grid cells.' |
|
|
STOP 1 |
|
|
ENDIF |
|
|
|
|
|
IF (decrois) xo(:nmax) = xo(nmax:1:- 1) |
|
|
xo(nmax + 1) = 90. |
|
|
|
|
|
END SUBROUTINE ord_coord |
|
|
|
|
|
!*********************************** |
|
|
|
|
|
function ord_coordm(xi) |
|
|
|
|
|
! From dyn3d/ord_coordm.F, version 1.1.1.1 2004/05/19 12:53:06 |
|
|
! Author : P. Le Van |
|
|
|
|
|
! This procedure converts to degrees, if necessary, and inverts the |
|
|
! order. |
|
|
|
|
|
use nr_util, only: pi |
|
|
|
|
|
IMPLICIT NONE |
|
|
|
|
|
REAL, intent(in):: xi(:) ! angle, in rad or degrees |
|
|
REAL ord_coordm(size(xi)) ! angle, in degrees |
|
|
|
|
|
!----------------------------- |
|
|
|
|
|
IF (xi(1) < 6.5) THEN |
|
|
! "xi" is in rad |
|
|
ord_coordm(:) = xi(size(xi):1:-1) * 180. / pi |
|
|
else |
|
|
! "xi" is in degrees |
|
|
ord_coordm(:) = xi(size(xi):1:-1) |
|
|
ENDIF |
|
|
|
|
|
END function ord_coordm |
|
|
|
|
102 |
end module inter_barxy_m |
end module inter_barxy_m |