- Timestamp:
- 2016-01-08T10:35:19+01:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4704_NOC5_MPP_BDY_UPDATE/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90
r4366 r6225 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) … … 105 108 REAL(wp) :: zlam1, zcos_alpha, zim1 , zjm1 , ze1, ze1deg 106 109 REAL(wp) :: zphi1, zsin_alpha, zim05, zjm05 110 INTEGER :: isrow ! index for ORCA1 starting row 111 INTEGER :: ie1e2u_v ! fag for u- & v-surface read in coordinate file or not 107 112 !!---------------------------------------------------------------------- 108 113 ! … … 120 125 WRITE(numout,*) ' meridional grid-spacing (meters) ppe2_m = ', ppe2_m 121 126 ENDIF 122 123 124 SELECT CASE( jphgr_msh ) ! type of horizontal mesh125 126 CASE ( 0 ) ! curvilinear coordinate on the sphere read in coordinate.nc file127 127 ! 128 ! 129 SELECT CASE( jphgr_msh ) ! type of horizontal mesh 130 ! 131 CASE ( 0 ) !== read in coordinate.nc file ==! 132 ! 128 133 IF(lwp) WRITE(numout,*) 129 134 IF(lwp) WRITE(numout,*) ' curvilinear coordinate on the sphere read in "coordinate" file' 130 131 CALL hgr_read ! Defaultl option : NetCDF file 132 133 ! ! ===================== 134 IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! ORCA R2 configuration 135 ! ! ===================== 136 IF( nn_cla == 0 ) THEN 137 ! 138 ii0 = 139 ; ii1 = 140 ! Gibraltar Strait (e2u = 20 km) 139 ij0 = 102 ; ij1 = 102 ; e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 20.e3 140 IF(lwp) WRITE(numout,*) 141 IF(lwp) WRITE(numout,*) ' orca_r2: Gibraltar : e2u reduced to 20 km' 142 ! 143 ii0 = 160 ; ii1 = 160 ! Bab el Mandeb (e2u = 18 km) 144 ij0 = 88 ; ij1 = 88 ; e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 18.e3 145 e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 30.e3 146 IF(lwp) WRITE(numout,*) 147 IF(lwp) WRITE(numout,*) ' orca_r2: Bab el Mandeb: e2u reduced to 30 km' 148 IF(lwp) WRITE(numout,*) ' e1v reduced to 18 km' 149 ENDIF 150 151 ii0 = 145 ; ii1 = 146 ! Danish Straits (e2u = 10 km) 152 ij0 = 116 ; ij1 = 116 ; e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 10.e3 153 IF(lwp) WRITE(numout,*) 154 IF(lwp) WRITE(numout,*) ' orca_r2: Danish Straits : e2u reduced to 10 km' 155 ! 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(:,:) 156 144 ENDIF 157 158 ! ! ===================== 159 IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN ! ORCA R1 configuration 160 ! ! ===================== 161 162 ii0 = 281 ; ii1 = 282 ! Gibraltar Strait (e2u = 20 km) 163 ij0 = 200 ; ij1 = 200 ; e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 20.e3 164 IF(lwp) WRITE(numout,*) 165 IF(lwp) WRITE(numout,*) ' orca_r1: Gibraltar : e2u reduced to 20 km' 166 167 ii0 = 314 ; ii1 = 315 ! Bhosporus Strait (e2u = 10 km) 168 ij0 = 208 ; ij1 = 208 ; e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 10.e3 169 IF(lwp) WRITE(numout,*) 170 IF(lwp) WRITE(numout,*) ' orca_r1: Bhosporus : e2u reduced to 10 km' 171 172 ii0 = 44 ; ii1 = 44 ! Lombok Strait (e1v = 13 km) 173 ij0 = 124 ; ij1 = 125 ; e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 13.e3 174 IF(lwp) WRITE(numout,*) 175 IF(lwp) WRITE(numout,*) ' orca_r1: Lombok : e1v reduced to 10 km' 176 177 ii0 = 48 ; ii1 = 48 ! Sumba Strait (e1v = 8 km) [closed from bathy_11 on] 178 ij0 = 124 ; ij1 = 125 ; e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 8.e3 179 IF(lwp) WRITE(numout,*) 180 IF(lwp) WRITE(numout,*) ' orca_r1: Sumba : e1v reduced to 8 km' 181 182 ii0 = 53 ; ii1 = 53 ! Ombai Strait (e1v = 13 km) 183 ij0 = 124 ; ij1 = 125 ; e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 13.e3 184 IF(lwp) WRITE(numout,*) 185 IF(lwp) WRITE(numout,*) ' orca_r1: Ombai : e1v reduced to 13 km' 186 187 ii0 = 56 ; ii1 = 56 ! Timor Passage (e1v = 20 km) 188 ij0 = 124 ; ij1 = 125 ; e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 20.e3 189 IF(lwp) WRITE(numout,*) 190 IF(lwp) WRITE(numout,*) ' orca_r1: Timor Passage : e1v reduced to 20 km' 191 192 ii0 = 55 ; ii1 = 55 ! West Halmahera Strait (e1v = 30 km) 193 ij0 = 141 ; ij1 = 142 ; e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 30.e3 194 IF(lwp) WRITE(numout,*) 195 IF(lwp) WRITE(numout,*) ' orca_r1: W Halmahera : e1v reduced to 30 km' 196 197 ii0 = 58 ; ii1 = 58 ! East Halmahera Strait (e1v = 50 km) 198 ij0 = 141 ; ij1 = 142 ; e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 50.e3 199 IF(lwp) WRITE(numout,*) 200 IF(lwp) WRITE(numout,*) ' orca_r1: E Halmahera : e1v reduced to 50 km' 201 202 ! 203 204 ! 205 ! 206 ! 207 ! 208 ENDIF 209 210 ! ! ====================== 211 IF( cp_cfg == "orca" .AND. jp_cfg == 05 ) THEN ! ORCA R05 configuration 212 ! ! ====================== 213 ii0 = 563 ; ii1 = 564 ! Gibraltar Strait (e2u = 20 km) 214 ij0 = 327 ; ij1 = 327 ; e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 20.e3 215 IF(lwp) WRITE(numout,*) 216 IF(lwp) WRITE(numout,*) ' orca_r05: Reduced e2u at the Gibraltar Strait' 217 ! 218 ii0 = 627 ; ii1 = 628 ! Bosphore Strait (e2u = 10 km) 219 ij0 = 343 ; ij1 = 343 ; e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 10.e3 220 IF(lwp) WRITE(numout,*) 221 IF(lwp) WRITE(numout,*) ' orca_r05: Reduced e2u at the Bosphore Strait' 222 ! 223 ii0 = 93 ; ii1 = 94 ! Sumba Strait (e2u = 40 km) 224 ij0 = 232 ; ij1 = 232 ; e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 40.e3 225 IF(lwp) WRITE(numout,*) 226 IF(lwp) WRITE(numout,*) ' orca_r05: Reduced e2u at the Sumba Strait' 227 ! 228 ii0 = 103 ; ii1 = 103 ! Ombai Strait (e2u = 15 km) 229 ij0 = 232 ; ij1 = 232 ; e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 15.e3 230 IF(lwp) WRITE(numout,*) 231 IF(lwp) WRITE(numout,*) ' orca_r05: Reduced e2u at the Ombai Strait' 232 ! 233 ii0 = 15 ; ii1 = 15 ! Palk Strait (e2u = 10 km) 234 ij0 = 270 ; ij1 = 270 ; e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 10.e3 235 IF(lwp) WRITE(numout,*) 236 IF(lwp) WRITE(numout,*) ' orca_r05: Reduced e2u at the Palk Strait' 237 ! 238 ii0 = 87 ; ii1 = 87 ! Lombok Strait (e1v = 10 km) 239 ij0 = 232 ; ij1 = 233 ; e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 10.e3 240 IF(lwp) WRITE(numout,*) 241 IF(lwp) WRITE(numout,*) ' orca_r05: Reduced e1v at the Lombok Strait' 242 ! 243 ! 244 ii0 = 662 ; ii1 = 662 ! Bab el Mandeb (e1v = 25 km) 245 ij0 = 276 ; ij1 = 276 ; e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 25.e3 246 IF(lwp) WRITE(numout,*) 247 IF(lwp) WRITE(numout,*) ' orca_r05: Reduced e1v at the Bab el Mandeb' 248 ! 249 ENDIF 250 251 252 ! N.B. : General case, lat and long function of both i and j indices: 253 ! e1t(ji,jj) = ra * rad * SQRT( ( cos( rad*gphit(ji,jj) ) * fsdila( zti, ztj ) )**2 & 254 ! + ( fsdiph( zti, ztj ) )**2 ) 255 ! e1u(ji,jj) = ra * rad * SQRT( ( cos( rad*gphiu(ji,jj) ) * fsdila( zui, zuj ) )**2 & 256 ! + ( fsdiph( zui, zuj ) )**2 ) 257 ! e1v(ji,jj) = ra * rad * SQRT( ( cos( rad*gphiv(ji,jj) ) * fsdila( zvi, zvj ) )**2 & 258 ! + ( fsdiph( zvi, zvj ) )**2 ) 259 ! e1f(ji,jj) = ra * rad * SQRT( ( cos( rad*gphif(ji,jj) ) * fsdila( zfi, zfj ) )**2 & 260 ! + ( fsdiph( zfi, zfj ) )**2 ) 261 ! 262 ! e2t(ji,jj) = ra * rad * SQRT( ( cos( rad*gphit(ji,jj) ) * fsdjla( zti, ztj ) )**2 & 263 ! + ( fsdjph( zti, ztj ) )**2 ) 264 ! e2u(ji,jj) = ra * rad * SQRT( ( cos( rad*gphiu(ji,jj) ) * fsdjla( zui, zuj ) )**2 & 265 ! + ( fsdjph( zui, zuj ) )**2 ) 266 ! e2v(ji,jj) = ra * rad * SQRT( ( cos( rad*gphiv(ji,jj) ) * fsdjla( zvi, zvj ) )**2 & 267 ! + ( fsdjph( zvi, zvj ) )**2 ) 268 ! e2f(ji,jj) = ra * rad * SQRT( ( cos( rad*gphif(ji,jj) ) * fsdjla( zfi, zfj ) )**2 & 269 ! + ( fsdjph( zfi, zfj ) )**2 ) 270 271 272 CASE ( 1 ) ! geographical mesh on the sphere with regular grid-spacing 273 145 ! 146 CASE ( 1 ) !== geographical mesh on the sphere with regular (in degree) grid-spacing ==! 147 ! 274 148 IF(lwp) WRITE(numout,*) 275 149 IF(lwp) WRITE(numout,*) ' geographical mesh on the sphere with regular grid-spacing' 276 150 IF(lwp) WRITE(numout,*) ' given by ppe1_deg and ppe2_deg' 277 151 ! 278 152 DO jj = 1, jpj 279 153 DO ji = 1, jpi 280 zti = FLOAT( ji - 1 + nimpp - 1 ) ; ztj = FLOAT( jj - 1 + njmpp - 1 )281 zui = FLOAT( ji - 1 + nimpp - 1 ) + 0.5 ; zuj = FLOAT( jj - 1 + njmpp - 1 )282 zvi = FLOAT( ji - 1 + nimpp - 1 ) ; zvj = FLOAT( jj - 1 + njmpp - 1 ) + 0.5283 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 284 158 ! Longitude 285 159 glamt(ji,jj) = ppglam0 + ppe1_deg * zti … … 304 178 END DO 305 179 END DO 306 307 308 CASE ( 2:3 ) ! f- or beta-plane with regular grid-spacing 309 180 ! 181 CASE ( 2:3 ) !== f- or beta-plane with regular grid-spacing ==! 182 ! 310 183 IF(lwp) WRITE(numout,*) 311 184 IF(lwp) WRITE(numout,*) ' f- or beta-plane with regular grid-spacing' 312 185 IF(lwp) WRITE(numout,*) ' given by ppe1_m and ppe2_m' 313 186 ! 314 187 ! Position coordinates (in kilometers) 315 188 ! ========== 316 glam0 = 0. e0189 glam0 = 0._wp 317 190 gphi0 = - ppe2_m * 1.e-3 318 191 ! 319 192 #if defined key_agrif 320 193 IF ( cp_cfg == 'eel' .AND. jp_cfg == 6 ) THEN ! for EEL6 configuration only … … 329 202 DO jj = 1, jpj 330 203 DO ji = 1, jpi 331 glamt(ji,jj) = glam0 + ppe1_m * 1.e-3 * ( FLOAT( ji - 1 + nimpp - 1 ) )332 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 ) 333 206 glamv(ji,jj) = glamt(ji,jj) 334 207 glamf(ji,jj) = glamu(ji,jj) 335 336 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 ) ) 337 210 gphiu(ji,jj) = gphit(ji,jj) 338 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 ) 339 212 gphif(ji,jj) = gphiv(ji,jj) 340 213 END DO 341 214 END DO 342 215 ! 343 216 ! Horizontal scale factors (in meters) 344 217 ! ====== … … 347 220 e1v(:,:) = ppe1_m ; e2v(:,:) = ppe2_m 348 221 e1f(:,:) = ppe1_m ; e2f(:,:) = ppe2_m 349 350 CASE ( 4 ) ! geographical mesh on the sphere, isotropic MERCATOR type351 222 ! 223 CASE ( 4 ) !== geographical mesh on the sphere, isotropic MERCATOR type ==! 224 ! 352 225 IF(lwp) WRITE(numout,*) 353 226 IF(lwp) WRITE(numout,*) ' geographical mesh on the sphere, MERCATOR type' 354 227 IF(lwp) WRITE(numout,*) ' longitudinal/latitudinal spacing given by ppe1_deg' 355 228 IF ( ppgphi0 == -90 ) CALL ctl_stop( ' Mercator grid cannot start at south pole !!!! ' ) 356 229 ! 357 230 ! Find index corresponding to the equator, given the grid spacing e1_deg 358 231 ! and the (approximate) southern latitude ppgphi0. … … 362 235 ijeq = ABS( 180./rpi * LOG( COS( zarg ) / SIN( zarg ) ) / ppe1_deg ) 363 236 IF( ppgphi0 > 0 ) ijeq = -ijeq 364 237 ! 365 238 IF(lwp) WRITE(numout,*) ' Index of the equator on the MERCATOR grid:', ijeq 366 239 ! 367 240 DO jj = 1, jpj 368 241 DO ji = 1, jpi 369 zti = FLOAT( ji - 1 + nimpp - 1 ) ; ztj = FLOAT( jj - ijeq + njmpp - 1 )370 zui = FLOAT( ji - 1 + nimpp - 1 ) + 0.5 ; zuj = FLOAT( jj - ijeq + njmpp - 1 )371 zvi = FLOAT( ji - 1 + nimpp - 1 ) ; zvj = FLOAT( jj - ijeq + njmpp - 1 ) + 0.5372 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 373 246 ! Longitude 374 247 glamt(ji,jj) = ppglam0 + ppe1_deg * zti … … 393 266 END DO 394 267 END DO 395 396 CASE ( 5 ) ! beta-plane with regular grid-spacing and rotated domain(GYRE configuration)397 268 ! 269 CASE ( 5 ) !== beta-plane with regular grid-spacing and rotated domain ==! (GYRE configuration) 270 ! 398 271 IF(lwp) WRITE(numout,*) 399 272 IF(lwp) WRITE(numout,*) ' beta-plane with regular grid-spacing and rotated domain (GYRE configuration)' 400 273 IF(lwp) WRITE(numout,*) ' given by ppe1_m and ppe2_m' 401 274 ! 402 275 ! Position coordinates (in kilometers) 403 276 ! ========== 404 277 ! 405 278 ! angle 45deg and ze1=106.e+3 / jp_cfg forced -> zlam1 = -85deg, zphi1 = 29degN 406 zlam1 = -85 407 zphi1 = 29279 zlam1 = -85._wp 280 zphi1 = 29._wp 408 281 ! resolution in meters 409 ze1 = 106000. / FLOAT(jp_cfg)282 ze1 = 106000. / REAL( jp_cfg , wp ) 410 283 ! benchmark: forced the resolution to be about 100 km 411 IF( nbench /= 0 ) ze1 = 106000. e0412 zsin_alpha = - SQRT( 2. ) / 2.413 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 414 287 ze1deg = ze1 / (ra * rad) 415 IF( nbench /= 0 ) ze1deg = ze1deg / FLOAT(jp_cfg)! benchmark: keep the lat/+lon416 ! ! at the right jp_cfg resolution417 glam0 = zlam1 + zcos_alpha * ze1deg * FLOAT( jpjglo-2)418 gphi0 = zphi1 + zsin_alpha * ze1deg * FLOAT( jpjglo-2)419 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 ! 420 293 IF( nprint==1 .AND. lwp ) THEN 421 294 WRITE(numout,*) ' ze1', ze1, 'cosalpha', zcos_alpha, 'sinalpha', zsin_alpha 422 295 WRITE(numout,*) ' ze1deg', ze1deg, 'glam0', glam0, 'gphi0', gphi0 423 296 ENDIF 424 297 ! 425 298 DO jj = 1, jpj 426 DO ji = 1, jpi427 zim1 = FLOAT( ji + nimpp - 1 ) - 1. ; zim05 = FLOAT( ji + nimpp - 1 ) - 1.5428 zjm1 = FLOAT( jj + njmpp - 1 ) - 1. ; zjm05 = FLOAT( jj + njmpp - 1 ) - 1.5429 430 glamf(ji,jj) = glam0 + zim1 * ze1deg * zcos_alpha + zjm1 * ze1deg * zsin_alpha431 gphif(ji,jj) = gphi0 - zim1 * ze1deg * zsin_alpha + zjm1 * ze1deg * zcos_alpha432 433 glamt(ji,jj) = glam0 + zim05 * ze1deg * zcos_alpha + zjm05 * ze1deg * zsin_alpha434 gphit(ji,jj) = gphi0 - zim05 * ze1deg * zsin_alpha + zjm05 * ze1deg * zcos_alpha435 436 glamu(ji,jj) = glam0 + zim1 * ze1deg * zcos_alpha + zjm05 * ze1deg * zsin_alpha437 gphiu(ji,jj) = gphi0 - zim1 * ze1deg * zsin_alpha + zjm05 * ze1deg * zcos_alpha438 439 glamv(ji,jj) = glam0 + zim05 * ze1deg * zcos_alpha + zjm1 * ze1deg * zsin_alpha440 gphiv(ji,jj) = gphi0 - zim05 * ze1deg * zsin_alpha + zjm1 * ze1deg * zcos_alpha441 END DO442 443 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 ! 444 317 ! Horizontal scale factors (in meters) 445 318 ! ====== … … 448 321 e1v(:,:) = ze1 ; e2v(:,:) = ze1 449 322 e1f(:,:) = ze1 ; e2f(:,:) = ze1 450 323 ! 451 324 CASE DEFAULT 452 325 WRITE(ctmp1,*) ' bad flag value for jphgr_msh = ', jphgr_msh 453 326 CALL ctl_stop( ctmp1 ) 454 327 ! 455 328 END SELECT 456 329 457 ! T-cell surface458 ! -------------- 459 e1e2t(:,:) = e1t(:,:) * e2t(:,:)460 461 ! Useful shortcuts (JC: note the duplicated e2e2t array ! Need some cleaning)462 ! ---------------------------------------------------------------------------463 e12t (:,:) = e1t(:,:) * e2t(:,:)464 e12u (:,:) = e1u(:,:) * e2u(:,:)465 e1 2v (:,:) = e1v(:,:) * e2v(:,:)466 e1 2f (:,:) = e1f(:,:) *e2f(:,:)467 r1_e12t (:,:) = 1._wp / e12t(:,:)468 r1_e12u (:,:) = 1._wp / e12u(:,:)469 r1_e12v (:,:) = 1._wp / e12v(:,:)470 r1_e12f (:,:) = 1._wp / e12f(:,:)471 r e2u_e1u(:,:) = e2u(:,:) / e1u(:,:)472 r e1v_e2v(:,:) = e1v(:,:) /e2v(:,:)473 474 ! Control printing : Grid informations (if not restart)475 ! ----------------476 477 IF( lwp .AND. .NOT.ln_rstart ) THEN330 ! 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. nn_print >=1 .AND. .NOT.ln_rstart ) THEN ! Control print : Grid informations (if not restart) 478 351 WRITE(numout,*) 479 352 WRITE(numout,*) ' longitude and e1 scale factors' … … 485 358 9300 FORMAT( 1x, i4, f8.2,1x, f8.2,1x, f8.2,1x, f8.2, 1x, & 486 359 f19.10, 1x, f19.10, 1x, f19.10, 1x, f19.10 ) 487 360 ! 488 361 WRITE(numout,*) 489 362 WRITE(numout,*) ' latitude and e2 scale factors' … … 495 368 ENDIF 496 369 497 498 IF( nprint == 1 .AND. lwp ) THEN499 WRITE(numout,*) ' e1u e2u '500 CALL prihre( e1u,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )501 CALL prihre( e2u,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )502 WRITE(numout,*) ' e1v e2v '503 CALL prihre( e1v,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )504 CALL prihre( e2v,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )505 WRITE(numout,*) ' e1f e2f '506 CALL prihre( e1f,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )507 CALL prihre( e2f,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )508 ENDIF509 510 370 511 371 ! ================= ! … … 514 374 515 375 SELECT CASE( jphgr_msh ) ! type of horizontal mesh 516 376 ! 517 377 CASE ( 0, 1, 4 ) ! mesh on the sphere 518 378 ! 519 379 ff(:,:) = 2. * omega * SIN( rad * gphif(:,:) ) 520 380 ! 521 381 CASE ( 2 ) ! f-plane at ppgphi0 522 382 ! 523 383 ff(:,:) = 2. * omega * SIN( rad * ppgphi0 ) 524 384 ! 525 385 IF(lwp) WRITE(numout,*) ' f-plane: Coriolis parameter = constant = ', ff(1,1) 526 386 ! 527 387 CASE ( 3 ) ! beta-plane 528 388 ! 529 389 zbeta = 2. * omega * COS( rad * ppgphi0 ) / ra ! beta at latitude ppgphi0 530 zphi0 = ppgphi0 - FLOAT( jpjglo/2) * ppe2_m / ( ra * rad ) ! latitude of the first row F-points531 390 zphi0 = ppgphi0 - REAL( jpjglo/2) * ppe2_m / ( ra * rad ) ! latitude of the first row F-points 391 ! 532 392 #if defined key_agrif 533 IF ( cp_cfg == 'eel' .AND. jp_cfg == 6 ) THEN! for EEL6 configuration only534 IF( .NOT. 535 zphi0 = ppgphi0 - FLOAT( Agrif_Parent(jpjglo)/2)*Agrif_Parent(ppe2_m) / (ra * rad)393 IF( cp_cfg == 'eel' .AND. jp_cfg == 6 ) THEN ! for EEL6 configuration only 394 IF( .NOT.Agrif_Root() ) THEN 395 zphi0 = ppgphi0 - REAL( Agrif_Parent(jpjglo)/2)*Agrif_Parent(ppe2_m) / (ra * rad) 536 396 ENDIF 537 397 ENDIF 538 398 #endif 539 399 zf0 = 2. * omega * SIN( rad * zphi0 ) ! compute f0 1st point south 540 400 ! 541 401 ff(:,:) = ( zf0 + zbeta * gphif(:,:) * 1.e+3 ) ! f = f0 +beta* y ( y=0 at south) 542 402 ! 543 403 IF(lwp) THEN 544 404 WRITE(numout,*) … … 553 413 IF(lwp) WRITE(numout,*) ' Coriolis parameter varies globally from ', zminff,' to ', zmaxff 554 414 END IF 555 415 ! 556 416 CASE ( 5 ) ! beta-plane and rotated domain (gyre configuration) 557 417 ! 558 418 zbeta = 2. * omega * COS( rad * ppgphi0 ) / ra ! beta at latitude ppgphi0 559 zphi0 = 15. e0! latitude of the first row F-points419 zphi0 = 15._wp ! latitude of the first row F-points 560 420 zf0 = 2. * omega * SIN( rad * zphi0 ) ! compute f0 1st point south 561 421 ! 562 422 ff(:,:) = ( zf0 + zbeta * ABS( gphif(:,:) - zphi0 ) * rad * ra ) ! f = f0 +beta* y ( y=0 at south) 563 423 ! 564 424 IF(lwp) THEN 565 425 WRITE(numout,*) … … 567 427 WRITE(numout,*) ' Coriolis parameter varies in this processor from ', ff(nldi,nldj),' to ', ff(nldi,nlej) 568 428 ENDIF 569 429 ! 570 430 IF( lk_mpp ) THEN 571 431 zminff=ff(nldi,nldj) … … 575 435 IF(lwp) WRITE(numout,*) ' Coriolis parameter varies globally from ', zminff,' to ', zmaxff 576 436 END IF 577 437 ! 578 438 END SELECT 579 439 … … 584 444 585 445 IF( nperio == 2 ) THEN 586 znorme = SQRT( SUM( gphiu(:,2) * gphiu(:,2) ) ) / FLOAT( jpi )446 znorme = SQRT( SUM( gphiu(:,2) * gphiu(:,2) ) ) / REAL( jpi ) 587 447 IF( znorme > 1.e-13 ) CALL ctl_stop( ' ===>>>> : symmetrical condition: rerun with good equator line' ) 588 448 ENDIF … … 593 453 594 454 595 SUBROUTINE hgr_read 455 SUBROUTINE hgr_read( ke1e2u_v ) 596 456 !!--------------------------------------------------------------------- 597 457 !! *** ROUTINE hgr_read *** 598 458 !! 599 !! ** Purpose : Read a coordinate file in NetCDF format 600 !! 601 !! ** Method : The mesh file has been defined trough a analytical 602 !! or semi-analytical method. It is read in a NetCDF file. 603 !! 459 !! ** Purpose : Read a coordinate file in NetCDF format using IOM 460 !! 604 461 !!---------------------------------------------------------------------- 605 462 USE iom 606 463 !! 464 INTEGER, INTENT( inout ) :: ke1e2u_v ! fag: e1e2u & e1e2v read in coordinate file (=1) or not (=0) 465 ! 607 466 INTEGER :: inum ! temporary logical unit 608 467 !!---------------------------------------------------------------------- 609 468 ! 610 469 IF(lwp) THEN 611 470 WRITE(numout,*) … … 613 472 WRITE(numout,*) '~~~~~~~~ jpiglo = ', jpiglo, ' jpjglo = ', jpjglo, ' jpk = ', jpk 614 473 ENDIF 615 474 ! 616 475 CALL iom_open( 'coordinates', inum ) 617 618 CALL iom_get( inum, jpdom_data, 'glamt', glamt ) 619 CALL iom_get( inum, jpdom_data, 'glamu', glamu ) 620 CALL iom_get( inum, jpdom_data, 'glamv', glamv ) 621 CALL iom_get( inum, jpdom_data, 'glamf', glamf ) 622 623 CALL iom_get( inum, jpdom_data, 'gphit', gphit ) 624 CALL iom_get( inum, jpdom_data, 'gphiu', gphiu ) 625 CALL iom_get( inum, jpdom_data, 'gphiv', gphiv ) 626 CALL iom_get( inum, jpdom_data, 'gphif', gphif ) 627 628 CALL iom_get( inum, jpdom_data, 'e1t', e1t ) 629 CALL iom_get( inum, jpdom_data, 'e1u', e1u ) 630 CALL iom_get( inum, jpdom_data, 'e1v', e1v ) 631 CALL iom_get( inum, jpdom_data, 'e1f', e1f ) 632 633 CALL iom_get( inum, jpdom_data, 'e2t', e2t ) 634 CALL iom_get( inum, jpdom_data, 'e2u', e2u ) 635 CALL iom_get( inum, jpdom_data, 'e2v', e2v ) 636 CALL iom_get( inum, jpdom_data, 'e2f', e2f ) 637 476 ! 477 CALL iom_get( inum, jpdom_data, 'glamt', glamt, lrowattr=ln_use_jattr ) 478 CALL iom_get( inum, jpdom_data, 'glamu', glamu, lrowattr=ln_use_jattr ) 479 CALL iom_get( inum, jpdom_data, 'glamv', glamv, lrowattr=ln_use_jattr ) 480 CALL iom_get( inum, jpdom_data, 'glamf', glamf, lrowattr=ln_use_jattr ) 481 ! 482 CALL iom_get( inum, jpdom_data, 'gphit', gphit, lrowattr=ln_use_jattr ) 483 CALL iom_get( inum, jpdom_data, 'gphiu', gphiu, lrowattr=ln_use_jattr ) 484 CALL iom_get( inum, jpdom_data, 'gphiv', gphiv, lrowattr=ln_use_jattr ) 485 CALL iom_get( inum, jpdom_data, 'gphif', gphif, lrowattr=ln_use_jattr ) 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 ! 638 506 CALL iom_close( inum ) 639 507
Note: See TracChangeset
for help on using the changeset viewer.
