Oasis3 4.0.2
|
00001 SUBROUTINE qgrho (prho, kto, kwg, 00002 $ px, py, 00003 $ px1, py1, kmsk1, kngx1, kngy1, 00004 $ psig, kvma1) 00005 C**** 00006 C ***************************** 00007 C * OASIS ROUTINE - LEVEL 3 * 00008 C * ------------- ------- * 00009 C ***************************** 00010 C 00011 C**** *qgrho* - Calculate weights and adresses at one location 00012 C 00013 C Purpose: 00014 C ------- 00015 C At one location xgr gives the kwg closest neighbours 00016 C adresses kto in grid 1 and their weight, which is function of 00017 C the distance and may be other grid dependant considerations. 00018 C 00019 C** Interface: 00020 C --------- 00021 C *CALL* *qgrho(prho, kto, kwg, 00022 C px, py, 00023 C px1, py1, kmsk1, kngx1, kngy1, 00024 C psig, kvma1)* 00025 C 00026 C Input: 00027 C ----- 00028 C kwg : maximum number of overlapped neighbors 00029 C px : longitude of the current point on target grid 00030 C py : latitude of the current point on source grid 00031 C px1 : longitudes for source grid (real 2D) 00032 C py1 : latitudes for source grid (real 2D) 00033 C kmsk1 : the mask for source grid (integer 2D) 00034 C kngx1 : number of longitudes for source grid 00035 C kngy1 : number of latitudes for source grid 00036 C kvma1 : the value of the mask for source grid 00037 C psig : variance of the gaussian 00038 C Output: 00039 C ------ 00040 C prho : the neighbors weights 00041 C kto : the neighbors adresses in the source grid 00042 C 00043 C Workspace: 00044 C --------- 00045 C None 00046 C 00047 C External: 00048 C -------- 00049 C qlsort, sqdis, qlins, qlgaus 00050 C 00051 C References: 00052 C ---------- 00053 C O. Thual, Simple ocean-atmosphere interpolation. 00054 C Part A: The method, EPICOA 0629 (1992) 00055 C Part B: Software implementation, EPICOA 0630 (1992) 00056 C See also OASIS manual (1995) 00057 C 00058 C History: 00059 C ------- 00060 C Version Programmer Date Description 00061 C ------- ---------- ---- ----------- 00062 C 1.1 O. Thual 93/04/15 created 00063 C 2.0 L. Terray 95/10/01 modified: new structure 00064 C 00065 C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 00066 C 00067 C* ---------------------------- Include files --------------------------- 00068 C 00069 USE mod_unit 00070 C 00071 C* ---------------------------- Argument declarations ------------------- 00072 C 00073 REAL (kind=ip_realwp_p) px1(kngx1,kngy1), py1(kngx1,kngy1) 00074 REAL (kind=ip_realwp_p) prho(kwg) 00075 INTEGER (kind=ip_intwp_p) kmsk1(kngx1,kngy1) 00076 INTEGER (kind=ip_intwp_p) kto(kwg) 00077 C 00078 C* ---------------------------- Poema verses ---------------------------- 00079 C 00080 C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 00081 C 00082 C* 1. Various initializations 00083 C ----------------------- 00084 C 00085 C* Value of the mask 00086 C 00087 ivma = kvma1 00088 ing = kngx1 * kngy1 00089 C 00090 C 00091 C* 2. Nearest neighbors search 00092 C ------------------------ 00093 C 00094 C* Checks 00095 C 00096 IF (kwg .GT. ing) THEN 00097 WRITE (UNIT = nulou,FMT = *) 'WARNING kwg = ',kwg, 00098 $ '.GT. kngx1*kngy1 = ',ing 00099 WRITE (UNIT = nulou,FMT = *) 00100 $ ' Number of searched neighbors greater ' 00101 WRITE (UNIT = nulou,FMT = *) 00102 $ ' than number of grid points ' 00103 iwg = ing 00104 WRITE (UNIT = nulou,FMT = *) 'We set iwg = ing = ',iwg 00105 ELSE 00106 iwg = kwg 00107 ENDIF 00108 C 00109 C 00110 C* Put first iwg unmasked points of source grid in kto and prho 00111 C 00112 iind = 0 00113 DO 210 j2 = 1, kngy1 00114 DO 220 j1 = 1, kngx1 00115 IF (kmsk1(j1,j2) .NE. ivma) THEN 00116 iind = iind + 1 00117 icum = (j2-1) * kngx1 + j1 00118 prho(iind) = sqdis (px, py, px1(j1,j2), py1(j1,j2)) 00119 kto(iind) = icum 00120 i1t = j1 00121 i2t = j2 00122 IF (iind .EQ. iwg) GO TO 225 00123 ENDIF 00124 220 CONTINUE 00125 210 CONTINUE 00126 225 CONTINUE 00127 C 00128 C* Check 00129 C 00130 IF (iind .NE. iwg) THEN 00131 WRITE(UNIT = nulou,FMT = *) 00132 $ ' WARNING: iind.ne.iwg ===> ',iind, iwg 00133 WRITE(UNIT = nulou,FMT = *) 00134 $ ' ******* ===> Program must stop ' 00135 CALL HALTE ('STOP in qgrho') 00136 ENDIF 00137 C 00138 C* Sorting of these first iwg points 00139 C 00140 CALL qlsort (prho, kto, iwg) 00141 C 00142 C* Loop on the remaining points of source grid 00143 C 00144 IF (iwg .NE. ing) THEN 00145 C 00146 C* Loop over the remaining longitudes for latitude i2t on source grid 00147 C if necessary 00148 C 00149 IF (i1t .LT. kngx1) THEN 00150 DO 230 j1 = i1t+1, kngx1 00151 IF (kmsk1(j1,i2t) .NE. ivma) THEN 00152 icum = (i2t-1) * kngx1 + j1 00153 C 00154 C* Calculate distance between the current source grid point and 00155 C the given target grid point 00156 C 00157 zrnew = sqdis (px, py, px1(j1,i2t), py1(j1,i2t)) 00158 C 00159 C* Insert in list at the correct rank 00160 C 00161 CALL qlins (prho, kto, iwg, zrnew, icum) 00162 ENDIF 00163 230 CONTINUE 00164 ENDIF 00165 C 00166 C* Loop over the remaining source grid points 00167 C 00168 DO 240 j2 = i2t+1, kngy1 00169 DO 250 j1 = 1, kngx1 00170 IF (kmsk1(j1,j2) .NE. ivma) THEN 00171 icum = (j2-1) * kngx1 + j1 00172 C 00173 C* Calculate distance between the current source grid point and 00174 C the given target grid point 00175 C 00176 zrnew = sqdis (px, py, px1(j1,j2), py1(j1,j2)) 00177 C 00178 C* Insert in list at the correct rank 00179 C 00180 CALL qlins (prho, kto, iwg, zrnew, icum) 00181 ENDIF 00182 250 CONTINUE 00183 240 CONTINUE 00184 ENDIF 00185 C 00186 C 00187 C* 3. Weight function 00188 C --------------- 00189 C 00190 C* Gaussian value for the weight function 00191 C 00192 zsum = 0. 00193 DO 310 jwg = 1, iwg 00194 zr = qlgaus (prho(jwg), psig) 00195 zsum = zsum + zr 00196 prho(jwg) = zr 00197 310 CONTINUE 00198 C 00199 C* Normalization 00200 C 00201 zthres = 1.e-6 00202 C 00203 C* If under threshold 00204 C 00205 IF (zsum .LT. zthres) THEN 00206 WRITE(UNIT=nulou,FMT=*) 'Gaussian sum under threshold value' 00207 DO 320 jwg = 1, iwg 00208 prho(jwg) = 1. / float(iwg) 00209 320 CONTINUE 00210 C 00211 ELSE 00212 C 00213 C* If above threshold 00214 C 00215 DO 330 jwg = 1, iwg 00216 prho(jwg) = prho(jwg) / zsum 00217 330 CONTINUE 00218 ENDIF 00219 C 00220 C* End of routine 00221 C 00222 RETURN 00223 END