Philippe 02/05/2016: moved all LIBTOOLS files in LIBTOOLS directory
[MNH-git_open_source-lfs.git] / LIBTOOLS / tools / diachro / src / EXTRACTDIA / ini2lalo.f90
1 !-----------------------------------------------------------------
2 !--------------- special set of characters for SCCS information
3 !-----------------------------------------------------------------
4 !      @(#) Lib:./s.ini2lalo.f90, Version:1.5, Date:03/06/05, Last modified:01/10/23
5 !-----------------------------------------------------------------
6 !     ######spl
7 MODULE MODI_INI2LALO
8 !###################
9 !
10 INTERFACE
11       SUBROUTINE INI2LALO(PLATLON,KNX,KNY,            &
12                           KIDEB,KIFIN,KJDEB,KJFIN,PDLON,PDLAT)
13 !
14 REAL, DIMENSION(4), INTENT(INOUT) :: PLATLON ! NSWE target domain bounds (deg)
15 INTEGER, INTENT(OUT) :: KNX,KNY        ! NUMBER OF TARGET POINTS IN X,Y
16 INTEGER, INTENT(IN), OPTIONAL :: KIDEB,KIFIN ! limites du
17 INTEGER, INTENT(IN), OPTIONAL :: KJDEB,KJFIN !zoom eventuel
18 REAL,    INTENT(OUT),OPTIONAL :: PDLON,PDLAT ! resolutions in LOn-LAt computed
19 !
20 END SUBROUTINE INI2LALO
21 END INTERFACE
22 END MODULE MODI_INI2LALO
23 !     ########################################################
24       SUBROUTINE INI2LALO(PLATLON,KNX,KNY,            &
25                           KIDEB,KIFIN,KJDEB,KJFIN,PDLON,PDLAT)
26 !     ########################################################
27 !
28 !!    PURPOSE
29 !!    -------
30 !
31 !-------------------------------------------------------------------------------
32 !
33 !*       0.    DECLARATIONS
34 !              ------------
35 !
36 USE MODD_CST,        ONLY : XRADIUS,XPI
37 USE MODD_PARAMETERS, ONLY : JPHEXT
38 USE MODD_GRID,       ONLY : XLAT0,XLATORI,XLONORI
39 USE MODD_DIM1,       ONLY : NIMAX,NJMAX
40 USE MODD_GRID1,      ONLY : XLAT,XMAP,XXHAT,XYHAT
41 !
42 USE MODE_GRIDPROJ
43 !
44 IMPLICIT NONE
45 !
46 !*       0.1   Arguments
47 !
48 REAL, DIMENSION(4), INTENT(INOUT) :: PLATLON ! NSWE target domain bounds (millideg)
49 INTEGER,            INTENT(OUT)   :: KNX,KNY ! nb of target points
50 INTEGER, INTENT(IN), OPTIONAL :: KIDEB,KIFIN ! limites du
51 INTEGER, INTENT(IN), OPTIONAL :: KJDEB,KJFIN !zoom eventuel
52 REAL   , INTENT(OUT),OPTIONAL :: PDLON,PDLAT ! resolutions in LOn-LAt computed
53 !
54 !*       0.2   Local variables
55 !
56 REAL :: ZLONW,ZLONE,ZLATN,ZLATS  ! LAT/LON rounded to nearest millidegree
57 REAL :: ZDX,ZDY                  ! increments in LAT/LON
58 REAL :: ZLATM                    ! extreme latitude of input domain
59 REAL :: ZLA,ZLO,ZI,ZJ,ZXHAT,ZYHAT
60 INTEGER, DIMENSION(2) :: IMAP
61 INTEGER :: II,IJ,IIN
62 INTEGER :: JX,JY
63 !
64 !------------------------------------------------------------------------------
65 !
66 !*       1.    CHECK Lat/Lon DOMAIN
67 !              --------------------
68 !
69 ZLATN=PLATLON(1) ; ZLATS=PLATLON(2)
70 ZLONW=PLATLON(3) ; ZLONE=PLATLON(4) 
71 !
72 ! round to nearest millidegree, longitudes in (0..360) interval
73 ZLATN=REAL(NINT(ZLATN))
74 ZLATS=REAL(NINT(ZLATS))
75 ZLONW=MOD(ZLONW,360000.) ; ZLONE=MOD(ZLONE,360000.)
76 ZLONW=REAL(NINT(ZLONW))
77 ZLONE=REAL(NINT(ZLONE))
78 PLATLON(1)=ZLATN
79 PLATLON(2)=ZLATS
80 PLATLON(3)=ZLONW
81 PLATLON(4)=ZLONE
82 !
83 ! check if domain is well-defined
84 IF(ABS(ZLATN)>90000.) THEN
85   PRINT*, 'INI2LALO: Bad N latitude - abort: ZLATN=',ZLATN
86   STOP
87 END IF
88 IF(ABS(ZLATS)>90000) THEN
89   PRINT*, 'INI2LALO: Bad S latitude - abort: ZLATS=',ZLATS
90   STOP
91 END IF
92 IF(ZLATN<=ZLATS) THEN
93   PRINT*, 'Bad latitude interval - abort'
94   STOP
95 END IF
96 !
97 ! compute optimum resolution
98 IF (PRESENT(KIDEB)) THEN
99   IMAP=MINLOC(XMAP(KIDEB:KIFIN,KJDEB:KJFIN)) 
100 ELSE
101   IMAP=MINLOC(XMAP(:,:)) 
102 END IF
103 ZLATM=XLAT(IMAP(1),IMAP(2))
104 ZDX=(XXHAT(3)-XXHAT(2))*180./XPI&
105                        /(XRADIUS*COS(ZLATM*XPI/180.))
106 ZDY=(XYHAT(3)-XYHAT(2))*180./XPI/XRADIUS
107 print*, 'INI2LALO: equivalent resolution in lat ',ZLATM,IMAP
108 print*, '        ',ZDX,ZDY
109 WRITE(6,'(A,I4,1X,I4,A,F6.1,A)')'INI2LALO: point where map scale is minimum ', &
110                                IMAP,' (lat ',ZLATM,')'
111 PRINT*,'equivalent resolution in lon. and lat.: ' ,ZDX,ZDY
112 !
113 ! compute number of points and
114 ! move E & S boundaries so that lon/lat increments are in millidegrees
115 ! (GRIB constraint)
116 KNX=NINT( (ZLONE-ZLONW)/(1000*ZDX) +1)
117 IF(KNX<0) KNX=NINT( (ZLONE-ZLONW+360000.)/(1000*ZDX) +1)
118 IF(ZDX/=REAL(NINT(ZDX*1000.))/1000.) THEN  ! need to fix longitude
119   ZDX=REAL(NINT(ZDX*1000.))
120   ZLONE=ZLONW+(KNX-1)*ZDX
121   IF(ZLONE>360000.) ZLONE=ZLONE-360000.
122   print*, 'INI2LALO: fixing E longitude to ',ZLONE
123   PLATLON(4)=ZLONE
124   ZDX=ZDX/1000.
125 ENDIF
126 !
127 KNY=NINT( (ZLATN-ZLATS)/(1000*ZDY) +1)
128 IF(ZDY/=REAL(NINT(ZDY*1000.))/1000.) THEN  ! need to fix latitude
129   ZDY=REAL(NINT(ZDY*1000.))
130   ZLATS=ZLATN-(KNY-1)*ZDY
131   IF(ABS(ZLATS)>90000.) THEN
132     STOP "TOO BIG DOMAIN in LATITUDE"
133   END IF
134   print*, 'INI2LALO: fixing S latitude to ',ZLATS
135   PLATLON(2)=ZLATS
136     ZDY=ZDY/1000.
137 ENDIF
138 !
139 IF(PRESENT(PDLON))THEN
140   PDLON=ZDX
141   PDLAT=ZDY
142 END IF
143 !
144 print*, 'INI2LALO: number of points in lon. and lat. domain:', KNX,KNY
145 IF (PRESENT(KIDEB)) THEN
146   PRINT*, 'number of points of input domain (i,j,i*j): ',(KIFIN-KIDEB+1),&
147          (KJFIN-KJDEB+1),(KIFIN-KIDEB+1)*(KJFIN-KJDEB+1)
148 ELSE
149   PRINT*, 'number of points of input domain (i,j,i*j): ',NIMAX,NJMAX,NIMAX*NJMAX
150 END IF
151 PRINT*, 'number of points of lon.-lat. domain(x,y,x*y):', KNX,KNY,KNX*KNY
152 !
153 ! check if target domain is inside file domain
154 IIN=0
155 DO JY=1,KNY ; DO JX=1,KNX
156   ZLO=MOD(ZLONW/1000.+ZDX*(JX-1),360.)
157   ZLA=ZLATN/1000.-ZDY*(JY-1)          ! output has N->S scanning
158   CALL SM_XYHAT(XLATORI,XLONORI,ZLA,ZLO,ZXHAT,ZYHAT)
159   II=MAX(MIN(COUNT(XXHAT(:)<ZXHAT),NIMAX+JPHEXT),1+JPHEXT)
160   IJ=MAX(MIN(COUNT(XYHAT(:)<ZYHAT),NJMAX+JPHEXT),1+JPHEXT)
161   ZI=(ZXHAT-XXHAT(II))/(XXHAT(II+1)-XXHAT(II))+FLOAT(II)-1
162   ZJ=(ZYHAT-XYHAT(IJ))/(XYHAT(IJ+1)-XYHAT(IJ))+FLOAT(IJ)-1
163   !
164   IF (      (ZI>=1.) .AND. (ZI<=NIMAX) &
165      .AND.  (ZJ>=1.) .AND. (ZJ<=NJMAX) ) THEN
166     IIN=IIN+1   ! points inside
167   ENDIF
168 END DO ; END DO
169 PRINT*, 'number of points of lon.-lat. domain inside input file one:', IIN
170 !
171 !
172 END SUBROUTINE INI2LALO