- Timestamp:
- 2015-09-13T09:42:41+02:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90
r5656 r5737 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) 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 mesh 127 128 CASE ( 0 ) ! curvilinear coordinate on the sphere read in coordinate.nc file 129 127 ! 128 ie1e2u_v = 0 ! set to unread e1e2u and e1e2v 129 ! 130 SELECT CASE( jphgr_msh ) ! type of horizontal mesh 131 ! 132 CASE ( 0 ) !== read in coordinate.nc file ==! 133 ! 130 134 IF(lwp) WRITE(numout,*) 131 135 IF(lwp) WRITE(numout,*) ' curvilinear coordinate on the sphere read in "coordinate" file' 132 133 CALL hgr_read ! Defaultl option : NetCDF file 134 136 ! 137 CALL hgr_read( ie1e2u_v ) 138 ! 139 IF( ie1e2u_v == 0 ) THEN ! e1e2u and e1e2v have not been read: compute them 140 ! ! e2u and e1v does not include a reduction in some strait: apply reduction 141 e1e2u (:,:) = e1u(:,:) * e2u(:,:) 142 e1e2v (:,:) = e1v(:,:) * e2v(:,:) 143 144 ! 135 145 ! ! ===================== 136 146 IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! ORCA R2 configuration … … 157 167 ! 158 168 ENDIF 159 160 !! =====================169 ! 170 ! ! ===================== 161 171 IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN ! ORCA R1 configuration 162 172 ! ! ===================== 163 173 ! This dirty section will be suppressed by simplification process: all this will come back in input files 164 174 ! Currently these hard-wired indices relate to configuration with 165 ! extend grid (jpjglo=332) 166 ! which had a grid-size of 362x292. 175 ! extend grid (jpjglo=332) which had a grid-size of 362x292. 167 176 ! 168 177 isrow = 332 - jpjglo … … 208 217 IF(lwp) WRITE(numout,*) ' orca_r1: E Halmahera : e1v reduced to 50 km' 209 218 ! 210 ! 211 ENDIF 212 219 ENDIF 220 ! 213 221 ! ! ====================== 214 222 IF( cp_cfg == "orca" .AND. jp_cfg == 05 ) THEN ! ORCA R05 configuration … … 251 259 ! 252 260 ENDIF 253 254 261 262 ! ! create 'coordinate_e1e2u_v.nc' file that contains 263 ! ! reduced scale factor in some strait but full e1e2u and e1e2v surfaces 264 IF( ie1e2u_v == 0 ) CALL dom_wri_coordinate 265 ! 266 ! 267 ENDIF 268 269 270 ! 255 271 ! N.B. : General case, lat and long function of both i and j indices: 256 272 ! e1t(ji,jj) = ra * rad * SQRT( ( cos( rad*gphit(ji,jj) ) * fsdila( zti, ztj ) )**2 & … … 271 287 ! e2f(ji,jj) = ra * rad * SQRT( ( cos( rad*gphif(ji,jj) ) * fsdjla( zfi, zfj ) )**2 & 272 288 ! + ( fsdjph( zfi, zfj ) )**2 ) 273 274 275 CASE ( 1 ) ! geographical mesh on the sphere with regular grid-spacing276 289 ! 290 ! 291 CASE ( 1 ) !== geographical mesh on the sphere with regular (in degree) grid-spacing ==! 292 ! 277 293 IF(lwp) WRITE(numout,*) 278 294 IF(lwp) WRITE(numout,*) ' geographical mesh on the sphere with regular grid-spacing' 279 295 IF(lwp) WRITE(numout,*) ' given by ppe1_deg and ppe2_deg' 280 296 ! 281 297 DO jj = 1, jpj 282 298 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.5299 zti = REAL( ji - 1 + nimpp - 1 ) ; ztj = REAL( jj - 1 + njmpp - 1 ) 300 zui = REAL( ji - 1 + nimpp - 1 ) + 0.5 ; zuj = REAL( jj - 1 + njmpp - 1 ) 301 zvi = REAL( ji - 1 + nimpp - 1 ) ; zvj = REAL( jj - 1 + njmpp - 1 ) + 0.5 302 zfi = REAL( ji - 1 + nimpp - 1 ) + 0.5 ; zfj = REAL( jj - 1 + njmpp - 1 ) + 0.5 287 303 ! Longitude 288 304 glamt(ji,jj) = ppglam0 + ppe1_deg * zti … … 307 323 END DO 308 324 END DO 309 310 311 CASE ( 2:3 ) ! f- or beta-plane with regular grid-spacing 312 325 ! 326 CASE ( 2:3 ) !== f- or beta-plane with regular grid-spacing ==! 327 ! 313 328 IF(lwp) WRITE(numout,*) 314 329 IF(lwp) WRITE(numout,*) ' f- or beta-plane with regular grid-spacing' 315 330 IF(lwp) WRITE(numout,*) ' given by ppe1_m and ppe2_m' 316 331 ! 317 332 ! Position coordinates (in kilometers) 318 333 ! ========== 319 334 glam0 = 0.e0 320 335 gphi0 = - ppe2_m * 1.e-3 321 336 ! 322 337 #if defined key_agrif 323 338 IF ( cp_cfg == 'eel' .AND. jp_cfg == 6 ) THEN ! for EEL6 configuration only … … 332 347 DO jj = 1, jpj 333 348 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 )349 glamt(ji,jj) = glam0 + ppe1_m * 1.e-3 * ( REAL( ji - 1 + nimpp - 1 ) ) 350 glamu(ji,jj) = glam0 + ppe1_m * 1.e-3 * ( REAL( ji - 1 + nimpp - 1 ) + 0.5 ) 336 351 glamv(ji,jj) = glamt(ji,jj) 337 352 glamf(ji,jj) = glamu(ji,jj) 338 339 gphit(ji,jj) = gphi0 + ppe2_m * 1.e-3 * ( FLOAT( jj - 1 + njmpp - 1 ) )353 ! 354 gphit(ji,jj) = gphi0 + ppe2_m * 1.e-3 * ( REAL( jj - 1 + njmpp - 1 ) ) 340 355 gphiu(ji,jj) = gphit(ji,jj) 341 gphiv(ji,jj) = gphi0 + ppe2_m * 1.e-3 * ( FLOAT( jj - 1 + njmpp - 1 ) + 0.5 )356 gphiv(ji,jj) = gphi0 + ppe2_m * 1.e-3 * ( REAL( jj - 1 + njmpp - 1 ) + 0.5 ) 342 357 gphif(ji,jj) = gphiv(ji,jj) 343 358 END DO 344 359 END DO 345 360 ! 346 361 ! Horizontal scale factors (in meters) 347 362 ! ====== … … 350 365 e1v(:,:) = ppe1_m ; e2v(:,:) = ppe2_m 351 366 e1f(:,:) = ppe1_m ; e2f(:,:) = ppe2_m 352 353 CASE ( 4 ) ! geographical mesh on the sphere, isotropic MERCATOR type354 367 ! 368 CASE ( 4 ) !== geographical mesh on the sphere, isotropic MERCATOR type ==! 369 ! 355 370 IF(lwp) WRITE(numout,*) 356 371 IF(lwp) WRITE(numout,*) ' geographical mesh on the sphere, MERCATOR type' 357 372 IF(lwp) WRITE(numout,*) ' longitudinal/latitudinal spacing given by ppe1_deg' 358 373 IF ( ppgphi0 == -90 ) CALL ctl_stop( ' Mercator grid cannot start at south pole !!!! ' ) 359 374 ! 360 375 ! Find index corresponding to the equator, given the grid spacing e1_deg 361 376 ! and the (approximate) southern latitude ppgphi0. … … 365 380 ijeq = ABS( 180./rpi * LOG( COS( zarg ) / SIN( zarg ) ) / ppe1_deg ) 366 381 IF( ppgphi0 > 0 ) ijeq = -ijeq 367 382 ! 368 383 IF(lwp) WRITE(numout,*) ' Index of the equator on the MERCATOR grid:', ijeq 369 384 ! 370 385 DO jj = 1, jpj 371 386 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.5387 zti = REAL( ji - 1 + nimpp - 1 ) ; ztj = REAL( jj - ijeq + njmpp - 1 ) 388 zui = REAL( ji - 1 + nimpp - 1 ) + 0.5 ; zuj = REAL( jj - ijeq + njmpp - 1 ) 389 zvi = REAL( ji - 1 + nimpp - 1 ) ; zvj = REAL( jj - ijeq + njmpp - 1 ) + 0.5 390 zfi = REAL( ji - 1 + nimpp - 1 ) + 0.5 ; zfj = REAL( jj - ijeq + njmpp - 1 ) + 0.5 376 391 ! Longitude 377 392 glamt(ji,jj) = ppglam0 + ppe1_deg * zti … … 396 411 END DO 397 412 END DO 398 399 CASE ( 5 ) ! beta-plane with regular grid-spacing and rotated domain(GYRE configuration)400 413 ! 414 CASE ( 5 ) !== beta-plane with regular grid-spacing and rotated domain ==! (GYRE configuration) 415 ! 401 416 IF(lwp) WRITE(numout,*) 402 417 IF(lwp) WRITE(numout,*) ' beta-plane with regular grid-spacing and rotated domain (GYRE configuration)' 403 418 IF(lwp) WRITE(numout,*) ' given by ppe1_m and ppe2_m' 404 419 ! 405 420 ! Position coordinates (in kilometers) 406 421 ! ========== 407 422 ! 408 423 ! angle 45deg and ze1=106.e+3 / jp_cfg forced -> zlam1 = -85deg, zphi1 = 29degN 409 zlam1 = -85 410 zphi1 = 29424 zlam1 = -85._wp 425 zphi1 = 29._wp 411 426 ! resolution in meters 412 ze1 = 106000. / FLOAT(jp_cfg)427 ze1 = 106000. / REAL( jp_cfg , wp ) 413 428 ! benchmark: forced the resolution to be about 100 km 414 429 IF( nbench /= 0 ) ze1 = 106000.e0 415 zsin_alpha = - SQRT( 2. ) / 2.416 zcos_alpha = SQRT( 2. ) / 2.430 zsin_alpha = - SQRT( 2._wp ) * 0.5_wp 431 zcos_alpha = SQRT( 2._wp ) * 0.5_wp 417 432 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 433 IF( nbench /= 0 ) ze1deg = ze1deg / REAL( jp_cfg , wp ) ! benchmark: keep the lat/+lon 434 ! ! at the right jp_cfg resolution 435 glam0 = zlam1 + zcos_alpha * ze1deg * REAL( jpjglo-2 , wp ) 436 gphi0 = zphi1 + zsin_alpha * ze1deg * REAL( jpjglo-2 , wp ) 437 ! 423 438 IF( nprint==1 .AND. lwp ) THEN 424 439 WRITE(numout,*) ' ze1', ze1, 'cosalpha', zcos_alpha, 'sinalpha', zsin_alpha 425 440 WRITE(numout,*) ' ze1deg', ze1deg, 'glam0', glam0, 'gphi0', gphi0 426 441 ENDIF 427 442 ! 428 443 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 444 DO ji = 1, jpi 445 zim1 = REAL( ji + nimpp - 1 ) - 1. ; zim05 = REAL( ji + nimpp - 1 ) - 1.5 446 zjm1 = REAL( jj + njmpp - 1 ) - 1. ; zjm05 = REAL( jj + njmpp - 1 ) - 1.5 447 ! 448 glamf(ji,jj) = glam0 + zim1 * ze1deg * zcos_alpha + zjm1 * ze1deg * zsin_alpha 449 gphif(ji,jj) = gphi0 - zim1 * ze1deg * zsin_alpha + zjm1 * ze1deg * zcos_alpha 450 ! 451 glamt(ji,jj) = glam0 + zim05 * ze1deg * zcos_alpha + zjm05 * ze1deg * zsin_alpha 452 gphit(ji,jj) = gphi0 - zim05 * ze1deg * zsin_alpha + zjm05 * ze1deg * zcos_alpha 453 ! 454 glamu(ji,jj) = glam0 + zim1 * ze1deg * zcos_alpha + zjm05 * ze1deg * zsin_alpha 455 gphiu(ji,jj) = gphi0 - zim1 * ze1deg * zsin_alpha + zjm05 * ze1deg * zcos_alpha 456 ! 457 glamv(ji,jj) = glam0 + zim05 * ze1deg * zcos_alpha + zjm1 * ze1deg * zsin_alpha 458 gphiv(ji,jj) = gphi0 - zim05 * ze1deg * zsin_alpha + zjm1 * ze1deg * zcos_alpha 459 END DO 460 END DO 461 ! 447 462 ! Horizontal scale factors (in meters) 448 463 ! ====== … … 451 466 e1v(:,:) = ze1 ; e2v(:,:) = ze1 452 467 e1f(:,:) = ze1 ; e2f(:,:) = ze1 453 468 ! 454 469 CASE DEFAULT 455 470 WRITE(ctmp1,*) ' bad flag value for jphgr_msh = ', jphgr_msh 456 471 CALL ctl_stop( ctmp1 ) 457 472 ! 458 473 END SELECT 459 474 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 475 ! associated horizontal metrics 476 ! ----------------------------- 477 ! 478 r1_e1t(:,:) = 1._wp / e1t(:,:) ; r1_e2t (:,:) = 1._wp / e2t(:,:) 479 r1_e1u(:,:) = 1._wp / e1u(:,:) ; r1_e2u (:,:) = 1._wp / e2u(:,:) 480 r1_e1v(:,:) = 1._wp / e1v(:,:) ; r1_e2v (:,:) = 1._wp / e2v(:,:) 481 r1_e1f(:,:) = 1._wp / e1f(:,:) ; r1_e2f (:,:) = 1._wp / e2f(:,:) 482 ! 483 e1e2t (:,:) = e1t(:,:) * e2t(:,:) ; r1_e1e2t(:,:) = 1._wp / e1e2t(:,:) 484 e1e2f (:,:) = e1f(:,:) * e2f(:,:) ; r1_e1e2f(:,:) = 1._wp / e1e2f(:,:) 485 IF( jphgr_msh /= 0 ) THEN ! e1e2u and e1e2v have not been set: compute them 486 e1e2u (:,:) = e1u(:,:) * e2u(:,:) 487 e1e2v (:,:) = e1v(:,:) * e2v(:,:) 488 ENDIF 489 r1_e1e2u(:,:) = 1._wp / e1e2u(:,:) ! compute their invert in both cases 490 r1_e1e2v(:,:) = 1._wp / e1e2v(:,:) 491 ! 492 e2_e1u(:,:) = e2u(:,:) / e1u(:,:) 493 e1_e2v(:,:) = e1v(:,:) / e2v(:,:) 494 495 IF( lwp .AND. .NOT.ln_rstart ) THEN ! Control print : Grid informations (if not restart) 489 496 WRITE(numout,*) 490 497 WRITE(numout,*) ' longitude and e1 scale factors' … … 496 503 9300 FORMAT( 1x, i4, f8.2,1x, f8.2,1x, f8.2,1x, f8.2, 1x, & 497 504 f19.10, 1x, f19.10, 1x, f19.10, 1x, f19.10 ) 498 505 ! 499 506 WRITE(numout,*) 500 507 WRITE(numout,*) ' latitude and e2 scale factors' … … 506 513 ENDIF 507 514 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 515 522 516 ! ================= ! … … 525 519 526 520 SELECT CASE( jphgr_msh ) ! type of horizontal mesh 527 521 ! 528 522 CASE ( 0, 1, 4 ) ! mesh on the sphere 529 523 ! 530 524 ff(:,:) = 2. * omega * SIN( rad * gphif(:,:) ) 531 525 ! 532 526 CASE ( 2 ) ! f-plane at ppgphi0 533 527 ! 534 528 ff(:,:) = 2. * omega * SIN( rad * ppgphi0 ) 535 529 ! 536 530 IF(lwp) WRITE(numout,*) ' f-plane: Coriolis parameter = constant = ', ff(1,1) 537 531 ! 538 532 CASE ( 3 ) ! beta-plane 539 533 ! 540 534 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 535 zphi0 = ppgphi0 - REAL( jpjglo/2) * ppe2_m / ( ra * rad ) ! latitude of the first row F-points 536 ! 543 537 #if defined key_agrif 544 538 IF ( cp_cfg == 'eel' .AND. jp_cfg == 6 ) THEN ! for EEL6 configuration only 545 539 IF( .NOT. Agrif_Root() ) THEN 546 zphi0 = ppgphi0 - FLOAT( Agrif_Parent(jpjglo)/2)*Agrif_Parent(ppe2_m) & 547 & / (ra * rad) 540 zphi0 = ppgphi0 - REAL( Agrif_Parent(jpjglo)/2)*Agrif_Parent(ppe2_m) / (ra * rad) 548 541 ENDIF 549 542 ENDIF 550 543 #endif 551 544 zf0 = 2. * omega * SIN( rad * zphi0 ) ! compute f0 1st point south 552 545 ! 553 546 ff(:,:) = ( zf0 + zbeta * gphif(:,:) * 1.e+3 ) ! f = f0 +beta* y ( y=0 at south) 554 547 ! 555 548 IF(lwp) THEN 556 549 WRITE(numout,*) … … 565 558 IF(lwp) WRITE(numout,*) ' Coriolis parameter varies globally from ', zminff,' to ', zmaxff 566 559 END IF 567 560 ! 568 561 CASE ( 5 ) ! beta-plane and rotated domain (gyre configuration) 569 562 ! 570 563 zbeta = 2. * omega * COS( rad * ppgphi0 ) / ra ! beta at latitude ppgphi0 571 564 zphi0 = 15.e0 ! latitude of the first row F-points 572 565 zf0 = 2. * omega * SIN( rad * zphi0 ) ! compute f0 1st point south 573 566 ! 574 567 ff(:,:) = ( zf0 + zbeta * ABS( gphif(:,:) - zphi0 ) * rad * ra ) ! f = f0 +beta* y ( y=0 at south) 575 568 ! 576 569 IF(lwp) THEN 577 570 WRITE(numout,*) … … 579 572 WRITE(numout,*) ' Coriolis parameter varies in this processor from ', ff(nldi,nldj),' to ', ff(nldi,nlej) 580 573 ENDIF 581 574 ! 582 575 IF( lk_mpp ) THEN 583 576 zminff=ff(nldi,nldj) … … 587 580 IF(lwp) WRITE(numout,*) ' Coriolis parameter varies globally from ', zminff,' to ', zmaxff 588 581 END IF 589 582 ! 590 583 END SELECT 591 584 … … 596 589 597 590 IF( nperio == 2 ) THEN 598 znorme = SQRT( SUM( gphiu(:,2) * gphiu(:,2) ) ) / FLOAT( jpi )591 znorme = SQRT( SUM( gphiu(:,2) * gphiu(:,2) ) ) / REAL( jpi ) 599 592 IF( znorme > 1.e-13 ) CALL ctl_stop( ' ===>>>> : symmetrical condition: rerun with good equator line' ) 600 593 ENDIF … … 605 598 606 599 607 SUBROUTINE hgr_read 600 SUBROUTINE hgr_read( ke1e2u_v ) 608 601 !!--------------------------------------------------------------------- 609 602 !! *** ROUTINE hgr_read *** 610 603 !! 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 !! 604 !! ** Purpose : Read a coordinate file in NetCDF format using IOM 605 !! 616 606 !!---------------------------------------------------------------------- 617 607 USE iom 618 608 !! 609 INTEGER, INTENT( inout ) :: ke1e2u_v ! fag: e1e2u & e1e2v read in coordinate file (=1) or not (=0) 610 ! 619 611 INTEGER :: inum ! temporary logical unit 620 612 !!---------------------------------------------------------------------- 621 613 ! 622 614 IF(lwp) THEN 623 615 WRITE(numout,*) … … 625 617 WRITE(numout,*) '~~~~~~~~ jpiglo = ', jpiglo, ' jpjglo = ', jpjglo, ' jpk = ', jpk 626 618 ENDIF 627 619 ! 628 620 CALL iom_open( 'coordinates', inum ) 629 621 ! 630 622 CALL iom_get( inum, jpdom_data, 'glamt', glamt, lrowattr=ln_use_jattr ) 631 623 CALL iom_get( inum, jpdom_data, 'glamu', glamu, lrowattr=ln_use_jattr ) 632 624 CALL iom_get( inum, jpdom_data, 'glamv', glamv, lrowattr=ln_use_jattr ) 633 625 CALL iom_get( inum, jpdom_data, 'glamf', glamf, lrowattr=ln_use_jattr ) 634 626 ! 635 627 CALL iom_get( inum, jpdom_data, 'gphit', gphit, lrowattr=ln_use_jattr ) 636 628 CALL iom_get( inum, jpdom_data, 'gphiu', gphiu, lrowattr=ln_use_jattr ) 637 629 CALL iom_get( inum, jpdom_data, 'gphiv', gphiv, lrowattr=ln_use_jattr ) 638 630 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 631 ! 632 CALL iom_get( inum, jpdom_data, 'e1t' , e1t , lrowattr=ln_use_jattr ) 633 CALL iom_get( inum, jpdom_data, 'e1u' , e1u , lrowattr=ln_use_jattr ) 634 CALL iom_get( inum, jpdom_data, 'e1v' , e1v , lrowattr=ln_use_jattr ) 635 CALL iom_get( inum, jpdom_data, 'e1f' , e1f , lrowattr=ln_use_jattr ) 636 ! 637 CALL iom_get( inum, jpdom_data, 'e2t' , e2t , lrowattr=ln_use_jattr ) 638 CALL iom_get( inum, jpdom_data, 'e2u' , e2u , lrowattr=ln_use_jattr ) 639 CALL iom_get( inum, jpdom_data, 'e2v' , e2v , lrowattr=ln_use_jattr ) 640 CALL iom_get( inum, jpdom_data, 'e2f' , e2f , lrowattr=ln_use_jattr ) 641 ! 642 IF( iom_varid( inum, 'e1e2u', ldstop = .FALSE. ) > 0 ) THEN 643 IF(lwp) WRITE(numout,*) 'hgr_read : e1e2u & e1e2v read in coordinates file' 644 CALL iom_get( inum, jpdom_data, 'e1e2u' , e1e2u , lrowattr=ln_use_jattr ) 645 CALL iom_get( inum, jpdom_data, 'e1e2v' , e1e2v , lrowattr=ln_use_jattr ) 646 ke1e2u_v = 1 647 ELSE 648 ke1e2u_v = 0 649 ENDIF 650 ! 650 651 CALL iom_close( inum ) 651 652 653 !!gm THIS is TO BE REMOVED !!!!!!! 654 652 655 ! need to be define for the extended grid south of -80S 653 656 ! some point are undefined but you need to have e1 and e2 .NE. 0 … … 676 679 e2f=1.0e2 677 680 END WHERE 681 !!gm end 678 682 679 683 END SUBROUTINE hgr_read
Note: See TracChangeset
for help on using the changeset viewer.
