/[lmdze]/trunk/dyn3d/Inter_barxy/inter_barxy.f
ViewVC logotype

Diff of /trunk/dyn3d/Inter_barxy/inter_barxy.f

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

trunk/dyn3d/inter_barxy.f revision 97 by guez, Fri Apr 25 14:58:31 2014 UTC trunk/dyn3d/Inter_barxy/inter_barxy.f revision 98 by guez, Tue May 13 17:23:16 2014 UTC
# Line 4  module inter_barxy_m Line 4  module inter_barxy_m
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)
# Line 29  contains Line 30  contains
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    
# Line 82  contains Line 83  contains
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)
# Line 98  contains Line 99  contains
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
Legend:
Removed from v.97  
changed lines
  Added in v.98
  ViewVC Help
Powered by ViewVC 1.1.21