Changeset 25
- Timestamp:
- 11/28/07 17:04:24 (17 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/procs/meshes/mesh_gaussian.pro
r2 r25 5 5 ; init grid, sf, masks for gaussian grid for SINTEX 6 6 ; 7 print,' T',trunc, ' gaussian grid inits.',format='(A5,I3,A20)'7 print,' T',trunc, ' gaussian grid inits.',format='(A5,I3,A20)' 8 8 9 ; use key_shiftg if grid and data file are not organized the same 10 ; key_shiftg(0) = longitude shift 11 ; t106 for SQs is jpi/2 12 ; key_shiftg(1) = -1 if latitude reversing needed 9 key_shift = 0 13 10 14 key_shift = 0 15 key_shiftg = [0, 0] 11 ; initialisation of character variables used in the execution of initncdf 12 ; (and computegrid as a consequence) 13 shift_txt = '' 14 plain_txt = '' 15 16 IF keyword_set(NO_SHIFT) THEN BEGIN 17 shift_txt = ', SHIFT = 0' 18 ENDIF 19 IF keyword_set(WHOLE_ARRAYS) THEN BEGIN 20 plain_txt = ', PLAIN = 1' 21 ;; Force YREVERSE = 0, ZREVERSE = 0, PERIODIC = 0, SHIFT = 0, STRIDE = [1, 1, 1] and 22 ;; suppress the automatic redefinition of the domain in case of x periodicity overlap, 23 ;; y periodicity overlap (ORCA type only) and mask border to 0. 24 ENDIF 25 26 IF debug_w THEN print, ' key_yreverse = ', key_yreverse 16 27 17 28 CASE trunc OF … … 43 54 44 55 glamt = 360.0*findgen(jpi)/float(jpi) 45 glamt = glamt#replicate(1, jpj)46 56 47 57 ; 2 Define latitudes … … 135 145 ENDCASE 136 146 137 gphit = replicate(1, jpi)#gphit138 139 147 ; 4. Define mask 140 148 … … 168 176 END 169 177 ENDCASE 170 171 IF key_shiftg(0) NE 0 THEN BEGIN172 glamt = shift(glamt, key_shiftg(0), 0)173 glamt(0:key_shiftg(0)-1, *) = glamt(0:key_shiftg(0)-1, *) - 360.0174 gphit = shift(gphit, key_shiftg(0), 0)175 tmask = shift(tmask, key_shiftg(0), 0)176 ENDIF177 178 IF key_shiftg(1) EQ -1 THEN BEGIN179 glamt = reverse(glamt, 2)180 gphit = reverse(gphit, 2)181 tmask = reverse(tmask, 2)182 ENDIF183 178 ; 184 ; definition of key_shift179 ; compute grid 185 180 ; 186 if keyword_set(glamboundary) AND jpi GT 1 then BEGIN 187 xaxis = glamt[*, 0] 188 case 1 OF 189 glamboundary[0] LT xaxis[0]:BEGIN 190 toobig = where(xaxis GT glamboundary[1]) 191 if toobig[0] NE -1 then begin 192 key_shift = n_elements(toobig) 193 toobig = where(glamt GT glamboundary[1]) 194 glamt[toobig] = glamt[toobig]-360 195 endif ELSE key_shift = 0 196 END 197 glamboundary[1] GT xaxis[jpi-1]:BEGIN 198 toosmall = where(xaxis LT glamboundary[0]) 199 if toosmall[0] NE -1 then begin 200 key_shift = -n_elements(toosmall) 201 toosmall = where(glamt LT glamboundary[0]) 202 glamt[toosmall] = glamt[toosmall]+360 203 ENDIF ELSE key_shift = 0 204 END 205 ELSE:key_shift = 0 206 endcase 207 ENDIF ELSE key_shift = 0 208 209 IF keyword_set(NO_SHIFT) THEN key_shift = 0 210 IF keyword_set(WHOLE_ARRAY) THEN key_shift = 0 211 212 print, ' key_shift =', key_shift 213 214 glamt = shift(glamt,key_shift, 0) 215 gphit = shift(gphit,key_shift, 0) 216 tmask = shift(tmask, key_shift, 0) 217 218 ; 3 Define scale factors 219 220 zrad=6371229.0 221 222 ; Exact method: integration on a sphere 223 224 e1t = abs(2*!pi*zrad*cos(gphit*!pi/180.0)/jpi) 225 226 IF gphit(1, 0) GT gphit(0, 0) THEN BEGIN 227 ytop = (shift(gphit, 0, -1)+gphit)*0.5*!pi/180.0 228 ybot = (shift(gphit, 0, +1)+gphit)*0.5*!pi/180.0 229 230 ytop(*, jpj-1) = !pi/2 231 ybot(*, 0) = -!pi/2 232 233 e2t = abs(zrad*(ytop-ybot)) 234 ENDIF ELSE BEGIN 235 ybot = (shift(gphit, 0, -1)+gphit)*0.5*!pi/180.0 236 ytop = (shift(gphit, 0, +1)+gphit)*0.5*!pi/180.0 237 238 ybot(*, jpj-1) = !pi/2 239 ytop(*, 0) = -!pi/2 240 241 e2t = abs(zrad*(ybot-ytop)) 242 ENDELSE 243 244 ; 5. vertical grid (hPa) 245 246 ; gdept = reverse(gdept) 247 gdepw = gdept 248 249 e3t = shift(gdept, 1)-gdept 250 e3t[0] = e3t[1] 251 252 253 key_periodique=1 254 255 glamu = glamt 256 gphiu = gphit 257 e1u = e1t 258 e2u = e2t 259 260 glamv = glamt 261 gphiv = gphit 262 e1v = e1t 263 e2v = e2t 264 265 ; used for coast plot 266 267 glamf = 0.5*(shift(glamt, -1, 0)+glamt) 268 glamf(jpi-1, *) = glamf(jpi-2, *) + (glamf(jpi-2, *)-glamf(jpi-3, *)) 269 270 IF gphit(1, 0) GT gphit(0, 0) THEN BEGIN 271 gphif = 0.5*(shift(gphit, 0, -1)+gphit) 272 gphif(*, jpj-1) = gphif(*, jpj-2) + (gphif(*, jpj-2)-gphif(*, jpj-3)) 273 ENDIF ELSE BEGIN 274 gphif = 0.5*(shift(gphit, 0, 1)+gphit) 275 gphif(*, 0) = gphif(*, 1) + (gphif(*, 2)-gphif(*, 1)) 276 ENDELSE 277 278 e1f = e1t 279 e2f = e2t 280 281 e3w = e3t 282 181 computegrid, xaxis = glamt, yaxis = gphit, mask = tmask, GLAMBOUNDARY = glamboundary, /periodic 182 283 183 key_offset = [0, 0, 0] 284 285 masked_data = 0 286 mesh_type = 'atm' 287 ; 288 ; indice i pour grille j moyenne zonale 289 ; 290 diaznl_idx = 1 291 292 tmask = reform(tmask, jpi*jpj) 293 tmask = reform(tmask#replicate(1, jpk), jpi, jpj, jpk) 184 185 print,' End of initialisation for gaussian grid T', trunc 294 186 295 187 return
Note: See TracChangeset
for help on using the changeset viewer.