Oasis3 4.0.2
|
00001 SUBROUTINE pminm (kflag, psurf, 00002 $ p2xi, p2xs, p2yi, p2ys, 00003 $ p1xi, p1xs, p1yi, p1ys) 00004 C**** 00005 C ***************************** 00006 C * OASIS ROUTINE - LEVEL 3 * 00007 C * ------------- ------- * 00008 C ***************************** 00009 C 00010 C**** *pminm* - calculate intersection between two meshes 00011 C 00012 C Purpose: 00013 C ------- 00014 C For two meshes defined by their inf. and sup. latitude and 00015 C longitude (p1xi,p1xs,p1yi,p1ys) and (p2xi,p2xs,p2yi,p2ys), checks 00016 C if they have any intersection and calculates the surface, psurf, 00017 C of that intersection. The 2 grids considered may be periodic 00018 C or non-periodic. 00019 C 00020 C N.B: Intersection surface is given by the following formula: 00021 C Si = r*r * (sin(lat2)-sin(lat1)) * (phi2-phi1) 00022 C with lat2, lat1 sup. and inf. latitudes 00023 C phi2, phi1 sup. and inf. longitudes 00024 C This results from the integration of dS=r*r*cos(lat)*d(lat)*d(phi) 00025 C 00026 C** Interface: 00027 C --------- 00028 C *CALL* *pminm(kflag, psurf, 00029 C p2xi, p2xs, p2yi, p2ys, 00030 C p1xi, p1xs, p1yi, p1ys)* 00031 C 00032 C Input: 00033 C ----- 00034 C p2xi : target grid square inferior longitude 00035 C p2xs : target grid square superior longitude 00036 C p2yi : target grid square inferior latitude 00037 C p2ys : target grid square superior latitude 00038 C p1xi : source grid square inferior longitude 00039 C p1xs : source grid square superior longitude 00040 C p1yi : source grid square inferior latitude 00041 C p1ys : source grid square superior latitude 00042 C 00043 C Output: 00044 C ------ 00045 C psurf : the intersection's surface 00046 C kflag : 1 if intersection is not empty, 0 otherwise 00047 C 00048 C Workspace: 00049 C --------- 00050 C None 00051 C 00052 C External: 00053 C -------- 00054 C None 00055 C 00056 C References: 00057 C ---------- 00058 C O. Thual, Simple ocean-atmosphere interpolation. 00059 C Part A: The method, EPICOA 0629 (1992) 00060 C Part B: Software implementation, EPICOA 0630 (1992) 00061 C See also OASIS manual (1995) 00062 C 00063 C History: 00064 C ------- 00065 C Version Programmer Date Description 00066 C ------- ---------- ---- ----------- 00067 C 1.1 O. Thual 93/04/15 created 00068 C 2.0 L. Terray 95/10/01 modified: new structure 00069 C 00070 C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 00071 C 00072 C* ---------------------------- Include files --------------------------- 00073 C 00074 USE mod_kinds_oasis 00075 USE mod_unit 00076 C 00077 C* ---------------------------- Poema verses ---------------------------- 00078 C 00079 C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 00080 C 00081 C* 1. Initialization 00082 C -------------- 00083 C 00084 kflag = 0 00085 psurf = 0. 00086 z2pi = 2. * acos(-1.) 00087 zconv = 1.74532925199432957692e-2 00088 C 00089 C* Assign local variables 00090 C 00091 zfi1i = p1xi 00092 zfi1s = p1xs 00093 zfi2i = p2xi 00094 zfi2s = p2xs 00095 C 00096 C* Put longitudes in the [0-360] interval to avoid problem in 00097 C surface calculation. 00098 C 00099 IF (p1xi .GE. 360.0) zfi1i = zfi1i - 360.0 00100 IF (p1xi .LT. 0. ) zfi1i = zfi1i + 360.0 00101 IF (p1xs .GT. 360.0) zfi1s = zfi1s - 360.0 00102 IF (p1xs .LE. 0. ) zfi1s = zfi1s + 360.0 00103 IF (p2xi .GE. 360.0) zfi2i = zfi2i - 360.0 00104 IF (p2xi .LT. 0. ) zfi2i = zfi2i + 360.0 00105 IF (p2xs .GT. 360.0) zfi2s = zfi2s - 360.0 00106 IF (p2xs .LE. 0. ) zfi2s = zfi2s + 360.0 00107 C 00108 C* Conversion degrees to radians 00109 C 00110 zfi1i = zfi1i * zconv 00111 zfi1s = zfi1s * zconv 00112 zth1i = p1yi * zconv 00113 zth1s = p1ys * zconv 00114 zfi2i = zfi2i * zconv 00115 zfi2s = zfi2s * zconv 00116 zth2i = p2yi * zconv 00117 zth2s = p2ys * zconv 00118 C 00119 C 00120 C* 2. Check and calculation of intersection 00121 C ------------------------------------- 00122 C 00123 C* Check of latitudinal intersection 00124 C 00125 IF (zth2s .GT. zth1i .AND. zth2i .LT. zth1s) THEN 00126 C 00127 C* Calculate the delta(latitude) 00128 C 00129 C* Latitude inf. 00130 C 00131 zlati = amax1(zth2i,zth1i) 00132 C 00133 C* Latitude sup. 00134 C 00135 zlats = amin1(zth2s,zth1s) 00136 C 00137 C* sin(lat_sup) - sin(lat_inf) 00138 C 00139 zdth = sin(zlats) - sin(zlati) 00140 C 00141 C* Check of longitudinal intersection 00142 C 00143 zfi1ir = zfi1i - z2pi 00144 zfi2ir = zfi2i - z2pi 00145 zfi2sr = zfi2s - z2pi 00146 C 00147 C* Case of a mesh 1 whose minimal longitude zfi1i is less than the 00148 C maximal one zfi1s (normal case). The opposite can occur only for 00149 C a periodic grid. 00150 C 00151 IF (zfi1i .LT. zfi1s) THEN 00152 C 00153 C - Check if it is similar for mesh 2 00154 C 00155 IF (zfi2i .LT. zfi2s) THEN 00156 C 00157 C - If yes, find the intersection 00158 C 00159 IF (zfi2s .GT. zfi1i .AND. zfi2i .LT. zfi1s) THEN 00160 kflag = 1 00161 zdfi = amin1(zfi2s,zfi1s) - amax1(zfi2i,zfi1i) 00162 psurf = zdth * zdfi 00163 ENDIF 00164 C 00165 C - If not, mesh 2 is periodic 00166 C 00167 ELSE IF (zfi2i .GT. zfi2s) THEN 00168 C 00169 C - Find intersection in that case 00170 C 00171 IF (zfi2s .GT. zfi1i .AND. zfi2ir .LT. zfi1s) THEN 00172 kflag = 1 00173 zdfi = amin1(zfi2s,zfi1s) - zfi1i 00174 psurf = zdth * zdfi 00175 ENDIF 00176 IF (zfi2s .GT. zfi1ir .AND. zfi2i .LT. zfi1s) THEN 00177 kflag = 1 00178 zdfi = amin1(zfi2i,zfi1s) - zfi1i 00179 psurf = zdth * zdfi 00180 ENDIF 00181 ENDIF 00182 ENDIF 00183 C 00184 C* Case of a mesh 1 whose minimal longitude zfi1i is located 00185 C on the right of the grid and whose maximum longitude zfi1s 00186 C is on the left of the grid. This occurs only for a periodic mesh 00187 C 00188 IF (zfi1i .GT. zfi1s) THEN 00189 C 00190 C - Check for mesh 2 00191 C 00192 IF (zfi2i .LT. zfi2s) THEN 00193 C 00194 C - Find intersection 00195 C 00196 IF (zfi2s .GT. zfi1ir .AND. zfi2i .LT. zfi1s) THEN 00197 kflag = 1 00198 zdfi = amin1(zfi2s,zfi1s) - zfi2i 00199 psurf = zdth * zdfi 00200 ELSE IF (zfi2sr .GT. zfi1ir .AND. 00201 $ zfi2ir .LT. zfi1s) THEN 00202 kflag = 1 00203 zdfi = zfi2s - amin1(zfi2i,zfi1i) 00204 psurf = zdth * zdfi 00205 ENDIF 00206 C 00207 C - Mesh 2 is periodic 00208 C 00209 ELSE IF (zfi2i .GT. zfi2s) THEN 00210 C 00211 C - Find intersection 00212 C 00213 IF (zfi2s .GT. zfi1ir .AND. zfi2ir .LT. zfi1s) THEN 00214 kflag = 1 00215 zdfi = amin1(zfi2s,zfi1s) - amax1(zfi2ir,zfi1ir) 00216 psurf = zdth * zdfi 00217 ENDIF 00218 ENDIF 00219 ENDIF 00220 ENDIF 00221 C 00222 C* End of routine 00223 C 00224 RETURN 00225 END 00226