- Timestamp:
- 2015-12-08T12:39:53+01:00 (8 years ago)
- Location:
- branches/UKMO/icebergs_restart_single_file/NEMOGCM/NEMO/OPA_SRC/DOM
- Files:
-
- 12 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/UKMO/icebergs_restart_single_file/NEMOGCM/NEMO/OPA_SRC/DOM/closea.F90
r6019 r6020 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 ! -
branches/UKMO/icebergs_restart_single_file/NEMOGCM/NEMO/OPA_SRC/DOM/daymod.F90
r6019 r6020 73 73 !!---------------------------------------------------------------------- 74 74 ! 75 ! max number of seconds between each restart 76 IF( REAL( nitend - nit000 + 1 ) * rdt > REAL( HUGE( nsec1jan000 ) ) ) THEN 77 CALL ctl_stop( 'The number of seconds between each restart exceeds the integer 4 max value: 2^31-1. ', & 78 & 'You must do a restart at higher frequency (or remove this stop and recompile the code in I8)' ) 79 ENDIF 75 80 ! all calendar staff is based on the fact that MOD( rday, rdttra(1) ) == 0 76 81 IF( MOD( rday , rdttra(1) ) /= 0. ) CALL ctl_stop( 'the time step must devide the number of second of in a day' ) … … 238 243 nday_year = 1 239 244 nsec_year = ndt05 240 IF( nsec1jan000 >= 2 * (2**30 - nsecd * nyear_len(1) / 2 ) ) THEN ! test integer 4 max value241 CALL ctl_stop( 'The number of seconds between Jan. 1st 00h of nit000 year and Jan. 1st 00h ', &242 & 'of the current year is exceeding the INTEGER 4 max VALUE: 2^31-1 -> 68.09 years in seconds', &243 & 'You must do a restart at higher frequency (or remove this STOP and recompile everything in I8)' )244 ENDIF245 245 nsec1jan000 = nsec1jan000 + nsecd * nyear_len(1) 246 246 IF( nleapy == 1 ) CALL day_mth -
branches/UKMO/icebergs_restart_single_file/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90
r6019 r6020 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 … … 45 45 LOGICAL , PUBLIC :: ln_crs !: Apply grid coarsening to dynamical model output or online passive tracers 46 46 47 !! Free surface parameters 48 !! ======================= 49 LOGICAL, PUBLIC :: ln_dynspg_exp !: Explicit free surface flag 50 LOGICAL, PUBLIC :: ln_dynspg_ts !: Split-Explicit free surface flag 51 47 52 !! Time splitting parameters 48 53 !! ========================= 49 54 LOGICAL, PUBLIC :: ln_bt_fw !: Forward integration of barotropic sub-stepping 50 55 LOGICAL, PUBLIC :: ln_bt_av !: Time averaging of barotropic variables 51 LOGICAL, PUBLIC :: ln_bt_ nn_auto!: Set number of barotropic iterations automatically56 LOGICAL, PUBLIC :: ln_bt_auto !: Set number of barotropic iterations automatically 52 57 INTEGER, PUBLIC :: nn_bt_flt !: Filter choice 53 58 INTEGER, PUBLIC :: nn_baro !: Number of barotropic iterations during one baroclinic step (rdt) 54 REAL(wp), PUBLIC :: rn_bt_cmax !: Maximum allowed courant number (used if ln_bt_ nn_auto=T)59 REAL(wp), PUBLIC :: rn_bt_cmax !: Maximum allowed courant number (used if ln_bt_auto=T) 55 60 56 61 !! Horizontal grid parameters for domhgr … … 107 112 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: r2dtra !: = 2*rdttra except at nit000 (=rdttra) if neuler=0 108 113 109 ! !!* Namelist namcla : cross land advection110 INTEGER, PUBLIC :: nn_cla !: =1 cross land advection for exchanges through some straits (ORCA2)111 112 114 !!---------------------------------------------------------------------- 113 115 !! space domain parameters … … 158 160 !! horizontal curvilinear coordinate and scale factors 159 161 !! --------------------------------------------------------------------- 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) 162 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: glamt , glamu, glamv , glamf !: longitude at t, u, v, f-points [degree] 163 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: gphit , gphiu, gphiv , gphif !: latitude at t, u, v, f-points [degree] 164 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) :: e1t , e2t , r1_e1t, r1_e2t !: t-point horizontal scale factors [m] 165 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) :: e1u , e2u , r1_e1u, r1_e2u !: horizontal scale factors at u-point [m] 166 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) :: e1v , e2v , r1_e1v, r1_e2v !: horizontal scale factors at v-point [m] 167 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) :: e1f , e2f , r1_e1f, r1_e2f !: horizontal scale factors at f-point [m] 168 ! 169 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: e1e2t , r1_e1e2t !: associated metrics at t-point 170 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: e1e2u , r1_e1e2u , e2_e1u !: associated metrics at u-point 171 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: e1e2v , r1_e1e2v , e1_e2v !: associated metrics at v-point 172 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: e1e2f , r1_e1e2f !: associated metrics at f-point 173 ! 174 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ff !: coriolis factor [1/s] 170 175 171 176 !!---------------------------------------------------------------------- … … 216 221 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ht_0 !: reference depth at t- points (meters) 217 222 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 223 225 224 INTEGER, PUBLIC :: nla10 !: deepest W level Above ~10m (nlb10 - 1) … … 265 264 266 265 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 266 273 267 !!---------------------------------------------------------------------- … … 333 327 INTEGER FUNCTION dom_oce_alloc() 334 328 !!---------------------------------------------------------------------- 335 INTEGER, DIMENSION(1 2) :: ierr329 INTEGER, DIMENSION(13) :: ierr 336 330 !!---------------------------------------------------------------------- 337 331 ierr(:) = 0 … … 346 340 & tpol(jpiglo) , fpol(jpiglo) , STAT=ierr(2) ) 347 341 ! 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) ) 342 ALLOCATE( glamt(jpi,jpj) , glamu(jpi,jpj) , glamv(jpi,jpj) , glamf(jpi,jpj) , & 343 & gphit(jpi,jpj) , gphiu(jpi,jpj) , gphiv(jpi,jpj) , gphif(jpi,jpj) , & 344 & e1t (jpi,jpj) , e2t (jpi,jpj) , r1_e1t(jpi,jpj) , r1_e2t(jpi,jpj) , & 345 & e1u (jpi,jpj) , e2u (jpi,jpj) , r1_e1u(jpi,jpj) , r1_e2u(jpi,jpj) , & 346 & e1v (jpi,jpj) , e2v (jpi,jpj) , r1_e1v(jpi,jpj) , r1_e2v(jpi,jpj) , & 347 & e1f (jpi,jpj) , e2f (jpi,jpj) , r1_e1f(jpi,jpj) , r1_e2f(jpi,jpj) , & 348 & e1e2t(jpi,jpj) , r1_e1e2t(jpi,jpj) , & 349 & e1e2u(jpi,jpj) , r1_e1e2u(jpi,jpj) , e2_e1u(jpi,jpj) , & 350 & e1e2v(jpi,jpj) , r1_e1e2v(jpi,jpj) , e1_e2v(jpi,jpj) , & 351 & e1e2f(jpi,jpj) , r1_e1e2f(jpi,jpj) , & 352 & ff (jpi,jpj) , STAT=ierr(3) ) 353 353 ! 354 354 ALLOCATE( gdep3w_0(jpi,jpj,jpk) , e3v_0(jpi,jpj,jpk) , e3f_0 (jpi,jpj,jpk) , & … … 364 364 & gdept_b (jpi,jpj,jpk) ,gdepw_b(jpi,jpj,jpk) , e3w_b (jpi,jpj,jpk) , & 365 365 & 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) )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) ) 370 370 #endif 371 371 ! 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) ) 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) , STAT=ierr(6) ) 379 374 ! 380 375 ALLOCATE( gdept_1d(jpk) , gdepw_1d(jpk) , & … … 387 382 & scosrf(jpi,jpj) , scobot(jpi,jpj) , & 388 383 & hifv (jpi,jpj) , hiff (jpi,jpj) , & 389 & hift (jpi,jpj) , hifu (jpi,jpj) , rx1 384 & hift (jpi,jpj) , hifu (jpi,jpj) , rx1(jpi,jpj) , STAT=ierr(8) ) 390 385 391 386 ALLOCATE( mbathy(jpi,jpj) , bathy(jpi,jpj) , & 392 387 & tmask_i(jpi,jpj) , umask_i(jpi,jpj), vmask_i(jpi,jpj), fmask_i(jpi,jpj), & 393 & bmask (jpi,jpj), &388 & bmask (jpi,jpj) , & 394 389 & mbkt (jpi,jpj) , mbku (jpi,jpj) , mbkv(jpi,jpj) , STAT=ierr(9) ) 395 390 396 391 ! (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) )392 ALLOCATE( misfdep(jpi,jpj) , risfdep(jpi,jpj), & 393 & mikt(jpi,jpj), miku(jpi,jpj), mikv(jpi,jpj) , & 394 & mikf(jpi,jpj), ssmask(jpi,jpj), STAT=ierr(10) ) 400 395 401 396 ALLOCATE( tmask(jpi,jpj,jpk) , umask(jpi,jpj,jpk), & … … 403 398 404 399 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 400 ! 410 401 dom_oce_alloc = MAXVAL(ierr) -
branches/UKMO/icebergs_restart_single_file/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90
r6019 r6020 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 !!====================================================================== -
branches/UKMO/icebergs_restart_single_file/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90
r6019 r6020 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 = 201 + 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 = 208 + 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 = 124 + 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 = 124 + 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 = 124 + 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 = 124 + 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 = 141 + 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 = 141 + 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) / (ra * rad)395 zphi0 = ppgphi0 - REAL( Agrif_Parent(jpjglo)/2)*Agrif_Parent(ppe2_m) / (ra * rad) 547 396 ENDIF 548 397 ENDIF 549 398 #endif 550 399 zf0 = 2. * omega * SIN( rad * zphi0 ) ! compute f0 1st point south 551 400 ! 552 401 ff(:,:) = ( zf0 + zbeta * gphif(:,:) * 1.e+3 ) ! f = f0 +beta* y ( y=0 at south) 553 402 ! 554 403 IF(lwp) THEN 555 404 WRITE(numout,*) … … 564 413 IF(lwp) WRITE(numout,*) ' Coriolis parameter varies globally from ', zminff,' to ', zmaxff 565 414 END IF 566 415 ! 567 416 CASE ( 5 ) ! beta-plane and rotated domain (gyre configuration) 568 417 ! 569 418 zbeta = 2. * omega * COS( rad * ppgphi0 ) / ra ! beta at latitude ppgphi0 570 zphi0 = 15. e0! latitude of the first row F-points419 zphi0 = 15._wp ! latitude of the first row F-points 571 420 zf0 = 2. * omega * SIN( rad * zphi0 ) ! compute f0 1st point south 572 421 ! 573 422 ff(:,:) = ( zf0 + zbeta * ABS( gphif(:,:) - zphi0 ) * rad * ra ) ! f = f0 +beta* y ( y=0 at south) 574 423 ! 575 424 IF(lwp) THEN 576 425 WRITE(numout,*) … … 578 427 WRITE(numout,*) ' Coriolis parameter varies in this processor from ', ff(nldi,nldj),' to ', ff(nldi,nlej) 579 428 ENDIF 580 429 ! 581 430 IF( lk_mpp ) THEN 582 431 zminff=ff(nldi,nldj) … … 586 435 IF(lwp) WRITE(numout,*) ' Coriolis parameter varies globally from ', zminff,' to ', zmaxff 587 436 END IF 588 437 ! 589 438 END SELECT 590 439 … … 595 444 596 445 IF( nperio == 2 ) THEN 597 znorme = SQRT( SUM( gphiu(:,2) * gphiu(:,2) ) ) / FLOAT( jpi )446 znorme = SQRT( SUM( gphiu(:,2) * gphiu(:,2) ) ) / REAL( jpi ) 598 447 IF( znorme > 1.e-13 ) CALL ctl_stop( ' ===>>>> : symmetrical condition: rerun with good equator line' ) 599 448 ENDIF … … 604 453 605 454 606 SUBROUTINE hgr_read 455 SUBROUTINE hgr_read( ke1e2u_v ) 607 456 !!--------------------------------------------------------------------- 608 457 !! *** ROUTINE hgr_read *** 609 458 !! 610 !! ** Purpose : Read a coordinate file in NetCDF format 611 !! 612 !! ** Method : The mesh file has been defined trough a analytical 613 !! or semi-analytical method. It is read in a NetCDF file. 614 !! 459 !! ** Purpose : Read a coordinate file in NetCDF format using IOM 460 !! 615 461 !!---------------------------------------------------------------------- 616 462 USE iom 617 463 !! 464 INTEGER, INTENT( inout ) :: ke1e2u_v ! fag: e1e2u & e1e2v read in coordinate file (=1) or not (=0) 465 ! 618 466 INTEGER :: inum ! temporary logical unit 619 467 !!---------------------------------------------------------------------- 620 468 ! 621 469 IF(lwp) THEN 622 470 WRITE(numout,*) … … 624 472 WRITE(numout,*) '~~~~~~~~ jpiglo = ', jpiglo, ' jpjglo = ', jpjglo, ' jpk = ', jpk 625 473 ENDIF 626 474 ! 627 475 CALL iom_open( 'coordinates', inum ) 628 476 ! 629 477 CALL iom_get( inum, jpdom_data, 'glamt', glamt, lrowattr=ln_use_jattr ) 630 478 CALL iom_get( inum, jpdom_data, 'glamu', glamu, lrowattr=ln_use_jattr ) 631 479 CALL iom_get( inum, jpdom_data, 'glamv', glamv, lrowattr=ln_use_jattr ) 632 480 CALL iom_get( inum, jpdom_data, 'glamf', glamf, lrowattr=ln_use_jattr ) 633 481 ! 634 482 CALL iom_get( inum, jpdom_data, 'gphit', gphit, lrowattr=ln_use_jattr ) 635 483 CALL iom_get( inum, jpdom_data, 'gphiu', gphiu, lrowattr=ln_use_jattr ) 636 484 CALL iom_get( inum, jpdom_data, 'gphiv', gphiv, lrowattr=ln_use_jattr ) 637 485 CALL iom_get( inum, jpdom_data, 'gphif', gphif, lrowattr=ln_use_jattr ) 638 639 CALL iom_get( inum, jpdom_data, 'e1t', e1t, lrowattr=ln_use_jattr ) 640 CALL iom_get( inum, jpdom_data, 'e1u', e1u, lrowattr=ln_use_jattr ) 641 CALL iom_get( inum, jpdom_data, 'e1v', e1v, lrowattr=ln_use_jattr ) 642 CALL iom_get( inum, jpdom_data, 'e1f', e1f, lrowattr=ln_use_jattr ) 643 644 CALL iom_get( inum, jpdom_data, 'e2t', e2t, lrowattr=ln_use_jattr ) 645 CALL iom_get( inum, jpdom_data, 'e2u', e2u, lrowattr=ln_use_jattr ) 646 CALL iom_get( inum, jpdom_data, 'e2v', e2v, lrowattr=ln_use_jattr ) 647 CALL iom_get( inum, jpdom_data, 'e2f', e2f, lrowattr=ln_use_jattr ) 648 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 ! 649 506 CALL iom_close( inum ) 650 507 508 !!gm THIS is TO BE REMOVED !!!!!!! 509 651 510 ! need to be define for the extended grid south of -80S 652 511 ! some point are undefined but you need to have e1 and e2 .NE. 0 … … 675 534 e2f=1.0e2 676 535 END WHERE 536 !!gm end 677 537 678 538 END SUBROUTINE hgr_read -
branches/UKMO/icebergs_restart_single_file/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90
r6019 r6020 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 … … 28 28 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 29 29 USE lib_mpp 30 USE dynspg_oce ! choice/control of key cpp for surface pressure gradient31 30 USE wrk_nemo ! Memory allocation 32 31 USE timing ! Timing … … 36 35 37 36 PUBLIC dom_msk ! routine called by inidom.F90 38 PUBLIC dom_msk_alloc ! routine called by nemogcm.F9039 37 40 38 ! !!* Namelist namlbc : lateral boundary condition * … … 42 40 LOGICAL, PUBLIC :: ln_vorlat ! consistency of vorticity boundary condition 43 41 ! with analytical eqs. 44 45 46 INTEGER, ALLOCATABLE, SAVE, DIMENSION(:,:) :: icoord ! Workspace for dom_msk_nsa()47 42 48 43 !! * Substitutions … … 54 49 !!---------------------------------------------------------------------- 55 50 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 51 70 52 SUBROUTINE dom_msk … … 129 111 !! tmask_i : interior ocean mask 130 112 !!---------------------------------------------------------------------- 131 ! 132 INTEGER :: ji, jj, jk ! dummy loop indices 113 INTEGER :: ji, jj, jk ! dummy loop indices 133 114 INTEGER :: iif, iil, ii0, ii1, ii ! local integers 134 115 INTEGER :: ijf, ijl, ij0, ij1 ! - - … … 199 180 END DO 200 181 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 182 ! Interior domain mask (used for global sum) 222 183 ! -------------------- … … 284 245 ! 3. Ocean/land mask at wu-, wv- and w points 285 246 !---------------------------------------------- 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)247 wmask (:,:,1) = tmask(:,:,1) ! surface 248 wumask(:,:,1) = umask(:,:,1) 249 wvmask(:,:,1) = vmask(:,:,1) 250 DO jk = 2, jpk ! interior values 251 wmask (:,:,jk) = tmask(:,:,jk) * tmask(:,:,jk-1) 252 wumask(:,:,jk) = umask(:,:,jk) * umask(:,:,jk-1) 253 wvmask(:,:,jk) = vmask(:,:,jk) * vmask(:,:,jk-1) 293 254 END DO 294 255 … … 339 300 ENDIF 340 301 341 342 ! mask for second order calculation of vorticity343 ! ----------------------------------------------344 CALL dom_msk_nsa345 346 347 302 ! Lateral boundary conditions on velocity (modify fmask) 348 303 ! --------------------------------------- … … 377 332 IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! ORCA_R2 configuration 378 333 ! ! 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 334 ! ! Gibraltar strait : partial slip (fmask=0.5) 335 ij0 = 101 ; ij1 = 101 336 ii0 = 139 ; ii1 = 140 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 0.5_wp 337 ij0 = 102 ; ij1 = 102 338 ii0 = 139 ; ii1 = 140 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 0.5_wp 339 ! 340 ! ! Bab el Mandeb : partial slip (fmask=1) 341 ij0 = 87 ; ij1 = 88 342 ii0 = 160 ; ii1 = 160 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 1._wp 343 ij0 = 88 ; ij1 = 88 344 ii0 = 159 ; ii1 = 159 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 1._wp 345 ! 393 346 ! ! Danish straits : strong slip (fmask > 2) 394 347 ! We keep this as an example but it is instable in this case … … 413 366 IF(lwp) WRITE(numout,*) ' Gibraltar ' 414 367 ii0 = 282 ; ii1 = 283 ! Gibraltar Strait 415 ij0 = 2 01 +isrow ; ij1 = 241 - isrow ; fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp368 ij0 = 241 - isrow ; ij1 = 241 - isrow ; fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp 416 369 417 370 IF(lwp) WRITE(numout,*) ' Bhosporus ' 418 371 ii0 = 314 ; ii1 = 315 ! Bhosporus Strait 419 ij0 = 2 08 +isrow ; ij1 = 248 - isrow ; fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp372 ij0 = 248 - isrow ; ij1 = 248 - isrow ; fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp 420 373 421 374 IF(lwp) WRITE(numout,*) ' Makassar (Top) ' 422 375 ii0 = 48 ; ii1 = 48 ! Makassar Strait (Top) 423 ij0 = 1 49 +isrow ; ij1 = 190 - isrow ; fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp376 ij0 = 189 - isrow ; ij1 = 190 - isrow ; fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp 424 377 425 378 IF(lwp) WRITE(numout,*) ' Lombok ' 426 379 ii0 = 44 ; ii1 = 44 ! Lombok Strait 427 ij0 = 1 24 +isrow ; ij1 = 165 - isrow ; fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp380 ij0 = 164 - isrow ; ij1 = 165 - isrow ; fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp 428 381 429 382 IF(lwp) WRITE(numout,*) ' Ombai ' 430 383 ii0 = 53 ; ii1 = 53 ! Ombai Strait 431 ij0 = 1 24 +isrow ; ij1 = 165 - isrow ; fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp384 ij0 = 164 - isrow ; ij1 = 165 - isrow ; fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp 432 385 433 386 IF(lwp) WRITE(numout,*) ' Timor Passage ' 434 387 ii0 = 56 ; ii1 = 56 ! Timor Passage 435 ij0 = 1 24 +isrow ; ij1 = 165 - isrow ; fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp388 ij0 = 164 - isrow ; ij1 = 165 - isrow ; fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp 436 389 437 390 IF(lwp) WRITE(numout,*) ' West Halmahera ' 438 391 ii0 = 58 ; ii1 = 58 ! West Halmahera Strait 439 ij0 = 1 41 +isrow ; ij1 = 182 - isrow ; fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp392 ij0 = 181 - isrow ; ij1 = 182 - isrow ; fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp 440 393 441 394 IF(lwp) WRITE(numout,*) ' East Halmahera ' 442 395 ii0 = 55 ; ii1 = 55 ! East Halmahera Strait 443 ij0 = 1 41 +isrow ; ij1 = 182 - isrow ; fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp396 ij0 = 181 - isrow ; ij1 = 182 - isrow ; fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp 444 397 ! 445 398 ENDIF … … 500 453 ! 501 454 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 455 696 456 !!====================================================================== -
branches/UKMO/icebergs_restart_single_file/NEMOGCM/NEMO/OPA_SRC/DOM/domngb.F90
r6019 r6020 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') -
branches/UKMO/icebergs_restart_single_file/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90
r6019 r6020 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 !!====================================================================== -
branches/UKMO/icebergs_restart_single_file/NEMOGCM/NEMO/OPA_SRC/DOM/domwri.F90
r6019 r6020 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 … … 215 279 CALL iom_rstput( 0, 0, inum4, 'gdept_1d' , gdept_1d ) ! ! stretched system 216 280 CALL iom_rstput( 0, 0, inum4, 'gdepw_1d' , gdepw_1d ) 281 CALL iom_rstput( 0, 0, inum4, 'gdept_0', gdept_0, ktype = jp_r4 ) 282 CALL iom_rstput( 0, 0, inum4, 'gdepw_0', gdepw_0, ktype = jp_r4 ) 217 283 ENDIF 218 284 -
branches/UKMO/icebergs_restart_single_file/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
r6019 r6020 219 219 & ppsur == pp_to_be_computed ) THEN 220 220 ! 221 #if defined key_agrif 222 za1 = ( ppdzmin - pphmax / FLOAT(jpkdta-1) ) & 223 & / ( TANH((1-ppkth)/ppacr) - ppacr/FLOAT(jpkdta-1) * ( LOG( COSH( (jpkdta - ppkth) / ppacr) )& 224 & - LOG( COSH( ( 1 - ppkth) / ppacr) ) ) ) 225 #else 221 226 za1 = ( ppdzmin - pphmax / FLOAT(jpkm1) ) & 222 227 & / ( TANH((1-ppkth)/ppacr) - ppacr/FLOAT(jpk-1) * ( LOG( COSH( (jpk - ppkth) / ppacr) ) & 223 228 & - LOG( COSH( ( 1 - ppkth) / ppacr) ) ) ) 229 #endif 224 230 za0 = ppdzmin - za1 * TANH( (1-ppkth) / ppacr ) 225 231 zsur = - za0 - za1 * ppacr * LOG( COSH( (1-ppkth) / ppacr ) ) … … 236 242 WRITE(numout,*) ' Uniform grid with ',jpk-1,' layers' 237 243 WRITE(numout,*) ' Total depth :', zhmax 244 #if defined key_agrif 245 WRITE(numout,*) ' Layer thickness:', zhmax/(jpkdta-1) 246 #else 238 247 WRITE(numout,*) ' Layer thickness:', zhmax/(jpk-1) 248 #endif 239 249 ELSE 240 250 IF( ppa1 == 0._wp .AND. ppa0 == 0._wp .AND. ppsur == 0._wp ) THEN … … 260 270 ! Reference z-coordinate (depth - scale factor at T- and W-points) 261 271 ! ====================== 262 IF( ppkth == 0._wp ) THEN ! uniform vertical grid 272 IF( ppkth == 0._wp ) THEN ! uniform vertical grid 273 #if defined key_agrif 274 za1 = zhmax / FLOAT(jpkdta-1) 275 #else 263 276 za1 = zhmax / FLOAT(jpk-1) 277 #endif 264 278 DO jk = 1, jpk 265 279 zw = FLOAT( jk ) … … 487 501 IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! ORCA R2 configuration 488 502 ! ! ===================== 489 IF( nn_cla == 0 ) THEN 490 ii0 = 140 ; ii1 = 140 ! Gibraltar Strait open 491 ij0 = 102 ; ij1 = 102 ! (Thomson, Ocean Modelling, 1995) 492 DO ji = mi0(ii0), mi1(ii1) 493 DO jj = mj0(ij0), mj1(ij1) 494 mbathy(ji,jj) = 15 495 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 496 509 END DO 497 IF(lwp) WRITE(numout,*)498 IF(lwp) WRITE(numout,*) ' orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0499 !500 ii0 = 160 ; ii1 = 160 ! Bab el mandeb Strait open501 ij0 = 88 ; ij1 = 88 ! (Thomson, Ocean Modelling, 1995)502 DO ji = mi0(ii0), mi1(ii1)503 DO jj = mj0(ij0), mj1(ij1)504 mbathy(ji,jj) = 12505 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 506 519 END DO 507 IF(lwp) WRITE(numout,*)508 IF(lwp) WRITE(numout,*) ' orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0509 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 510 523 ! 511 524 ENDIF … … 531 544 ! 532 545 IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! ORCA R2 configuration 533 ! 534 IF( nn_cla == 0 ) THEN 535 ii0 = 140 ; ii1 = 140 ! Gibraltar Strait open 536 ij0 = 102 ; ij1 = 102 ! (Thomson, Ocean Modelling, 1995) 537 DO ji = mi0(ii0), mi1(ii1) 538 DO jj = mj0(ij0), mj1(ij1) 539 bathy(ji,jj) = 284._wp 540 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 541 552 END DO 542 IF(lwp) WRITE(numout,*)543 IF(lwp) WRITE(numout,*) ' orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0544 !545 ii0 = 160 ; ii1 = 160 ! Bab el mandeb Strait open546 ij0 = 88 ; ij1 = 88 ! (Thomson, Ocean Modelling, 1995)547 DO ji = mi0(ii0), mi1(ii1)548 DO jj = mj0(ij0), mj1(ij1)549 bathy(ji,jj) = 137._wp550 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 551 562 END DO 552 IF(lwp) WRITE(numout,*)553 IF(lwp) WRITE(numout,*) ' orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0554 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 555 566 ! 556 567 ENDIF -
branches/UKMO/icebergs_restart_single_file/NEMOGCM/NEMO/OPA_SRC/DOM/dtatsd.F90
r6019 r6020 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 ! -
branches/UKMO/icebergs_restart_single_file/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90
r6019 r6020 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 USE dynspg_oce ! pressure gradient schemes40 USE dynspg_flt ! filtered free surface41 USE sol_oce ! ocean solver variables42 37 ! 43 38 USE in_out_manager ! I/O manager … … 76 71 ! 77 72 78 IF(lwp) WRITE(numout,*) ' '73 IF(lwp) WRITE(numout,*) 79 74 IF(lwp) WRITE(numout,*) 'istate_ini : Initialization of the dynamics and tracers' 80 75 IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 81 76 82 CALL dta_tsd_init! Initialisation of T & S input data83 IF( lk_c1d ) CALL dta_uvd_init! Initialization of U & V input data77 CALL dta_tsd_init ! Initialisation of T & S input data 78 IF( lk_c1d ) CALL dta_uvd_init ! Initialization of U & V input data 84 79 85 80 rhd (:,:,: ) = 0._wp ; rhop (:,:,: ) = 0._wp ! set one for all to 0 at level jpk … … 103 98 ub (:,:,:) = 0._wp ; un (:,:,:) = 0._wp 104 99 vb (:,:,:) = 0._wp ; vn (:,:,:) = 0._wp 105 rotb (:,:,:) = 0._wp ; rotn (:,:,:) = 0._wp 106 hdivb(:,:,:) = 0._wp ; hdivn(:,:,:) = 0._wp 100 hdivn(:,:,:) = 0._wp 107 101 ! 108 102 IF( cp_cfg == 'eel' ) THEN … … 119 113 ENDIF 120 114 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 )115 CALL wrk_alloc( jpi,jpj,jpk,2, zuvd ) 122 116 CALL dta_uvd( nit000, zuvd ) 123 117 ub(:,:,:) = zuvd(:,:,:,1) ; un(:,:,:) = ub(:,:,:) 124 118 vb(:,:,:) = zuvd(:,:,:,2) ; vn(:,:,:) = vb(:,:,:) 125 CALL wrk_dealloc( jpi, jpj, jpk, 2,zuvd )119 CALL wrk_dealloc( jpi,jpj,jpk,2, zuvd ) 126 120 ENDIF 127 121 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 122 ! 140 123 ! - ML - sshn could be modified by istate_eel, so that initialization of fse3t_b is done here … … 142 125 DO jk = 1, jpk 143 126 fse3t_b(:,:,jk) = fse3t_n(:,:,jk) 144 END DO127 END DO 145 128 ENDIF 146 129 ! 147 ENDIF148 !149 IF( lk_agrif ) THEN ! read free surface arrays in restart file150 IF( ln_rstart ) THEN151 IF( lk_dynspg_flt ) THEN ! read or initialize the following fields152 ! ! gcx, gcxb for agrif_opa_init153 IF( sol_oce_alloc() > 0 ) CALL ctl_stop('agrif sol_oce_alloc: allocation of arrays failed')154 CALL flt_rst( nit000, 'READ' )155 ENDIF156 ENDIF ! explicit case not coded yet with AGRIF157 130 ENDIF 158 131 ! … … 163 136 ! 164 137 ! 165 un_b(:,:) = 0._wp ;vn_b(:,:) = 0._wp166 ub_b(:,:) = 0._wp ;vb_b(:,:) = 0._wp138 un_b(:,:) = 0._wp ; vn_b(:,:) = 0._wp 139 ub_b(:,:) = 0._wp ; vb_b(:,:) = 0._wp 167 140 ! 168 141 DO jk = 1, jpkm1 … … 202 175 !! References : Philander ??? 203 176 !!---------------------------------------------------------------------- 204 INTEGER :: ji, jj, jk205 REAL(wp) :: zsal = 35.50 177 INTEGER :: ji, jj, jk 178 REAL(wp) :: zsal = 35.50_wp 206 179 !!---------------------------------------------------------------------- 207 180 ! … … 233 206 !! and relative vorticity fields 234 207 !!---------------------------------------------------------------------- 235 USE div cur ! hor. divergence & rel. vorticity (div_cur routine)208 USE divhor ! hor. divergence (div_hor routine) 236 209 USE iom 237 210 ! … … 282 255 tsb(:,:,:,jp_sal) = tsn(:,:,:,jp_sal) 283 256 ! 284 ! set the dynamics: U,V, hdiv , rot(and ssh if necessary)257 ! set the dynamics: U,V, hdiv (and ssh if necessary) 285 258 ! ---------------- 286 259 ! Start EEL5 configuration with barotropic geostrophic velocities … … 318 291 ENDIF 319 292 ! 320 CALL div_cur( nit000 ) ! horizontal divergence and relative vorticity (curl) 293 !!gm Check here call to div_hor should not be necessary 294 !!gm div_hor call runoffs not sure they are defined at that level 295 CALL div_hor( nit000 ) ! horizontal divergence and relative vorticity (curl) 321 296 ! N.B. the vertical velocity will be computed from the horizontal divergence field 322 297 ! in istate by a call to wzv routine … … 371 346 !! 372 347 !! ** Method : - set temprature field 373 !! - set salinity field348 !! - set salinity field 374 349 !!---------------------------------------------------------------------- 375 350 INTEGER :: ji, jj, jk ! dummy loop indices … … 457 432 !! p=integral [ rau*g dz ] 458 433 !!---------------------------------------------------------------------- 459 USE dynspg ! surface pressure gradient (dyn_spg routine) 460 USE divcur ! hor. divergence & rel. vorticity (div_cur routine) 434 USE divhor ! hor. divergence (div_hor routine) 461 435 USE lbclnk ! ocean lateral boundary condition (or mpp link) 462 436 ! 463 437 INTEGER :: ji, jj, jk ! dummy loop indices 464 INTEGER :: indic ! ???465 438 REAL(wp) :: zmsv, zphv, zmsu, zphu, zalfg ! temporary scalars 466 439 REAL(wp), POINTER, DIMENSION(:,:,:) :: zprn 467 440 !!---------------------------------------------------------------------- 468 441 ! 469 CALL wrk_alloc( jpi, jpj, jpk,zprn)442 CALL wrk_alloc( jpi,jpj,jpk, zprn) 470 443 ! 471 444 IF(lwp) WRITE(numout,*) … … 529 502 vb(:,:,:) = vn(:,:,:) 530 503 531 ! WARNING !!!!! 532 ! after initializing u and v, we need to calculate the initial streamfunction bsf. 533 ! Otherwise, only the trend will be computed and the model will blow up (inconsistency). 534 ! to do that, we call dyn_spg with a special trick: 535 ! we fill ua and va with the velocities divided by dt, and the streamfunction will be brought to the 536 ! right value assuming the velocities have been set up in one time step. 537 ! we then set bsfd to zero (first guess for next step is d(psi)/dt = 0.) 538 ! sets up s false trend to calculate the barotropic streamfunction. 539 540 ua(:,:,:) = ub(:,:,:) / rdt 541 va(:,:,:) = vb(:,:,:) / rdt 542 543 ! calls dyn_spg. we assume euler time step, starting from rest. 544 indic = 0 545 CALL dyn_spg( nit000, indic ) ! surface pressure gradient 546 547 ! the new velocity is ua*rdt 548 549 CALL lbc_lnk( ua, 'U', -1. ) 550 CALL lbc_lnk( va, 'V', -1. ) 551 552 ub(:,:,:) = ua(:,:,:) * rdt 553 vb(:,:,:) = va(:,:,:) * rdt 554 ua(:,:,:) = 0.e0 555 va(:,:,:) = 0.e0 556 un(:,:,:) = ub(:,:,:) 557 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) 504 ! 505 !!gm Check here call to div_hor should not be necessary 506 !!gm div_hor call runoffs not sure they are defined at that level 507 CALL div_hor( nit000 ) ! now horizontal divergence 508 ! 509 CALL wrk_dealloc( jpi,jpj,jpk, zprn) 567 510 ! 568 511 END SUBROUTINE istate_uvg
Note: See TracChangeset
for help on using the changeset viewer.