Philippe 23/02/2016: lfi2cdf: modif for JPHEXT/=1
[MNH-git_open_source-lfs.git] / LIBTOOLS / tools / diachro / src / EXTRACTDIA / int2lalo.f90
1 !----------------------------------------------------------------
2 !--------------- special set of characters for SCCS information
3 !-----------------------------------------------------------------
4 !      @(#) Lib:./s.int2lalo.f90, Version:1.8, Date:03/06/05, Last modified:01/10/15
5 !-----------------------------------------------------------------
6 !     ######spl
7 MODULE MODI_INT2LALO
8 !###################
9 !
10 INTERFACE
11       SUBROUTINE INT2LALO(HHORTYPE,P3D,PLATLON,PSVAL,PLALO)
12 !
13 CHARACTER(LEN=4),     INTENT(IN) :: HHORTYPE ! type of horizontal interpolation
14 REAL,DIMENSION(:,:,:),INTENT(IN) :: P3D  ! input 3d array s->n, w->e
15 REAL,DIMENSION(4),    INTENT(IN) :: PLATLON  ! NSWE target domain bounds (milliDEGS)
16 REAL,                 INTENT(IN) :: PSVAL    ! value for missing data
17 REAL,DIMENSION(:,:,:),INTENT(OUT):: PLALO    ! output interpolated LAT/LON field
18 !
19 END SUBROUTINE INT2LALO
20 END INTERFACE
21 END MODULE MODI_INT2LALO
22 !
23 !####################################################
24 SUBROUTINE INT2LALO(HHORTYPE,P3D,PLATLON,PSVAL,PLALO)
25 !####################################################
26 !
27 !!    PURPOSE
28 !!    -------
29 !       Interpolates data from a conformal grid to a lat/lon grid
30 !
31 !!**  METHOD
32 !!    ------
33 !!       Input is the conformal data (S->N scanning) and lat/lon domain 
34 !!      definition.
35 !!       Output is the lat/lon data in N->S scanning (required).
36 !-------------------------------------------------------------------------------
37 !
38 !*       0.    DECLARATIONS
39 !              ------------
40 !
41 USE MODD_CST,        ONLY : XRADIUS,XPI
42 USE MODD_GRID,       ONLY : XLONORI,XLATORI
43 USE MODD_PARAMETERS, ONLY : XUNDEF,JPHEXT
44 USE MODD_DIM1,       ONLY : NIMAX,NJMAX,NKMAX
45 USE MODD_GRID1,      ONLY : XXHAT,XYHAT
46 !
47 USE MODE_GRIDPROJ
48 !
49 IMPLICIT NONE
50 !
51 !*       0.1   Arguments
52 !
53 CHARACTER(LEN=4),     INTENT(IN) :: HHORTYPE ! type of horizontal interpolation
54 REAL,DIMENSION(:,:,:),INTENT(IN) :: P3D      ! input 3d array s->n, w->e
55 REAL,DIMENSION(4),    INTENT(IN) :: PLATLON  ! NSWE target domain bounds (milliDEGS)
56 REAL,                 INTENT(IN) :: PSVAL    ! value for missing data
57 REAL,DIMENSION(:,:,:),INTENT(OUT):: PLALO    ! output interpolated LAT/LON field
58                                              !with N->S scanning
59 !
60 !*       0.2   Local variables
61 !
62 REAL :: ZLONW,ZLONE,ZLATN,ZLATS   !(degres)
63 REAL :: ZXHAT,ZYHAT
64 REAL :: ZDX,ZDY                  ! TARGET INCREMENTS IN LAT/LON
65 REAL :: ZLA,ZLO,ZAB,ZCD,ZI,ZJ,ZXR,ZYR
66 REAL :: ZEPS
67 INTEGER :: JX,JY,JK,INX,INY,II,IJ,IK,IM,IN
68 !
69 !------------------------------------------------------------------------------
70 !
71 !*       1.    INITIALISATION
72 !              --------------
73 !
74 ZEPS=1.E-10
75 PLALO(:,:,:)= PSVAL
76 INX=SIZE(PLALO,1) ; INY=SIZE(PLALO,2) ; IK=SIZE(PLALO,3)
77 ZLONW=PLATLON(3)/1000.  ; ZLONE=PLATLON(4)/1000. ; ZLATN=PLATLON(1)/1000. ; ZLATS=PLATLON(2)/1000.
78 !
79 ZDX=(ZLONE-ZLONW)/(INX-1)
80 IF (ZDX<0) ZDX=(ZLONE-ZLONW+360.)/(INX-1)
81 ZDY=(ZLATN-ZLATS)/(INY-1)
82 print*, 'INT2LALO: target increments:',ZDX,ZDY
83 PRINT*, 'INT2LALO: target increments:',ZDX,ZDY
84 !
85 !------------------------------------------------------------------------------
86 !
87 !*       2.    INTERPOLATION
88 !              -------------
89 !
90 !print*,'av interp.: ',minval(P3D),minloc(p3d),maxval(p3d),maxloc(p3d)
91 DO JK=1,IK
92   DO JY=1,INY ; DO JX=1,INX
93     ZLO=MOD(ZLONW+ZDX*(JX-1),360.)
94     ZLA=ZLATN-ZDY*(JY-1)          ! output has N->S scanning
95  !   print*,ZLO,ZLA,JX,JY
96     CALL SM_XYHAT(XLATORI,XLONORI,ZLA,ZLO,ZXHAT,ZYHAT)
97     II=MAX(MIN(COUNT(XXHAT(:)<ZXHAT),NIMAX+JPHEXT),1+JPHEXT)
98     IJ=MAX(MIN(COUNT(XYHAT(:)<ZYHAT),NJMAX+JPHEXT),1+JPHEXT)
99     ZI=(ZXHAT-XXHAT(II))/(XXHAT(II+1)-XXHAT(II))+FLOAT(II)
100     ZJ=(ZYHAT-XYHAT(IJ))/(XYHAT(IJ+1)-XYHAT(IJ))+FLOAT(IJ)
101     !
102 !!!!!!!!!!!!!!! PLUSE DE DECALAGE D INDICES ENTRE P3D ET LE TABLEAUX MNH!!!!!!!!! 
103     IM=II
104     IN=IJ
105 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
106     !
107     IF (HHORTYPE=='NEAR') THEN 
108     ! NEARest neighbour method on conformal plane
109       IF (      (II>=1+JPHEXT) .AND. (II<=NIMAX+JPHEXT) &
110          .AND.  (IJ>=1+JPHEXT) .AND. (IJ<=NJMAX+JPHEXT) ) THEN
111         PLALO(JX,JY,JK)=P3D(IM,IN,JK) ! take nearest-neighbour value
112       ENDIF
113     !
114     ELSEIF (HHORTYPE == 'BILI') THEN
115     ! LInear interpolation method on conformal plane
116       IF (      (ZI>=1+JPHEXT) .AND. (ZI<=NIMAX+JPHEXT) &
117          .AND.  (ZJ>=1+JPHEXT) .AND. (ZJ<=NJMAX+JPHEXT) ) THEN
118         IF (ALL(ABS(P3D(IM:IM+1,IN:IN+1,JK)-XUNDEF)>=ZEPS) ) THEN
119         ! take the 4 surrounding values and apply bilinear interpolation
120           ZXR=ZI-REAL(II) ; ZYR=ZJ-REAL(IJ)  ! coordinates inside rectangle
121           ZAB= (1.-ZXR)*P3D(IM,  IN,  JK)  &
122               +    ZXR *P3D(IM+1,IN,  JK)
123           ZCD= (1.-ZXR)*P3D(IM,  IN+1,JK)  &
124               +    ZXR *P3D(IM+1,IN+1,JK)
125           PLALO(JX,JY,JK)= (1.-ZYR)*ZAB + ZYR*ZCD
126         ENDIF
127       ENDIF
128     ELSE
129       print*, 'Horizontal type interpolation unknown ',HHORTYPE
130     ENDIF
131   ENDDO ; ENDDO 
132 ENDDO
133 !
134 !------------------------------------------------------------------------------
135 !
136 !*       3.    EXTENSION
137 !              ---------
138 !
139 !IF (IK/=1) CALL EXTENDLAM
140 !print*,'ap interp.: ',minval(Plalo),minloc(plalo),maxval(plalo),maxloc(plalo)
141 !
142 RETURN
143 CONTAINS 
144 !-----------------------------
145 SUBROUTINE EXTENDLAM
146 !  PURPOSE: EXTEND AN INTERPOLATED LAT/LON FIELD OUTSIDE THE kAl MODEL
147 !           DOMAIN BY REMOVING ALL ITS UNDEFINED VALUES.
148 !  METHOD: REPLACED ALL UNDEFINED VALUES BY AVERAGE OF DEFINED VALUES.
149 !
150 REAL ZS
151 INTEGER IPOP,JI,JJ,JK
152
153 ! COMPUTE AVERAGE OF DEFINED VALUES
154 ZS=0. ; IPOP=0
155 DO JK=1,IK
156   DO JJ=1,INY ; DO JI=1,INX
157     IF(PLALO(JI,JJ,JK)/=PSVAL)THEN
158       ZS=ZS+PLALO(JI,JJ,JK)
159       IPOP=IPOP+1
160     ENDIF
161   ENDDO ; ENDDO
162 ENDDO
163 ZS=ZS/(FLOAT(IPOP)+TINY(ZS))
164 !
165 ! Replace ALL UNDEFINED VALUES BY THE AVERAGE
166 DO JK=1,IK
167   DO JJ=1,INY ; DO JI=1,INX
168     IF(PLALO(JI,JJ,JK)==PSVAL) PLALO(JI,JJ,JK)=ZS
169   ENDDO ; ENDDO
170 ENDDO
171 !  
172 RETURN
173 END SUBROUTINE EXTENDLAM
174
175 END SUBROUTINE INT2LALO