Changeset 5836 for trunk/NEMOGCM/NEMO/OPA_SRC/DOM
- Timestamp:
- 2015-10-26T15:49:40+01:00 (9 years ago)
- Location:
- trunk/NEMOGCM/NEMO/OPA_SRC/DOM
- Files:
-
- 11 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/OPA_SRC/DOM/closea.F90
r5506 r5836 71 71 !! =2 put at location runoff 72 72 !!---------------------------------------------------------------------- 73 INTEGER :: jc 74 INTEGER :: isrow! local index75 !!---------------------------------------------------------------------- 76 73 INTEGER :: jc ! dummy loop indices 74 INTEGER :: isrow ! local index 75 !!---------------------------------------------------------------------- 76 ! 77 77 IF(lwp) WRITE(numout,*) 78 78 IF(lwp) WRITE(numout,*)'dom_clo : closed seas ' 79 79 IF(lwp) WRITE(numout,*)'~~~~~~~' 80 80 ! 81 81 ! initial values 82 82 ncsnr(:) = 1 ; ncsi1(:) = 1 ; ncsi2(:) = 1 ; ncsir(:,:) = 1 83 83 ncstt(:) = 0 ; ncsj1(:) = 1 ; ncsj2(:) = 1 ; ncsjr(:,:) = 1 84 84 ! 85 85 ! set the closed seas (in data domain indices) 86 86 ! ------------------- 87 87 ! 88 88 IF( cp_cfg == "orca" ) THEN 89 89 ! -
trunk/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90
r5123 r5836 7 7 !! History : 1.0 ! 2005-10 (A. Beckmann, G. Madec) reactivate s-coordinate 8 8 !! 3.3 ! 2010-11 (G. Madec) add mbk. arrays associated to the deepest ocean level 9 !! 4.0! 2011-01 (A. R. Porter, STFC Daresbury) dynamical allocation9 !! 3.4 ! 2011-01 (A. R. Porter, STFC Daresbury) dynamical allocation 10 10 !! 3.5 ! 2012 (S. Mocavero, I. Epicoco) Add arrays associated 11 11 !! to the optimization of BDY communications 12 !! 3.7 ! 2015-11 (G. Madec) introduce surface and scale factor ratio 12 13 !!---------------------------------------------------------------------- 13 14 … … 20 21 21 22 IMPLICIT NONE 22 PUBLIC ! allows the acces to par_oce when dom_oce is used 23 ! ! exception to coding rules... to be suppressed ??? 23 PUBLIC ! allows the acces to par_oce when dom_oce is used (exception to coding rules) 24 24 25 25 PUBLIC dom_oce_alloc ! Called from nemogcm.F90 … … 107 107 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: r2dtra !: = 2*rdttra except at nit000 (=rdttra) if neuler=0 108 108 109 ! !!* Namelist namcla : cross land advection110 INTEGER, PUBLIC :: nn_cla !: =1 cross land advection for exchanges through some straits (ORCA2)111 112 109 !!---------------------------------------------------------------------- 113 110 !! space domain parameters … … 158 155 !! horizontal curvilinear coordinate and scale factors 159 156 !! --------------------------------------------------------------------- 160 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: glamt, glamu !: longitude of t-, u-, v- and f-points (degre) 161 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: glamv, glamf !: 162 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: gphit, gphiu !: latitude of t-, u-, v- and f-points (degre) 163 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: gphiv, gphif !: 164 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) :: e1t, e2t, r1_e1t, r1_e2t !: horizontal scale factors and inverse at t-point (m) 165 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) :: e1u, e2u, r1_e1u, r1_e2u !: horizontal scale factors and inverse at u-point (m) 166 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) :: e1v, e2v, r1_e1v, r1_e2v !: horizontal scale factors and inverse at v-point (m) 167 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) :: e1f, e2f, r1_e1f, r1_e2f !: horizontal scale factors and inverse at f-point (m) 168 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: e1e2t !: surface at t-point (m2) 169 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ff !: coriolis factor (2.*omega*sin(yphi) ) (s-1) 157 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: glamt , glamu, glamv , glamf !: longitude at t, u, v, f-points [degree] 158 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: gphit , gphiu, gphiv , gphif !: latitude at t, u, v, f-points [degree] 159 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) :: e1t , e2t , r1_e1t, r1_e2t !: t-point horizontal scale factors [m] 160 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) :: e1u , e2u , r1_e1u, r1_e2u !: horizontal scale factors at u-point [m] 161 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) :: e1v , e2v , r1_e1v, r1_e2v !: horizontal scale factors at v-point [m] 162 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) :: e1f , e2f , r1_e1f, r1_e2f !: horizontal scale factors at f-point [m] 163 ! 164 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: e1e2t , r1_e1e2t !: associated metrics at t-point 165 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: e1e2u , r1_e1e2u , e2_e1u !: associated metrics at u-point 166 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: e1e2v , r1_e1e2v , e1_e2v !: associated metrics at v-point 167 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: e1e2f , r1_e1e2f !: associated metrics at f-point 168 ! 169 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ff !: coriolis factor [1/s] 170 170 171 171 !!---------------------------------------------------------------------- … … 216 216 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ht_0 !: reference depth at t- points (meters) 217 217 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hu_0 , hv_0 !: reference depth at u- and v-points (meters) 218 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: re2u_e1u !: scale factor coeffs at u points (e2u/e1u)219 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: re1v_e2v !: scale factor coeffs at v points (e1v/e2v)220 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: e12t , r1_e12t !: horizontal cell surface and inverse at t points221 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: e12u , r1_e12u !: horizontal cell surface and inverse at u points222 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: e12v , r1_e12v !: horizontal cell surface and inverse at v points223 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: e12f , r1_e12f !: horizontal cell surface and inverse at f points224 218 225 219 INTEGER, PUBLIC :: nla10 !: deepest W level Above ~10m (nlb10 - 1) … … 265 259 266 260 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: tpol, fpol !: north fold mask (jperio= 3 or 4) 267 268 #if defined key_noslip_accurate269 INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,: ) :: npcoa !: ???270 INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: nicoa, njcoa !: ???271 #endif272 261 273 262 !!---------------------------------------------------------------------- … … 333 322 INTEGER FUNCTION dom_oce_alloc() 334 323 !!---------------------------------------------------------------------- 335 INTEGER, DIMENSION(1 2) :: ierr324 INTEGER, DIMENSION(13) :: ierr 336 325 !!---------------------------------------------------------------------- 337 326 ierr(:) = 0 … … 346 335 & tpol(jpiglo) , fpol(jpiglo) , STAT=ierr(2) ) 347 336 ! 348 ALLOCATE( glamt(jpi,jpj) , gphit(jpi,jpj) , e1t(jpi,jpj) , e2t(jpi,jpj) , r1_e1t(jpi,jpj) , r1_e2t(jpi,jpj) , & 349 & glamu(jpi,jpj) , gphiu(jpi,jpj) , e1u(jpi,jpj) , e2u(jpi,jpj) , r1_e1u(jpi,jpj) , r1_e2u(jpi,jpj) , & 350 & glamv(jpi,jpj) , gphiv(jpi,jpj) , e1v(jpi,jpj) , e2v(jpi,jpj) , r1_e1v(jpi,jpj) , r1_e2v(jpi,jpj) , & 351 & glamf(jpi,jpj) , gphif(jpi,jpj) , e1f(jpi,jpj) , e2f(jpi,jpj) , r1_e1f(jpi,jpj) , r1_e2f(jpi,jpj) , & 352 & e1e2t(jpi,jpj) , ff (jpi,jpj) , STAT=ierr(3) ) 337 ALLOCATE( glamt(jpi,jpj) , glamu(jpi,jpj) , glamv(jpi,jpj) , glamf(jpi,jpj) , & 338 & gphit(jpi,jpj) , gphiu(jpi,jpj) , gphiv(jpi,jpj) , gphif(jpi,jpj) , & 339 & e1t (jpi,jpj) , e2t (jpi,jpj) , r1_e1t(jpi,jpj) , r1_e2t(jpi,jpj) , & 340 & e1u (jpi,jpj) , e2u (jpi,jpj) , r1_e1u(jpi,jpj) , r1_e2u(jpi,jpj) , & 341 & e1v (jpi,jpj) , e2v (jpi,jpj) , r1_e1v(jpi,jpj) , r1_e2v(jpi,jpj) , & 342 & e1f (jpi,jpj) , e2f (jpi,jpj) , r1_e1f(jpi,jpj) , r1_e2f(jpi,jpj) , & 343 & e1e2t(jpi,jpj) , r1_e1e2t(jpi,jpj) , & 344 & e1e2u(jpi,jpj) , r1_e1e2u(jpi,jpj) , e2_e1u(jpi,jpj) , & 345 & e1e2v(jpi,jpj) , r1_e1e2v(jpi,jpj) , e1_e2v(jpi,jpj) , & 346 & e1e2f(jpi,jpj) , r1_e1e2f(jpi,jpj) , & 347 & ff (jpi,jpj) , STAT=ierr(3) ) 353 348 ! 354 349 ALLOCATE( gdep3w_0(jpi,jpj,jpk) , e3v_0(jpi,jpj,jpk) , e3f_0 (jpi,jpj,jpk) , & … … 364 359 & gdept_b (jpi,jpj,jpk) ,gdepw_b(jpi,jpj,jpk) , e3w_b (jpi,jpj,jpk) , & 365 360 & e3t_a (jpi,jpj,jpk) , e3u_a (jpi,jpj,jpk) , e3v_a (jpi,jpj,jpk) , & 366 & ehu_a (jpi,jpj) , ehv_a(jpi,jpj), &367 & ehur_a (jpi,jpj) , ehvr_a(jpi,jpj), &368 & ehu_b (jpi,jpj) , ehv_b(jpi,jpj), &369 & ehur_b (jpi,jpj) , ehvr_b(jpi,jpj), STAT=ierr(5) )361 & ehu_a (jpi,jpj) , ehv_a (jpi,jpj), & 362 & ehur_a (jpi,jpj) , ehvr_a(jpi,jpj), & 363 & ehu_b (jpi,jpj) , ehv_b (jpi,jpj), & 364 & ehur_b (jpi,jpj) , ehvr_b(jpi,jpj), STAT=ierr(5) ) 370 365 #endif 371 366 ! 372 ALLOCATE( hu (jpi,jpj) , hur (jpi,jpj) , hu_0(jpi,jpj) , ht_0 (jpi,jpj) , & 373 & hv (jpi,jpj) , hvr (jpi,jpj) , hv_0(jpi,jpj) , ht (jpi,jpj) , & 374 & re2u_e1u(jpi,jpj) , re1v_e2v(jpi,jpj) , & 375 & e12t (jpi,jpj) , r1_e12t (jpi,jpj) , & 376 & e12u (jpi,jpj) , r1_e12u (jpi,jpj) , & 377 & e12v (jpi,jpj) , r1_e12v (jpi,jpj) , & 378 & e12f (jpi,jpj) , r1_e12f (jpi,jpj) , STAT=ierr(6) ) 367 ALLOCATE( hu(jpi,jpj) , hur(jpi,jpj) , hu_0(jpi,jpj) , ht_0(jpi,jpj) , & 368 & hv(jpi,jpj) , hvr(jpi,jpj) , hv_0(jpi,jpj) , ht (jpi,jpj) , STAT=ierr(6) ) 379 369 ! 380 370 ALLOCATE( gdept_1d(jpk) , gdepw_1d(jpk) , & … … 387 377 & scosrf(jpi,jpj) , scobot(jpi,jpj) , & 388 378 & hifv (jpi,jpj) , hiff (jpi,jpj) , & 389 & hift (jpi,jpj) , hifu (jpi,jpj) , rx1 379 & hift (jpi,jpj) , hifu (jpi,jpj) , rx1(jpi,jpj) , STAT=ierr(8) ) 390 380 391 381 ALLOCATE( mbathy(jpi,jpj) , bathy(jpi,jpj) , & 392 382 & tmask_i(jpi,jpj) , umask_i(jpi,jpj), vmask_i(jpi,jpj), fmask_i(jpi,jpj), & 393 & bmask (jpi,jpj), &383 & bmask (jpi,jpj) , & 394 384 & mbkt (jpi,jpj) , mbku (jpi,jpj) , mbkv(jpi,jpj) , STAT=ierr(9) ) 395 385 396 386 ! (ISF) Allocation of basic array 397 ALLOCATE( misfdep(jpi,jpj) , risfdep(jpi,jpj), &398 & mikt(jpi,jpj), miku(jpi,jpj), mikv(jpi,jpj) , &399 & mikf(jpi,jpj), ssmask(jpi,jpj), STAT=ierr(10) )387 ALLOCATE( misfdep(jpi,jpj) , risfdep(jpi,jpj), & 388 & mikt(jpi,jpj), miku(jpi,jpj), mikv(jpi,jpj) , & 389 & mikf(jpi,jpj), ssmask(jpi,jpj), STAT=ierr(10) ) 400 390 401 391 ALLOCATE( tmask(jpi,jpj,jpk) , umask(jpi,jpj,jpk), & … … 403 393 404 394 ALLOCATE( wmask(jpi,jpj,jpk) , wumask(jpi,jpj,jpk), wvmask(jpi,jpj,jpk) , STAT=ierr(12) ) 405 406 #if defined key_noslip_accurate407 ALLOCATE( npcoa(4,jpk), nicoa(2*(jpi+jpj),4,jpk), njcoa(2*(jpi+jpj),4,jpk), STAT=ierr(12) )408 #endif409 395 ! 410 396 dom_oce_alloc = MAXVAL(ierr) -
trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90
r5363 r5836 19 19 !! dom_nam : read and contral domain namelists 20 20 !! dom_ctl : control print for the ocean domain 21 !! dom_stiff : diagnose maximum grid stiffness/hydrostatic consistency (s-coordinate) 21 22 !!---------------------------------------------------------------------- 22 23 USE oce ! ocean variables … … 25 26 USE phycst ! physical constants 26 27 USE closea ! closed seas 27 USE in_out_manager ! I/O manager28 USE lib_mpp ! distributed memory computing library29 30 28 USE domhgr ! domain: set the horizontal mesh 31 29 USE domzgr ! domain: set the vertical mesh … … 36 34 USE c1d ! 1D vertical configuration 37 35 USE dyncor_c1d ! Coriolis term (c1d case) (cor_c1d routine) 36 ! 37 USE in_out_manager ! I/O manager 38 USE lib_mpp ! distributed memory computing library 39 USE lbclnk ! ocean lateral boundary condition (or mpp link) 38 40 USE timing ! Timing 39 USE lbclnk ! ocean lateral boundary condition (or mpp link)40 41 41 42 IMPLICIT NONE … … 81 82 ENDIF 82 83 ! 83 CALL dom_nam ! read namelist ( namrun, namdom , namcla)84 CALL dom_nam ! read namelist ( namrun, namdom ) 84 85 CALL dom_clo ! Closed seas and lake 85 86 CALL dom_hgr ! Horizontal mesh … … 88 89 IF( ln_sco ) CALL dom_stiff ! Maximum stiffness ratio/hydrostatic consistency 89 90 ! 90 ht_0(:,:) = 0. 0_wp! Reference ocean depth at T-points91 hu_0(:,:) = 0. 0_wp! Reference ocean depth at U-points92 hv_0(:,:) = 0. 0_wp! Reference ocean depth at V-points91 ht_0(:,:) = 0._wp ! Reference ocean depth at T-points 92 hu_0(:,:) = 0._wp ! Reference ocean depth at U-points 93 hv_0(:,:) = 0._wp ! Reference ocean depth at V-points 93 94 DO jk = 1, jpk 94 95 ht_0(:,:) = ht_0(:,:) + e3t_0(:,:,jk) * tmask(:,:,jk) … … 97 98 END DO 98 99 ! 99 IF( lk_vvl )CALL dom_vvl_init ! Vertical variable mesh100 IF( lk_vvl ) CALL dom_vvl_init ! Vertical variable mesh 100 101 ! 101 102 IF( lk_c1d ) CALL cor_c1d ! 1D configuration: Coriolis set at T-point … … 131 132 !! ** input : - namrun namelist 132 133 !! - namdom namelist 133 !! - namcla namelist134 134 !! - namnc4 namelist ! "key_netcdf4" only 135 135 !!---------------------------------------------------------------------- … … 146 146 & ppsur, ppa0, ppa1, ppkth, ppacr, ppdzmin, pphmax, ldbletanh, & 147 147 & ppa2, ppkth2, ppacr2 148 NAMELIST/namcla/ nn_cla149 148 #if defined key_netcdf4 150 149 NAMELIST/namnc4/ nn_nchunks_i, nn_nchunks_j, nn_nchunks_k, ln_nc4zip … … 155 154 REWIND( numnam_ref ) ! Namelist namrun in reference namelist : Parameters of the run 156 155 READ ( numnam_ref, namrun, IOSTAT = ios, ERR = 901) 157 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namrun in reference namelist', lwp )156 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namrun in reference namelist', lwp ) 158 157 159 158 REWIND( numnam_cfg ) ! Namelist namrun in configuration namelist : Parameters of the run 160 159 READ ( numnam_cfg, namrun, IOSTAT = ios, ERR = 902 ) 161 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namrun in configuration namelist', lwp )160 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namrun in configuration namelist', lwp ) 162 161 IF(lwm) WRITE ( numond, namrun ) 163 162 ! … … 251 250 904 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdom in configuration namelist', lwp ) 252 251 IF(lwm) WRITE ( numond, namdom ) 253 252 ! 254 253 IF(lwp) THEN 255 254 WRITE(numout,*) … … 293 292 WRITE(numout,*) ' ppacr2 = ', ppacr2 294 293 ENDIF 295 294 ! 296 295 ntopo = nn_bathy ! conversion DOCTOR names into model names (this should disappear soon) 297 296 e3zps_min = rn_e3zps_min … … 304 303 rdtmax = rn_rdtmin 305 304 rdth = rn_rdth 306 307 REWIND( numnam_ref ) ! Namelist namcla in reference namelist : Cross land advection308 READ ( numnam_ref, namcla, IOSTAT = ios, ERR = 905)309 905 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namcla in reference namelist', lwp )310 311 REWIND( numnam_cfg ) ! Namelist namcla in configuration namelist : Cross land advection312 READ ( numnam_cfg, namcla, IOSTAT = ios, ERR = 906 )313 906 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namcla in configuration namelist', lwp )314 IF(lwm) WRITE( numond, namcla )315 316 IF(lwp) THEN317 WRITE(numout,*)318 WRITE(numout,*) ' Namelist namcla'319 WRITE(numout,*) ' cross land advection nn_cla = ', nn_cla320 ENDIF321 IF ( nn_cla .EQ. 1 ) THEN322 IF ( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! ORCA R2323 CONTINUE324 ELSE325 CALL ctl_stop( 'STOP', 'Cross land advation iplemented only for ORCA2 configuration: cp_cfg = "orca" and jp_cfg = 2 ' )326 ENDIF327 ENDIF328 305 329 306 #if defined key_netcdf4 … … 409 386 END SUBROUTINE dom_ctl 410 387 388 411 389 SUBROUTINE dom_stiff 412 390 !!---------------------------------------------------------------------- … … 427 405 REAL(wp), DIMENSION(4) :: zr1 428 406 !!---------------------------------------------------------------------- 429 rx1(:,:) = 0. e0430 zrxmax = 0. e0431 zr1(:) = 0. e0432 407 rx1(:,:) = 0._wp 408 zrxmax = 0._wp 409 zr1(:) = 0._wp 410 ! 433 411 DO ji = 2, jpim1 434 412 DO jj = 2, jpjm1 … … 455 433 END DO 456 434 END DO 457 458 435 CALL lbc_lnk( rx1, 'T', 1. ) 459 460 zrxmax = MAXVAL( rx1)461 436 ! 437 zrxmax = MAXVAL( rx1 ) 438 ! 462 439 IF( lk_mpp ) CALL mpp_max( zrxmax ) ! max over the global domain 463 440 ! 464 441 IF(lwp) THEN 465 442 WRITE(numout,*) … … 467 444 WRITE(numout,*) '~~~~~~~~~' 468 445 ENDIF 469 446 ! 470 447 END SUBROUTINE dom_stiff 471 472 473 448 474 449 !!====================================================================== -
trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90
r5656 r5836 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 -
trunk/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90
r5552 r5836 17 17 !! - ! 2005-11 (V. Garnier) Surface pressure gradient organization 18 18 !! 3.2 ! 2009-07 (R. Benshila) Suppression of rigid-lid option 19 !! 3.6 ! 2015-05 (P. Mathiot) ISF: add wmask,wumask and wvmask 19 20 !!---------------------------------------------------------------------- 20 21 21 22 !!---------------------------------------------------------------------- 22 23 !! dom_msk : compute land/ocean mask 23 !! dom_msk_nsa : update land/ocean mask when no-slip accurate option is used.24 24 !!---------------------------------------------------------------------- 25 25 USE oce ! ocean dynamics and tracers … … 36 36 37 37 PUBLIC dom_msk ! routine called by inidom.F90 38 PUBLIC dom_msk_alloc ! routine called by nemogcm.F9039 38 40 39 ! !!* Namelist namlbc : lateral boundary condition * … … 42 41 LOGICAL, PUBLIC :: ln_vorlat ! consistency of vorticity boundary condition 43 42 ! with analytical eqs. 44 45 46 INTEGER, ALLOCATABLE, SAVE, DIMENSION(:,:) :: icoord ! Workspace for dom_msk_nsa()47 43 48 44 !! * Substitutions … … 54 50 !!---------------------------------------------------------------------- 55 51 CONTAINS 56 57 INTEGER FUNCTION dom_msk_alloc()58 !!---------------------------------------------------------------------59 !! *** FUNCTION dom_msk_alloc ***60 !!---------------------------------------------------------------------61 dom_msk_alloc = 062 #if defined key_noslip_accurate63 ALLOCATE(icoord(jpi*jpj*jpk,3), STAT=dom_msk_alloc)64 #endif65 IF( dom_msk_alloc /= 0 ) CALL ctl_warn('dom_msk_alloc: failed to allocate icoord array')66 !67 END FUNCTION dom_msk_alloc68 69 52 70 53 SUBROUTINE dom_msk … … 129 112 !! tmask_i : interior ocean mask 130 113 !!---------------------------------------------------------------------- 131 ! 132 INTEGER :: ji, jj, jk ! dummy loop indices 114 INTEGER :: ji, jj, jk ! dummy loop indices 133 115 INTEGER :: iif, iil, ii0, ii1, ii ! local integers 134 116 INTEGER :: ijf, ijl, ij0, ij1 ! - - … … 199 181 END DO 200 182 201 !!gm ????202 #if defined key_zdfkpp203 IF( cp_cfg == 'orca' ) THEN204 IF( jp_cfg == 2 ) THEN ! land point on Bab el Mandeb zonal section205 ij0 = 87 ; ij1 = 88206 ii0 = 160 ; ii1 = 161207 tmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 0._wp208 ELSE209 IF(lwp) WRITE(numout,*)210 IF(lwp) WRITE(numout,cform_war)211 IF(lwp) WRITE(numout,*)212 IF(lwp) WRITE(numout,*)' A mask must be applied on Bab el Mandeb strait'213 IF(lwp) WRITE(numout,*)' in case of ORCAs configurations'214 IF(lwp) WRITE(numout,*)' This is a problem which is not yet solved'215 IF(lwp) WRITE(numout,*)216 ENDIF217 ENDIF218 #endif219 !!gm end220 221 183 ! Interior domain mask (used for global sum) 222 184 ! -------------------- … … 284 246 ! 3. Ocean/land mask at wu-, wv- and w points 285 247 !---------------------------------------------- 286 wmask (:,:,1) = tmask(:,:,1) ! ????????287 wumask(:,:,1) = umask(:,:,1) ! ????????288 wvmask(:,:,1) = vmask(:,:,1) ! ????????289 DO jk =2,jpk290 wmask (:,:,jk) =tmask(:,:,jk) * tmask(:,:,jk-1)291 wumask(:,:,jk) =umask(:,:,jk) * umask(:,:,jk-1)292 wvmask(:,:,jk) =vmask(:,:,jk) * vmask(:,:,jk-1)248 wmask (:,:,1) = tmask(:,:,1) ! surface 249 wumask(:,:,1) = umask(:,:,1) 250 wvmask(:,:,1) = vmask(:,:,1) 251 DO jk = 2, jpk ! interior values 252 wmask (:,:,jk) = tmask(:,:,jk) * tmask(:,:,jk-1) 253 wumask(:,:,jk) = umask(:,:,jk) * umask(:,:,jk-1) 254 wvmask(:,:,jk) = vmask(:,:,jk) * vmask(:,:,jk-1) 293 255 END DO 294 256 … … 339 301 ENDIF 340 302 341 342 ! mask for second order calculation of vorticity343 ! ----------------------------------------------344 CALL dom_msk_nsa345 346 347 303 ! Lateral boundary conditions on velocity (modify fmask) 348 304 ! --------------------------------------- … … 377 333 IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! ORCA_R2 configuration 378 334 ! ! Increased lateral friction near of some straits 379 IF( nn_cla == 0 ) THEN 380 ! ! Gibraltar strait : partial slip (fmask=0.5) 381 ij0 = 101 ; ij1 = 101 382 ii0 = 139 ; ii1 = 140 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 0.5_wp 383 ij0 = 102 ; ij1 = 102 384 ii0 = 139 ; ii1 = 140 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 0.5_wp 385 ! 386 ! ! Bab el Mandeb : partial slip (fmask=1) 387 ij0 = 87 ; ij1 = 88 388 ii0 = 160 ; ii1 = 160 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 1._wp 389 ij0 = 88 ; ij1 = 88 390 ii0 = 159 ; ii1 = 159 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 1._wp 391 ! 392 ENDIF 335 ! ! Gibraltar strait : partial slip (fmask=0.5) 336 ij0 = 101 ; ij1 = 101 337 ii0 = 139 ; ii1 = 140 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 0.5_wp 338 ij0 = 102 ; ij1 = 102 339 ii0 = 139 ; ii1 = 140 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 0.5_wp 340 ! 341 ! ! Bab el Mandeb : partial slip (fmask=1) 342 ij0 = 87 ; ij1 = 88 343 ii0 = 160 ; ii1 = 160 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 1._wp 344 ij0 = 88 ; ij1 = 88 345 ii0 = 159 ; ii1 = 159 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 1._wp 346 ! 393 347 ! ! Danish straits : strong slip (fmask > 2) 394 348 ! We keep this as an example but it is instable in this case … … 500 454 ! 501 455 END SUBROUTINE dom_msk 502 503 #if defined key_noslip_accurate504 !!----------------------------------------------------------------------505 !! 'key_noslip_accurate' : accurate no-slip boundary condition506 !!----------------------------------------------------------------------507 508 SUBROUTINE dom_msk_nsa509 !!---------------------------------------------------------------------510 !! *** ROUTINE dom_msk_nsa ***511 !!512 !! ** Purpose :513 !!514 !! ** Method :515 !!516 !! ** Action :517 !!----------------------------------------------------------------------518 INTEGER :: ji, jj, jk, jl ! dummy loop indices519 INTEGER :: ine, inw, ins, inn, itest, ierror, iind, ijnd520 REAL(wp) :: zaa521 !!---------------------------------------------------------------------522 !523 IF( nn_timing == 1 ) CALL timing_start('dom_msk_nsa')524 !525 IF(lwp) WRITE(numout,*)526 IF(lwp) WRITE(numout,*) 'dom_msk_nsa : noslip accurate boundary condition'527 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ using Schchepetkin and O Brian scheme'528 IF( lk_mpp ) CALL ctl_stop( ' mpp version is not yet implemented' )529 530 ! mask for second order calculation of vorticity531 ! ----------------------------------------------532 ! noslip boundary condition: fmask=1 at convex corner, store533 ! index of straight coast meshes ( 'west', refering to a coast,534 ! means west of the ocean, aso)535 536 DO jk = 1, jpk537 DO jl = 1, 4538 npcoa(jl,jk) = 0539 DO ji = 1, 2*(jpi+jpj)540 nicoa(ji,jl,jk) = 0541 njcoa(ji,jl,jk) = 0542 END DO543 END DO544 END DO545 546 IF( jperio == 2 ) THEN547 WRITE(numout,*) ' '548 WRITE(numout,*) ' symetric boundary conditions need special'549 WRITE(numout,*) ' treatment not implemented. we stop.'550 STOP551 ENDIF552 553 ! convex corners554 555 DO jk = 1, jpkm1556 DO jj = 1, jpjm1557 DO ji = 1, jpim1558 zaa = tmask(ji ,jj,jk) + tmask(ji ,jj+1,jk) &559 &+ tmask(ji+1,jj,jk) + tmask(ji+1,jj+1,jk)560 IF( ABS(zaa-3._wp) <= 0.1_wp ) fmask(ji,jj,jk) = 1._wp561 END DO562 END DO563 END DO564 565 ! north-south straight coast566 567 DO jk = 1, jpkm1568 inw = 0569 ine = 0570 DO jj = 2, jpjm1571 DO ji = 2, jpim1572 zaa = tmask(ji+1,jj,jk) + tmask(ji+1,jj+1,jk)573 IF( ABS(zaa-2._wp) <= 0.1_wp .AND. fmask(ji,jj,jk) == 0._wp ) THEN574 inw = inw + 1575 nicoa(inw,1,jk) = ji576 njcoa(inw,1,jk) = jj577 IF( nprint == 1 ) WRITE(numout,*) ' west : ', jk, inw, ji, jj578 ENDIF579 zaa = tmask(ji,jj,jk) + tmask(ji,jj+1,jk)580 IF( ABS(zaa-2._wp) <= 0.1_wp .AND. fmask(ji,jj,jk) == 0._wp ) THEN581 ine = ine + 1582 nicoa(ine,2,jk) = ji583 njcoa(ine,2,jk) = jj584 IF( nprint == 1 ) WRITE(numout,*) ' east : ', jk, ine, ji, jj585 ENDIF586 END DO587 END DO588 npcoa(1,jk) = inw589 npcoa(2,jk) = ine590 END DO591 592 ! west-east straight coast593 594 DO jk = 1, jpkm1595 ins = 0596 inn = 0597 DO jj = 2, jpjm1598 DO ji =2, jpim1599 zaa = tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk)600 IF( ABS(zaa-2._wp) <= 0.1_wp .AND. fmask(ji,jj,jk) == 0._wp ) THEN601 ins = ins + 1602 nicoa(ins,3,jk) = ji603 njcoa(ins,3,jk) = jj604 IF( nprint == 1 ) WRITE(numout,*) ' south : ', jk, ins, ji, jj605 ENDIF606 zaa = tmask(ji+1,jj,jk) + tmask(ji,jj,jk)607 IF( ABS(zaa-2._wp) <= 0.1_wp .AND. fmask(ji,jj,jk) == 0._wp ) THEN608 inn = inn + 1609 nicoa(inn,4,jk) = ji610 njcoa(inn,4,jk) = jj611 IF( nprint == 1 ) WRITE(numout,*) ' north : ', jk, inn, ji, jj612 ENDIF613 END DO614 END DO615 npcoa(3,jk) = ins616 npcoa(4,jk) = inn617 END DO618 619 itest = 2 * ( jpi + jpj )620 DO jk = 1, jpk621 IF( npcoa(1,jk) > itest .OR. npcoa(2,jk) > itest .OR. &622 npcoa(3,jk) > itest .OR. npcoa(4,jk) > itest ) THEN623 624 WRITE(ctmp1,*) ' level jk = ',jk625 WRITE(ctmp2,*) ' straight coast index arraies are too small.:'626 WRITE(ctmp3,*) ' npe, npw, nps, npn = ', npcoa(1,jk), npcoa(2,jk), &627 & npcoa(3,jk), npcoa(4,jk)628 WRITE(ctmp4,*) ' 2*(jpi+jpj) = ',itest,'. we stop.'629 CALL ctl_stop( ctmp1, ctmp2, ctmp3, ctmp4 )630 ENDIF631 END DO632 633 ierror = 0634 iind = 0635 ijnd = 0636 IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) iind = 2637 IF( nperio == 3 .OR. nperio == 4 .OR. nperio == 5 .OR. nperio == 6 ) ijnd = 2638 DO jk = 1, jpk639 DO jl = 1, npcoa(1,jk)640 IF( nicoa(jl,1,jk)+3 > jpi+iind ) THEN641 ierror = ierror+1642 icoord(ierror,1) = nicoa(jl,1,jk)643 icoord(ierror,2) = njcoa(jl,1,jk)644 icoord(ierror,3) = jk645 ENDIF646 END DO647 DO jl = 1, npcoa(2,jk)648 IF(nicoa(jl,2,jk)-2 < 1-iind ) THEN649 ierror = ierror + 1650 icoord(ierror,1) = nicoa(jl,2,jk)651 icoord(ierror,2) = njcoa(jl,2,jk)652 icoord(ierror,3) = jk653 ENDIF654 END DO655 DO jl = 1, npcoa(3,jk)656 IF( njcoa(jl,3,jk)+3 > jpj+ijnd ) THEN657 ierror = ierror + 1658 icoord(ierror,1) = nicoa(jl,3,jk)659 icoord(ierror,2) = njcoa(jl,3,jk)660 icoord(ierror,3) = jk661 ENDIF662 END DO663 DO jl = 1, npcoa(4,jk)664 IF( njcoa(jl,4,jk)-2 < 1) THEN665 ierror=ierror + 1666 icoord(ierror,1) = nicoa(jl,4,jk)667 icoord(ierror,2) = njcoa(jl,4,jk)668 icoord(ierror,3) = jk669 ENDIF670 END DO671 END DO672 673 IF( ierror > 0 ) THEN674 IF(lwp) WRITE(numout,*)675 IF(lwp) WRITE(numout,*) ' Problem on lateral conditions'676 IF(lwp) WRITE(numout,*) ' Bad marking off at points:'677 DO jl = 1, ierror678 IF(lwp) WRITE(numout,*) 'Level:',icoord(jl,3), &679 & ' Point(',icoord(jl,1),',',icoord(jl,2),')'680 END DO681 CALL ctl_stop( 'We stop...' )682 ENDIF683 !684 IF( nn_timing == 1 ) CALL timing_stop('dom_msk_nsa')685 !686 END SUBROUTINE dom_msk_nsa687 688 #else689 !!----------------------------------------------------------------------690 !! Default option : Empty routine691 !!----------------------------------------------------------------------692 SUBROUTINE dom_msk_nsa693 END SUBROUTINE dom_msk_nsa694 #endif695 456 696 457 !!====================================================================== -
trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domngb.F90
r3294 r5836 37 37 !! -> not good if located at too high latitude... 38 38 !!---------------------------------------------------------------------- 39 !40 39 REAL(wp) , INTENT(in ) :: plon, plat ! longitude,latitude of the point 41 40 INTEGER , INTENT( out) :: kii, kjj ! i-,j-index of the closes grid point … … 49 48 IF( nn_timing == 1 ) CALL timing_start('dom_ngb') 50 49 ! 51 CALL wrk_alloc( jpi, jpj,zglam, zgphi, zmask, zdist )50 CALL wrk_alloc( jpi,jpj, zglam, zgphi, zmask, zdist ) 52 51 ! 53 52 zmask(:,:) = 0._wp … … 76 75 ENDIF 77 76 ! 78 CALL wrk_dealloc( jpi, jpj,zglam, zgphi, zmask, zdist )77 CALL wrk_dealloc( jpi,jpj, zglam, zgphi, zmask, zdist ) 79 78 ! 80 79 IF( nn_timing == 1 ) CALL timing_stop('dom_ngb') -
trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90
r5506 r5836 10 10 !! 3.6 ! 2014-11 (P. Mathiot) add ice shelf capability 11 11 !!---------------------------------------------------------------------- 12 !! 'key_vvl' variable volume 13 !!---------------------------------------------------------------------- 12 14 13 !!---------------------------------------------------------------------- 15 14 !! dom_vvl_init : define initial vertical scale factors, depths and column thickness … … 19 18 !! dom_vvl_rst : read/write restart file 20 19 !! dom_vvl_ctl : Check the vvl options 21 !! dom_vvl_orca_fix : Recompute some area-weighted interpolations of vertical scale factors22 !! : to account for manual changes to e[1,2][u,v] in some Straits23 20 !!---------------------------------------------------------------------- 24 !! * Modules used25 21 USE oce ! ocean dynamics and tracers 26 22 USE dom_oce ! ocean space and time domain … … 37 33 PRIVATE 38 34 39 !! * Routine accessibility40 35 PUBLIC dom_vvl_init ! called by domain.F90 41 36 PUBLIC dom_vvl_sf_nxt ! called by step.F90 42 37 PUBLIC dom_vvl_sf_swp ! called by step.F90 43 38 PUBLIC dom_vvl_interpol ! called by dynnxt.F90 44 PRIVATE dom_vvl_orca_fix ! called by dom_vvl_interpol 45 46 !!* Namelist nam_vvl 47 LOGICAL , PUBLIC :: ln_vvl_zstar = .FALSE. ! zstar vertical coordinate 48 LOGICAL , PUBLIC :: ln_vvl_ztilde = .FALSE. ! ztilde vertical coordinate 49 LOGICAL , PUBLIC :: ln_vvl_layer = .FALSE. ! level vertical coordinate 50 LOGICAL , PUBLIC :: ln_vvl_ztilde_as_zstar = .FALSE. ! ztilde vertical coordinate 51 LOGICAL , PUBLIC :: ln_vvl_zstar_at_eqtor = .FALSE. ! ztilde vertical coordinate 52 LOGICAL , PUBLIC :: ln_vvl_kepe = .FALSE. ! kinetic/potential energy transfer 53 ! ! conservation: not used yet 54 REAL(wp) :: rn_ahe3 ! thickness diffusion coefficient 55 REAL(wp) :: rn_rst_e3t ! ztilde to zstar restoration timescale [days] 56 REAL(wp) :: rn_lf_cutoff ! cutoff frequency for low-pass filter [days] 57 REAL(wp) :: rn_zdef_max ! maximum fractional e3t deformation 58 LOGICAL , PUBLIC :: ln_vvl_dbg = .FALSE. ! debug control prints 59 60 !! * Module variables 61 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: un_td, vn_td ! thickness diffusion transport 62 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hdiv_lf ! low frequency part of hz divergence 63 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_b, tilde_e3t_n ! baroclinic scale factors 64 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_a, dtilde_e3t_a ! baroclinic scale factors 65 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: frq_rst_e3t ! retoring period for scale factors 66 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: frq_rst_hdv ! retoring period for low freq. divergence 39 40 ! !!* Namelist nam_vvl 41 LOGICAL , PUBLIC :: ln_vvl_zstar = .FALSE. ! zstar vertical coordinate 42 LOGICAL , PUBLIC :: ln_vvl_ztilde = .FALSE. ! ztilde vertical coordinate 43 LOGICAL , PUBLIC :: ln_vvl_layer = .FALSE. ! level vertical coordinate 44 LOGICAL , PUBLIC :: ln_vvl_ztilde_as_zstar = .FALSE. ! ztilde vertical coordinate 45 LOGICAL , PUBLIC :: ln_vvl_zstar_at_eqtor = .FALSE. ! ztilde vertical coordinate 46 LOGICAL , PUBLIC :: ln_vvl_kepe = .FALSE. ! kinetic/potential energy transfer 47 ! ! conservation: not used yet 48 REAL(wp) :: rn_ahe3 ! thickness diffusion coefficient 49 REAL(wp) :: rn_rst_e3t ! ztilde to zstar restoration timescale [days] 50 REAL(wp) :: rn_lf_cutoff ! cutoff frequency for low-pass filter [days] 51 REAL(wp) :: rn_zdef_max ! maximum fractional e3t deformation 52 LOGICAL , PUBLIC :: ln_vvl_dbg = .FALSE. ! debug control prints 53 54 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: un_td, vn_td ! thickness diffusion transport 55 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hdiv_lf ! low frequency part of hz divergence 56 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_b, tilde_e3t_n ! baroclinic scale factors 57 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_a, dtilde_e3t_a ! baroclinic scale factors 58 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: frq_rst_e3t ! retoring period for scale factors 59 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: frq_rst_hdv ! retoring period for low freq. divergence 67 60 68 61 !! * Substitutions … … 372 365 DO jj = 1, jpjm1 373 366 DO ji = 1, fs_jpim1 ! vector opt. 374 un_td(ji,jj,jk) = rn_ahe3 * umask(ji,jj,jk) * re2u_e1u(ji,jj)&375 &* ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj ,jk) )376 vn_td(ji,jj,jk) = rn_ahe3 * vmask(ji,jj,jk) * re1v_e2v(ji,jj)&377 &* ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji ,jj+1,jk) )367 un_td(ji,jj,jk) = rn_ahe3 * umask(ji,jj,jk) * e2_e1u(ji,jj) & 368 & * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj ,jk) ) 369 vn_td(ji,jj,jk) = rn_ahe3 * vmask(ji,jj,jk) * e1_e2v(ji,jj) & 370 & * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji ,jj+1,jk) ) 378 371 zwu(ji,jj) = zwu(ji,jj) + un_td(ji,jj,jk) 379 372 zwv(ji,jj) = zwv(ji,jj) + vn_td(ji,jj,jk) … … 394 387 tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) + ( un_td(ji-1,jj ,jk) - un_td(ji,jj,jk) & 395 388 & + vn_td(ji ,jj-1,jk) - vn_td(ji,jj,jk) & 396 & ) * r1_e1 2t(ji,jj)389 & ) * r1_e1e2t(ji,jj) 397 390 END DO 398 391 END DO … … 695 688 !! - vertical interpolation: simple averaging 696 689 !!---------------------------------------------------------------------- 697 !! * Arguments 698 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in ) :: pe3_in ! input e3 to be interpolated 699 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) :: pe3_out ! output interpolated e3 700 CHARACTER(LEN=*), INTENT( in ) :: pout ! grid point of out scale factors 701 ! ! = 'U', 'V', 'W, 'F', 'UW' or 'VW' 702 !! * Local declarations 690 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pe3_in ! input e3 to be interpolated 691 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pe3_out ! output interpolated e3 692 CHARACTER(LEN=*) , INTENT(in ) :: pout ! grid point of out scale factors 693 ! ! = 'U', 'V', 'W, 'F', 'UW' or 'VW' 694 ! 703 695 INTEGER :: ji, jj, jk ! dummy loop indices 704 LOGICAL :: l_is_orca ! local logical705 ! !----------------------------------------------------------------------696 !!---------------------------------------------------------------------- 697 ! 706 698 IF( nn_timing == 1 ) CALL timing_start('dom_vvl_interpol') 707 ! 708 l_is_orca = .FALSE. 709 IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) l_is_orca = .TRUE. ! ORCA R2 configuration - will need to correct some locations 710 711 SELECT CASE ( pout ) 712 ! ! ------------------------------------- ! 713 CASE( 'U' ) ! interpolation from T-point to U-point ! 714 ! ! ------------------------------------- ! 715 ! horizontal surface weighted interpolation 699 ! 700 SELECT CASE ( pout ) !== type of interpolation ==! 701 ! 702 CASE( 'U' ) !* from T- to U-point : hor. surface weighted mean 716 703 DO jk = 1, jpk 717 704 DO jj = 1, jpjm1 718 705 DO ji = 1, fs_jpim1 ! vector opt. 719 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * r1_e1 2u(ji,jj) &720 & * ( e1 2t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) &721 & + e1 2t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) )706 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * r1_e1e2u(ji,jj) & 707 & * ( e1e2t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) & 708 & + e1e2t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) ) 722 709 END DO 723 710 END DO 724 711 END DO 725 !726 IF( l_is_orca ) CALL dom_vvl_orca_fix( pe3_in, pe3_out, pout )727 ! boundary conditions728 712 CALL lbc_lnk( pe3_out(:,:,:), 'U', 1._wp ) 729 713 pe3_out(:,:,:) = pe3_out(:,:,:) + e3u_0(:,:,:) 730 ! ! ------------------------------------- ! 731 CASE( 'V' ) ! interpolation from T-point to V-point ! 732 ! ! ------------------------------------- ! 733 ! horizontal surface weighted interpolation 714 ! 715 CASE( 'V' ) !* from T- to V-point : hor. surface weighted mean 734 716 DO jk = 1, jpk 735 717 DO jj = 1, jpjm1 736 718 DO ji = 1, fs_jpim1 ! vector opt. 737 pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) * r1_e1 2v(ji,jj) &738 & * ( e1 2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) &739 & + e1 2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) )719 pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) * r1_e1e2v(ji,jj) & 720 & * ( e1e2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) & 721 & + e1e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) ) 740 722 END DO 741 723 END DO 742 724 END DO 743 !744 IF( l_is_orca ) CALL dom_vvl_orca_fix( pe3_in, pe3_out, pout )745 ! boundary conditions746 725 CALL lbc_lnk( pe3_out(:,:,:), 'V', 1._wp ) 747 726 pe3_out(:,:,:) = pe3_out(:,:,:) + e3v_0(:,:,:) 748 ! ! ------------------------------------- ! 749 CASE( 'F' ) ! interpolation from U-point to F-point ! 750 ! ! ------------------------------------- ! 751 ! horizontal surface weighted interpolation 727 ! 728 CASE( 'F' ) !* from U-point to F-point : hor. surface weighted mean 752 729 DO jk = 1, jpk 753 730 DO jj = 1, jpjm1 754 731 DO ji = 1, fs_jpim1 ! vector opt. 755 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) * r1_e1 2f(ji,jj) &756 & * ( e1 2u(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3u_0(ji,jj ,jk) ) &757 & + e1 2u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) )732 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) * r1_e1e2f(ji,jj) & 733 & * ( e1e2u(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3u_0(ji,jj ,jk) ) & 734 & + e1e2u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) ) 758 735 END DO 759 736 END DO 760 737 END DO 761 !762 IF( l_is_orca ) CALL dom_vvl_orca_fix( pe3_in, pe3_out, pout )763 ! boundary conditions764 738 CALL lbc_lnk( pe3_out(:,:,:), 'F', 1._wp ) 765 739 pe3_out(:,:,:) = pe3_out(:,:,:) + e3f_0(:,:,:) 766 ! ! ------------------------------------- ! 767 CASE( 'W' ) ! interpolation from T-point to W-point ! 768 ! ! ------------------------------------- ! 769 ! vertical simple interpolation 740 ! 741 CASE( 'W' ) !* from T- to W-point : vertical simple mean 742 ! 770 743 pe3_out(:,:,1) = e3w_0(:,:,1) + pe3_in(:,:,1) - e3t_0(:,:,1) 771 ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing 744 ! - ML - The use of mask in this formulea enables the special treatment of the last w-point without indirect adressing 745 !!gm BUG? use here wmask in case of ISF ? to be checked 772 746 DO jk = 2, jpk 773 747 pe3_out(:,:,jk) = e3w_0(:,:,jk) + ( 1.0_wp - 0.5_wp * tmask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3t_0(:,:,jk-1) ) & 774 748 & + 0.5_wp * tmask(:,:,jk) * ( pe3_in(:,:,jk ) - e3t_0(:,:,jk ) ) 775 749 END DO 776 ! ! -------------------------------------- ! 777 CASE( 'UW' ) ! interpolation from U-point to UW-point ! 778 ! ! -------------------------------------- ! 779 ! vertical simple interpolation 750 ! 751 CASE( 'UW' ) !* from U- to UW-point : vertical simple mean 752 ! 780 753 pe3_out(:,:,1) = e3uw_0(:,:,1) + pe3_in(:,:,1) - e3u_0(:,:,1) 781 754 ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing 755 !!gm BUG? use here wumask in case of ISF ? to be checked 782 756 DO jk = 2, jpk 783 757 pe3_out(:,:,jk) = e3uw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * umask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3u_0(:,:,jk-1) ) & 784 758 & + 0.5_wp * umask(:,:,jk) * ( pe3_in(:,:,jk ) - e3u_0(:,:,jk ) ) 785 759 END DO 786 ! ! -------------------------------------- ! 787 CASE( 'VW' ) ! interpolation from V-point to VW-point ! 788 ! ! -------------------------------------- ! 789 ! vertical simple interpolation 760 ! 761 CASE( 'VW' ) !* from V- to VW-point : vertical simple mean 762 ! 790 763 pe3_out(:,:,1) = e3vw_0(:,:,1) + pe3_in(:,:,1) - e3v_0(:,:,1) 791 764 ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing 765 !!gm BUG? use here wvmask in case of ISF ? to be checked 792 766 DO jk = 2, jpk 793 767 pe3_out(:,:,jk) = e3vw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * vmask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3v_0(:,:,jk-1) ) & … … 796 770 END SELECT 797 771 ! 798 799 772 IF( nn_timing == 1 ) CALL timing_stop('dom_vvl_interpol') 800 773 ! 801 774 END SUBROUTINE dom_vvl_interpol 775 802 776 803 777 SUBROUTINE dom_vvl_rst( kt, cdrw ) … … 813 787 !! they are set to 0. 814 788 !!---------------------------------------------------------------------- 815 !! * Arguments816 789 INTEGER , INTENT(in) :: kt ! ocean time-step 817 790 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag 818 ! ! * Local declarations791 ! 819 792 INTEGER :: jk 820 793 INTEGER :: id1, id2, id3, id4, id5 ! local integers … … 911 884 END IF 912 885 ENDIF 913 886 ! 914 887 ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN ! Create restart file 915 888 ! ! =================== … … 931 904 CALL iom_rstput( kt, nitrst, numrow, 'hdiv_lf', hdiv_lf(:,:,:) ) 932 905 ENDIF 933 934 ENDIF 906 ! 907 ENDIF 908 ! 935 909 IF( nn_timing == 1 ) CALL timing_stop('dom_vvl_rst') 936 910 ! 937 911 END SUBROUTINE dom_vvl_rst 938 912 … … 945 919 !! for vertical coordinate 946 920 !!---------------------------------------------------------------------- 947 INTEGER :: ioptio 948 INTEGER :: ios 949 921 INTEGER :: ioptio, ios 922 !! 950 923 NAMELIST/nam_vvl/ ln_vvl_zstar, ln_vvl_ztilde, ln_vvl_layer, ln_vvl_ztilde_as_zstar, & 951 &ln_vvl_zstar_at_eqtor , rn_ahe3 , rn_rst_e3t , &952 &rn_lf_cutoff , rn_zdef_max , ln_vvl_dbg ! not yet implemented: ln_vvl_kepe924 & ln_vvl_zstar_at_eqtor , rn_ahe3 , rn_rst_e3t , & 925 & rn_lf_cutoff , rn_zdef_max , ln_vvl_dbg ! not yet implemented: ln_vvl_kepe 953 926 !!---------------------------------------------------------------------- 954 927 ! 955 928 REWIND( numnam_ref ) ! Namelist nam_vvl in reference namelist : 956 929 READ ( numnam_ref, nam_vvl, IOSTAT = ios, ERR = 901) 957 930 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_vvl in reference namelist', lwp ) 958 931 ! 959 932 REWIND( numnam_cfg ) ! Namelist nam_vvl in configuration namelist : Parameters of the run 960 933 READ ( numnam_cfg, nam_vvl, IOSTAT = ios, ERR = 902 ) 961 934 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_vvl in configuration namelist', lwp ) 962 935 IF(lwm) WRITE ( numond, nam_vvl ) 963 936 ! 964 937 IF(lwp) THEN ! Namelist print 965 938 WRITE(numout,*) … … 994 967 WRITE(numout,*) ' ln_vvl_dbg = ', ln_vvl_dbg 995 968 ENDIF 996 969 ! 997 970 ioptio = 0 ! Parameter control 998 IF( ln_vvl_ztilde_as_zstar ) ln_vvl_ztilde = .true.999 IF( ln_vvl_zstar ) 1000 IF( ln_vvl_ztilde ) 1001 IF( ln_vvl_layer ) 1002 971 IF( ln_vvl_ztilde_as_zstar ) ln_vvl_ztilde = .true. 972 IF( ln_vvl_zstar ) ioptio = ioptio + 1 973 IF( ln_vvl_ztilde ) ioptio = ioptio + 1 974 IF( ln_vvl_layer ) ioptio = ioptio + 1 975 ! 1003 976 IF( ioptio /= 1 ) CALL ctl_stop( 'Choose ONE vertical coordinate in namelist nam_vvl' ) 1004 977 IF( .NOT. ln_vvl_zstar .AND. nn_isf .NE. 0) CALL ctl_stop( 'Only vvl_zstar has been tested with ice shelf cavity' ) 1005 978 ! 1006 979 IF(lwp) THEN ! Print the choice 1007 980 WRITE(numout,*) … … 1014 987 ! IF( .NOT. ln_vvl_kepe ) WRITE(numout,*) ' kinetic to potential energy transfer : option not used' 1015 988 ENDIF 1016 989 ! 1017 990 #if defined key_agrif 1018 991 IF (.NOT.Agrif_Root()) CALL ctl_stop( 'AGRIF not implemented with non-linear free surface (key_vvl)' ) 1019 992 #endif 1020 993 ! 1021 994 END SUBROUTINE dom_vvl_ctl 1022 1023 SUBROUTINE dom_vvl_orca_fix( pe3_in, pe3_out, pout )1024 !!---------------------------------------------------------------------1025 !! *** ROUTINE dom_vvl_orca_fix ***1026 !!1027 !! ** Purpose : Correct surface weighted, horizontally interpolated,1028 !! scale factors at locations that have been individually1029 !! modified in domhgr. Such modifications break the1030 !! relationship between e12t and e1u*e2u etc.1031 !! Recompute some scale factors ignoring the modified metric.1032 !!----------------------------------------------------------------------1033 !! * Arguments1034 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in ) :: pe3_in ! input e3 to be interpolated1035 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) :: pe3_out ! output interpolated e31036 CHARACTER(LEN=*), INTENT( in ) :: pout ! grid point of out scale factors1037 ! ! = 'U', 'V', 'W, 'F', 'UW' or 'VW'1038 !! * Local declarations1039 INTEGER :: ji, jj, jk ! dummy loop indices1040 INTEGER :: ij0, ij1, ii0, ii1 ! dummy loop indices1041 INTEGER :: isrow ! index for ORCA1 starting row1042 !! acc1043 !! Hmm with the time splitting these "fixes" seem to do more harm than good. Temporarily disabled for1044 !! the ORCA2 tests (by changing jp_cfg test from 2 to 3) pending further investigations1045 !!1046 ! ! =====================1047 IF( cp_cfg == "orca" .AND. jp_cfg == 3 ) THEN ! ORCA R2 configuration1048 ! ! =====================1049 !! acc1050 IF( nn_cla == 0 ) THEN1051 !1052 ii0 = 139 ; ii1 = 140 ! Gibraltar Strait (e2u was modified)1053 ij0 = 102 ; ij1 = 1021054 DO jk = 1, jpkm11055 DO jj = mj0(ij0), mj1(ij1)1056 DO ji = mi0(ii0), mi1(ii1)1057 SELECT CASE ( pout )1058 CASE( 'U' )1059 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) &1060 & * ( e1t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) &1061 & + e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &1062 & ) / e1u(ji,jj) + e3u_0(ji,jj,jk)1063 CASE( 'F' )1064 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) &1065 & * ( e1u(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3u_0(ji ,jj,jk) ) &1066 & + e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &1067 & ) / e1f(ji,jj) + e3f_0(ji,jj,jk)1068 END SELECT1069 END DO1070 END DO1071 END DO1072 !1073 ii0 = 160 ; ii1 = 160 ! Bab el Mandeb (e2u and e1v were modified)1074 ij0 = 88 ; ij1 = 881075 DO jk = 1, jpkm11076 DO jj = mj0(ij0), mj1(ij1)1077 DO ji = mi0(ii0), mi1(ii1)1078 SELECT CASE ( pout )1079 CASE( 'U' )1080 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) &1081 & * ( e1t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) &1082 & + e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &1083 & ) / e1u(ji,jj) + e3u_0(ji,jj,jk)1084 CASE( 'V' )1085 pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) &1086 & * ( e2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) &1087 & + e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &1088 & ) / e2v(ji,jj) + e3v_0(ji,jj,jk)1089 CASE( 'F' )1090 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) &1091 & * ( e1u(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3u_0(ji ,jj,jk) ) &1092 & + e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &1093 & ) / e1f(ji,jj) + e3f_0(ji,jj,jk)1094 END SELECT1095 END DO1096 END DO1097 END DO1098 ENDIF1099 1100 ii0 = 145 ; ii1 = 146 ! Danish Straits (e2u was modified)1101 ij0 = 116 ; ij1 = 1161102 DO jk = 1, jpkm11103 DO jj = mj0(ij0), mj1(ij1)1104 DO ji = mi0(ii0), mi1(ii1)1105 SELECT CASE ( pout )1106 CASE( 'U' )1107 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) &1108 & * ( e1t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) &1109 & + e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &1110 & ) / e1u(ji,jj) + e3u_0(ji,jj,jk)1111 CASE( 'F' )1112 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) &1113 & * ( e1u(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3u_0(ji ,jj,jk) ) &1114 & + e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &1115 & ) / e1f(ji,jj) + e3f_0(ji,jj,jk)1116 END SELECT1117 END DO1118 END DO1119 END DO1120 ENDIF1121 !1122 ! ! =====================1123 IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN ! ORCA R1 configuration1124 ! ! =====================1125 ! This dirty section will be suppressed by simplification process:1126 ! all this will come back in input files1127 ! Currently these hard-wired indices relate to configuration with1128 ! extend grid (jpjglo=332)1129 ! which had a grid-size of 362x292.1130 isrow = 332 - jpjglo1131 !1132 ii0 = 282 ; ii1 = 283 ! Gibraltar Strait (e2u was modified)1133 ij0 = 241 - isrow ; ij1 = 241 - isrow1134 DO jk = 1, jpkm11135 DO jj = mj0(ij0), mj1(ij1)1136 DO ji = mi0(ii0), mi1(ii1)1137 SELECT CASE ( pout )1138 CASE( 'U' )1139 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) &1140 & * ( e1t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) &1141 & + e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &1142 & ) / e1u(ji,jj) + e3u_0(ji,jj,jk)1143 CASE( 'F' )1144 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) &1145 & * ( e1u(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3u_0(ji ,jj,jk) ) &1146 & + e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &1147 & ) / e1f(ji,jj) + e3f_0(ji,jj,jk)1148 END SELECT1149 END DO1150 END DO1151 END DO1152 !1153 ii0 = 314 ; ii1 = 315 ! Bhosporus Strait (e2u was modified)1154 ij0 = 248 - isrow ; ij1 = 248 - isrow1155 DO jk = 1, jpkm11156 DO jj = mj0(ij0), mj1(ij1)1157 DO ji = mi0(ii0), mi1(ii1)1158 SELECT CASE ( pout )1159 CASE( 'U' )1160 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) &1161 & * ( e1t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) &1162 & + e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &1163 & ) / e1u(ji,jj) + e3u_0(ji,jj,jk)1164 CASE( 'F' )1165 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) &1166 & * ( e1u(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3u_0(ji ,jj,jk) ) &1167 & + e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &1168 & ) / e1f(ji,jj) + e3f_0(ji,jj,jk)1169 END SELECT1170 END DO1171 END DO1172 END DO1173 !1174 ii0 = 44 ; ii1 = 44 ! Lombok Strait (e1v was modified)1175 ij0 = 164 - isrow ; ij1 = 165 - isrow1176 DO jk = 1, jpkm11177 DO jj = mj0(ij0), mj1(ij1)1178 DO ji = mi0(ii0), mi1(ii1)1179 SELECT CASE ( pout )1180 CASE( 'V' )1181 pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) &1182 & * ( e2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) &1183 & + e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &1184 & ) / e2v(ji,jj) + e3v_0(ji,jj,jk)1185 END SELECT1186 END DO1187 END DO1188 END DO1189 !1190 ii0 = 48 ; ii1 = 48 ! Sumba Strait (e1v was modified) [closed from bathy_11 on]1191 ij0 = 164 - isrow ; ij1 = 165 - isrow1192 DO jk = 1, jpkm11193 DO jj = mj0(ij0), mj1(ij1)1194 DO ji = mi0(ii0), mi1(ii1)1195 SELECT CASE ( pout )1196 CASE( 'V' )1197 pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) &1198 & * ( e2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) &1199 & + e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &1200 & ) / e2v(ji,jj) + e3v_0(ji,jj,jk)1201 END SELECT1202 END DO1203 END DO1204 END DO1205 !1206 ii0 = 53 ; ii1 = 53 ! Ombai Strait (e1v was modified)1207 ij0 = 164 - isrow ; ij1 = 165 - isrow1208 DO jk = 1, jpkm11209 DO jj = mj0(ij0), mj1(ij1)1210 DO ji = mi0(ii0), mi1(ii1)1211 SELECT CASE ( pout )1212 CASE( 'V' )1213 pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) &1214 & * ( e2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) &1215 & + e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &1216 & ) / e2v(ji,jj) + e3v_0(ji,jj,jk)1217 END SELECT1218 END DO1219 END DO1220 END DO1221 !1222 ii0 = 56 ; ii1 = 56 ! Timor Passage (e1v was modified)1223 ij0 = 164 - isrow ; ij1 = 165 - isrow1224 DO jk = 1, jpkm11225 DO jj = mj0(ij0), mj1(ij1)1226 DO ji = mi0(ii0), mi1(ii1)1227 SELECT CASE ( pout )1228 CASE( 'V' )1229 pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) &1230 & * ( e2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) &1231 & + e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &1232 & ) / e2v(ji,jj) + e3v_0(ji,jj,jk)1233 END SELECT1234 END DO1235 END DO1236 END DO1237 !1238 ii0 = 55 ; ii1 = 55 ! West Halmahera Strait (e1v was modified)1239 ij0 = 181 - isrow ; ij1 = 182 - isrow1240 DO jk = 1, jpkm11241 DO jj = mj0(ij0), mj1(ij1)1242 DO ji = mi0(ii0), mi1(ii1)1243 SELECT CASE ( pout )1244 CASE( 'V' )1245 pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) &1246 & * ( e2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) &1247 & + e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &1248 & ) / e2v(ji,jj) + e3v_0(ji,jj,jk)1249 END SELECT1250 END DO1251 END DO1252 END DO1253 !1254 ii0 = 58 ; ii1 = 58 ! East Halmahera Strait (e1v was modified)1255 ij0 = 181 - isrow ; ij1 = 182 - isrow1256 DO jk = 1, jpkm11257 DO jj = mj0(ij0), mj1(ij1)1258 DO ji = mi0(ii0), mi1(ii1)1259 SELECT CASE ( pout )1260 CASE( 'V' )1261 pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) &1262 & * ( e2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) &1263 & + e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &1264 & ) / e2v(ji,jj) + e3v_0(ji,jj,jk)1265 END SELECT1266 END DO1267 END DO1268 END DO1269 ENDIF1270 ! ! =====================1271 IF( cp_cfg == "orca" .AND. jp_cfg == 05 ) THEN ! ORCA R05 configuration1272 ! ! =====================1273 !1274 ii0 = 563 ; ii1 = 564 ! Gibraltar Strait (e2u was modified)1275 ij0 = 327 ; ij1 = 3271276 DO jk = 1, jpkm11277 DO jj = mj0(ij0), mj1(ij1)1278 DO ji = mi0(ii0), mi1(ii1)1279 SELECT CASE ( pout )1280 CASE( 'U' )1281 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) &1282 & * ( e1t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) &1283 & + e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &1284 & ) / e1u(ji,jj) + e3u_0(ji,jj,jk)1285 CASE( 'F' )1286 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) &1287 & * ( e1u(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3u_0(ji ,jj,jk) ) &1288 & + e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &1289 & ) / e1f(ji,jj) + e3f_0(ji,jj,jk)1290 END SELECT1291 END DO1292 END DO1293 END DO1294 !1295 ii0 = 627 ; ii1 = 628 ! Bosphorus Strait (e2u was modified)1296 ij0 = 343 ; ij1 = 3431297 DO jk = 1, jpkm11298 DO jj = mj0(ij0), mj1(ij1)1299 DO ji = mi0(ii0), mi1(ii1)1300 SELECT CASE ( pout )1301 CASE( 'U' )1302 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) &1303 & * ( e1t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) &1304 & + e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &1305 & ) / e1u(ji,jj) + e3u_0(ji,jj,jk)1306 CASE( 'F' )1307 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) &1308 & * ( e1u(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3u_0(ji ,jj,jk) ) &1309 & + e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &1310 & ) / e1f(ji,jj) + e3f_0(ji,jj,jk)1311 END SELECT1312 END DO1313 END DO1314 END DO1315 !1316 ii0 = 93 ; ii1 = 94 ! Sumba Strait (e2u was modified)1317 ij0 = 232 ; ij1 = 2321318 DO jk = 1, jpkm11319 DO jj = mj0(ij0), mj1(ij1)1320 DO ji = mi0(ii0), mi1(ii1)1321 SELECT CASE ( pout )1322 CASE( 'U' )1323 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) &1324 & * ( e1t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) &1325 & + e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &1326 & ) / e1u(ji,jj) + e3u_0(ji,jj,jk)1327 CASE( 'F' )1328 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) &1329 & * ( e1u(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3u_0(ji ,jj,jk) ) &1330 & + e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &1331 & ) / e1f(ji,jj) + e3f_0(ji,jj,jk)1332 END SELECT1333 END DO1334 END DO1335 END DO1336 !1337 ii0 = 103 ; ii1 = 103 ! Ombai Strait (e2u was modified)1338 ij0 = 232 ; ij1 = 2321339 DO jk = 1, jpkm11340 DO jj = mj0(ij0), mj1(ij1)1341 DO ji = mi0(ii0), mi1(ii1)1342 SELECT CASE ( pout )1343 CASE( 'U' )1344 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) &1345 & * ( e1t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) &1346 & + e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &1347 & ) / e1u(ji,jj) + e3u_0(ji,jj,jk)1348 CASE( 'F' )1349 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) &1350 & * ( e1u(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3u_0(ji ,jj,jk) ) &1351 & + e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &1352 & ) / e1f(ji,jj) + e3f_0(ji,jj,jk)1353 END SELECT1354 END DO1355 END DO1356 END DO1357 !1358 ii0 = 15 ; ii1 = 15 ! Palk Strait (e2u was modified)1359 ij0 = 270 ; ij1 = 2701360 DO jk = 1, jpkm11361 DO jj = mj0(ij0), mj1(ij1)1362 DO ji = mi0(ii0), mi1(ii1)1363 SELECT CASE ( pout )1364 CASE( 'U' )1365 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) &1366 & * ( e1t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) &1367 & + e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &1368 & ) / e1u(ji,jj) + e3u_0(ji,jj,jk)1369 CASE( 'F' )1370 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) &1371 & * ( e1u(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3u_0(ji ,jj,jk) ) &1372 & + e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &1373 & ) / e1f(ji,jj) + e3f_0(ji,jj,jk)1374 END SELECT1375 END DO1376 END DO1377 END DO1378 !1379 ii0 = 87 ; ii1 = 87 ! Lombok Strait (e1v was modified)1380 ij0 = 232 ; ij1 = 2331381 DO jk = 1, jpkm11382 DO jj = mj0(ij0), mj1(ij1)1383 DO ji = mi0(ii0), mi1(ii1)1384 SELECT CASE ( pout )1385 CASE( 'V' )1386 pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) &1387 & * ( e2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) &1388 & + e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &1389 & ) / e2v(ji,jj) + e3v_0(ji,jj,jk)1390 END SELECT1391 END DO1392 END DO1393 END DO1394 !1395 ii0 = 662 ; ii1 = 662 ! Bab el Mandeb (e1v was modified)1396 ij0 = 276 ; ij1 = 2761397 DO jk = 1, jpkm11398 DO jj = mj0(ij0), mj1(ij1)1399 DO ji = mi0(ii0), mi1(ii1)1400 SELECT CASE ( pout )1401 CASE( 'V' )1402 pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) &1403 & * ( e2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) &1404 & + e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &1405 & ) / e2v(ji,jj) + e3v_0(ji,jj,jk)1406 END SELECT1407 END DO1408 END DO1409 END DO1410 ENDIF1411 END SUBROUTINE dom_vvl_orca_fix1412 995 1413 996 !!====================================================================== -
trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domwri.F90
r5603 r5836 7 7 !! 8.1 ! 1999-11 (M. Imbard) NetCDF FORMAT with IOIPSL 8 8 !! NEMO 1.0 ! 2002-08 (G. Madec) F90 and several file 9 !! 3.0 ! 2008-01 (S. Masson) add dom_uniq 10 !! 4.0 ! 2011-01 (A. R. Porter, STFC Daresbury) dynamical allocation 9 !! 3.0 ! 2008-01 (S. Masson) add dom_uniq 11 10 !!---------------------------------------------------------------------- 12 11 13 12 !!---------------------------------------------------------------------- 14 13 !! dom_wri : create and write mesh and mask file(s) 15 !! dom_uniq : 14 !! dom_uniq : identify unique point of a grid (TUVF) 16 15 !!---------------------------------------------------------------------- 17 16 USE dom_oce ! ocean space and time domain … … 26 25 PRIVATE 27 26 28 PUBLIC dom_wri! routine called by inidom.F9029 27 PUBLIC dom_wri ! routine called by inidom.F90 28 PUBLIC dom_wri_coordinate ! routine called by domhgr.F90 30 29 !! * Substitutions 31 30 # include "vectopt_loop_substitute.h90" … … 36 35 !!---------------------------------------------------------------------- 37 36 CONTAINS 37 38 SUBROUTINE dom_wri_coordinate 39 !!---------------------------------------------------------------------- 40 !! *** ROUTINE dom_wri_coordinate *** 41 !! 42 !! ** Purpose : Create the NetCDF file which contains all the 43 !! standard coordinate information plus the surface, 44 !! e1e2u and e1e2v. By doing so, those surface will 45 !! not be changed by the reduction of e1u or e2v scale 46 !! factors in some straits. 47 !! NB: call just after the read of standard coordinate 48 !! and the reduction of scale factors in some straits 49 !! 50 !! ** output file : coordinate_e1e2u_v.nc 51 !!---------------------------------------------------------------------- 52 INTEGER :: inum0 ! temprary units for 'coordinate_e1e2u_v.nc' file 53 CHARACTER(len=21) :: clnam0 ! filename (mesh and mask informations) 54 ! ! workspaces 55 REAL(wp), POINTER, DIMENSION(:,: ) :: zprt, zprw 56 REAL(wp), POINTER, DIMENSION(:,:,:) :: zdepu, zdepv 57 !!---------------------------------------------------------------------- 58 ! 59 IF( nn_timing == 1 ) CALL timing_start('dom_wri_coordinate') 60 ! 61 IF(lwp) WRITE(numout,*) 62 IF(lwp) WRITE(numout,*) 'dom_wri_coordinate : create NetCDF coordinate file' 63 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~~~' 64 65 clnam0 = 'coordinate_e1e2u_v' ! filename (mesh and mask informations) 66 67 ! create 'coordinate_e1e2u_v.nc' file 68 ! ============================ 69 ! 70 CALL iom_open( TRIM(clnam0), inum0, ldwrt = .TRUE., kiolib = jprstlib ) 71 ! 72 ! ! horizontal mesh (inum3) 73 CALL iom_rstput( 0, 0, inum0, 'glamt', glamt, ktype = jp_r4 ) ! ! latitude 74 CALL iom_rstput( 0, 0, inum0, 'glamu', glamu, ktype = jp_r4 ) 75 CALL iom_rstput( 0, 0, inum0, 'glamv', glamv, ktype = jp_r4 ) 76 CALL iom_rstput( 0, 0, inum0, 'glamf', glamf, ktype = jp_r4 ) 77 78 CALL iom_rstput( 0, 0, inum0, 'gphit', gphit, ktype = jp_r4 ) ! ! longitude 79 CALL iom_rstput( 0, 0, inum0, 'gphiu', gphiu, ktype = jp_r4 ) 80 CALL iom_rstput( 0, 0, inum0, 'gphiv', gphiv, ktype = jp_r4 ) 81 CALL iom_rstput( 0, 0, inum0, 'gphif', gphif, ktype = jp_r4 ) 82 83 CALL iom_rstput( 0, 0, inum0, 'e1t', e1t, ktype = jp_r8 ) ! ! e1 scale factors 84 CALL iom_rstput( 0, 0, inum0, 'e1u', e1u, ktype = jp_r8 ) 85 CALL iom_rstput( 0, 0, inum0, 'e1v', e1v, ktype = jp_r8 ) 86 CALL iom_rstput( 0, 0, inum0, 'e1f', e1f, ktype = jp_r8 ) 87 88 CALL iom_rstput( 0, 0, inum0, 'e2t', e2t, ktype = jp_r8 ) ! ! e2 scale factors 89 CALL iom_rstput( 0, 0, inum0, 'e2u', e2u, ktype = jp_r8 ) 90 CALL iom_rstput( 0, 0, inum0, 'e2v', e2v, ktype = jp_r8 ) 91 CALL iom_rstput( 0, 0, inum0, 'e2f', e2f, ktype = jp_r8 ) 92 93 CALL iom_rstput( 0, 0, inum0, 'e1e2u', e1e2u, ktype = jp_r8 ) 94 CALL iom_rstput( 0, 0, inum0, 'e1e2v', e1e2v, ktype = jp_r8 ) 95 96 CALL iom_close( inum0 ) 97 ! 98 IF( nn_timing == 1 ) CALL timing_stop('dom_wri_coordinate') 99 ! 100 END SUBROUTINE dom_wri_coordinate 101 38 102 39 103 SUBROUTINE dom_wri -
trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
r5656 r5836 501 501 IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! ORCA R2 configuration 502 502 ! ! ===================== 503 IF( nn_cla == 0 ) THEN 504 ii0 = 140 ; ii1 = 140 ! Gibraltar Strait open 505 ij0 = 102 ; ij1 = 102 ! (Thomson, Ocean Modelling, 1995) 506 DO ji = mi0(ii0), mi1(ii1) 507 DO jj = mj0(ij0), mj1(ij1) 508 mbathy(ji,jj) = 15 509 END DO 503 ! 504 ii0 = 140 ; ii1 = 140 ! Gibraltar Strait open 505 ij0 = 102 ; ij1 = 102 ! (Thomson, Ocean Modelling, 1995) 506 DO ji = mi0(ii0), mi1(ii1) 507 DO jj = mj0(ij0), mj1(ij1) 508 mbathy(ji,jj) = 15 510 509 END DO 511 IF(lwp) WRITE(numout,*)512 IF(lwp) WRITE(numout,*) ' orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0513 !514 ii0 = 160 ; ii1 = 160 ! Bab el mandeb Strait open515 ij0 = 88 ; ij1 = 88 ! (Thomson, Ocean Modelling, 1995)516 DO ji = mi0(ii0), mi1(ii1)517 DO jj = mj0(ij0), mj1(ij1)518 mbathy(ji,jj) = 12519 END DO510 END DO 511 IF(lwp) WRITE(numout,*) 512 IF(lwp) WRITE(numout,*) ' orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0 513 ! 514 ii0 = 160 ; ii1 = 160 ! Bab el mandeb Strait open 515 ij0 = 88 ; ij1 = 88 ! (Thomson, Ocean Modelling, 1995) 516 DO ji = mi0(ii0), mi1(ii1) 517 DO jj = mj0(ij0), mj1(ij1) 518 mbathy(ji,jj) = 12 520 519 END DO 521 IF(lwp) WRITE(numout,*)522 IF(lwp) WRITE(numout,*) ' orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0523 ENDIF520 END DO 521 IF(lwp) WRITE(numout,*) 522 IF(lwp) WRITE(numout,*) ' orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0 524 523 ! 525 524 ENDIF … … 545 544 ! 546 545 IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! ORCA R2 configuration 547 ! 548 IF( nn_cla == 0 ) THEN 549 ii0 = 140 ; ii1 = 140 ! Gibraltar Strait open 550 ij0 = 102 ; ij1 = 102 ! (Thomson, Ocean Modelling, 1995) 551 DO ji = mi0(ii0), mi1(ii1) 552 DO jj = mj0(ij0), mj1(ij1) 553 bathy(ji,jj) = 284._wp 554 END DO 546 ! 547 ii0 = 140 ; ii1 = 140 ! Gibraltar Strait open 548 ij0 = 102 ; ij1 = 102 ! (Thomson, Ocean Modelling, 1995) 549 DO ji = mi0(ii0), mi1(ii1) 550 DO jj = mj0(ij0), mj1(ij1) 551 bathy(ji,jj) = 284._wp 555 552 END DO 556 IF(lwp) WRITE(numout,*)557 IF(lwp) WRITE(numout,*) ' orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0558 !559 ii0 = 160 ; ii1 = 160 ! Bab el mandeb Strait open560 ij0 = 88 ; ij1 = 88 ! (Thomson, Ocean Modelling, 1995)561 DO ji = mi0(ii0), mi1(ii1)562 DO jj = mj0(ij0), mj1(ij1)563 bathy(ji,jj) = 137._wp564 END DO553 END DO 554 IF(lwp) WRITE(numout,*) 555 IF(lwp) WRITE(numout,*) ' orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0 556 ! 557 ii0 = 160 ; ii1 = 160 ! Bab el mandeb Strait open 558 ij0 = 88 ; ij1 = 88 ! (Thomson, Ocean Modelling, 1995) 559 DO ji = mi0(ii0), mi1(ii1) 560 DO jj = mj0(ij0), mj1(ij1) 561 bathy(ji,jj) = 137._wp 565 562 END DO 566 IF(lwp) WRITE(numout,*)567 IF(lwp) WRITE(numout,*) ' orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0568 ENDIF563 END DO 564 IF(lwp) WRITE(numout,*) 565 IF(lwp) WRITE(numout,*) ' orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0 569 566 ! 570 567 ENDIF -
trunk/NEMOGCM/NEMO/OPA_SRC/DOM/dtatsd.F90
r5215 r5836 174 174 END DO 175 175 END DO 176 IF( nn_cla == 1 ) THEN ! Cross Land advection 177 il0 = 138 ; il1 = 138 ! set T & S profile at Gibraltar Strait 178 ij0 = 101 ; ij1 = 102 179 ii0 = 139 ; ii1 = 139 180 DO jl = mi0(il0), mi1(il1) 181 DO jj = mj0(ij0), mj1(ij1) 182 DO ji = mi0(ii0), mi1(ii1) 183 sf_tsd(jp_tem)%fnow(ji,jj,:) = sf_tsd(jp_tem)%fnow(jl,jj,:) 184 sf_tsd(jp_sal)%fnow(ji,jj,:) = sf_tsd(jp_sal)%fnow(jl,jj,:) 185 END DO 186 END DO 187 END DO 188 il0 = 164 ; il1 = 164 ! set T & S profile at Bab el Mandeb Strait 189 ij0 = 87 ; ij1 = 88 190 ii0 = 161 ; ii1 = 163 191 DO jl = mi0(il0), mi1(il1) 192 DO jj = mj0(ij0), mj1(ij1) 193 DO ji = mi0(ii0), mi1(ii1) 194 sf_tsd(jp_tem)%fnow(ji,jj,:) = sf_tsd(jp_tem)%fnow(jl,jj,:) 195 sf_tsd(jp_sal)%fnow(ji,jj,:) = sf_tsd(jp_sal)%fnow(jl,jj,:) 196 END DO 197 END DO 198 END DO 199 ELSE ! No Cross Land advection 200 ij0 = 87 ; ij1 = 96 ! Reduced temperature in Red Sea 201 ii0 = 148 ; ii1 = 160 202 sf_tsd(jp_tem)%fnow( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 4:10 ) = 7.0_wp 203 sf_tsd(jp_tem)%fnow( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 11:13 ) = 6.5_wp 204 sf_tsd(jp_tem)%fnow( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 14:20 ) = 6.0_wp 205 ENDIF 176 ij0 = 87 ; ij1 = 96 ! Reduced temperature in Red Sea 177 ii0 = 148 ; ii1 = 160 178 sf_tsd(jp_tem)%fnow( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 4:10 ) = 7.0_wp 179 sf_tsd(jp_tem)%fnow( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 11:13 ) = 6.5_wp 180 sf_tsd(jp_tem)%fnow( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 14:20 ) = 6.0_wp 206 181 ENDIF 207 182 ! -
trunk/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90
r5332 r5836 29 29 USE daymod ! calendar 30 30 USE eosbn2 ! eq. of state, Brunt Vaisala frequency (eos routine) 31 USE ldftra _oce ! ocean active tracers: lateral physics31 USE ldftra ! lateral physics: ocean active tracers 32 32 USE zdf_oce ! ocean vertical physics 33 33 USE phycst ! physical constants 34 34 USE dtatsd ! data temperature and salinity (dta_tsd routine) 35 35 USE dtauvd ! data: U & V current (dta_uvd routine) 36 USE zpshde ! partial step: hor. derivative (zps_hde routine)37 USE eosbn2 ! equation of state (eos bn2 routine)38 36 USE domvvl ! varying vertical mesh 39 37 USE dynspg_oce ! pressure gradient schemes … … 76 74 ! 77 75 78 IF(lwp) WRITE(numout,*) ' '76 IF(lwp) WRITE(numout,*) 79 77 IF(lwp) WRITE(numout,*) 'istate_ini : Initialization of the dynamics and tracers' 80 78 IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 81 79 82 CALL dta_tsd_init! Initialisation of T & S input data83 IF( lk_c1d ) CALL dta_uvd_init! Initialization of U & V input data80 CALL dta_tsd_init ! Initialisation of T & S input data 81 IF( lk_c1d ) CALL dta_uvd_init ! Initialization of U & V input data 84 82 85 83 rhd (:,:,: ) = 0._wp ; rhop (:,:,: ) = 0._wp ! set one for all to 0 at level jpk … … 103 101 ub (:,:,:) = 0._wp ; un (:,:,:) = 0._wp 104 102 vb (:,:,:) = 0._wp ; vn (:,:,:) = 0._wp 105 rotb (:,:,:) = 0._wp ; rotn (:,:,:) = 0._wp 106 hdivb(:,:,:) = 0._wp ; hdivn(:,:,:) = 0._wp 103 hdivn(:,:,:) = 0._wp 107 104 ! 108 105 IF( cp_cfg == 'eel' ) THEN … … 119 116 ENDIF 120 117 IF ( ln_uvd_init .AND. lk_c1d ) THEN ! read 3D U and V data at nit000 121 CALL wrk_alloc( jpi, jpj, jpk, 2,zuvd )118 CALL wrk_alloc( jpi,jpj,jpk,2, zuvd ) 122 119 CALL dta_uvd( nit000, zuvd ) 123 120 ub(:,:,:) = zuvd(:,:,:,1) ; un(:,:,:) = ub(:,:,:) 124 121 vb(:,:,:) = zuvd(:,:,:,2) ; vn(:,:,:) = vb(:,:,:) 125 CALL wrk_dealloc( jpi, jpj, jpk, 2,zuvd )122 CALL wrk_dealloc( jpi,jpj,jpk,2, zuvd ) 126 123 ENDIF 127 124 ENDIF 128 !129 CALL eos( tsb, rhd, rhop, gdept_0(:,:,:) ) ! before potential and in situ densities130 #if ! defined key_c1d131 IF( ln_zps .AND. .NOT. ln_isfcav) &132 & CALL zps_hde ( nit000, jpts, tsb, gtsu, gtsv, & ! Partial steps: before horizontal gradient133 & rhd, gru , grv ) ! of t, s, rd at the last ocean level134 IF( ln_zps .AND. ln_isfcav) &135 & CALL zps_hde_isf( nit000, jpts, tsb, gtsu, gtsv, & ! Partial steps for top cell (ISF)136 & rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv , &137 & gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi ) ! of t, s, rd at the last ocean level138 #endif139 125 ! 140 126 ! - ML - sshn could be modified by istate_eel, so that initialization of fse3t_b is done here … … 142 128 DO jk = 1, jpk 143 129 fse3t_b(:,:,jk) = fse3t_n(:,:,jk) 144 END DO130 END DO 145 131 ENDIF 146 132 ! … … 163 149 ! 164 150 ! 165 un_b(:,:) = 0._wp ;vn_b(:,:) = 0._wp166 ub_b(:,:) = 0._wp ;vb_b(:,:) = 0._wp151 un_b(:,:) = 0._wp ; vn_b(:,:) = 0._wp 152 ub_b(:,:) = 0._wp ; vb_b(:,:) = 0._wp 167 153 ! 168 154 DO jk = 1, jpkm1 … … 202 188 !! References : Philander ??? 203 189 !!---------------------------------------------------------------------- 204 INTEGER :: ji, jj, jk205 REAL(wp) :: zsal = 35.50 190 INTEGER :: ji, jj, jk 191 REAL(wp) :: zsal = 35.50_wp 206 192 !!---------------------------------------------------------------------- 207 193 ! … … 233 219 !! and relative vorticity fields 234 220 !!---------------------------------------------------------------------- 235 USE div cur ! hor. divergence & rel. vorticity (div_cur routine)221 USE divhor ! hor. divergence (div_hor routine) 236 222 USE iom 237 223 ! … … 282 268 tsb(:,:,:,jp_sal) = tsn(:,:,:,jp_sal) 283 269 ! 284 ! set the dynamics: U,V, hdiv , rot(and ssh if necessary)270 ! set the dynamics: U,V, hdiv (and ssh if necessary) 285 271 ! ---------------- 286 272 ! Start EEL5 configuration with barotropic geostrophic velocities … … 318 304 ENDIF 319 305 ! 320 CALL div_cur( nit000 ) ! horizontal divergence and relative vorticity (curl) 306 !!gm Check here call to div_hor should not be necessary 307 !!gm div_hor call runoffs not sure they are defined at that level 308 CALL div_hor( nit000 ) ! horizontal divergence and relative vorticity (curl) 321 309 ! N.B. the vertical velocity will be computed from the horizontal divergence field 322 310 ! in istate by a call to wzv routine … … 371 359 !! 372 360 !! ** Method : - set temprature field 373 !! - set salinity field361 !! - set salinity field 374 362 !!---------------------------------------------------------------------- 375 363 INTEGER :: ji, jj, jk ! dummy loop indices … … 458 446 !!---------------------------------------------------------------------- 459 447 USE dynspg ! surface pressure gradient (dyn_spg routine) 460 USE div cur ! hor. divergence & rel. vorticity (div_cur routine)448 USE divhor ! hor. divergence (div_hor routine) 461 449 USE lbclnk ! ocean lateral boundary condition (or mpp link) 462 450 ! … … 467 455 !!---------------------------------------------------------------------- 468 456 ! 469 CALL wrk_alloc( jpi, jpj, jpk,zprn)457 CALL wrk_alloc( jpi,jpj,jpk, zprn) 470 458 ! 471 459 IF(lwp) WRITE(numout,*) … … 544 532 indic = 0 545 533 CALL dyn_spg( nit000, indic ) ! surface pressure gradient 546 534 ! 547 535 ! the new velocity is ua*rdt 548 536 ! 549 537 CALL lbc_lnk( ua, 'U', -1. ) 550 538 CALL lbc_lnk( va, 'V', -1. ) … … 556 544 un(:,:,:) = ub(:,:,:) 557 545 vn(:,:,:) = vb(:,:,:) 558 559 ! Compute the divergence and curl 560 561 CALL div_cur( nit000 ) ! now horizontal divergence and curl 562 563 hdivb(:,:,:) = hdivn(:,:,:) ! set the before to the now value 564 rotb (:,:,:) = rotn (:,:,:) ! set the before to the now value 565 ! 566 CALL wrk_dealloc( jpi, jpj, jpk, zprn) 546 ! 547 !!gm Check here call to div_hor should not be necessary 548 !!gm div_hor call runoffs not sure they are defined at that level 549 CALL div_hor( nit000 ) ! now horizontal divergence 550 ! 551 CALL wrk_dealloc( jpi,jpj,jpk, zprn) 567 552 ! 568 553 END SUBROUTINE istate_uvg
Note: See TracChangeset
for help on using the changeset viewer.