Philippe 02/05/2016: moved all LIBTOOLS files in LIBTOOLS directory
[MNH-git_open_source-lfs.git] / LIBTOOLS / tools / diachro / src / DIAPRO / coupe_fordiachro.f90
1 !     ######spl
2       SUBROUTINE COUPE_FORDIACHRO(PTABI,PTABO,K)
3 !     ##########################################
4 !
5 !!****  *COUPE_FORDIACHRO* - Vertical cross-section interpolation
6 !!
7 !!    PURPOSE
8 !!    -------
9 !         Interpolates 2D vertical cross-sections within the Meso-NH 3D
10 !       arrays. Model fields, iapprpriate gridpoint altitudes as well as 
11 !       appropriate topography height are interpolated. 
12 !
13 !!**  METHOD
14 !!    ------
15 !!        The general case of a vertical cross-section along any oblique 
16 !!      direction with respect to the x-y model axes is considered. Simple
17 !!      linear interpolation is done.
18 !!
19 !!    EXTERNAL
20 !!    --------
21 !!      None
22 !!
23 !!    IMPLICIT ARGUMENTS
24 !!    ------------------
25 !!      Module MODD_COORD  : declares gridpoint coordinates (TRACE use)
26 !!       XXX,XXY,XXZ    : coordinate values for all the MESO-NH grids
27 !!       XXZS           : topography values for all the MESO_NH grids
28 !!       XDSX, XDSY     : projections on the MESO-NH cartesian axes of the XDS
29 !!                        oblique abscissa (meters), for all grid locations
30 !!
31 !!      Module MODD_GRID1  : declares grid variables (Model module)
32 !!       XZZ            : true z altitude for the current NMGRID grid location
33 !!
34 !!      Module MODN_NCAR   : defines NAM_DIRTRA_POS namelist (form. NCAR common)
35 !!       XSPVAL         : Special value
36 !!
37 !!      Module MODD_CVERT  :  Declares work arrays for vertical cross-sections
38 !!       XWORKZ         : working array for true altitude storage (all grids)
39 !!       XWZ            : working array for topography (all grids)
40 !!
41 !!      Module MODN_PARA   : Defines NAM_DOMAIN_POS namelist (form. PARA common)
42 !!       NLMAX          :  Number of points horizontally along
43 !!                         the vertical section
44 !!       Module MODD_DIM1 : contains dimensions of data arrays
45 !!        NKMAX  : z array dimension
46 !!
47 !!      Module MODD_PARAMETERS : Contains array border depths
48 !!       JPHEXT         : Horizontal external points number
49 !!       JPVEXT         : Vertical external points number
50 !!
51 !!
52 !!      Module MODD_NMGRID  : declares global variable  NMGRID
53 !!       NMGRID         : Current MESO-NH grid indicator
54 !!
55 !!    REFERENCE
56 !!    ---------
57 !!
58 !!      MESO-NH User's Manual, TRACE Post Processing sections, Version 1.0:
59 !!       + Book1: Concepts and Fundamentals, to appear in 1994;
60 !!       + Book2: Technical Reference and Flowcharts, to appear in 1994;
61 !!       + Book3: Tutorial, November 1994.
62 !!
63 !!     NCAR Graphics Technical documentation, UNIX version 3.2,
64 !!     Scientific computing division, NCAR/UCAR, Boulder, USA.
65 !!      Volume 1: Fundamentals, Vers. 1, May 1993
66 !!      Volume 2: Contouring and mapping tutorial, Vers. 2, May 1993
67 !!
68 !!
69 !!    AUTHOR
70 !!    ------
71 !!      J. Duron    * Laboratoire d'Aerologie *
72 !!
73 !!    MODIFICATIONS
74 !!    -------------
75 !!      Original       06/06/94
76 !!      Updated   PM   06/01/59
77 !-------------------------------------------------------------------------
78 !
79 !*     0.   DECLARATIONS
80 !           ------------
81 !
82 USE MODD_COORD
83 USE MODD_GRID1
84 USE MODN_NCAR
85 USE MODD_CVERT
86 USE MODD_MEMCV
87 USE MODN_PARA
88 USE MODD_PARAMETERS
89 USE MODD_NMGRID 
90 USE MODD_RESOLVCAR
91 USE MODD_TYPE_AND_LH 
92 !
93 IMPLICIT NONE
94 !
95 !*     0.1  Dummy arguments and results
96 !
97 REAL,DIMENSION(:,:)      :: PTABI    ! Input data array to be interpolated
98 REAL,DIMENSION(:)        :: PTABO    ! Returned interpolated 2D array
99 INTEGER                  ::  K       ! Model level where interpolation is done
100 !
101 !*     0.2  Local variables
102 !
103 REAL    :: ZCIINF,ZCISUP,ZCJINF,ZCJSUP,ZXX,ZYY
104 INTEGER :: IMIM1,IMI,IMJM1,IMJ,JILOOP, JI, JJ
105 INTEGER :: IIU, IJU, IKU, IKB, IKE
106 !
107 !------------------------------------------------------------------------------
108 !
109 !*     1.   PERFORMING VERTICAL INTERPOLATIONS
110 !           ----------------------------------
111 !
112 !*     1.0  Presetting array extends
113 !
114 !print *,' ++++coupe NMGRID DIRCUR ',NMGRID,CDIRCUR(1:LEN_TRIM(CDIRCUR))
115 IIU=NIMAX+2*JPHEXT
116 IJU=NJMAX+2*JPHEXT
117 IKU=NKMAX+2*JPVEXT
118 IKE=IKU-JPVEXT
119 IKB=1+JPVEXT
120 ! Oct 2000 prise en compte du 2D horizontal -> PH=CH+CV
121 ! Ajout LKCP Avril 2001 -> prise en compte du 3D compresse en K
122 IF(LCHXY .OR. LKCP)THEN
123   IKU=1; IKB=1; IKE=1
124   if(nverbia > 0)then
125     print *,' **coupe NKL NKH LCHXY ',NKL,NKH,LCHXY
126   endif
127 ENDIF
128 !
129 !*     1.1  Scans along-section  oblique x axis 
130 !
131 DO JILOOP=1,NLMAX   ! Start of general X- scanning loop
132 !
133 !*     1.2  Locates the current gridpoint along the cross-section 
134 !*          oblique x-axis within the Meso-NH grid
135 !
136   ZXX=XDSX(JILOOP,NMGRID)   ! Collects the model X- and Y- axes projections
137   ZYY=XDSY(JILOOP,NMGRID)   ! onto the oblique vertical section plane
138                             ! for the current (new) section point
139 !
140 !*     1.3   The current section point is located along
141 !*           the x- axis
142 !
143     DO JI=2,IIU
144       IF(ZXX.LE.XXX(JI,NMGRID).AND.ZXX.GE.XXX(JI-1,NMGRID))GO TO 1
145     ENDDO
146
147 1 CONTINUE
148
149   IMIM1=MAX(1,JI-1)   
150   IMI=JI             ! JI is the index of the first model bin to the left
151 !
152 !*    1.4    Then, it is located along
153 !*           the Y- axis
154 !
155     DO JJ=2,IJU
156       IF(ZYY.LE.XXY(JJ,NMGRID).AND.ZYY.GE.XXY(JJ-1,NMGRID))GO TO 2
157     ENDDO
158
159 2 CONTINUE
160
161   IMJM1=MAX(1,JJ-1)
162   IMJ=JJ             ! JJ is the index of the first model bin below
163 !
164 !*   1.5     Finally the X- and Y- distances between the current section 
165 !*           point and closest model box to the left-bottom are calculated
166 !
167   IF(IMI==IMIM1)THEN       ! Left wall special case
168     ZCIINF=0.
169     ZCISUP=0.
170   ELSE
171     !print *,'XX(IMI IMIM1) ZXX',XX(IMI),XX(IMIM1),ZXX
172     ZCIINF=(XXX(IMI,NMGRID)-ZXX)/MAX(1.E-10,(XXX(IMI,NMGRID)-XXX(IMIM1,NMGRID)))
173     ZCISUP=(ZXX-XXX(IMIM1,NMGRID))/MAX(1.E-10,(XXX(IMI,NMGRID)-XXX(IMIM1,NMGRID)))
174   END IF
175   !
176   IF(IMJ==IMJM1)THEN      ! Bottom wall special case
177     ZCJINF=0.
178     ZCJSUP=0.
179   ELSE
180     !PRINT *,'XY(IMJ IMJM1) ZXY',XY(IMJ),XY(IMJM1),ZYY
181     ZCJINF=(XXY(IMJ,NMGRID)-ZYY)/MAX(1.E-10,(XXY(IMJ,NMGRID)-XXY(IMJM1,NMGRID)))
182     ZCJSUP=(ZYY-XXY(IMJM1,NMGRID))/MAX(1.E-10,(XXY(IMJ,NMGRID)-XXY(IMJM1,NMGRID)))
183   END IF
184 !
185 !*   1.6     Computes the interpolated altitude of the 
186 !*           current section point
187 !
188 if(nverbia > 1)then
189 print *,' ** coupe   AV XWORKZ  K= ',K,' size XWORKZ ',size(XWORKZ,1),size(XWORKZ,2),size(XWORKZ,3), &
190 ' IMIM1,IMJM1,IMJ,IMI ',IMIM1,IMJM1,IMJ,IMI
191 print *,' ** coupe   AV XWORKZ  ZCIINF,ZCJINF,ZCISUP,ZCJSUP,NMGRID ',ZCIINF,ZCJINF,ZCISUP,ZCJSUP,NMGRID
192 print *,' ** coupe   AV XWORKZ JILOOP XZZ(IMIM1,IMJM1,K),XZZ(IMIM1,IMJ,K),XZZ(IMI,IMJM1,K),XZZ(IMI,IMJ,K) ',&
193 JILOOP,XZZ(IMIM1,IMJM1,K),XZZ(IMIM1,IMJ,K),XZZ(IMI,IMJM1,K),XZZ(IMI,IMJ,K)
194 endif
195
196   XWORKZ(JILOOP,K,NMGRID)=ZCIINF*ZCJINF*XZZ(IMIM1,IMJM1,K)+  &
197          ZCIINF*ZCJSUP*XZZ(IMIM1,IMJ,K)+                     &
198          ZCISUP*ZCJINF*XZZ(IMI,IMJM1,K)+                     &
199          ZCISUP*ZCJSUP*XZZ(IMI,IMJ,K)    
200
201 if(nverbia > 1)then
202 print *,' ** coupe   AP XWORKZ  K= ',K,' size XZZ ',size(XZZ,1),size(XZZ,2),size(XZZ,3),' IMIM1,IMJM1,IMJ,IMI ',IMIM1,IMJM1,IMJ,IMI
203 endif
204 !
205 !*   1.7     Computes the interpolated value of the field for
206 !*           current section point
207
208 !  Modifs for diachro
209 ! Avril 2001 Ajout LKCP -> prise en compte 3D compresse sur K
210   IF((K.LT.MAX(NKL,IKB).OR.K.GT.MIN(NKH,IKE)) .AND. .NOT.LKCP)THEN
211     PTABO(JILOOP)=XSPVAL
212   ELSE
213 ! Ajout pour les PH definis avec _CV__K_  ou _Z_ etc... le 10/3/99
214 !   IF(LCV .AND. LCH)THEN
215 ! idem (02/04/04) pour les _CV_ classiques (cas obs2mesonh)
216 ! Ds ce cas on ne travaille pas necessairemnet sur les niveaux du modele
217 ! mais sur un plan Z ou PR ou TK ou EV qui peut contenir des valeurs speciales
218 ! XSPVAL. Il faut donc en tenir compte en limite de relief
219 ! A PEAUFINER
220 ! Le calcul des altitudes et du relief est fait mais je ne m'en sers pas
221 ! Il peut etre aberrant . PENSER a les eliminer avec LPRINT et LPRINTXY
222       IF((PTABI(IMIM1,IMJM1)==XSPVAL .AND. PTABI(IMIM1,IMJ)==XSPVAL).OR.&
223          (PTABI(IMIM1,IMJM1)==XSPVAL .AND. PTABI(IMI,IMJM1)==XSPVAL).OR.&
224          (PTABI(IMIM1,IMJM1)==XSPVAL .AND. PTABI(IMI,IMJ)==XSPVAL).OR. &
225          (PTABI(IMI,IMJM1)==XSPVAL .AND. PTABI(IMI,IMJ)==XSPVAL).OR. &
226          (PTABI(IMIM1,IMJ)==XSPVAL .AND. PTABI(IMI,IMJ)==XSPVAL).OR. &
227          (PTABI(IMIM1,IMJ)==XSPVAL .AND. PTABI(IMI,IMJM1)==XSPVAL))THEN
228          PTABO(JILOOP)=XSPVAL
229       ELSE IF(PTABI(IMIM1,IMJM1)==XSPVAL .AND. PTABI(IMIM1,IMJ)/=XSPVAL.AND.&
230               PTABI(IMI,IMJM1)/=XSPVAL .AND. PTABI(IMI,IMJ)/=XSPVAL)THEN
231         PTABO(JILOOP)=                                   &
232           ZCIINF*ZCJSUP*PTABI(IMIM1,IMJ)+            &
233           ZCISUP*ZCJINF*PTABI(IMI,IMJM1)+            &
234           ZCISUP*ZCJSUP*PTABI(IMI,IMJ)    
235       ELSE IF(PTABI(IMIM1,IMJM1)/=XSPVAL .AND. PTABI(IMIM1,IMJ)==XSPVAL.AND.&
236               PTABI(IMI,IMJM1)/=XSPVAL .AND. PTABI(IMI,IMJ)/=XSPVAL)THEN
237         PTABO(JILOOP)=ZCIINF*ZCJINF*PTABI(IMIM1,IMJM1)+  &
238           ZCISUP*ZCJINF*PTABI(IMI,IMJM1)+            &
239           ZCISUP*ZCJSUP*PTABI(IMI,IMJ)    
240       ELSE IF(PTABI(IMIM1,IMJM1)/=XSPVAL .AND. PTABI(IMIM1,IMJ)/=XSPVAL.AND.&
241               PTABI(IMI,IMJM1)==XSPVAL .AND. PTABI(IMI,IMJ)/=XSPVAL)THEN
242         PTABO(JILOOP)=ZCIINF*ZCJINF*PTABI(IMIM1,IMJM1)+  &
243           ZCIINF*ZCJSUP*PTABI(IMIM1,IMJ)+            &
244           ZCISUP*ZCJSUP*PTABI(IMI,IMJ)    
245       ELSE IF(PTABI(IMIM1,IMJM1)/=XSPVAL .AND. PTABI(IMIM1,IMJ)/=XSPVAL.AND.&
246               PTABI(IMI,IMJM1)/=XSPVAL .AND. PTABI(IMI,IMJ)==XSPVAL)THEN
247         PTABO(JILOOP)=ZCIINF*ZCJINF*PTABI(IMIM1,IMJM1)+  &
248           ZCIINF*ZCJSUP*PTABI(IMIM1,IMJ)+            &
249           ZCISUP*ZCJINF*PTABI(IMI,IMJM1)
250       ELSE
251         PTABO(JILOOP)=ZCIINF*ZCJINF*PTABI(IMIM1,IMJM1)+  &
252           ZCIINF*ZCJSUP*PTABI(IMIM1,IMJ)+            &
253           ZCISUP*ZCJINF*PTABI(IMI,IMJM1)+            &
254           ZCISUP*ZCJSUP*PTABI(IMI,IMJ)    
255       ENDIF
256
257 !   ELSE
258 ! Cas habituel
259 !   PTABO(JILOOP)=ZCIINF*ZCJINF*PTABI(IMIM1,IMJM1)+  &
260 !      ZCIINF*ZCJSUP*PTABI(IMIM1,IMJ)+            &
261 !      ZCISUP*ZCJINF*PTABI(IMI,IMJM1)+            &
262 !      ZCISUP*ZCJSUP*PTABI(IMI,IMJ)    
263 !   ENDIF
264   END IF
265 !
266 !*   1.8     Computes the interpolated topography height for
267 !*           current section point
268 !
269 if(nverbia > 1)then
270   print *,' ** coupe   AV XWZ '
271 endif
272
273 XWZ(JILOOP,NMGRID)=ZCIINF*ZCJINF*XXZS(IMIM1,IMJM1,NMGRID)+  &
274           ZCIINF*ZCJSUP*XXZS(IMIM1,IMJ,NMGRID)+             &
275           ZCISUP*ZCJINF*XXZS(IMI,IMJM1,NMGRID)+             &
276           ZCISUP*ZCJSUP*XXZS(IMI,IMJ,NMGRID)    
277 !
278 ENDDO                     ! End of the general X- scanning loop
279
280 if(nverbia > 0)then
281 print *,' >>>> SORTIE COUPE  NMGRID= ',NMGRID,' size(XWZ)' ,size(XWZ,1)
282 endif
283 if(nverbia > 1)then
284 print *,' XWZ ',XWZ(:,NMGRID)
285 endif
286 !
287 RETURN
288 !------------------------------------------------------------------------
289 !
290 !*   2.     EXIT
291 !           ----
292 !
293 END SUBROUTINE COUPE_FORDIACHRO