Changeset 6041 for branches/2015/dev_r5776_UKMO2_OBS_efficiency_improvs/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90
- Timestamp:
- 2015-12-14T10:06:06+01:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/dev_r5776_UKMO2_OBS_efficiency_improvs/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90
r5656 r6041 14 14 !! use of parameters in par_CONFIG-Rxx.h90, not in namelist 15 15 !! - ! 2004-05 (A. Koch-Larrouy) Add Gyre configuration 16 !! 4.0 ! 2011-02 (G. Madec) add cell surface (e1e2t) 16 !! 3.7 ! 2015-09 (G. Madec, S. Flavoni) add cell surface and their inverse 17 !! add optional read of e1e2u & e1e2v 17 18 !!---------------------------------------------------------------------- 18 19 … … 23 24 USE dom_oce ! ocean space and time domain 24 25 USE phycst ! physical constants 26 USE domwri ! write 'meshmask.nc' & 'coordinate_e1e2u_v.nc' files 27 ! 25 28 USE in_out_manager ! I/O manager 26 29 USE lib_mpp ! MPP library … … 35 38 36 39 !!---------------------------------------------------------------------- 37 !! NEMO/OPA 4.0 , NEMO Consortium (2011)40 !! NEMO/OPA 3.7 , NEMO Consortium (2014) 38 41 !! $Id$ 39 42 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 106 109 REAL(wp) :: zphi1, zsin_alpha, zim05, zjm05 107 110 INTEGER :: isrow ! index for ORCA1 starting row 108 111 INTEGER :: ie1e2u_v ! fag for u- & v-surface read in coordinate file or not 109 112 !!---------------------------------------------------------------------- 110 113 ! … … 122 125 WRITE(numout,*) ' meridional grid-spacing (meters) ppe2_m = ', ppe2_m 123 126 ENDIF 124 125 126 SELECT CASE( jphgr_msh ) ! type of horizontal mesh127 128 CASE ( 0 ) ! curvilinear coordinate on the sphere read in coordinate.nc file129 127 ! 128 ! 129 SELECT CASE( jphgr_msh ) ! type of horizontal mesh 130 ! 131 CASE ( 0 ) !== read in coordinate.nc file ==! 132 ! 130 133 IF(lwp) WRITE(numout,*) 131 134 IF(lwp) WRITE(numout,*) ' curvilinear coordinate on the sphere read in "coordinate" file' 132 133 CALL hgr_read ! Defaultl option : NetCDF file 134 135 ! ! ===================== 136 IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! ORCA R2 configuration 137 ! ! ===================== 138 IF( nn_cla == 0 ) THEN 139 ! 140 ii0 = 139 ; ii1 = 140 ! Gibraltar Strait (e2u = 20 km) 141 ij0 = 102 ; ij1 = 102 ; e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 20.e3 142 IF(lwp) WRITE(numout,*) 143 IF(lwp) WRITE(numout,*) ' orca_r2: Gibraltar : e2u reduced to 20 km' 144 ! 145 ii0 = 160 ; ii1 = 160 ! Bab el Mandeb (e2u = 18 km) 146 ij0 = 88 ; ij1 = 88 ; e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 18.e3 147 e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 30.e3 148 IF(lwp) WRITE(numout,*) 149 IF(lwp) WRITE(numout,*) ' orca_r2: Bab el Mandeb: e2u reduced to 30 km' 150 IF(lwp) WRITE(numout,*) ' e1v reduced to 18 km' 151 ENDIF 152 153 ii0 = 145 ; ii1 = 146 ! Danish Straits (e2u = 10 km) 154 ij0 = 116 ; ij1 = 116 ; e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 10.e3 155 IF(lwp) WRITE(numout,*) 156 IF(lwp) WRITE(numout,*) ' orca_r2: Danish Straits : e2u reduced to 10 km' 157 ! 135 ! 136 ie1e2u_v = 0 ! set to unread e1e2u and e1e2v 137 ! 138 CALL hgr_read( ie1e2u_v ) ! read the coordinate.nc file 139 ! 140 IF( ie1e2u_v == 0 ) THEN ! e1e2u and e1e2v have not been read: compute them 141 ! ! e2u and e1v does not include a reduction in some strait: apply reduction 142 e1e2u (:,:) = e1u(:,:) * e2u(:,:) 143 e1e2v (:,:) = e1v(:,:) * e2v(:,:) 158 144 ENDIF 159 160 ! ! ===================== 161 IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN ! ORCA R1 configuration 162 ! ! ===================== 163 ! This dirty section will be suppressed by simplification process: all this will come back in input files 164 ! Currently these hard-wired indices relate to configuration with 165 ! extend grid (jpjglo=332) 166 ! which had a grid-size of 362x292. 167 ! 168 isrow = 332 - jpjglo 169 ! 170 ii0 = 282 ; ii1 = 283 ! Gibraltar Strait (e2u = 20 km) 171 ij0 = 241 - isrow ; ij1 = 241 - isrow ; e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 20.e3 172 IF(lwp) WRITE(numout,*) 173 IF(lwp) WRITE(numout,*) ' orca_r1: Gibraltar : e2u reduced to 20 km' 174 175 ii0 = 314 ; ii1 = 315 ! Bhosporus Strait (e2u = 10 km) 176 ij0 = 248 - isrow ; ij1 = 248 - isrow ; e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 10.e3 177 IF(lwp) WRITE(numout,*) 178 IF(lwp) WRITE(numout,*) ' orca_r1: Bhosporus : e2u reduced to 10 km' 179 180 ii0 = 44 ; ii1 = 44 ! Lombok Strait (e1v = 13 km) 181 ij0 = 164 - isrow ; ij1 = 165 - isrow ; e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 13.e3 182 IF(lwp) WRITE(numout,*) 183 IF(lwp) WRITE(numout,*) ' orca_r1: Lombok : e1v reduced to 10 km' 184 185 ii0 = 48 ; ii1 = 48 ! Sumba Strait (e1v = 8 km) [closed from bathy_11 on] 186 ij0 = 164 - isrow ; ij1 = 165 - isrow ; e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 8.e3 187 IF(lwp) WRITE(numout,*) 188 IF(lwp) WRITE(numout,*) ' orca_r1: Sumba : e1v reduced to 8 km' 189 190 ii0 = 53 ; ii1 = 53 ! Ombai Strait (e1v = 13 km) 191 ij0 = 164 - isrow ; ij1 = 165 - isrow ; e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 13.e3 192 IF(lwp) WRITE(numout,*) 193 IF(lwp) WRITE(numout,*) ' orca_r1: Ombai : e1v reduced to 13 km' 194 195 ii0 = 56 ; ii1 = 56 ! Timor Passage (e1v = 20 km) 196 ij0 = 164 - isrow ; ij1 = 145 - isrow ; e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 20.e3 197 IF(lwp) WRITE(numout,*) 198 IF(lwp) WRITE(numout,*) ' orca_r1: Timor Passage : e1v reduced to 20 km' 199 200 ii0 = 55 ; ii1 = 55 ! West Halmahera Strait (e1v = 30 km) 201 ij0 = 181 - isrow ; ij1 = 182 - isrow ; e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 30.e3 202 IF(lwp) WRITE(numout,*) 203 IF(lwp) WRITE(numout,*) ' orca_r1: W Halmahera : e1v reduced to 30 km' 204 205 ii0 = 58 ; ii1 = 58 ! East Halmahera Strait (e1v = 50 km) 206 ij0 = 181 - isrow ; ij1 = 182 - isrow ; e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 50.e3 207 IF(lwp) WRITE(numout,*) 208 IF(lwp) WRITE(numout,*) ' orca_r1: E Halmahera : e1v reduced to 50 km' 209 ! 210 ! 211 ENDIF 212 213 ! ! ====================== 214 IF( cp_cfg == "orca" .AND. jp_cfg == 05 ) THEN ! ORCA R05 configuration 215 ! ! ====================== 216 ii0 = 563 ; ii1 = 564 ! Gibraltar Strait (e2u = 20 km) 217 ij0 = 327 ; ij1 = 327 ; e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 20.e3 218 IF(lwp) WRITE(numout,*) 219 IF(lwp) WRITE(numout,*) ' orca_r05: Reduced e2u at the Gibraltar Strait' 220 ! 221 ii0 = 627 ; ii1 = 628 ! Bosphore Strait (e2u = 10 km) 222 ij0 = 343 ; ij1 = 343 ; e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 10.e3 223 IF(lwp) WRITE(numout,*) 224 IF(lwp) WRITE(numout,*) ' orca_r05: Reduced e2u at the Bosphore Strait' 225 ! 226 ii0 = 93 ; ii1 = 94 ! Sumba Strait (e2u = 40 km) 227 ij0 = 232 ; ij1 = 232 ; e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 40.e3 228 IF(lwp) WRITE(numout,*) 229 IF(lwp) WRITE(numout,*) ' orca_r05: Reduced e2u at the Sumba Strait' 230 ! 231 ii0 = 103 ; ii1 = 103 ! Ombai Strait (e2u = 15 km) 232 ij0 = 232 ; ij1 = 232 ; e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 15.e3 233 IF(lwp) WRITE(numout,*) 234 IF(lwp) WRITE(numout,*) ' orca_r05: Reduced e2u at the Ombai Strait' 235 ! 236 ii0 = 15 ; ii1 = 15 ! Palk Strait (e2u = 10 km) 237 ij0 = 270 ; ij1 = 270 ; e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 10.e3 238 IF(lwp) WRITE(numout,*) 239 IF(lwp) WRITE(numout,*) ' orca_r05: Reduced e2u at the Palk Strait' 240 ! 241 ii0 = 87 ; ii1 = 87 ! Lombok Strait (e1v = 10 km) 242 ij0 = 232 ; ij1 = 233 ; e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 10.e3 243 IF(lwp) WRITE(numout,*) 244 IF(lwp) WRITE(numout,*) ' orca_r05: Reduced e1v at the Lombok Strait' 245 ! 246 ! 247 ii0 = 662 ; ii1 = 662 ! Bab el Mandeb (e1v = 25 km) 248 ij0 = 276 ; ij1 = 276 ; e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 25.e3 249 IF(lwp) WRITE(numout,*) 250 IF(lwp) WRITE(numout,*) ' orca_r05: Reduced e1v at the Bab el Mandeb' 251 ! 252 ENDIF 253 254 255 ! N.B. : General case, lat and long function of both i and j indices: 256 ! e1t(ji,jj) = ra * rad * SQRT( ( cos( rad*gphit(ji,jj) ) * fsdila( zti, ztj ) )**2 & 257 ! + ( fsdiph( zti, ztj ) )**2 ) 258 ! e1u(ji,jj) = ra * rad * SQRT( ( cos( rad*gphiu(ji,jj) ) * fsdila( zui, zuj ) )**2 & 259 ! + ( fsdiph( zui, zuj ) )**2 ) 260 ! e1v(ji,jj) = ra * rad * SQRT( ( cos( rad*gphiv(ji,jj) ) * fsdila( zvi, zvj ) )**2 & 261 ! + ( fsdiph( zvi, zvj ) )**2 ) 262 ! e1f(ji,jj) = ra * rad * SQRT( ( cos( rad*gphif(ji,jj) ) * fsdila( zfi, zfj ) )**2 & 263 ! + ( fsdiph( zfi, zfj ) )**2 ) 264 ! 265 ! e2t(ji,jj) = ra * rad * SQRT( ( cos( rad*gphit(ji,jj) ) * fsdjla( zti, ztj ) )**2 & 266 ! + ( fsdjph( zti, ztj ) )**2 ) 267 ! e2u(ji,jj) = ra * rad * SQRT( ( cos( rad*gphiu(ji,jj) ) * fsdjla( zui, zuj ) )**2 & 268 ! + ( fsdjph( zui, zuj ) )**2 ) 269 ! e2v(ji,jj) = ra * rad * SQRT( ( cos( rad*gphiv(ji,jj) ) * fsdjla( zvi, zvj ) )**2 & 270 ! + ( fsdjph( zvi, zvj ) )**2 ) 271 ! e2f(ji,jj) = ra * rad * SQRT( ( cos( rad*gphif(ji,jj) ) * fsdjla( zfi, zfj ) )**2 & 272 ! + ( fsdjph( zfi, zfj ) )**2 ) 273 274 275 CASE ( 1 ) ! geographical mesh on the sphere with regular grid-spacing 276 145 ! 146 CASE ( 1 ) !== geographical mesh on the sphere with regular (in degree) grid-spacing ==! 147 ! 277 148 IF(lwp) WRITE(numout,*) 278 149 IF(lwp) WRITE(numout,*) ' geographical mesh on the sphere with regular grid-spacing' 279 150 IF(lwp) WRITE(numout,*) ' given by ppe1_deg and ppe2_deg' 280 151 ! 281 152 DO jj = 1, jpj 282 153 DO ji = 1, jpi 283 zti = FLOAT( ji - 1 + nimpp - 1 ) ; ztj = FLOAT( jj - 1 + njmpp - 1 )284 zui = FLOAT( ji - 1 + nimpp - 1 ) + 0.5 ; zuj = FLOAT( jj - 1 + njmpp - 1 )285 zvi = FLOAT( ji - 1 + nimpp - 1 ) ; zvj = FLOAT( jj - 1 + njmpp - 1 ) + 0.5286 zfi = FLOAT( ji - 1 + nimpp - 1 ) + 0.5 ; zfj = FLOAT( jj - 1 + njmpp - 1 ) + 0.5154 zti = REAL( ji - 1 + nimpp - 1 ) ; ztj = REAL( jj - 1 + njmpp - 1 ) 155 zui = REAL( ji - 1 + nimpp - 1 ) + 0.5 ; zuj = REAL( jj - 1 + njmpp - 1 ) 156 zvi = REAL( ji - 1 + nimpp - 1 ) ; zvj = REAL( jj - 1 + njmpp - 1 ) + 0.5 157 zfi = REAL( ji - 1 + nimpp - 1 ) + 0.5 ; zfj = REAL( jj - 1 + njmpp - 1 ) + 0.5 287 158 ! Longitude 288 159 glamt(ji,jj) = ppglam0 + ppe1_deg * zti … … 307 178 END DO 308 179 END DO 309 310 311 CASE ( 2:3 ) ! f- or beta-plane with regular grid-spacing 312 180 ! 181 CASE ( 2:3 ) !== f- or beta-plane with regular grid-spacing ==! 182 ! 313 183 IF(lwp) WRITE(numout,*) 314 184 IF(lwp) WRITE(numout,*) ' f- or beta-plane with regular grid-spacing' 315 185 IF(lwp) WRITE(numout,*) ' given by ppe1_m and ppe2_m' 316 186 ! 317 187 ! Position coordinates (in kilometers) 318 188 ! ========== 319 glam0 = 0. e0189 glam0 = 0._wp 320 190 gphi0 = - ppe2_m * 1.e-3 321 191 ! 322 192 #if defined key_agrif 323 193 IF ( cp_cfg == 'eel' .AND. jp_cfg == 6 ) THEN ! for EEL6 configuration only … … 332 202 DO jj = 1, jpj 333 203 DO ji = 1, jpi 334 glamt(ji,jj) = glam0 + ppe1_m * 1.e-3 * ( FLOAT( ji - 1 + nimpp - 1 ) )335 glamu(ji,jj) = glam0 + ppe1_m * 1.e-3 * ( FLOAT( ji - 1 + nimpp - 1 ) + 0.5 )204 glamt(ji,jj) = glam0 + ppe1_m * 1.e-3 * ( REAL( ji - 1 + nimpp - 1 ) ) 205 glamu(ji,jj) = glam0 + ppe1_m * 1.e-3 * ( REAL( ji - 1 + nimpp - 1 ) + 0.5 ) 336 206 glamv(ji,jj) = glamt(ji,jj) 337 207 glamf(ji,jj) = glamu(ji,jj) 338 339 gphit(ji,jj) = gphi0 + ppe2_m * 1.e-3 * ( FLOAT( jj - 1 + njmpp - 1 ) )208 ! 209 gphit(ji,jj) = gphi0 + ppe2_m * 1.e-3 * ( REAL( jj - 1 + njmpp - 1 ) ) 340 210 gphiu(ji,jj) = gphit(ji,jj) 341 gphiv(ji,jj) = gphi0 + ppe2_m * 1.e-3 * ( FLOAT( jj - 1 + njmpp - 1 ) + 0.5 )211 gphiv(ji,jj) = gphi0 + ppe2_m * 1.e-3 * ( REAL( jj - 1 + njmpp - 1 ) + 0.5 ) 342 212 gphif(ji,jj) = gphiv(ji,jj) 343 213 END DO 344 214 END DO 345 215 ! 346 216 ! Horizontal scale factors (in meters) 347 217 ! ====== … … 350 220 e1v(:,:) = ppe1_m ; e2v(:,:) = ppe2_m 351 221 e1f(:,:) = ppe1_m ; e2f(:,:) = ppe2_m 352 353 CASE ( 4 ) ! geographical mesh on the sphere, isotropic MERCATOR type354 222 ! 223 CASE ( 4 ) !== geographical mesh on the sphere, isotropic MERCATOR type ==! 224 ! 355 225 IF(lwp) WRITE(numout,*) 356 226 IF(lwp) WRITE(numout,*) ' geographical mesh on the sphere, MERCATOR type' 357 227 IF(lwp) WRITE(numout,*) ' longitudinal/latitudinal spacing given by ppe1_deg' 358 228 IF ( ppgphi0 == -90 ) CALL ctl_stop( ' Mercator grid cannot start at south pole !!!! ' ) 359 229 ! 360 230 ! Find index corresponding to the equator, given the grid spacing e1_deg 361 231 ! and the (approximate) southern latitude ppgphi0. … … 365 235 ijeq = ABS( 180./rpi * LOG( COS( zarg ) / SIN( zarg ) ) / ppe1_deg ) 366 236 IF( ppgphi0 > 0 ) ijeq = -ijeq 367 237 ! 368 238 IF(lwp) WRITE(numout,*) ' Index of the equator on the MERCATOR grid:', ijeq 369 239 ! 370 240 DO jj = 1, jpj 371 241 DO ji = 1, jpi 372 zti = FLOAT( ji - 1 + nimpp - 1 ) ; ztj = FLOAT( jj - ijeq + njmpp - 1 )373 zui = FLOAT( ji - 1 + nimpp - 1 ) + 0.5 ; zuj = FLOAT( jj - ijeq + njmpp - 1 )374 zvi = FLOAT( ji - 1 + nimpp - 1 ) ; zvj = FLOAT( jj - ijeq + njmpp - 1 ) + 0.5375 zfi = FLOAT( ji - 1 + nimpp - 1 ) + 0.5 ; zfj = FLOAT( jj - ijeq + njmpp - 1 ) + 0.5242 zti = REAL( ji - 1 + nimpp - 1 ) ; ztj = REAL( jj - ijeq + njmpp - 1 ) 243 zui = REAL( ji - 1 + nimpp - 1 ) + 0.5 ; zuj = REAL( jj - ijeq + njmpp - 1 ) 244 zvi = REAL( ji - 1 + nimpp - 1 ) ; zvj = REAL( jj - ijeq + njmpp - 1 ) + 0.5 245 zfi = REAL( ji - 1 + nimpp - 1 ) + 0.5 ; zfj = REAL( jj - ijeq + njmpp - 1 ) + 0.5 376 246 ! Longitude 377 247 glamt(ji,jj) = ppglam0 + ppe1_deg * zti … … 396 266 END DO 397 267 END DO 398 399 CASE ( 5 ) ! beta-plane with regular grid-spacing and rotated domain(GYRE configuration)400 268 ! 269 CASE ( 5 ) !== beta-plane with regular grid-spacing and rotated domain ==! (GYRE configuration) 270 ! 401 271 IF(lwp) WRITE(numout,*) 402 272 IF(lwp) WRITE(numout,*) ' beta-plane with regular grid-spacing and rotated domain (GYRE configuration)' 403 273 IF(lwp) WRITE(numout,*) ' given by ppe1_m and ppe2_m' 404 274 ! 405 275 ! Position coordinates (in kilometers) 406 276 ! ========== 407 277 ! 408 278 ! angle 45deg and ze1=106.e+3 / jp_cfg forced -> zlam1 = -85deg, zphi1 = 29degN 409 zlam1 = -85 410 zphi1 = 29279 zlam1 = -85._wp 280 zphi1 = 29._wp 411 281 ! resolution in meters 412 ze1 = 106000. / FLOAT(jp_cfg)282 ze1 = 106000. / REAL( jp_cfg , wp ) 413 283 ! benchmark: forced the resolution to be about 100 km 414 IF( nbench /= 0 ) ze1 = 106000. e0415 zsin_alpha = - SQRT( 2. ) / 2.416 zcos_alpha = SQRT( 2. ) / 2.284 IF( nbench /= 0 ) ze1 = 106000._wp 285 zsin_alpha = - SQRT( 2._wp ) * 0.5_wp 286 zcos_alpha = SQRT( 2._wp ) * 0.5_wp 417 287 ze1deg = ze1 / (ra * rad) 418 IF( nbench /= 0 ) ze1deg = ze1deg / FLOAT(jp_cfg)! benchmark: keep the lat/+lon419 ! ! at the right jp_cfg resolution420 glam0 = zlam1 + zcos_alpha * ze1deg * FLOAT( jpjglo-2)421 gphi0 = zphi1 + zsin_alpha * ze1deg * FLOAT( jpjglo-2)422 288 IF( nbench /= 0 ) ze1deg = ze1deg / REAL( jp_cfg , wp ) ! benchmark: keep the lat/+lon 289 ! ! at the right jp_cfg resolution 290 glam0 = zlam1 + zcos_alpha * ze1deg * REAL( jpjglo-2 , wp ) 291 gphi0 = zphi1 + zsin_alpha * ze1deg * REAL( jpjglo-2 , wp ) 292 ! 423 293 IF( nprint==1 .AND. lwp ) THEN 424 294 WRITE(numout,*) ' ze1', ze1, 'cosalpha', zcos_alpha, 'sinalpha', zsin_alpha 425 295 WRITE(numout,*) ' ze1deg', ze1deg, 'glam0', glam0, 'gphi0', gphi0 426 296 ENDIF 427 297 ! 428 298 DO jj = 1, jpj 429 DO ji = 1, jpi430 zim1 = FLOAT( ji + nimpp - 1 ) - 1. ; zim05 = FLOAT( ji + nimpp - 1 ) - 1.5431 zjm1 = FLOAT( jj + njmpp - 1 ) - 1. ; zjm05 = FLOAT( jj + njmpp - 1 ) - 1.5432 433 glamf(ji,jj) = glam0 + zim1 * ze1deg * zcos_alpha + zjm1 * ze1deg * zsin_alpha434 gphif(ji,jj) = gphi0 - zim1 * ze1deg * zsin_alpha + zjm1 * ze1deg * zcos_alpha435 436 glamt(ji,jj) = glam0 + zim05 * ze1deg * zcos_alpha + zjm05 * ze1deg * zsin_alpha437 gphit(ji,jj) = gphi0 - zim05 * ze1deg * zsin_alpha + zjm05 * ze1deg * zcos_alpha438 439 glamu(ji,jj) = glam0 + zim1 * ze1deg * zcos_alpha + zjm05 * ze1deg * zsin_alpha440 gphiu(ji,jj) = gphi0 - zim1 * ze1deg * zsin_alpha + zjm05 * ze1deg * zcos_alpha441 442 glamv(ji,jj) = glam0 + zim05 * ze1deg * zcos_alpha + zjm1 * ze1deg * zsin_alpha443 gphiv(ji,jj) = gphi0 - zim05 * ze1deg * zsin_alpha + zjm1 * ze1deg * zcos_alpha444 END DO445 446 299 DO ji = 1, jpi 300 zim1 = REAL( ji + nimpp - 1 ) - 1. ; zim05 = REAL( ji + nimpp - 1 ) - 1.5 301 zjm1 = REAL( jj + njmpp - 1 ) - 1. ; zjm05 = REAL( jj + njmpp - 1 ) - 1.5 302 ! 303 glamf(ji,jj) = glam0 + zim1 * ze1deg * zcos_alpha + zjm1 * ze1deg * zsin_alpha 304 gphif(ji,jj) = gphi0 - zim1 * ze1deg * zsin_alpha + zjm1 * ze1deg * zcos_alpha 305 ! 306 glamt(ji,jj) = glam0 + zim05 * ze1deg * zcos_alpha + zjm05 * ze1deg * zsin_alpha 307 gphit(ji,jj) = gphi0 - zim05 * ze1deg * zsin_alpha + zjm05 * ze1deg * zcos_alpha 308 ! 309 glamu(ji,jj) = glam0 + zim1 * ze1deg * zcos_alpha + zjm05 * ze1deg * zsin_alpha 310 gphiu(ji,jj) = gphi0 - zim1 * ze1deg * zsin_alpha + zjm05 * ze1deg * zcos_alpha 311 ! 312 glamv(ji,jj) = glam0 + zim05 * ze1deg * zcos_alpha + zjm1 * ze1deg * zsin_alpha 313 gphiv(ji,jj) = gphi0 - zim05 * ze1deg * zsin_alpha + zjm1 * ze1deg * zcos_alpha 314 END DO 315 END DO 316 ! 447 317 ! Horizontal scale factors (in meters) 448 318 ! ====== … … 451 321 e1v(:,:) = ze1 ; e2v(:,:) = ze1 452 322 e1f(:,:) = ze1 ; e2f(:,:) = ze1 453 323 ! 454 324 CASE DEFAULT 455 325 WRITE(ctmp1,*) ' bad flag value for jphgr_msh = ', jphgr_msh 456 326 CALL ctl_stop( ctmp1 ) 457 327 ! 458 328 END SELECT 459 329 460 ! T-cell surface 461 ! -------------- 462 e1e2t(:,:) = e1t(:,:) * e2t(:,:) 463 464 ! Useful shortcuts (JC: note the duplicated e2e2t array ! Need some cleaning) 465 ! --------------------------------------------------------------------------- 466 e12t (:,:) = e1t(:,:) * e2t(:,:) 467 e12u (:,:) = e1u(:,:) * e2u(:,:) 468 e12v (:,:) = e1v(:,:) * e2v(:,:) 469 e12f (:,:) = e1f(:,:) * e2f(:,:) 470 r1_e12t (:,:) = 1._wp / e12t(:,:) 471 r1_e12u (:,:) = 1._wp / e12u(:,:) 472 r1_e12v (:,:) = 1._wp / e12v(:,:) 473 r1_e12f (:,:) = 1._wp / e12f(:,:) 474 re2u_e1u(:,:) = e2u(:,:) / e1u(:,:) 475 re1v_e2v(:,:) = e1v(:,:) / e2v(:,:) 476 r1_e1t (:,:) = 1._wp / e1t(:,:) 477 r1_e1u (:,:) = 1._wp / e1u(:,:) 478 r1_e1v (:,:) = 1._wp / e1v(:,:) 479 r1_e1f (:,:) = 1._wp / e1f(:,:) 480 r1_e2t (:,:) = 1._wp / e2t(:,:) 481 r1_e2u (:,:) = 1._wp / e2u(:,:) 482 r1_e2v (:,:) = 1._wp / e2v(:,:) 483 r1_e2f (:,:) = 1._wp / e2f(:,:) 484 485 ! Control printing : Grid informations (if not restart) 486 ! ---------------- 487 488 IF( lwp .AND. .NOT.ln_rstart ) THEN 330 ! associated horizontal metrics 331 ! ----------------------------- 332 ! 333 r1_e1t(:,:) = 1._wp / e1t(:,:) ; r1_e2t (:,:) = 1._wp / e2t(:,:) 334 r1_e1u(:,:) = 1._wp / e1u(:,:) ; r1_e2u (:,:) = 1._wp / e2u(:,:) 335 r1_e1v(:,:) = 1._wp / e1v(:,:) ; r1_e2v (:,:) = 1._wp / e2v(:,:) 336 r1_e1f(:,:) = 1._wp / e1f(:,:) ; r1_e2f (:,:) = 1._wp / e2f(:,:) 337 ! 338 e1e2t (:,:) = e1t(:,:) * e2t(:,:) ; r1_e1e2t(:,:) = 1._wp / e1e2t(:,:) 339 e1e2f (:,:) = e1f(:,:) * e2f(:,:) ; r1_e1e2f(:,:) = 1._wp / e1e2f(:,:) 340 IF( jphgr_msh /= 0 ) THEN ! e1e2u and e1e2v have not been set: compute them 341 e1e2u (:,:) = e1u(:,:) * e2u(:,:) 342 e1e2v (:,:) = e1v(:,:) * e2v(:,:) 343 ENDIF 344 r1_e1e2u(:,:) = 1._wp / e1e2u(:,:) ! compute their invert in both cases 345 r1_e1e2v(:,:) = 1._wp / e1e2v(:,:) 346 ! 347 e2_e1u(:,:) = e2u(:,:) / e1u(:,:) 348 e1_e2v(:,:) = e1v(:,:) / e2v(:,:) 349 350 IF( lwp .AND. .NOT.ln_rstart ) THEN ! Control print : Grid informations (if not restart) 489 351 WRITE(numout,*) 490 352 WRITE(numout,*) ' longitude and e1 scale factors' … … 496 358 9300 FORMAT( 1x, i4, f8.2,1x, f8.2,1x, f8.2,1x, f8.2, 1x, & 497 359 f19.10, 1x, f19.10, 1x, f19.10, 1x, f19.10 ) 498 360 ! 499 361 WRITE(numout,*) 500 362 WRITE(numout,*) ' latitude and e2 scale factors' … … 506 368 ENDIF 507 369 508 509 IF( nprint == 1 .AND. lwp ) THEN510 WRITE(numout,*) ' e1u e2u '511 CALL prihre( e1u,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )512 CALL prihre( e2u,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )513 WRITE(numout,*) ' e1v e2v '514 CALL prihre( e1v,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )515 CALL prihre( e2v,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )516 WRITE(numout,*) ' e1f e2f '517 CALL prihre( e1f,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )518 CALL prihre( e2f,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )519 ENDIF520 521 370 522 371 ! ================= ! … … 525 374 526 375 SELECT CASE( jphgr_msh ) ! type of horizontal mesh 527 376 ! 528 377 CASE ( 0, 1, 4 ) ! mesh on the sphere 529 378 ! 530 379 ff(:,:) = 2. * omega * SIN( rad * gphif(:,:) ) 531 380 ! 532 381 CASE ( 2 ) ! f-plane at ppgphi0 533 382 ! 534 383 ff(:,:) = 2. * omega * SIN( rad * ppgphi0 ) 535 384 ! 536 385 IF(lwp) WRITE(numout,*) ' f-plane: Coriolis parameter = constant = ', ff(1,1) 537 386 ! 538 387 CASE ( 3 ) ! beta-plane 539 388 ! 540 389 zbeta = 2. * omega * COS( rad * ppgphi0 ) / ra ! beta at latitude ppgphi0 541 zphi0 = ppgphi0 - FLOAT( jpjglo/2) * ppe2_m / ( ra * rad ) ! latitude of the first row F-points542 390 zphi0 = ppgphi0 - REAL( jpjglo/2) * ppe2_m / ( ra * rad ) ! latitude of the first row F-points 391 ! 543 392 #if defined key_agrif 544 393 IF ( cp_cfg == 'eel' .AND. jp_cfg == 6 ) THEN ! for EEL6 configuration only 545 394 IF( .NOT. Agrif_Root() ) THEN 546 zphi0 = ppgphi0 - FLOAT( Agrif_Parent(jpjglo)/2)*Agrif_Parent(ppe2_m) & 547 & / (ra * rad) 395 zphi0 = ppgphi0 - REAL( Agrif_Parent(jpjglo)/2)*Agrif_Parent(ppe2_m) / (ra * rad) 548 396 ENDIF 549 397 ENDIF 550 398 #endif 551 399 zf0 = 2. * omega * SIN( rad * zphi0 ) ! compute f0 1st point south 552 400 ! 553 401 ff(:,:) = ( zf0 + zbeta * gphif(:,:) * 1.e+3 ) ! f = f0 +beta* y ( y=0 at south) 554 402 ! 555 403 IF(lwp) THEN 556 404 WRITE(numout,*) … … 565 413 IF(lwp) WRITE(numout,*) ' Coriolis parameter varies globally from ', zminff,' to ', zmaxff 566 414 END IF 567 415 ! 568 416 CASE ( 5 ) ! beta-plane and rotated domain (gyre configuration) 569 417 ! 570 418 zbeta = 2. * omega * COS( rad * ppgphi0 ) / ra ! beta at latitude ppgphi0 571 zphi0 = 15. e0! latitude of the first row F-points419 zphi0 = 15._wp ! latitude of the first row F-points 572 420 zf0 = 2. * omega * SIN( rad * zphi0 ) ! compute f0 1st point south 573 421 ! 574 422 ff(:,:) = ( zf0 + zbeta * ABS( gphif(:,:) - zphi0 ) * rad * ra ) ! f = f0 +beta* y ( y=0 at south) 575 423 ! 576 424 IF(lwp) THEN 577 425 WRITE(numout,*) … … 579 427 WRITE(numout,*) ' Coriolis parameter varies in this processor from ', ff(nldi,nldj),' to ', ff(nldi,nlej) 580 428 ENDIF 581 429 ! 582 430 IF( lk_mpp ) THEN 583 431 zminff=ff(nldi,nldj) … … 587 435 IF(lwp) WRITE(numout,*) ' Coriolis parameter varies globally from ', zminff,' to ', zmaxff 588 436 END IF 589 437 ! 590 438 END SELECT 591 439 … … 596 444 597 445 IF( nperio == 2 ) THEN 598 znorme = SQRT( SUM( gphiu(:,2) * gphiu(:,2) ) ) / FLOAT( jpi )446 znorme = SQRT( SUM( gphiu(:,2) * gphiu(:,2) ) ) / REAL( jpi ) 599 447 IF( znorme > 1.e-13 ) CALL ctl_stop( ' ===>>>> : symmetrical condition: rerun with good equator line' ) 600 448 ENDIF … … 605 453 606 454 607 SUBROUTINE hgr_read 455 SUBROUTINE hgr_read( ke1e2u_v ) 608 456 !!--------------------------------------------------------------------- 609 457 !! *** ROUTINE hgr_read *** 610 458 !! 611 !! ** Purpose : Read a coordinate file in NetCDF format 612 !! 613 !! ** Method : The mesh file has been defined trough a analytical 614 !! or semi-analytical method. It is read in a NetCDF file. 615 !! 459 !! ** Purpose : Read a coordinate file in NetCDF format using IOM 460 !! 616 461 !!---------------------------------------------------------------------- 617 462 USE iom 618 463 !! 464 INTEGER, INTENT( inout ) :: ke1e2u_v ! fag: e1e2u & e1e2v read in coordinate file (=1) or not (=0) 465 ! 619 466 INTEGER :: inum ! temporary logical unit 620 467 !!---------------------------------------------------------------------- 621 468 ! 622 469 IF(lwp) THEN 623 470 WRITE(numout,*) … … 625 472 WRITE(numout,*) '~~~~~~~~ jpiglo = ', jpiglo, ' jpjglo = ', jpjglo, ' jpk = ', jpk 626 473 ENDIF 627 474 ! 628 475 CALL iom_open( 'coordinates', inum ) 629 476 ! 630 477 CALL iom_get( inum, jpdom_data, 'glamt', glamt, lrowattr=ln_use_jattr ) 631 478 CALL iom_get( inum, jpdom_data, 'glamu', glamu, lrowattr=ln_use_jattr ) 632 479 CALL iom_get( inum, jpdom_data, 'glamv', glamv, lrowattr=ln_use_jattr ) 633 480 CALL iom_get( inum, jpdom_data, 'glamf', glamf, lrowattr=ln_use_jattr ) 634 481 ! 635 482 CALL iom_get( inum, jpdom_data, 'gphit', gphit, lrowattr=ln_use_jattr ) 636 483 CALL iom_get( inum, jpdom_data, 'gphiu', gphiu, lrowattr=ln_use_jattr ) 637 484 CALL iom_get( inum, jpdom_data, 'gphiv', gphiv, lrowattr=ln_use_jattr ) 638 485 CALL iom_get( inum, jpdom_data, 'gphif', gphif, lrowattr=ln_use_jattr ) 639 640 CALL iom_get( inum, jpdom_data, 'e1t', e1t, lrowattr=ln_use_jattr ) 641 CALL iom_get( inum, jpdom_data, 'e1u', e1u, lrowattr=ln_use_jattr ) 642 CALL iom_get( inum, jpdom_data, 'e1v', e1v, lrowattr=ln_use_jattr ) 643 CALL iom_get( inum, jpdom_data, 'e1f', e1f, lrowattr=ln_use_jattr ) 644 645 CALL iom_get( inum, jpdom_data, 'e2t', e2t, lrowattr=ln_use_jattr ) 646 CALL iom_get( inum, jpdom_data, 'e2u', e2u, lrowattr=ln_use_jattr ) 647 CALL iom_get( inum, jpdom_data, 'e2v', e2v, lrowattr=ln_use_jattr ) 648 CALL iom_get( inum, jpdom_data, 'e2f', e2f, lrowattr=ln_use_jattr ) 649 486 ! 487 CALL iom_get( inum, jpdom_data, 'e1t' , e1t , lrowattr=ln_use_jattr ) 488 CALL iom_get( inum, jpdom_data, 'e1u' , e1u , lrowattr=ln_use_jattr ) 489 CALL iom_get( inum, jpdom_data, 'e1v' , e1v , lrowattr=ln_use_jattr ) 490 CALL iom_get( inum, jpdom_data, 'e1f' , e1f , lrowattr=ln_use_jattr ) 491 ! 492 CALL iom_get( inum, jpdom_data, 'e2t' , e2t , lrowattr=ln_use_jattr ) 493 CALL iom_get( inum, jpdom_data, 'e2u' , e2u , lrowattr=ln_use_jattr ) 494 CALL iom_get( inum, jpdom_data, 'e2v' , e2v , lrowattr=ln_use_jattr ) 495 CALL iom_get( inum, jpdom_data, 'e2f' , e2f , lrowattr=ln_use_jattr ) 496 ! 497 IF( iom_varid( inum, 'e1e2u', ldstop = .FALSE. ) > 0 ) THEN 498 IF(lwp) WRITE(numout,*) 'hgr_read : e1e2u & e1e2v read in coordinates file' 499 CALL iom_get( inum, jpdom_data, 'e1e2u' , e1e2u , lrowattr=ln_use_jattr ) 500 CALL iom_get( inum, jpdom_data, 'e1e2v' , e1e2v , lrowattr=ln_use_jattr ) 501 ke1e2u_v = 1 502 ELSE 503 ke1e2u_v = 0 504 ENDIF 505 ! 650 506 CALL iom_close( inum ) 651 507 508 !!gm THIS is TO BE REMOVED !!!!!!! 509 652 510 ! need to be define for the extended grid south of -80S 653 511 ! some point are undefined but you need to have e1 and e2 .NE. 0 … … 676 534 e2f=1.0e2 677 535 END WHERE 536 !!gm end 678 537 679 538 END SUBROUTINE hgr_read
Note: See TracChangeset
for help on using the changeset viewer.