Changeset 7397
- Timestamp:
- 2016-11-30T15:13:05+01:00 (7 years ago)
- Location:
- branches/2016/dev_NOC_CMCC_merge_2016/NEMOGCM
- Files:
-
- 16 added
- 23 edited
- 1 copied
-
ARCH/arch-THALASSA.fcm (added)
-
CONFIG/SHARED/namelist_ref (modified) (2 diffs)
-
CONFIG/WAD_TEST_CASES (added)
-
CONFIG/WAD_TEST_CASES/EXP00 (added)
-
CONFIG/WAD_TEST_CASES/EXP00/README.py (added)
-
CONFIG/WAD_TEST_CASES/EXP00/domain_def.xml (added)
-
CONFIG/WAD_TEST_CASES/EXP00/field_def.xml (added)
-
CONFIG/WAD_TEST_CASES/EXP00/iodef.xml (added)
-
CONFIG/WAD_TEST_CASES/EXP00/matpoly2.py (added)
-
CONFIG/WAD_TEST_CASES/EXP00/namelist_cfg (added)
-
CONFIG/WAD_TEST_CASES/EXP00/namelist_ref (added)
-
CONFIG/WAD_TEST_CASES/MY_SRC (added)
-
CONFIG/WAD_TEST_CASES/MY_SRC/domwri.F90 (added)
-
CONFIG/WAD_TEST_CASES/MY_SRC/dynadv_cen2.F90 (added)
-
CONFIG/WAD_TEST_CASES/MY_SRC/dynadv_ubs.F90 (added)
-
CONFIG/WAD_TEST_CASES/MY_SRC/dynzdf_imp.F90 (added)
-
CONFIG/WAD_TEST_CASES/MY_SRC/iom.F90 (copied) (copied from branches/2016/dev_NOC_CMCC_merge_2016/NEMOGCM/NEMO/OPA_SRC/IOM/iom.F90) (3 diffs)
-
CONFIG/WAD_TEST_CASES/cpp_WAD_TEST_CASES.fcm (added)
-
CONFIG/cfg.txt (modified) (1 diff)
-
NEMO/OPA_SRC/BDY/bdy_oce.F90 (modified) (3 diffs)
-
NEMO/OPA_SRC/BDY/bdydta.F90 (modified) (3 diffs)
-
NEMO/OPA_SRC/BDY/bdydyn3d.F90 (modified) (3 diffs)
-
NEMO/OPA_SRC/BDY/bdyini.F90 (modified) (6 diffs)
-
NEMO/OPA_SRC/BDY/bdylib.F90 (modified) (2 diffs)
-
NEMO/OPA_SRC/BDY/bdytra.F90 (modified) (2 diffs)
-
NEMO/OPA_SRC/DOM/dommsk.F90 (modified) (4 diffs)
-
NEMO/OPA_SRC/DOM/domvvl.F90 (modified) (2 diffs)
-
NEMO/OPA_SRC/DOM/domzgr.F90 (modified) (4 diffs)
-
NEMO/OPA_SRC/DOM/istate.F90 (modified) (2 diffs)
-
NEMO/OPA_SRC/DYN/dynhpg.F90 (modified) (17 diffs)
-
NEMO/OPA_SRC/DYN/dynkeg.F90 (modified) (4 diffs)
-
NEMO/OPA_SRC/DYN/dynspg_ts.F90 (modified) (13 diffs)
-
NEMO/OPA_SRC/DYN/sshwzv.F90 (modified) (2 diffs)
-
NEMO/OPA_SRC/DYN/wet_dry.F90 (modified) (18 diffs)
-
NEMO/OPA_SRC/IOM/iom.F90 (modified) (4 diffs)
-
NEMO/OPA_SRC/IOM/iom_nf90.F90 (modified) (5 diffs)
-
NEMO/OPA_SRC/LBC/mppini_2.h90 (modified) (3 diffs)
-
NEMO/OPA_SRC/LDF/ldfdyn.F90 (modified) (9 diffs)
-
NEMO/OPA_SRC/SBC/fldread.F90 (modified) (17 diffs)
-
NEMO/TOP_SRC/trcbdy.F90 (modified) (3 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/2016/dev_NOC_CMCC_merge_2016/NEMOGCM/CONFIG/SHARED/namelist_ref
r7299 r7397 653 653 ln_vol = .false. ! total volume correction (see nn_volctl parameter) 654 654 nn_volctl = 1 ! = 0, the total water flux across open boundaries is zero 655 nb_jpk_bdy = -1 ! number of levels in the bdy data (set < 0 if consistent with planned run) 655 656 / 656 657 !----------------------------------------------------------------------- … … 942 943 ! ! = 30 F(i,j,k)=c2d*c1d 943 944 ! ! = 31 F(i,j,k)=F(grid spacing and local velocity) 945 ! ! = 32 F(i,j,k)=F(local gridscale and deformation rate) 946 ! Caution in 20 and 30 cases the coefficient have to be given for a 1 degree grid (~111km) 944 947 rn_ahm_0 = 40000. ! horizontal laplacian eddy viscosity [m2/s] 945 948 rn_ahm_b = 0. ! background eddy viscosity for ldf_iso [m2/s] 946 949 rn_bhm_0 = 1.e+12 ! horizontal bilaplacian eddy viscosity [m4/s] 947 ! 948 ! Caution in 20 and 30 cases the coefficient have to be given for a 1 degree grid (~111km) 950 ! ! Smagorinsky settings (nn_ahm_ijk_t = 32) : 951 rn_csmc = 3.5 ! Smagorinsky constant of proportionality 952 rn_minfac = 1.0 ! multiplier of theorectical lower limit 953 rn_maxfac = 1.0 ! multiplier of theorectical upper limit 949 954 / 950 955 -
branches/2016/dev_NOC_CMCC_merge_2016/NEMOGCM/CONFIG/WAD_TEST_CASES/MY_SRC/iom.F90
r7385 r7397 78 78 !!---------------------------------------------------------------------- 79 79 !! NEMO/OPA 3.3 , NEMO Consortium (2010) 80 !! $Id $80 !! $Id: iom.F90 6140 2015-12-21 11:35:23Z timgraham $ 81 81 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 82 82 !!---------------------------------------------------------------------- … … 114 114 CASE (30) ; CALL xios_set_context_attr(TRIM(clname), calendar_type= "D360") 115 115 END SELECT 116 WRITE(cldate,"(i4.4,'-',i2.2,'-',i2.2,' ',i2.2,':',i2.2,':00')") nyear,nmonth,nday,nhour,nminute116 WRITE(cldate,"(i4.4,'-',i2.2,'-',i2.2,' 00:00:00')") nyear,nmonth,nday 117 117 CALL xios_set_context_attr(TRIM(clname), start_date=cldate ) 118 118 … … 172 172 z_bnds(1:jpkm1,2) = gdepw_1d(2:jpk) 173 173 z_bnds(jpk: ,2) = gdepw_1d(jpk) + e3t_1d(jpk) 174 CALL iom_set_axis_attr( "deptht", bounds=z_bnds )175 CALL iom_set_axis_attr( "depthu", bounds=z_bnds )176 CALL iom_set_axis_attr( "depthv", bounds=z_bnds )174 !CALL iom_set_axis_attr( "deptht", bounds=z_bnds ) 175 !CALL iom_set_axis_attr( "depthu", bounds=z_bnds ) 176 !CALL iom_set_axis_attr( "depthv", bounds=z_bnds ) 177 177 z_bnds(: ,2) = gdept_1d(:) 178 178 z_bnds(2:jpk,1) = gdept_1d(1:jpkm1) 179 179 z_bnds(1 ,1) = gdept_1d(1) - e3w_1d(1) 180 CALL iom_set_axis_attr( "depthw", bounds=z_bnds )180 !CALL iom_set_axis_attr( "depthw", bounds=z_bnds ) 181 181 182 182 # if defined key_floats -
branches/2016/dev_NOC_CMCC_merge_2016/NEMOGCM/CONFIG/cfg.txt
r6403 r7397 11 11 GYRE OPA_SRC 12 12 ORCA2_LIM_PISCES OPA_SRC LIM_SRC_2 NST_SRC TOP_SRC 13 WAD_TEST_CASES OPA_SRC -
branches/2016/dev_NOC_CMCC_merge_2016/NEMOGCM/NEMO/OPA_SRC/BDY/bdy_oce.F90
r7299 r7397 47 47 LOGICAL :: ll_tem 48 48 LOGICAL :: ll_sal 49 LOGICAL :: ll_fvl 49 50 REAL(wp), POINTER, DIMENSION(:) :: ssh 50 51 REAL(wp), POINTER, DIMENSION(:) :: u2d … … 91 92 ! 92 93 INTEGER :: nb_bdy !: number of open boundary sets 94 INTEGER :: nb_jpk_bdy !: number of levels in the bdy data (set < 0 if consistent with planned run) 93 95 INTEGER, DIMENSION(jp_bdy) :: nn_rimwidth !: boundary rim width for Flow Relaxation Scheme 94 96 INTEGER :: nn_volctl !: = 0 the total volume will have the variability of the surface Flux E-P … … 134 136 !: =1 => some data to be read in from data files 135 137 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET :: dta_global !: workspace for reading in global data arrays (unstr. bdy) 138 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET :: dta_global_z !: workspace for reading in global depth arrays (unstr. bdy) 139 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET :: dta_global_dz !: workspace for reading in global depth arrays (unstr. bdy) 136 140 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET :: dta_global2 !: workspace for reading in global data arrays (struct. bdy) 141 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET :: dta_global2_z !: workspace for reading in global depth arrays (struct. bdy) 142 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET :: dta_global2_dz !: workspace for reading in global depth arrays (struct. bdy) 137 143 !$AGRIF_DO_NOT_TREAT 138 144 TYPE(OBC_INDEX), DIMENSION(jp_bdy), TARGET :: idx_bdy !: bdy indices (local process) -
branches/2016/dev_NOC_CMCC_merge_2016/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90
r7299 r7397 264 264 265 265 jend = jstart + dta%nread(2) - 1 266 CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), map=nbmap_ptr(jstart:jend), & 267 & kit=jit, kt_offset=time_offset ) 266 IF( ln_full_vel_array(ib_bdy) ) THEN 267 CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), map=nbmap_ptr(jstart:jend), & 268 & kit=jit, kt_offset=time_offset , jpk_bdy=nb_jpk_bdy, fvl=ln_full_vel_array(ib_bdy) ) 269 ELSE 270 CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), map=nbmap_ptr(jstart:jend), & 271 & kit=jit, kt_offset=time_offset ) 272 ENDIF 268 273 269 274 ! If full velocities in boundary data then extract barotropic velocities from 3D fields … … 330 335 jend = jstart + dta%nread(1) - 1 331 336 CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), & 332 & map=nbmap_ptr(jstart:jend), kt_offset=time_offset )337 & map=nbmap_ptr(jstart:jend), kt_offset=time_offset, jpk_bdy=nb_jpk_bdy, fvl=ln_full_vel_array(ib_bdy) ) 333 338 ENDIF 334 339 ! If full velocities in boundary data then split into barotropic and baroclinic data … … 456 461 NAMELIST/nambdy_dta/ bn_a_i, bn_ht_i, bn_ht_s 457 462 #endif 458 NAMELIST/nambdy_dta/ ln_full_vel 463 NAMELIST/nambdy_dta/ ln_full_vel, nb_jpk_bdy 459 464 !!--------------------------------------------------------------------------- 460 465 ! -
branches/2016/dev_NOC_CMCC_merge_2016/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn3d.F90
r7299 r7397 53 53 CASE('orlanski' ) ; CALL bdy_dyn3d_orlanski( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy, ll_npo=.false. ) 54 54 CASE('orlanski_npo'); CALL bdy_dyn3d_orlanski( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy, ll_npo=.true. ) 55 CASE('zerograd') ; CALL bdy_dyn3d_zgrad( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt, ib_bdy ) 56 CASE('neumann') ; CALL bdy_dyn3d_nmn( idx_bdy(ib_bdy), ib_bdy ) 55 57 CASE DEFAULT ; CALL ctl_stop( 'bdy_dyn3d : unrecognised option for open boundaries for baroclinic velocities' ) 56 58 END SELECT … … 106 108 END SUBROUTINE bdy_dyn3d_spe 107 109 110 SUBROUTINE bdy_dyn3d_zgrad( idx, dta, kt , ib_bdy ) 111 !!---------------------------------------------------------------------- 112 !! *** SUBROUTINE bdy_dyn3d_zgrad *** 113 !! 114 !! ** Purpose : - Enforce a zero gradient of normal velocity 115 !! 116 !!---------------------------------------------------------------------- 117 INTEGER :: kt 118 TYPE(OBC_INDEX), INTENT(in) :: idx ! OBC indices 119 TYPE(OBC_DATA), INTENT(in) :: dta ! OBC external data 120 INTEGER, INTENT(in) :: ib_bdy ! BDY set index 121 !! 122 INTEGER :: jb, jk ! dummy loop indices 123 INTEGER :: ii, ij, igrd ! local integers 124 REAL(wp) :: zwgt ! boundary weight 125 INTEGER :: fu, fv 126 !!---------------------------------------------------------------------- 127 ! 128 IF( nn_timing == 1 ) CALL timing_start('bdy_dyn3d_zgrad') 129 ! 130 igrd = 2 ! Copying tangential velocity into bdy points 131 DO jb = 1, idx%nblenrim(igrd) 132 DO jk = 1, jpkm1 133 ii = idx%nbi(jb,igrd) 134 ij = idx%nbj(jb,igrd) 135 fu = ABS( ABS (NINT( idx%flagu(jb,igrd) ) ) - 1 ) 136 ua(ii,ij,jk) = ua(ii,ij,jk) * REAL( 1 - fu) + ( ua(ii,ij+fu,jk) * umask(ii,ij+fu,jk) & 137 &+ ua(ii,ij-fu,jk) * umask(ii,ij-fu,jk) ) * umask(ii,ij,jk) * REAL( fu ) 138 END DO 139 END DO 140 ! 141 igrd = 3 ! Copying tangential velocity into bdy points 142 DO jb = 1, idx%nblenrim(igrd) 143 DO jk = 1, jpkm1 144 ii = idx%nbi(jb,igrd) 145 ij = idx%nbj(jb,igrd) 146 fv = ABS( ABS (NINT( idx%flagv(jb,igrd) ) ) - 1 ) 147 va(ii,ij,jk) = va(ii,ij,jk) * REAL( 1 - fv ) + ( va(ii+fv,ij,jk) * vmask(ii+fv,ij,jk) & 148 &+ va(ii-fv,ij,jk) * vmask(ii-fv,ij,jk) ) * vmask(ii,ij,jk) * REAL( fv ) 149 END DO 150 END DO 151 CALL lbc_bdy_lnk( ua, 'U', -1., ib_bdy ) ! Boundary points should be updated 152 CALL lbc_bdy_lnk( va, 'V', -1., ib_bdy ) 153 ! 154 IF( kt .eq. nit000 ) CLOSE( unit = 102 ) 155 156 IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn3d_zgrad') 157 158 END SUBROUTINE bdy_dyn3d_zgrad 108 159 109 160 SUBROUTINE bdy_dyn3d_zro( idx, dta, kt, ib_bdy ) … … 292 343 END SUBROUTINE bdy_dyn3d_dmp 293 344 345 SUBROUTINE bdy_dyn3d_nmn( idx, ib_bdy ) 346 !!---------------------------------------------------------------------- 347 !! *** SUBROUTINE bdy_dyn3d_nmn *** 348 !! 349 !! - Apply Neumann condition to baroclinic velocities. 350 !! - Wrapper routine for bdy_nmn 351 !! 352 !! 353 !!---------------------------------------------------------------------- 354 TYPE(OBC_INDEX), INTENT(in) :: idx ! OBC indices 355 INTEGER, INTENT(in) :: ib_bdy ! BDY set index 356 357 INTEGER :: jb, igrd ! dummy loop indices 358 !!---------------------------------------------------------------------- 359 360 IF( nn_timing == 1 ) CALL timing_start('bdy_dyn3d_nmn') 361 ! 362 !! Note that at this stage the ub and ua arrays contain the baroclinic velocities. 363 ! 364 igrd = 2 ! Neumann bc on u-velocity; 365 ! 366 CALL bdy_nmn( idx, igrd, ua ) 367 368 igrd = 3 ! Neumann bc on v-velocity 369 ! 370 CALL bdy_nmn( idx, igrd, va ) 371 ! 372 CALL lbc_bdy_lnk( ua, 'U', -1., ib_bdy ) ! Boundary points should be updated 373 CALL lbc_bdy_lnk( va, 'V', -1., ib_bdy ) 374 ! 375 IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn3d_nmn') 376 ! 377 END SUBROUTINE bdy_dyn3d_nmn 378 294 379 !!====================================================================== 295 380 END MODULE bdydyn3d -
branches/2016/dev_NOC_CMCC_merge_2016/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90
r7299 r7397 71 71 & cn_ice_lim, nn_ice_lim_dta, & 72 72 & rn_ice_tem, rn_ice_sal, rn_ice_age, & 73 & ln_vol, nn_volctl, nn_rimwidth 73 & ln_vol, nn_volctl, nn_rimwidth, nb_jpk_bdy 74 74 ! 75 75 INTEGER :: ios ! Local integer output status for namelist read … … 244 244 dta_bdy(ib_bdy)%ll_u3d = .true. 245 245 dta_bdy(ib_bdy)%ll_v3d = .true. 246 CASE('neumann') 247 IF(lwp) WRITE(numout,*) ' Neumann conditions' 248 dta_bdy(ib_bdy)%ll_u3d = .false. 249 dta_bdy(ib_bdy)%ll_v3d = .false. 250 CASE('zerograd') 251 IF(lwp) WRITE(numout,*) ' Zero gradient for baroclinic velocities' 252 dta_bdy(ib_bdy)%ll_u3d = .false. 253 dta_bdy(ib_bdy)%ll_v3d = .false. 246 254 CASE('zero') 247 255 IF(lwp) WRITE(numout,*) ' Zero baroclinic velocities (runoff case)' … … 412 420 IF(lwp) WRITE(numout,*) 'No volume correction applied at open boundaries' 413 421 IF(lwp) WRITE(numout,*) 422 ENDIF 423 IF( nb_jpk_bdy > 0 ) THEN 424 IF(lwp) WRITE(numout,*) '*** open boundary will be interpolate in the vertical onto the native grid ***' 425 ELSE 426 IF(lwp) WRITE(numout,*) '*** open boundary will be read straight onto the native grid without vertical interpolation ***' 414 427 ENDIF 415 428 ENDIF … … 534 547 & nbrdta(jpbdta, jpbgrd, nb_bdy) ) 535 548 536 ALLOCATE( dta_global(jpbdtau, 1, jpk) ) 537 IF ( icount>0 ) ALLOCATE( dta_global2(jpbdtas, nrimmax, jpk) ) 549 IF( nb_jpk_bdy>0 ) THEN 550 ALLOCATE( dta_global(jpbdtau, 1, nb_jpk_bdy) ) 551 ALLOCATE( dta_global_z(jpbdtau, 1, nb_jpk_bdy) ) 552 ALLOCATE( dta_global_dz(jpbdtau, 1, nb_jpk_bdy) ) 553 ELSE 554 ALLOCATE( dta_global(jpbdtau, 1, jpk) ) 555 ALLOCATE( dta_global_z(jpbdtau, 1, jpk) ) ! needed ?? TODO 556 ALLOCATE( dta_global_dz(jpbdtau, 1, jpk) )! needed ?? TODO 557 ENDIF 558 559 IF ( icount>0 ) THEN 560 IF( nb_jpk_bdy>0 ) THEN 561 ALLOCATE( dta_global2(jpbdtas, nrimmax, nb_jpk_bdy) ) 562 ALLOCATE( dta_global2_z(jpbdtas, nrimmax, nb_jpk_bdy) ) 563 ALLOCATE( dta_global2_dz(jpbdtas, nrimmax, nb_jpk_bdy) ) 564 ELSE 565 ALLOCATE( dta_global2(jpbdtas, nrimmax, jpk) ) 566 ALLOCATE( dta_global2_z(jpbdtas, nrimmax, jpk) ) ! needed ?? TODO 567 ALLOCATE( dta_global2_dz(jpbdtas, nrimmax, jpk) )! needed ?? TODO 568 ENDIF 569 ENDIF 538 570 ! 539 571 ENDIF … … 1127 1159 ! = 0 elsewhere 1128 1160 1161 bdytmask(:,:) = ssmask(:,:) 1162 1129 1163 IF( ln_mask_file ) THEN 1130 1164 CALL iom_open( cn_mask_file, inum ) … … 1143 1177 CALL lbc_lnk( bdyumask(:,:), 'U', 1. ) ; CALL lbc_lnk( bdyvmask(:,:), 'V', 1. ) ! Lateral boundary cond. 1144 1178 1145 1146 ! Mask corrections1147 ! ----------------1148 DO ik = 1, jpkm11149 DO ij = 1, jpj1150 DO ii = 1, jpi1151 tmask(ii,ij,ik) = tmask(ii,ij,ik) * bdytmask(ii,ij)1152 umask(ii,ij,ik) = umask(ii,ij,ik) * bdyumask(ii,ij)1153 vmask(ii,ij,ik) = vmask(ii,ij,ik) * bdyvmask(ii,ij)1154 END DO1155 END DO1156 DO ij = 2, jpjm11157 DO ii = 2, jpim11158 fmask(ii,ij,ik) = fmask(ii,ij,ik) * bdytmask(ii,ij ) * bdytmask(ii+1,ij ) &1159 & * bdytmask(ii,ij+1) * bdytmask(ii+1,ij+1)1160 END DO1161 END DO1162 END DO1163 tmask_i (:,:) = ssmask(:,:) * tmask_i(:,:)1164 !1165 1179 ENDIF ! ln_mask_file=.TRUE. 1166 1180 1167 bdytmask(:,:) = ssmask(:,:)1168 1181 IF( .NOT.ln_mask_file ) THEN 1169 1182 ! If .not. ln_mask_file then we need to derive mask on U and V grid from mask on T grid here. -
branches/2016/dev_NOC_CMCC_merge_2016/NEMOGCM/NEMO/OPA_SRC/BDY/bdylib.F90
r7299 r7397 97 97 ! 98 98 END SUBROUTINE bdy_spe 99 100 SUBROUTINE bdy_nmn( idx, pta )101 !!----------------------------------------------------------------------102 !! *** SUBROUTINE bdy_nmn ***103 !!104 !! ** Purpose : Duplicate the value for tracers at open boundaries.105 !!106 !!----------------------------------------------------------------------107 TYPE(OBC_INDEX), INTENT(in) :: idx ! OBC indices108 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pta ! tracer trend109 !!110 REAL(wp) :: zwgt ! boundary weight111 INTEGER :: ib, ik, igrd ! dummy loop indices112 INTEGER :: ii, ij, zcoef, zcoef1, zcoef2, ip, jp ! 2D addresses113 !!----------------------------------------------------------------------114 !115 IF( nn_timing == 1 ) CALL timing_start('bdy_nmn')116 !117 igrd = 1 ! Everything is at T-points here118 DO ib = 1, idx%nblenrim(igrd)119 ii = idx%nbi(ib,igrd)120 ij = idx%nbj(ib,igrd)121 DO ik = 1, jpkm1122 ! search the sense of the gradient123 zcoef1 = bdytmask(ii-1,ij ) + bdytmask(ii+1,ij )124 zcoef2 = bdytmask(ii ,ij-1) + bdytmask(ii ,ij+1)125 IF ( zcoef1+zcoef2 == 0) THEN126 ! corner127 zcoef = tmask(ii-1,ij,ik) + tmask(ii+1,ij,ik) + tmask(ii,ij-1,ik) + tmask(ii,ij+1,ik)128 pta(ii,ij,ik) = pta(ii-1,ij ,ik) * tmask(ii-1,ij ,ik) + &129 & pta(ii+1,ij ,ik) * tmask(ii+1,ij ,ik) + &130 & pta(ii ,ij-1,ik) * tmask(ii ,ij-1,ik) + &131 & pta(ii ,ij+1,ik) * tmask(ii ,ij+1,ik)132 pta(ii,ij,ik) = ( pta(ii,ij,ik) / MAX( 1, zcoef) ) * tmask(ii,ij,ik)133 ELSE134 ip = bdytmask(ii+1,ij ) - bdytmask(ii-1,ij )135 jp = bdytmask(ii ,ij+1) - bdytmask(ii ,ij-1)136 pta(ii,ij,ik) = pta(ii+ip,ij+jp,ik) * tmask(ii+ip,ij+jp,ik)137 ENDIF138 END DO139 END DO140 !141 IF( nn_timing == 1 ) CALL timing_stop('bdy_nmn')142 !143 END SUBROUTINE bdy_nmn144 99 145 100 SUBROUTINE bdy_orl( idx, ptb, pta, dta, ll_npo ) … … 490 445 END SUBROUTINE bdy_orlanski_3d 491 446 447 SUBROUTINE bdy_nmn( idx, igrd, phia ) 448 !!---------------------------------------------------------------------- 449 !! *** SUBROUTINE bdy_nmn *** 450 !! 451 !! ** Purpose : Duplicate the value at open boundaries, zero gradient. 452 !! 453 !!---------------------------------------------------------------------- 454 INTEGER, INTENT(in) :: igrd ! grid index 455 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: phia ! model after 3D field (to be updated) 456 TYPE(OBC_INDEX), INTENT(in) :: idx ! OBC indices 457 !! 458 REAL(wp) :: zcoef, zcoef1, zcoef2 459 REAL(wp), POINTER, DIMENSION(:,:,:) :: pmask ! land/sea mask for field 460 REAL(wp), POINTER, DIMENSION(:,:) :: bdypmask ! land/sea mask for field 461 INTEGER :: ib, ik ! dummy loop indices 462 INTEGER :: ii, ij, ip, jp ! 2D addresses 463 !!---------------------------------------------------------------------- 464 !! 465 IF( nn_timing == 1 ) CALL timing_start('bdy_nmn') 466 ! 467 SELECT CASE(igrd) 468 CASE(1) 469 pmask => tmask(:,:,:) 470 bdypmask => bdytmask(:,:) 471 CASE(2) 472 pmask => umask(:,:,:) 473 bdypmask => bdyumask(:,:) 474 CASE(3) 475 pmask => vmask(:,:,:) 476 bdypmask => bdyvmask(:,:) 477 CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for igrd in bdy_nmn' ) 478 END SELECT 479 DO ib = 1, idx%nblenrim(igrd) 480 ii = idx%nbi(ib,igrd) 481 ij = idx%nbj(ib,igrd) 482 DO ik = 1, jpkm1 483 ! search the sense of the gradient 484 zcoef1 = bdypmask(ii-1,ij )*pmask(ii-1,ij,ik) + bdypmask(ii+1,ij )*pmask(ii+1,ij,ik) 485 zcoef2 = bdypmask(ii ,ij-1)*pmask(ii,ij-1,ik) + bdypmask(ii ,ij+1)*pmask(ii,ij+1,ik) 486 IF ( nint(zcoef1+zcoef2) == 0) THEN 487 ! corner **** we probably only want to set the tangentail component for the dynamics here 488 zcoef = pmask(ii-1,ij,ik) + pmask(ii+1,ij,ik) + pmask(ii,ij-1,ik) + pmask(ii,ij+1,ik) 489 IF (zcoef > .5_wp) THEN ! Only set none isolated points. 490 phia(ii,ij,ik) = phia(ii-1,ij ,ik) * pmask(ii-1,ij ,ik) + & 491 & phia(ii+1,ij ,ik) * pmask(ii+1,ij ,ik) + & 492 & phia(ii ,ij-1,ik) * pmask(ii ,ij-1,ik) + & 493 & phia(ii ,ij+1,ik) * pmask(ii ,ij+1,ik) 494 phia(ii,ij,ik) = ( phia(ii,ij,ik) / zcoef ) * pmask(ii,ij,ik) 495 ELSE 496 phia(ii,ij,ik) = phia(ii,ij ,ik) * pmask(ii,ij ,ik) 497 ENDIF 498 ELSEIF ( nint(zcoef1+zcoef2) == 2) THEN 499 ! oblique corner **** we probably only want to set the normal component for the dynamics here 500 zcoef = pmask(ii-1,ij,ik)*bdypmask(ii-1,ij ) + pmask(ii+1,ij,ik)*bdypmask(ii+1,ij ) + & 501 & pmask(ii,ij-1,ik)*bdypmask(ii,ij -1 ) + pmask(ii,ij+1,ik)*bdypmask(ii,ij+1 ) 502 phia(ii,ij,ik) = phia(ii-1,ij ,ik) * pmask(ii-1,ij ,ik)*bdypmask(ii-1,ij ) + & 503 & phia(ii+1,ij ,ik) * pmask(ii+1,ij ,ik)*bdypmask(ii+1,ij ) + & 504 & phia(ii ,ij-1,ik) * pmask(ii ,ij-1,ik)*bdypmask(ii,ij -1 ) + & 505 & phia(ii ,ij+1,ik) * pmask(ii ,ij+1,ik)*bdypmask(ii,ij+1 ) 506 507 phia(ii,ij,ik) = ( phia(ii,ij,ik) / MAX(1._wp, zcoef) ) * pmask(ii,ij,ik) 508 ELSE 509 ip = nint(bdypmask(ii+1,ij )*pmask(ii+1,ij,ik) - bdypmask(ii-1,ij )*pmask(ii-1,ij,ik)) 510 jp = nint(bdypmask(ii ,ij+1)*pmask(ii,ij+1,ik) - bdypmask(ii ,ij-1)*pmask(ii,ij-1,ik)) 511 phia(ii,ij,ik) = phia(ii+ip,ij+jp,ik) * pmask(ii+ip,ij+jp,ik) 512 ENDIF 513 END DO 514 END DO 515 ! 516 IF( nn_timing == 1 ) CALL timing_stop('bdy_nmn') 517 ! 518 END SUBROUTINE bdy_nmn 519 492 520 !!====================================================================== 493 521 END MODULE bdylib -
branches/2016/dev_NOC_CMCC_merge_2016/NEMOGCM/NEMO/OPA_SRC/BDY/bdytra.F90
r7299 r7397 48 48 INTEGER, INTENT(in) :: kt ! Main time step counter 49 49 ! 50 INTEGER :: ib_bdy, jn ! Loop indeces51 TYPE(ztrabdy), DIMENSION(jpts) :: zdta ! Temporary data structure50 INTEGER :: ib_bdy, jn, igrd ! Loop indeces 51 TYPE(ztrabdy), DIMENSION(jpts) :: zdta ! Temporary data structure 52 52 !!---------------------------------------------------------------------- 53 igrd = 1 53 54 54 55 DO ib_bdy=1, nb_bdy … … 63 64 CASE('frs' ) ; CALL bdy_frs ( idx_bdy(ib_bdy), tsa(:,:,:,jn), zdta(jn)%tra ) 64 65 CASE('specified' ) ; CALL bdy_spe ( idx_bdy(ib_bdy), tsa(:,:,:,jn), zdta(jn)%tra ) 65 CASE('neumann' ) ; CALL bdy_nmn ( idx_bdy(ib_bdy), 66 CASE('neumann' ) ; CALL bdy_nmn ( idx_bdy(ib_bdy), igrd , tsa(:,:,:,jn) ) 66 67 CASE('orlanski' ) ; CALL bdy_orl ( idx_bdy(ib_bdy), tsb(:,:,:,jn), tsa(:,:,:,jn), zdta(jn)%tra, ll_npo=.false. ) 67 68 CASE('orlanski_npo') ; CALL bdy_orl ( idx_bdy(ib_bdy), tsb(:,:,:,jn), tsa(:,:,:,jn), zdta(jn)%tra, ll_npo=.true. ) -
branches/2016/dev_NOC_CMCC_merge_2016/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90
r7299 r7397 24 24 USE oce ! ocean dynamics and tracers 25 25 USE dom_oce ! ocean space and time domain 26 !26 USE bdy_oce 27 27 USE in_out_manager ! I/O manager 28 28 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 29 29 USE lib_mpp ! 30 USE iom 30 31 USE wrk_nemo ! Memory allocation 31 32 USE timing ! Timing … … 102 103 INTEGER :: iif, iil, ii0, ii1, ii ! local integers 103 104 INTEGER :: ijf, ijl, ij0, ij1 ! - - 104 INTEGER :: ios 105 INTEGER :: ios, inum 105 106 INTEGER :: isrow ! index for ORCA1 starting row 106 107 INTEGER , POINTER, DIMENSION(:,:) :: imsk … … 108 109 !! 109 110 NAMELIST/namlbc/ rn_shlat, ln_vorlat 111 #if defined key_bdy 112 NAMELIST/nambdy/ nb_bdy, ln_coords_file, cn_coords_file, & 113 & ln_mask_file, cn_mask_file, cn_dyn2d, nn_dyn2d_dta, & 114 & cn_dyn3d, nn_dyn3d_dta, cn_tra, nn_tra_dta, & 115 & ln_tra_dmp, ln_dyn3d_dmp, rn_time_dmp, rn_time_dmp_out, & 116 & cn_ice_lim, nn_ice_lim_dta, & 117 & rn_ice_tem, rn_ice_sal, rn_ice_age, & 118 & ln_vol, nn_volctl, nn_rimwidth, nb_jpk_bdy 119 #endif 110 120 !!--------------------------------------------------------------------- 111 121 ! … … 155 165 END DO 156 166 167 #if defined key_bdy 168 REWIND( numnam_ref ) ! Namelist nambdy in reference namelist :Unstructured open boundaries 169 READ ( numnam_ref, nambdy, IOSTAT = ios, ERR = 903) 170 903 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambdy in reference namelist', lwp ) 171 172 REWIND( numnam_cfg ) ! Namelist nambdy in configuration namelist :Unstructured open boundaries 173 READ ( numnam_cfg, nambdy, IOSTAT = ios, ERR = 904 ) 174 904 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambdy in configuration namelist', lwp ) 175 IF(lwm) WRITE ( numond, nambdy ) 176 177 IF( ln_mask_file ) THEN ! correct for bdy mask 178 CALL iom_open( cn_mask_file, inum ) 179 CALL iom_get ( inum, jpdom_data, 'bdy_msk', bdytmask(:,:) ) 180 CALL iom_close( inum ) 181 182 ! Mask corrections 183 ! ---------------- 184 DO jk = 1, jpkm1 185 DO jj = 1, jpj 186 DO ji = 1, jpi 187 tmask(ji,jj,jk) = tmask(ji,jj,jk) * bdytmask(ji,jj) 188 END DO 189 END DO 190 END DO 191 ENDIF 192 #endif 157 193 ! (ISF) define barotropic mask and mask the ice shelf point 158 194 ssmask(:,:)=tmask(:,:,1) ! at this stage ice shelf is not masked -
branches/2016/dev_NOC_CMCC_merge_2016/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90
r6351 r7397 874 874 ! 875 875 ELSE !* Initialize at "rest" 876 e3t_b(:,:,:) = e3t_0(:,:,:) 877 e3t_n(:,:,:) = e3t_0(:,:,:) 878 sshn(:,:) = 0.0_wp 879 880 IF( ln_wd ) THEN 876 ! 877 IF( ln_wd .AND. ( cp_cfg == 'wad' ) ) THEN 878 ! 879 CALL wad_istate ! WAD test configuration : start from 880 ! uniform T-S fields and initial ssh slope 881 ! needs to be called here and in istate which is called later. 882 ! Adjust vertical metrics 883 DO jk = 1, jpk 884 e3t_n(:,:,jk) = e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) & 885 & / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk) & 886 & + e3t_0(:,:,jk) * (1._wp -tmask(:,:,jk)) 887 END DO 888 e3t_b(:,:,:) = e3t_n(:,:,:) 889 ! 890 ELSEIF( ln_wd ) THEN 891 ! 881 892 DO jj = 1, jpj 882 893 DO ji = 1, jpi 883 894 IF( e3t_0(ji,jj,1) <= 0.5_wp * rn_wdmin1 ) THEN 884 e3t_b(ji,jj,:) = 0.5_wp * rn_wdmin1 885 e3t_n(ji,jj,:) = 0.5_wp * rn_wdmin1 886 e3t_a(ji,jj,:) = 0.5_wp * rn_wdmin1 895 e3t_b(ji,jj,:) = 0.5_wp * rn_wdmin1 896 e3t_n(ji,jj,:) = 0.5_wp * rn_wdmin1 897 e3t_a(ji,jj,:) = 0.5_wp * rn_wdmin1 887 898 sshb(ji,jj) = rn_wdmin1 - bathy(ji,jj) 888 899 sshn(ji,jj) = rn_wdmin1 - bathy(ji,jj) … … 891 902 ENDDO 892 903 ENDDO 904 ! 905 ELSE 906 ! 907 e3t_b(:,:,:) = e3t_0(:,:,:) 908 e3t_n(:,:,:) = e3t_0(:,:,:) 909 sshn(:,:) = 0.0_wp 910 ! 893 911 END IF 894 912 -
branches/2016/dev_NOC_CMCC_merge_2016/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
r6492 r7397 421 421 IF(lwp) WRITE(numout,*) ' Depth = rn_bathy read in namelist' 422 422 zdta(:,:) = rn_bathy 423 ! 424 IF( cp_cfg == 'wad' ) THEN 425 SELECT CASE ( jp_cfg ) 426 ! ! ==================== 427 CASE ( 1 ) ! WAD 1 configuration 428 ! ! ==================== 429 ! 430 IF(lwp) WRITE(numout,*) 431 IF(lwp) WRITE(numout,*) 'zgr_bat : Closed box with EW linear bottom slope' 432 IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 433 ! 434 zdta = 1.5_wp 435 DO ji = 10, jpidta 436 zi = MIN(FLOAT(ji - 10)/FLOAT(jpidta - 10), 1.0 ) 437 zdta(ji,:) = MAX(rn_bathy*zi, 1.5) 438 IF(lwp)write(numout,*) 'ZDTA ',ji,zi,zdta(ji,1) 439 END DO 440 !!DO ji = 1, jpidta 441 !! zi = 1.0-EXP(-0.045*(ji-25.0)**2) 442 !! zdta(ji,:) = MAX(rn_bathy*zi, 1.5) 443 !! IF(lwp)write(numout,*) 'ZDTA ',ji,zi,zdta(ji,1) 444 !!END DO 445 zdta(1:2,:) = -2._wp 446 zdta(jpidta-1:jpidta,:) = -2._wp 447 zdta(:,1) = -2._wp 448 zdta(:,jpjdta) = -2._wp 449 zdta(:,1:3) = -2._wp 450 zdta(:,jpjdta-2:jpjdta) = -2._wp 451 ! ! ==================== 452 CASE ( 2, 3 ) ! WAD 2 or 3 configuration 453 ! ! ==================== 454 ! 455 IF(lwp) WRITE(numout,*) 456 IF(lwp) WRITE(numout,*) 'zgr_bat : Parobolic EW channel' 457 IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 458 ! 459 DO ji = 1, jpidta 460 zi = MAX(1.0-FLOAT((ji-25)**2)/484.0, 0.0 ) 461 zi = MAX(1.0-FLOAT((ji-25)**2)/484.0, -2.0 ) 462 zdta(ji,:) = MAX(rn_bathy*zi, -20.0) 463 IF(lwp)write(numout,*) 'ZDTA ',ji,zi,zdta(ji,1) 464 END DO 465 zdta(1:2,:) = -2._wp 466 zdta(jpidta-1:jpidta,:) = -2._wp 467 zdta(:,1) = -2._wp 468 zdta(:,jpjdta) = -2._wp 469 zdta(:,1:3) = -2._wp 470 zdta(:,jpjdta-2:jpjdta) = -2._wp 471 ! ! ==================== 472 CASE ( 4 ) ! WAD 4 configuration 473 ! ! ==================== 474 ! 475 IF(lwp) WRITE(numout,*) 476 IF(lwp) WRITE(numout,*) 'zgr_bat : Parobolic bowl' 477 IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 478 ! 479 DO ji = 1, jpidta 480 zi = MAX(1.0-FLOAT((ji-25)**2)/484.0, -2.0 ) 481 DO jj = 1, jpjdta 482 zj = MAX(1.0-FLOAT((jj-17)**2)/196.0, -2.0 ) 483 zdta(ji,jj) = MAX(rn_bathy*zi*zj, -2.0) 484 END DO 485 IF(lwp)write(numout,*) 'ZDTA ',ji,zi,zdta(ji,1) 486 END DO 487 zdta(1:2,:) = -2._wp 488 zdta(jpidta-1:jpidta,:) = -2._wp 489 zdta(:,1) = -2._wp 490 zdta(:,jpjdta) = -2._wp 491 zdta(:,1:3) = -2._wp 492 zdta(:,jpjdta-2:jpjdta) = -2._wp 493 ! ! =========================== 494 CASE ( 5 ) ! WAD 5 configuration 495 ! ! ==================== 496 ! 497 IF(lwp) WRITE(numout,*) 498 IF(lwp) WRITE(numout,*) 'zgr_bat : Double slope with shelf' 499 IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 500 ! 501 DO ji = 1, jpidta 502 zi = MIN(FLOAT(ji)/FLOAT(jpidta - 5), 1.0 ) 503 zdta(ji,:) = MAX(rn_bathy*zi, 0.5) 504 IF(lwp)write(numout,*) 'ZDTA ',ji,zi,zdta(ji,1) 505 END DO 506 DO ji = jpidta,46,-1 507 zdta(ji,:) = 10.0 508 IF(lwp)write(numout,*) 'ZDTA ',ji,zi,zdta(ji,1) 509 END DO 510 DO ji = 46,20,-1 511 zi = 7.5/25. 512 zdta(ji,:) = MAX(10. - zi*(47.-ji),2.5) 513 IF(lwp)write(numout,*) 'ZDTA ',ji,zi,zdta(ji,1) 514 END DO 515 DO ji = 19,15,-1 516 zdta(ji,:) = 2.5 517 IF(lwp)write(numout,*) 'ZDTA ',ji,zi,zdta(ji,1) 518 END DO 519 DO ji = 15,4,-1 520 zi = 2.0/11.0 521 zdta(ji,:) = MAX(2.5 - zi*(16-ji), 0.5) 522 IF(lwp)write(numout,*) 'ZDTA ',ji,zi,zdta(ji,1) 523 END DO 524 DO ji = 4,1,-1 525 zdta(ji,:) = 0.5 526 IF(lwp)write(numout,*) 'ZDTA ',ji,zi,zdta(ji,1) 527 END DO 528 ! ! =========================== 529 zdta(1:2,:) = -4._wp 530 zdta(jpidta-1:jpidta,:) = -4._wp 531 zdta(:,1) = -4._wp 532 zdta(:,jpjdta) = -4._wp 533 zdta(:,1:3) = -4._wp 534 zdta(:,jpjdta-2:jpjdta) = -4._wp 535 ! ! =========================== 536 CASE ( 6 ) ! WAD 6 configuration 537 ! ! ==================== 538 ! 539 IF(lwp) WRITE(numout,*) 540 IF(lwp) WRITE(numout,*) 'zgr_bat : Parabolic channel with gaussian ridge' 541 IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 542 ! 543 DO ji = 1, jpidta 544 zi = MAX(1.0-FLOAT((ji-25)**2)/484.0, -2.0 ) 545 zj = 0.95*MAX(EXP(-1.0*FLOAT((ji-25)**2)/32.0) , 0.0 ) 546 zdta(ji,:) = MAX(rn_bathy*(zi-zj), -2.0) 547 IF(lwp)write(numout,*) 'ZDTA ',ji,zi,zdta(ji,1) 548 END DO 549 zdta(1:2,:) = -4._wp 550 zdta(jpidta-1:jpidta,:) = -4._wp 551 zdta(:,1) = -4._wp 552 zdta(:,jpjdta) = -4._wp 553 zdta(:,1:3) = -4._wp 554 zdta(:,jpjdta-2:jpjdta) = -4._wp 555 ! ! =========================== 556 CASE DEFAULT 557 ! ! =========================== 558 WRITE(ctmp1,*) 'WAD test with a ', jp_cfg,' option is not coded' 559 ! 560 CALL ctl_stop( ctmp1 ) 561 ! 562 END SELECT 563 END IF 564 ! 423 565 IF( ln_sco ) THEN ! s-coordinate (zsc ): idta()=jpk 424 566 idta(:,:) = jpkm1 … … 2193 2335 CALL lbc_lnk( e3vw_0, 'V', 1._wp ) 2194 2336 ! 2195 IF( .NOT.ln_wd ) THEN 2196 WHERE( e3t_0 (:,:,:) == 0._wp ) e3t_0 (:,:,:) = 1._wp 2197 WHERE( e3u_0 (:,:,:) == 0._wp ) e3u_0 (:,:,:) = 1._wp 2198 WHERE( e3v_0 (:,:,:) == 0._wp ) e3v_0 (:,:,:) = 1._wp 2199 WHERE( e3f_0 (:,:,:) == 0._wp ) e3f_0 (:,:,:) = 1._wp 2200 WHERE( e3w_0 (:,:,:) == 0._wp ) e3w_0 (:,:,:) = 1._wp 2201 WHERE( e3uw_0(:,:,:) == 0._wp ) e3uw_0(:,:,:) = 1._wp 2202 WHERE( e3vw_0(:,:,:) == 0._wp ) e3vw_0(:,:,:) = 1._wp 2203 END IF 2337 WHERE( e3t_0 (:,:,:) == 0._wp ) e3t_0 (:,:,:) = 1._wp 2338 WHERE( e3u_0 (:,:,:) == 0._wp ) e3u_0 (:,:,:) = 1._wp 2339 WHERE( e3v_0 (:,:,:) == 0._wp ) e3v_0 (:,:,:) = 1._wp 2340 WHERE( e3f_0 (:,:,:) == 0._wp ) e3f_0 (:,:,:) = 1._wp 2341 WHERE( e3w_0 (:,:,:) == 0._wp ) e3w_0 (:,:,:) = 1._wp 2342 WHERE( e3uw_0(:,:,:) == 0._wp ) e3uw_0(:,:,:) = 1._wp 2343 WHERE( e3vw_0(:,:,:) == 0._wp ) e3vw_0(:,:,:) = 1._wp 2204 2344 2205 2345 #if defined key_agrif … … 2303 2443 DO jk = 1, mbathy(ji,jj) 2304 2444 ! check coordinate is monotonically increasing 2305 IF (e3w_ n(ji,jj,jk) <= 0._wp .OR. e3t_n(ji,jj,jk) <= 0._wp ) THEN2445 IF (e3w_0(ji,jj,jk) <= 0._wp .OR. e3t_0(ji,jj,jk) <= 0._wp ) THEN 2306 2446 WRITE(ctmp1,*) 'ERROR zgr_sco : e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk 2307 2447 WRITE(numout,*) 'ERROR zgr_sco : e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk 2308 WRITE(numout,*) 'e3w',e3w_ n(ji,jj,:)2309 WRITE(numout,*) 'e3t',e3t_ n(ji,jj,:)2448 WRITE(numout,*) 'e3w',e3w_0(ji,jj,:) 2449 WRITE(numout,*) 'e3t',e3t_0(ji,jj,:) 2310 2450 CALL ctl_stop( ctmp1 ) 2311 2451 ENDIF 2312 2452 ! and check it has never gone negative 2313 IF( gdepw_ n(ji,jj,jk) < 0._wp .OR. gdept_n(ji,jj,jk) < 0._wp ) THEN2453 IF( gdepw_0(ji,jj,jk) < 0._wp .OR. gdept_0(ji,jj,jk) < 0._wp ) THEN 2314 2454 WRITE(ctmp1,*) 'ERROR zgr_sco : gdepw or gdept =< 0 at point (i,j,k)= ', ji, jj, jk 2315 2455 WRITE(numout,*) 'ERROR zgr_sco : gdepw or gdept =< 0 at point (i,j,k)= ', ji, jj, jk 2316 WRITE(numout,*) 'gdepw',gdepw_ n(ji,jj,:)2317 WRITE(numout,*) 'gdept',gdept_ n(ji,jj,:)2456 WRITE(numout,*) 'gdepw',gdepw_0(ji,jj,:) 2457 WRITE(numout,*) 'gdept',gdept_0(ji,jj,:) 2318 2458 CALL ctl_stop( ctmp1 ) 2319 2459 ENDIF 2320 2460 ! and check it never exceeds the total depth 2321 IF( gdepw_ n(ji,jj,jk) > hbatt(ji,jj) ) THEN2461 IF( gdepw_0(ji,jj,jk) > hbatt(ji,jj) ) THEN 2322 2462 WRITE(ctmp1,*) 'ERROR zgr_sco : gdepw > hbatt at point (i,j,k)= ', ji, jj, jk 2323 2463 WRITE(numout,*) 'ERROR zgr_sco : gdepw > hbatt at point (i,j,k)= ', ji, jj, jk 2324 WRITE(numout,*) 'gdepw',gdepw_ n(ji,jj,:)2464 WRITE(numout,*) 'gdepw',gdepw_0(ji,jj,:) 2325 2465 CALL ctl_stop( ctmp1 ) 2326 2466 ENDIF … … 2329 2469 DO jk = 1, mbathy(ji,jj)-1 2330 2470 ! and check it never exceeds the total depth 2331 IF( gdept_ n(ji,jj,jk) > hbatt(ji,jj) ) THEN2471 IF( gdept_0(ji,jj,jk) > hbatt(ji,jj) ) THEN 2332 2472 WRITE(ctmp1,*) 'ERROR zgr_sco : gdept > hbatt at point (i,j,k)= ', ji, jj, jk 2333 2473 WRITE(numout,*) 'ERROR zgr_sco : gdept > hbatt at point (i,j,k)= ', ji, jj, jk 2334 WRITE(numout,*) 'gdept',gdept_ n(ji,jj,:)2474 WRITE(numout,*) 'gdept',gdept_0(ji,jj,:) 2335 2475 CALL ctl_stop( ctmp1 ) 2336 2476 ENDIF -
branches/2016/dev_NOC_CMCC_merge_2016/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90
r6140 r7397 36 36 USE domvvl ! varying vertical mesh 37 37 USE iscplrst ! ice sheet coupling 38 USE wet_dry ! wetting and drying (needed for wad_istate) 38 39 ! 39 40 USE in_out_manager ! I/O manager … … 105 106 ELSEIF( cp_cfg == 'gyre' ) THEN 106 107 CALL istate_gyre ! GYRE configuration : start from pre-defined T-S fields 108 ELSEIF( cp_cfg == 'wad' ) THEN 109 CALL wad_istate ! WAD test configuration : start from pre-defined T-S fields and initial ssh slope 107 110 ELSE ! Initial T-S, U-V fields read in files 108 111 IF ( ln_tsd_init ) THEN ! read 3D T and S data at nit000 -
branches/2016/dev_NOC_CMCC_merge_2016/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
r6152 r7397 432 432 INTEGER :: ji, jj, jk, jii, jjj ! dummy loop indices 433 433 REAL(wp) :: zcoef0, zuap, zvap, znad, ztmp ! temporary scalars 434 LOGICAL :: ll_tmp1, ll_tmp2 , ll_tmp3! local logical variables434 LOGICAL :: ll_tmp1, ll_tmp2 ! local logical variables 435 435 REAL(wp), POINTER, DIMENSION(:,:,:) :: zhpi, zhpj 436 436 REAL(wp), POINTER, DIMENSION(:,:) :: zcpx, zcpy !W/D pressure filter … … 438 438 ! 439 439 CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj ) 440 IF( ln_wd) CALL wrk_alloc( jpi,jpj, zcpx, zcpy )440 IF( ln_wd ) CALL wrk_alloc( jpi,jpj, zcpx, zcpy ) 441 441 ! 442 442 IF( kt == nit000 ) THEN … … 451 451 ENDIF 452 452 ! 453 IF( ln_wd) THEN453 IF( ln_wd ) THEN 454 454 DO jj = 2, jpjm1 455 455 DO ji = 2, jpim1 456 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) 457 ll_tmp2 = MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji+1,jj) + bathy(ji+1,jj)) > rn_wdmin1 + rn_wdmin2 458 ll_tmp3 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) + & 459 & rn_wdmin1 + rn_wdmin2 460 461 IF(ll_tmp1.AND.ll_tmp2) THEN 456 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji+1,jj) ) > & 457 & MAX( -bathy(ji,jj) , -bathy(ji+1,jj) ) .AND. & 458 & MAX( sshn(ji,jj) + bathy(ji,jj), sshn(ji+1,jj) + bathy(ji+1,jj) ) & 459 & > rn_wdmin1 + rn_wdmin2 460 ll_tmp2 = MAX( sshn(ji,jj) , sshn(ji+1,jj) ) > & 461 & MAX( -bathy(ji,jj) , -bathy(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 462 463 IF(ll_tmp1) THEN 462 464 zcpx(ji,jj) = 1.0_wp 463 wduflt(ji,jj) = 1.0_wp 464 ELSE IF(ll_tmp3) THEN 465 ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen ! here 466 zcpx(ji,jj) = ABS((sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) / & 467 & (sshn(ji+1,jj) - sshn(ji,jj))) 468 wduflt(ji,jj) = 1.0_wp 465 ELSE IF(ll_tmp2) THEN 466 ! no worries about sshn(ji+1,jj) - sshn(ji ,jj) = 0, it won't happen ! here 467 zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) & 468 & / (sshn(ji+1,jj) - sshn(ji ,jj)) ) 469 469 ELSE 470 470 zcpx(ji,jj) = 0._wp 471 wduflt(ji,jj) = 0.0_wp472 471 END IF 473 472 474 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) 475 ll_tmp2 = MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji,jj+1) + bathy(ji,jj+1)) > rn_wdmin1 + rn_wdmin2 476 ll_tmp3 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) + & 477 & rn_wdmin1 + rn_wdmin2 478 479 IF(ll_tmp1.AND.ll_tmp2) THEN 473 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji,jj+1) ) > & 474 & MAX( -bathy(ji,jj) , -bathy(ji,jj+1) ) .AND. & 475 & MAX( sshn(ji,jj) + bathy(ji,jj), sshn(ji,jj+1) + bathy(ji,jj+1) ) & 476 & > rn_wdmin1 + rn_wdmin2 477 ll_tmp2 = MAX( sshn(ji,jj) , sshn(ji,jj+1) ) > & 478 & MAX( -bathy(ji,jj) , -bathy(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 479 480 IF(ll_tmp1) THEN 480 481 zcpy(ji,jj) = 1.0_wp 481 wdvflt(ji,jj) = 1.0_wp 482 ELSE IF(ll_tmp3) THEN 483 ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen ! here 484 zcpy(ji,jj) = ABS((sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) / & 485 & (sshn(ji,jj+1) - sshn(ji,jj))) 486 wdvflt(ji,jj) = 1.0_wp 482 ELSE IF(ll_tmp2) THEN 483 ! no worries about sshn(ji,jj+1) - sshn(ji,jj ) = 0, it won't happen ! here 484 zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) & 485 & / (sshn(ji,jj+1) - sshn(ji,jj )) ) 487 486 ELSE 488 487 zcpy(ji,jj) = 0._wp 489 wdvflt(ji,jj) = 0.0_wp490 488 END IF 491 489 END DO 492 490 END DO 493 491 CALL lbc_lnk( zcpx, 'U', 1._wp ) ; CALL lbc_lnk( zcpy, 'V', 1._wp ) 494 ENDIF 495 492 END IF 496 493 497 494 ! Surface value … … 510 507 511 508 512 IF( ln_wd) THEN509 IF( ln_wd ) THEN 513 510 514 511 zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj) … … 541 538 & * ( gde3w_n(ji ,jj+1,jk) - gde3w_n(ji,jj,jk) ) * r1_e2v(ji,jj) 542 539 543 IF( ln_wd) THEN540 IF( ln_wd ) THEN 544 541 zhpi(ji,jj,jk) = zhpi(ji,jj,jk) * zcpx(ji,jj) 545 542 zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj) … … 556 553 ! 557 554 CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj ) 558 IF( ln_wd) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy )555 IF( ln_wd ) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy ) 559 556 ! 560 557 END SUBROUTINE hpg_sco … … 701 698 CALL wrk_alloc( jpi, jpj, jpk, drhox, drhoy, drhoz, drhou, drhov, drhow ) 702 699 CALL wrk_alloc( jpi, jpj, jpk, rho_i, rho_j, rho_k, zhpi, zhpj ) 703 IF( ln_wd) CALL wrk_alloc( jpi,jpj, zcpx, zcpy )704 ! 705 ! 706 IF( ln_wd) THEN700 IF( ln_wd ) CALL wrk_alloc( jpi,jpj, zcpx, zcpy ) 701 ! 702 ! 703 IF( ln_wd ) THEN 707 704 DO jj = 2, jpjm1 708 705 DO ji = 2, jpim1 709 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) & 710 & .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji+1,jj) + bathy(ji+1,jj)) & 711 & > rn_wdmin1 + rn_wdmin2 712 ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) +& 713 & rn_wdmin1 + rn_wdmin2 706 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji+1,jj) ) > & 707 & MAX( -bathy(ji,jj) , -bathy(ji+1,jj) ) .AND. & 708 & MAX( sshn(ji,jj) + bathy(ji,jj), sshn(ji+1,jj) + bathy(ji+1,jj) ) & 709 & > rn_wdmin1 + rn_wdmin2 710 ll_tmp2 = MAX( sshn(ji,jj) , sshn(ji+1,jj) ) > & 711 & MAX( -bathy(ji,jj) , -bathy(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 714 712 715 713 IF(ll_tmp1) THEN 716 714 zcpx(ji,jj) = 1.0_wp 717 715 ELSE IF(ll_tmp2) THEN 718 ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen ! here719 zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) /&720 & (sshn(ji+1,jj) - sshn(ji,jj)))716 ! no worries about sshn(ji+1,jj) - sshn(ji ,jj) = 0, it won't happen ! here 717 zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) & 718 & / (sshn(ji+1,jj) - sshn(ji ,jj)) ) 721 719 ELSE 722 720 zcpx(ji,jj) = 0._wp 723 721 END IF 724 722 725 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) & 726 & .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji,jj+1) + bathy(ji,jj+1)) & 727 & > rn_wdmin1 + rn_wdmin2 728 ll_tmp2 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) +& 729 & rn_wdmin1 + rn_wdmin2 723 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji,jj+1) ) > & 724 & MAX( -bathy(ji,jj) , -bathy(ji,jj+1) ) .AND. & 725 & MAX( sshn(ji,jj) + bathy(ji,jj), sshn(ji,jj+1) + bathy(ji,jj+1) ) & 726 & > rn_wdmin1 + rn_wdmin2 727 ll_tmp2 = MAX( sshn(ji,jj) , sshn(ji,jj+1) ) > & 728 & MAX( -bathy(ji,jj) , -bathy(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 730 729 731 730 IF(ll_tmp1) THEN 732 731 zcpy(ji,jj) = 1.0_wp 733 732 ELSE IF(ll_tmp2) THEN 734 ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen ! here735 zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) /&736 & (sshn(ji,jj+1) - sshn(ji,jj)))733 ! no worries about sshn(ji,jj+1) - sshn(ji,jj ) = 0, it won't happen ! here 734 zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) & 735 & / (sshn(ji,jj+1) - sshn(ji,jj )) ) 737 736 ELSE 738 737 zcpy(ji,jj) = 0._wp … … 741 740 END DO 742 741 CALL lbc_lnk( zcpx, 'U', 1._wp ) ; CALL lbc_lnk( zcpy, 'V', 1._wp ) 743 ENDIF 744 742 END IF 745 743 746 744 IF( kt == nit000 ) THEN … … 913 911 zhpi(ji,jj,1) = ( rho_k(ji+1,jj ,1) - rho_k(ji,jj,1) - rho_i(ji,jj,1) ) * r1_e1u(ji,jj) 914 912 zhpj(ji,jj,1) = ( rho_k(ji ,jj+1,1) - rho_k(ji,jj,1) - rho_j(ji,jj,1) ) * r1_e2v(ji,jj) 915 IF( ln_wd) THEN913 IF( ln_wd ) THEN 916 914 zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj) 917 915 zhpj(ji,jj,1) = zhpj(ji,jj,1) * zcpy(ji,jj) … … 936 934 & + ( ( rho_k(ji,jj+1,jk) - rho_k(ji,jj,jk ) ) & 937 935 & -( rho_j(ji,jj ,jk) - rho_j(ji,jj,jk-1) ) ) * r1_e2v(ji,jj) 938 IF( ln_wd) THEN936 IF( ln_wd ) THEN 939 937 zhpi(ji,jj,jk) = zhpi(ji,jj,jk) * zcpx(ji,jj) 940 938 zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj) … … 950 948 CALL wrk_dealloc( jpi, jpj, jpk, drhox, drhoy, drhoz, drhou, drhov, drhow ) 951 949 CALL wrk_dealloc( jpi, jpj, jpk, rho_i, rho_j, rho_k, zhpi, zhpj ) 952 IF( ln_wd) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy )950 IF( ln_wd ) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy ) 953 951 ! 954 952 END SUBROUTINE hpg_djc … … 987 985 CALL wrk_alloc( jpi,jpj,jpk, zdept, zrhh ) 988 986 CALL wrk_alloc( jpi,jpj, zsshu_n, zsshv_n ) 989 IF( ln_wd) CALL wrk_alloc( jpi,jpj, zcpx, zcpy )987 IF( ln_wd ) CALL wrk_alloc( jpi,jpj, zcpx, zcpy ) 990 988 ! 991 989 IF( kt == nit000 ) THEN … … 1000 998 IF( ln_linssh ) znad = 0._wp 1001 999 1002 IF( ln_wd) THEN1000 IF( ln_wd ) THEN 1003 1001 DO jj = 2, jpjm1 1004 1002 DO ji = 2, jpim1 1005 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) & 1006 & .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji+1,jj) + bathy(ji+1,jj)) & 1007 & > rn_wdmin1 + rn_wdmin2 1008 ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) +& 1009 & rn_wdmin1 + rn_wdmin2 1003 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji+1,jj) ) > & 1004 & MAX( -bathy(ji,jj) , -bathy(ji+1,jj) ) .AND. & 1005 & MAX( sshn(ji,jj) + bathy(ji,jj), sshn(ji+1,jj) + bathy(ji+1,jj) ) & 1006 & > rn_wdmin1 + rn_wdmin2 1007 ll_tmp2 = MAX( sshn(ji,jj) , sshn(ji+1,jj) ) > & 1008 & MAX( -bathy(ji,jj) , -bathy(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 1010 1009 1011 1010 IF(ll_tmp1) THEN 1012 1011 zcpx(ji,jj) = 1.0_wp 1013 1012 ELSE IF(ll_tmp2) THEN 1014 ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen ! here1015 zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) /&1016 & (sshn(ji+1,jj) - sshn(ji,jj)))1013 ! no worries about sshn(ji+1,jj) - sshn(ji ,jj) = 0, it won't happen ! here 1014 zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) & 1015 & / (sshn(ji+1,jj) - sshn(ji ,jj)) ) 1017 1016 ELSE 1018 1017 zcpx(ji,jj) = 0._wp 1019 1018 END IF 1020 1019 1021 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) & 1022 & .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji,jj+1) + bathy(ji,jj+1)) & 1023 & > rn_wdmin1 + rn_wdmin2 1024 ll_tmp2 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) +& 1025 & rn_wdmin1 + rn_wdmin2 1026 1027 IF(ll_tmp1.OR.ll_tmp2) THEN 1020 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji,jj+1) ) > & 1021 & MAX( -bathy(ji,jj) , -bathy(ji,jj+1) ) .AND. & 1022 & MAX( sshn(ji,jj) + bathy(ji,jj), sshn(ji,jj+1) + bathy(ji,jj+1) ) & 1023 & > rn_wdmin1 + rn_wdmin2 1024 ll_tmp2 = MAX( sshn(ji,jj) , sshn(ji,jj+1) ) > & 1025 & MAX( -bathy(ji,jj) , -bathy(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 1026 1027 IF(ll_tmp1) THEN 1028 1028 zcpy(ji,jj) = 1.0_wp 1029 1029 ELSE IF(ll_tmp2) THEN 1030 ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen ! here1031 zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) /&1032 & (sshn(ji,jj+1) - sshn(ji,jj)))1030 ! no worries about sshn(ji,jj+1) - sshn(ji,jj ) = 0, it won't happen ! here 1031 zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) & 1032 & / (sshn(ji,jj+1) - sshn(ji,jj )) ) 1033 1033 ELSE 1034 1034 zcpy(ji,jj) = 0._wp … … 1037 1037 END DO 1038 1038 CALL lbc_lnk( zcpx, 'U', 1._wp ) ; CALL lbc_lnk( zcpy, 'V', 1._wp ) 1039 END IF1039 END IF 1040 1040 1041 1041 ! Clean 3-D work arrays … … 1221 1221 zdpdx2 = zcoef0 * r1_e1u(ji,jj) * REAL(jis-jid, wp) * (zpwes + zpwed) 1222 1222 ENDIF 1223 IF( ln_wd) THEN1223 IF( ln_wd ) THEN 1224 1224 zdpdx1 = zdpdx1 * zcpx(ji,jj) 1225 1225 zdpdx2 = zdpdx2 * zcpx(ji,jj) … … 1280 1280 zdpdy2 = zcoef0 * r1_e2v(ji,jj) * REAL(jjs-jjd, wp) * (zpnss + zpnsd ) 1281 1281 ENDIF 1282 IF( ln_wd) THEN1282 IF( ln_wd ) THEN 1283 1283 zdpdy1 = zdpdy1 * zcpy(ji,jj) 1284 1284 zdpdy2 = zdpdy2 * zcpy(ji,jj) … … 1295 1295 CALL wrk_dealloc( jpi,jpj,jpk, zdept, zrhh ) 1296 1296 CALL wrk_dealloc( jpi,jpj, zsshu_n, zsshv_n ) 1297 IF( ln_wd) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy )1297 IF( ln_wd ) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy ) 1298 1298 ! 1299 1299 END SUBROUTINE hpg_prj -
branches/2016/dev_NOC_CMCC_merge_2016/NEMOGCM/NEMO/OPA_SRC/DYN/dynkeg.F90
r5328 r7397 24 24 USE wrk_nemo ! Memory Allocation 25 25 USE timing ! Timing 26 #if defined key_bdy 27 USE bdy_oce ! ocean open boundary conditions 28 #endif 26 29 27 30 IMPLICIT NONE … … 78 81 REAL(wp), POINTER, DIMENSION(:,:,:) :: zhke 79 82 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdu, ztrdv 83 #if defined key_bdy 84 INTEGER :: jb ! dummy loop indices 85 INTEGER :: ii, ij, igrd, ib_bdy ! local integers 86 INTEGER :: fu, fv 87 #endif 80 88 !!---------------------------------------------------------------------- 81 89 ! … … 98 106 zhke(:,:,jpk) = 0._wp 99 107 108 #if defined key_bdy 109 ! Maria Luneva & Fred Wobus: July-2016 110 ! compensate for lack of turbulent kinetic energy on liquid bdy points 111 DO ib_bdy = 1, nb_bdy 112 IF( cn_dyn3d(ib_bdy) /= 'none' ) THEN 113 igrd = 2 ! Copying normal velocity into points outside bdy 114 DO jb = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 115 DO jk = 1, jpkm1 116 ii = idx_bdy(ib_bdy)%nbi(jb,igrd) 117 ij = idx_bdy(ib_bdy)%nbj(jb,igrd) 118 fu = NINT( idx_bdy(ib_bdy)%flagu(jb,igrd) ) 119 un(ii-fu,ij,jk) = un(ii,ij,jk) * umask(ii,ij,jk) 120 END DO 121 END DO 122 ! 123 igrd = 3 ! Copying normal velocity into points outside bdy 124 DO jb = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 125 DO jk = 1, jpkm1 126 ii = idx_bdy(ib_bdy)%nbi(jb,igrd) 127 ij = idx_bdy(ib_bdy)%nbj(jb,igrd) 128 fv = NINT( idx_bdy(ib_bdy)%flagv(jb,igrd) ) 129 vn(ii,ij-fv,jk) = vn(ii,ij,jk) * vmask(ii,ij,jk) 130 END DO 131 END DO 132 ENDIF 133 ENDDO 134 #endif 135 100 136 SELECT CASE ( kscheme ) !== Horizontal kinetic energy at T-point ==! 101 137 ! … … 133 169 ! 134 170 END SELECT 171 172 #if defined key_bdy 173 ! restore velocity masks at points outside boundary 174 un(:,:,:) = un(:,:,:) * umask(:,:,:) 175 vn(:,:,:) = vn(:,:,:) * vmask(:,:,:) 176 #endif 177 178 135 179 ! 136 180 DO jk = 1, jpkm1 !== grad( KE ) added to the general momentum trends ==! -
branches/2016/dev_NOC_CMCC_merge_2016/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
r7299 r7397 156 156 REAL(wp), POINTER, DIMENSION(:,:) :: zhf 157 157 REAL(wp), POINTER, DIMENSION(:,:) :: zcpx, zcpy ! Wetting/Dying gravity filter coef. 158 REAL(wp), POINTER, DIMENSION(:,:) :: wduflt1, wdvflt1 ! Wetting/Dying velocity filter coef.159 158 !!---------------------------------------------------------------------- 160 159 ! … … 168 167 CALL wrk_alloc( jpi,jpj, zsshu_a, zsshv_a ) 169 168 CALL wrk_alloc( jpi,jpj, zhf ) 170 IF( ln_wd ) CALL wrk_alloc( jpi, jpj, zcpx, zcpy , wduflt1, wdvflt1)169 IF( ln_wd ) CALL wrk_alloc( jpi, jpj, zcpx, zcpy ) 171 170 ! 172 171 zmdi=1.e+20 ! missing data indicator for masking … … 374 373 IF( .NOT.ln_linssh ) THEN ! Variable volume : remove surface pressure gradient 375 374 IF( ln_wd ) THEN ! Calculating and applying W/D gravity filters 376 wduflt1(:,:) = 1.0_wp377 wdvflt1(:,:) = 1.0_wp378 DO jj = 2, jpjm1379 DO ji = 2, jpim1380 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj))&381 & .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji+1,jj) + bathy(ji+1,jj)) &382 & > rn_wdmin1 + rn_wdmin2383 ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) &384 & + rn_wdmin1 + rn_wdmin2375 DO jj = 2, jpjm1 376 DO ji = 2, jpim1 377 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji+1,jj) ) > & 378 & MAX( -bathy(ji,jj) , -bathy(ji+1,jj) ) .AND. & 379 & MAX( sshn(ji,jj) + bathy(ji,jj), sshn(ji+1,jj) + bathy(ji+1,jj) ) & 380 & > rn_wdmin1 + rn_wdmin2 381 ll_tmp2 = MAX( sshn(ji,jj) , sshn(ji+1,jj) ) > & 382 & MAX( -bathy(ji,jj) , -bathy(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 383 385 384 IF(ll_tmp1) THEN 386 zcpx(ji,jj) 387 ELSE IF(ll_tmp2) THEN388 ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happenhere389 zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) &390 & /(sshn(ji+1,jj) - sshn(ji,jj)))385 zcpx(ji,jj) = 1.0_wp 386 ELSE IF(ll_tmp2) THEN 387 ! no worries about sshn(ji+1,jj) - sshn(ji ,jj) = 0, it won't happen ! here 388 zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) & 389 & / (sshn(ji+1,jj) - sshn(ji ,jj)) ) 391 390 ELSE 392 zcpx(ji,jj) = 0._wp 393 wduflt1(ji,jj) = 0.0_wp 391 zcpx(ji,jj) = 0._wp 394 392 END IF 395 396 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) & 397 & .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji,jj+1) + bathy(ji,jj+1)) & 398 & > rn_wdmin1 + rn_wdmin2 399 ll_tmp2 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) & 400 & + rn_wdmin1 + rn_wdmin2 393 394 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji,jj+1) ) > & 395 & MAX( -bathy(ji,jj) , -bathy(ji,jj+1) ) .AND. & 396 & MAX( sshn(ji,jj) + bathy(ji,jj), sshn(ji,jj+1) + bathy(ji,jj+1) ) & 397 & > rn_wdmin1 + rn_wdmin2 398 ll_tmp2 = MAX( sshn(ji,jj) , sshn(ji,jj+1) ) > & 399 & MAX( -bathy(ji,jj) , -bathy(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 400 401 401 IF(ll_tmp1) THEN 402 zcpy(ji,jj)= 1.0_wp403 ELSE IF(ll_tmp2) THEN404 ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happenhere405 zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) &406 & /(sshn(ji,jj+1) - sshn(ji,jj)))402 zcpy(ji,jj) = 1.0_wp 403 ELSE IF(ll_tmp2) THEN 404 ! no worries about sshn(ji,jj+1) - sshn(ji,jj ) = 0, it won't happen ! here 405 zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) & 406 & / (sshn(ji,jj+1) - sshn(ji,jj )) ) 407 407 ELSE 408 zcpy(ji,jj) = 0._wp 409 wdvflt1(ji,jj) = 0.0_wp 410 ENDIF 411 412 END DO 408 zcpy(ji,jj) = 0._wp 409 END IF 410 END DO 413 411 END DO 414 415 CALL lbc_lnk( zcpx, 'U', 1._wp ) ; CALL lbc_lnk( zcpy, 'V', 1._wp ) 416 412 417 413 DO jj = 2, jpjm1 418 414 DO ji = 2, jpim1 419 zu_trd(ji,jj) = (zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) &420 & * r1_e1u(ji,jj) ) * zcpx(ji,jj) * wduflt1(ji,jj)421 zv_trd(ji,jj) = (zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) &422 & * r1_e2v(ji,jj) ) * zcpy(ji,jj) * wdvflt1(ji,jj)415 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) & 416 & * r1_e1u(ji,jj) * zcpx(ji,jj) 417 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) & 418 & * r1_e2v(ji,jj) * zcpy(ji,jj) 423 419 END DO 424 420 END DO … … 567 563 ENDIF 568 564 569 IF( ln_wd ) THEN !preserve the positivity of water depth570 !ssh[b,n,a] should have already been processed for this571 sshbb_e(:,:) = MAX(sshbb_e(:,:), rn_wdmin1 - bathy(:,:))572 sshb_e(:,:) = MAX(sshb_e(:,:) , rn_wdmin1 - bathy(:,:))573 ENDIF574 565 ! 575 566 IF (ln_bt_fw) THEN ! FORWARD integration: start from NOW fields … … 644 635 zhup2_e (:,:) = hu_0(:,:) + zwx(:,:) ! Ocean depth at U- and V-points 645 636 zhvp2_e (:,:) = hv_0(:,:) + zwy(:,:) 646 IF( ln_wd ) THEN647 zhup2_e(:,:) = MAX(zhup2_e (:,:), rn_wdmin1)648 zhvp2_e(:,:) = MAX(zhvp2_e (:,:), rn_wdmin1)649 END IF650 637 ELSE 651 638 zhup2_e (:,:) = hu_n(:,:) … … 699 686 END DO 700 687 END DO 688 701 689 ssha_e(:,:) = ( sshn_e(:,:) - rdtbt * ( zssh_frc(:,:) + zhdiv(:,:) ) ) * ssmask(:,:) 702 IF( ln_wd ) ssha_e(:,:) = MAX(ssha_e(:,:), rn_wdmin1 - bathy(:,:))690 703 691 CALL lbc_lnk( ssha_e, 'T', 1._wp ) 704 692 … … 746 734 zsshp2_e(:,:) = za0 * ssha_e(:,:) + za1 * sshn_e (:,:) & 747 735 & + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:) 736 748 737 IF( ln_wd ) THEN ! Calculating and applying W/D gravity filters 749 wduflt1(:,:) = 1._wp750 wdvflt1(:,:) = 1._wp751 738 DO jj = 2, jpjm1 752 DO ji = 2, jpim1 753 ll_tmp1 = MIN( zsshp2_e(ji,jj), zsshp2_e(ji+1,jj) ) > MAX( -bathy(ji,jj), -bathy(ji+1,jj) ) & 754 & .AND. MAX( zsshp2_e(ji,jj) + bathy(ji,jj), zsshp2_e(ji+1,jj) + bathy(ji+1,jj) ) & 755 & > rn_wdmin1 + rn_wdmin2 756 ll_tmp2 = MAX( zsshp2_e(ji,jj), zsshp2_e(ji+1,jj) ) > MAX( -bathy(ji,jj), -bathy(ji+1,jj) ) & 757 & + rn_wdmin1 + rn_wdmin2 758 IF(ll_tmp1) THEN 759 zcpx(ji,jj) = 1._wp 760 ELSE IF(ll_tmp2) THEN 761 ! no worries about zsshp2_e(ji+1,jj)-zsshp2_e(ji,jj) = 0, it won't happen here 762 zcpx(ji,jj) = ABS( (zsshp2_e(ji+1,jj) + bathy(ji+1,jj) - zsshp2_e(ji,jj) - bathy(ji,jj)) & 763 & / (zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj)) ) 764 ELSE 765 zcpx(ji,jj) = 0._wp 766 wduflt1(ji,jj) = 0._wp 767 END IF 768 769 ll_tmp1 = MIN( zsshp2_e(ji,jj), zsshp2_e(ji,jj+1) ) > MAX( -bathy(ji,jj), -bathy(ji,jj+1) ) & 770 & .AND. MAX( zsshp2_e(ji,jj) + bathy(ji,jj), zsshp2_e(ji,jj+1) + bathy(ji,jj+1) ) & 771 & > rn_wdmin1 + rn_wdmin2 772 ll_tmp2 = MAX( zsshp2_e(ji,jj), zsshp2_e(ji,jj+1) ) > MAX( -bathy(ji,jj), -bathy(ji,jj+1) ) & 773 & + rn_wdmin1 + rn_wdmin2 774 IF(ll_tmp1) THEN 775 zcpy(ji,jj) = 1._wp 776 ELSE IF(ll_tmp2) THEN 777 ! no worries about zsshp2_e(ji,jj+1)-zsshp2_e(ji,jj) = 0, it won't happen here 778 zcpy(ji,jj) = ABS( (zsshp2_e(ji,jj+1) + bathy(ji,jj+1) - zsshp2_e(ji,jj) - bathy(ji,jj)) & 779 & / (zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj)) ) 780 ELSE 781 zcpy(ji,jj) = 0._wp 782 wdvflt1(ji,jj) = 0._wp 783 END IF 739 DO ji = 2, jpim1 740 ll_tmp1 = MIN( zsshp2_e(ji,jj) , zsshp2_e(ji+1,jj) ) > & 741 & MAX( -bathy(ji,jj) , -bathy(ji+1,jj) ) .AND. & 742 & MAX( zsshp2_e(ji,jj) + bathy(ji,jj), zsshp2_e(ji+1,jj) + bathy(ji+1,jj) ) & 743 & > rn_wdmin1 + rn_wdmin2 744 ll_tmp2 = MAX( zsshp2_e(ji,jj) , zsshp2_e(ji+1,jj) ) > & 745 & MAX( -bathy(ji,jj) , -bathy(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 746 747 IF(ll_tmp1) THEN 748 zcpx(ji,jj) = 1.0_wp 749 ELSE IF(ll_tmp2) THEN 750 ! no worries about zsshp2_e(ji+1,jj) - zsshp2_e(ji ,jj) = 0, it won't happen ! here 751 zcpx(ji,jj) = ABS( (zsshp2_e(ji+1,jj) + bathy(ji+1,jj) - zsshp2_e(ji,jj) - bathy(ji,jj)) & 752 & / (zsshp2_e(ji+1,jj) - zsshp2_e(ji ,jj)) ) 753 ELSE 754 zcpx(ji,jj) = 0._wp 755 END IF 756 757 ll_tmp1 = MIN( zsshp2_e(ji,jj) , zsshp2_e(ji,jj+1) ) > & 758 & MAX( -bathy(ji,jj) , -bathy(ji,jj+1) ) .AND. & 759 & MAX( zsshp2_e(ji,jj) + bathy(ji,jj), zsshp2_e(ji,jj+1) + bathy(ji,jj+1) ) & 760 & > rn_wdmin1 + rn_wdmin2 761 ll_tmp2 = MAX( zsshp2_e(ji,jj) , zsshp2_e(ji,jj+1) ) > & 762 & MAX( -bathy(ji,jj) , -bathy(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 763 764 IF(ll_tmp1) THEN 765 zcpy(ji,jj) = 1.0_wp 766 ELSE IF(ll_tmp2) THEN 767 ! no worries about zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj ) = 0, it won't happen ! here 768 zcpy(ji,jj) = ABS( (zsshp2_e(ji,jj+1) + bathy(ji,jj+1) - zsshp2_e(ji,jj) - bathy(ji,jj)) & 769 & / (zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj )) ) 770 ELSE 771 zcpy(ji,jj) = 0._wp 772 END IF 784 773 END DO 785 END DO 786 CALL lbc_lnk( zcpx, 'U', 1._wp ) ; CALL lbc_lnk( zcpy, 'V', 1._wp ) 787 ENDIF 774 END DO 775 END IF 788 776 ! 789 777 ! Compute associated depths at U and V points: … … 803 791 END DO 804 792 805 IF( ln_wd ) THEN806 zhust_e(:,:) = MAX(zhust_e (:,:), rn_wdmin1 )807 zhvst_e(:,:) = MAX(zhvst_e (:,:), rn_wdmin1 )808 END IF809 810 793 ENDIF 811 794 ! … … 885 868 zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj) 886 869 zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj) 887 zwx(ji,jj) = zu_spg * zcpx(ji,jj) 888 zwy(ji,jj) = zv_spg * zcpy(ji,jj) 870 zwx(ji,jj) = zu_spg * zcpx(ji,jj) * wdmask(ji,jj) * wdmask(ji+1, jj) 871 zwy(ji,jj) = zv_spg * zcpy(ji,jj) * wdmask(ji,jj) * wdmask(ji, jj+1) 889 872 END DO 890 873 END DO … … 924 907 DO ji = fs_2, fs_jpim1 ! vector opt. 925 908 926 IF( ln_wd ) THEN 927 zhura = MAX(hu_0(ji,jj) + zsshu_a(ji,jj), rn_wdmin1) 928 zhvra = MAX(hv_0(ji,jj) + zsshv_a(ji,jj), rn_wdmin1) 929 ELSE 930 zhura = hu_0(ji,jj) + zsshu_a(ji,jj) 931 zhvra = hv_0(ji,jj) + zsshv_a(ji,jj) 932 END IF 909 zhura = hu_0(ji,jj) + zsshu_a(ji,jj) 910 zhvra = hv_0(ji,jj) + zsshv_a(ji,jj) 933 911 zhura = ssumask(ji,jj)/(zhura + 1._wp - ssumask(ji,jj)) 934 912 zhvra = ssvmask(ji,jj)/(zhvra + 1._wp - ssvmask(ji,jj)) … … 950 928 ! 951 929 IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only) 952 IF( ln_wd ) THEN 953 hu_e (:,:) = MAX(hu_0(:,:) + zsshu_a(:,:), rn_wdmin1) 954 hv_e (:,:) = MAX(hv_0(:,:) + zsshv_a(:,:), rn_wdmin1) 955 ELSE 956 hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:) 957 hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:) 958 END IF 930 hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:) 931 hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:) 959 932 hur_e(:,:) = ssumask(:,:) / ( hu_e(:,:) + 1._wp - ssumask(:,:) ) 960 933 hvr_e(:,:) = ssvmask(:,:) / ( hv_e(:,:) + 1._wp - ssvmask(:,:) ) … … 1020 993 ! 1021 994 ! Update barotropic trend: 1022 IF( ln_dynadv_vec .OR. ln_linssh ) THEN 1023 DO jk=1,jpkm1 1024 ua(:,:,jk) = ua(:,:,jk) + ( ua_b(:,:) - ub_b(:,:) ) * z1_2dt_b 1025 va(:,:,jk) = va(:,:,jk) + ( va_b(:,:) - vb_b(:,:) ) * z1_2dt_b 1026 END DO 995 IF(ln_wd) THEN 996 IF( ln_dynadv_vec .OR. ln_linssh ) THEN 997 DO jk=1,jpkm1 998 ua(:,:,jk) = ua(:,:,jk) + ( ua_b(:,:) - ub_b(:,:) ) * z1_2dt_b * wdmask(:,:) 999 va(:,:,jk) = va(:,:,jk) + ( va_b(:,:) - vb_b(:,:) ) * z1_2dt_b * wdmask(:,:) 1000 END DO 1001 ELSE 1002 ! At this stage, ssha has been corrected: compute new depths at velocity points 1003 DO jj = 1, jpjm1 1004 DO ji = 1, jpim1 ! NO Vector Opt. 1005 zsshu_a(ji,jj) = z1_2 * umask(ji,jj,1) * r1_e1e2u(ji,jj) & 1006 & * ( e1e2t(ji ,jj) * ssha(ji ,jj) & 1007 & + e1e2t(ji+1,jj) * ssha(ji+1,jj) ) 1008 zsshv_a(ji,jj) = z1_2 * vmask(ji,jj,1) * r1_e1e2v(ji,jj) & 1009 & * ( e1e2t(ji,jj ) * ssha(ji,jj ) & 1010 & + e1e2t(ji,jj+1) * ssha(ji,jj+1) ) 1011 END DO 1012 END DO 1013 CALL lbc_lnk_multi( zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions 1014 ! 1015 DO jk=1,jpkm1 1016 ua(:,:,jk) = ua(:,:,jk) + r1_hu_n(:,:) * ( ua_b(:,:) - ub_b(:,:) * hu_b(:,:) ) * z1_2dt_b * wdmask(:,:) 1017 va(:,:,jk) = va(:,:,jk) + r1_hv_n(:,:) * ( va_b(:,:) - vb_b(:,:) * hv_b(:,:) ) * z1_2dt_b * wdmask(:,:) 1018 END DO 1019 ! Save barotropic velocities not transport: 1020 ua_b(:,:) = ua_b(:,:) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - ssumask(:,:) ) 1021 va_b(:,:) = va_b(:,:) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - ssvmask(:,:) ) 1022 ENDIF 1027 1023 ELSE 1028 ! At this stage, ssha has been corrected: compute new depths at velocity points 1029 DO jj = 1, jpjm1 1030 DO ji = 1, jpim1 ! NO Vector Opt. 1031 zsshu_a(ji,jj) = z1_2 * umask(ji,jj,1) * r1_e1e2u(ji,jj) & 1032 & * ( e1e2t(ji ,jj) * ssha(ji ,jj) & 1033 & + e1e2t(ji+1,jj) * ssha(ji+1,jj) ) 1034 zsshv_a(ji,jj) = z1_2 * vmask(ji,jj,1) * r1_e1e2v(ji,jj) & 1035 & * ( e1e2t(ji,jj ) * ssha(ji,jj ) & 1036 & + e1e2t(ji,jj+1) * ssha(ji,jj+1) ) 1037 END DO 1038 END DO 1039 CALL lbc_lnk_multi( zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions 1040 ! 1041 DO jk=1,jpkm1 1042 ua(:,:,jk) = ua(:,:,jk) + r1_hu_n(:,:) * ( ua_b(:,:) - ub_b(:,:) * hu_b(:,:) ) * z1_2dt_b 1043 va(:,:,jk) = va(:,:,jk) + r1_hv_n(:,:) * ( va_b(:,:) - vb_b(:,:) * hv_b(:,:) ) * z1_2dt_b 1044 END DO 1045 ! Save barotropic velocities not transport: 1046 ua_b(:,:) = ua_b(:,:) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - ssumask(:,:) ) 1047 va_b(:,:) = va_b(:,:) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - ssvmask(:,:) ) 1048 ENDIF 1024 IF( ln_dynadv_vec .OR. ln_linssh ) THEN 1025 DO jk=1,jpkm1 1026 ua(:,:,jk) = ua(:,:,jk) + ( ua_b(:,:) - ub_b(:,:) ) * z1_2dt_b 1027 va(:,:,jk) = va(:,:,jk) + ( va_b(:,:) - vb_b(:,:) ) * z1_2dt_b 1028 END DO 1029 ELSE 1030 ! At this stage, ssha has been corrected: compute new depths at velocity points 1031 DO jj = 1, jpjm1 1032 DO ji = 1, jpim1 ! NO Vector Opt. 1033 zsshu_a(ji,jj) = z1_2 * umask(ji,jj,1) * r1_e1e2u(ji,jj) & 1034 & * ( e1e2t(ji ,jj) * ssha(ji ,jj) & 1035 & + e1e2t(ji+1,jj) * ssha(ji+1,jj) ) 1036 zsshv_a(ji,jj) = z1_2 * vmask(ji,jj,1) * r1_e1e2v(ji,jj) & 1037 & * ( e1e2t(ji,jj ) * ssha(ji,jj ) & 1038 & + e1e2t(ji,jj+1) * ssha(ji,jj+1) ) 1039 END DO 1040 END DO 1041 CALL lbc_lnk_multi( zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions 1042 ! 1043 DO jk=1,jpkm1 1044 ua(:,:,jk) = ua(:,:,jk) + r1_hu_n(:,:) * ( ua_b(:,:) - ub_b(:,:) * hu_b(:,:) ) * z1_2dt_b 1045 va(:,:,jk) = va(:,:,jk) + r1_hv_n(:,:) * ( va_b(:,:) - vb_b(:,:) * hv_b(:,:) ) * z1_2dt_b 1046 END DO 1047 ! Save barotropic velocities not transport: 1048 ua_b(:,:) = ua_b(:,:) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - ssumask(:,:) ) 1049 va_b(:,:) = va_b(:,:) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - ssvmask(:,:) ) 1050 ENDIF 1051 1052 END IF 1049 1053 ! 1050 1054 DO jk = 1, jpkm1 … … 1082 1086 CALL wrk_dealloc( jpi,jpj, zsshu_a, zsshv_a ) 1083 1087 CALL wrk_dealloc( jpi,jpj, zhf ) 1084 IF( ln_wd ) CALL wrk_dealloc( jpi, jpj, zcpx, zcpy , wduflt1, wdvflt1)1088 IF( ln_wd ) CALL wrk_dealloc( jpi, jpj, zcpx, zcpy ) 1085 1089 ! 1086 1090 IF ( ln_diatmb ) THEN -
branches/2016/dev_NOC_CMCC_merge_2016/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90
r7299 r7397 87 87 ENDIF 88 88 ! 89 CALL div_hor( kt ) ! Horizontal divergence 90 ! 91 z2dt = 2._wp * rdt ! set time step size (Euler/Leapfrog) 89 z2dt = 2._wp * rdt ! set time step size (Euler/Leapfrog) 92 90 IF( neuler == 0 .AND. kt == nit000 ) z2dt = rdt 91 zcoef = 0.5_wp * r1_rau0 93 92 94 93 ! !------------------------------! 95 94 ! ! After Sea Surface Height ! 96 95 ! !------------------------------! 96 IF(ln_wd) THEN 97 CALL wad_lmt(sshb, zcoef * (emp_b(:,:) + emp(:,:)), z2dt) 98 ENDIF 99 100 CALL div_hor( kt ) ! Horizontal divergence 101 ! 97 102 zhdiv(:,:) = 0._wp 98 103 DO jk = 1, jpkm1 ! Horizontal divergence of barotropic transports … … 103 108 ! compute the vertical velocity which can be used to compute the non-linear terms of the momentum equations. 104 109 ! 105 zcoef = 0.5_wp * r1_rau0106 107 IF(ln_wd) CALL wad_lmt(sshb, zcoef * (emp_b(:,:) + emp(:,:)), z2dt)108 109 110 ssha(:,:) = ( sshb(:,:) - z2dt * ( zcoef * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) ) ) * ssmask(:,:) 110 111 -
branches/2016/dev_NOC_CMCC_merge_2016/NEMOGCM/NEMO/OPA_SRC/DYN/wet_dry.F90
r6152 r7397 33 33 !! --------------------------------------------------------------------- 34 34 35 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: wduflt, wdvflt !: u- and v- filter36 35 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: wdmask !: u- and v- limiter 37 36 … … 46 45 PUBLIC wad_lmt ! routine called by sshwzv.F90 47 46 PUBLIC wad_lmt_bt ! routine called by dynspg_ts.F90 47 PUBLIC wad_istate ! routine called by istate.F90 and domvvl.F90 48 48 49 49 !! * Substitutions … … 87 87 88 88 IF(ln_wd) THEN 89 ALLOCATE( wd uflt(jpi,jpj), wdvflt(jpi,jpj), wdmask(jpi,jpj), STAT=ierr )89 ALLOCATE( wdmask(jpi,jpj), STAT=ierr ) 90 90 IF( ierr /= 0 ) CALL ctl_stop('STOP', 'wad_init : Array allocation error') 91 91 ENDIF … … 145 145 ! Horizontal Flux in u and v direction 146 146 DO jk = 1, jpkm1 147 DO jj = 1, jpj m1148 DO ji = 1, jpi m1147 DO jj = 1, jpj 148 DO ji = 1, jpi 149 149 zflxu(ji,jj) = zflxu(ji,jj) + e3u_n(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk) 150 150 zflxv(ji,jj) = zflxv(ji,jj) + e3v_n(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk) … … 156 156 zflxv(:,:) = zflxv(:,:) * e1v(:,:) 157 157 158 DO jj = 2, jpjm1 159 DO ji = 2, jpim1 158 wdmask(:,:) = 1 159 DO jj = 2, jpj 160 DO ji = 2, jpi 160 161 161 162 IF(tmask(ji, jj, 1) < 0.5_wp) CYCLE ! we don't care about land cells … … 168 169 169 170 zdep2 = bathy(ji,jj) + sshb1(ji,jj) - rn_wdmin1 170 IF(zdep2 <0._wp) THEN !add more safty, but not necessary171 IF(zdep2 .le. 0._wp) THEN !add more safty, but not necessary 171 172 !zdep2 = 0._wp 172 173 sshb1(ji,jj) = rn_wdmin1 - bathy(ji,jj) 174 wdmask(ji,jj) = 0._wp 173 175 END IF 174 176 ENDDO … … 183 185 zflxv1(:,:) = zflxv(:,:) * zwdlmtv(:,:) 184 186 185 DO jj = 2, jpj m1186 DO ji = 2, jpi m1187 DO jj = 2, jpj 188 DO ji = 2, jpi 187 189 188 wdmask(ji,jj) = 0189 190 IF(tmask(ji, jj, 1) < 0.5_wp) CYCLE 190 191 IF(bathy(ji,jj) > zdepwd) CYCLE … … 202 203 IF(zdep1 > zdep2) THEN 203 204 zflag = 1 204 wdmask(ji, jj) = 1205 wdmask(ji, jj) = 0 205 206 zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zflxp(ji,jj) * z2dt ) 206 207 zcoef = max(zcoef, 0._wp) … … 209 210 IF(zflxu1(ji-1,jj) < 0._wp) zwdlmtu(ji-1,jj) = zcoef 210 211 IF(zflxv1(ji, jj) > 0._wp) zwdlmtv(ji ,jj) = zcoef 211 IF(zflxv1(ji,jj-1) < 0._wp) zwdlmtv(ji -1,jj) = zcoef212 IF(zflxv1(ji,jj-1) < 0._wp) zwdlmtv(ji,jj-1) = zcoef 212 213 END IF 213 214 END DO ! ji loop … … 231 232 CALL lbc_lnk( un, 'U', -1. ) 232 233 CALL lbc_lnk( vn, 'V', -1. ) 234 ! 235 un_b(:,:) = un_b(:,:) * zwdlmtu(:, :) 236 vn_b(:,:) = vn_b(:,:) * zwdlmtv(:, :) 237 CALL lbc_lnk( un_b, 'U', -1. ) 238 CALL lbc_lnk( vn_b, 'V', -1. ) 233 239 234 240 IF(zflag == 1 .AND. lwp) WRITE(numout,*) 'Need more iterations in wad_lmt!!!' … … 291 297 zflxp(:,:) = 0._wp 292 298 zflxn(:,:) = 0._wp 293 !zflxu(:,:) = 0._wp294 !zflxv(:,:) = 0._wp295 299 296 300 zwdlmtu(:,:) = 1._wp … … 299 303 ! Horizontal Flux in u and v direction 300 304 301 !zflxu(:,:) = zflxu(:,:) * e2u(:,:) 302 !zflxv(:,:) = zflxv(:,:) * e1v(:,:) 303 304 DO jj = 2, jpjm1 305 DO ji = 2, jpim1 305 DO jj = 2, jpj 306 DO ji = 2, jpi 306 307 307 308 IF(tmask(ji, jj, 1) < 0.5_wp) CYCLE ! we don't care about land cells … … 314 315 315 316 zdep2 = bathy(ji,jj) + sshn_e(ji,jj) - rn_wdmin1 316 IF(zdep2 < 0._wp) THEN !add more safty, but not necessary317 !zdep2 = 0._wp318 sshn_e(ji,jj) = rn_wdmin1 - bathy(ji,jj)319 END IF320 317 ENDDO 321 318 END DO … … 329 326 zflxv1(:,:) = zflxv(:,:) * zwdlmtv(:,:) 330 327 331 DO jj = 2, jpj m1332 DO ji = 2, jpi m1328 DO jj = 2, jpj 329 DO ji = 2, jpi 333 330 334 wdmask(ji,jj) = 0335 331 IF(tmask(ji, jj, 1) < 0.5_wp) CYCLE 336 332 IF(bathy(ji,jj) > zdepwd) CYCLE … … 349 345 IF(zdep1 > zdep2) THEN 350 346 zflag = 1 351 !wdmask(ji, jj) = 1352 347 zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zflxp(ji,jj) * z2dt ) 353 348 zcoef = max(zcoef, 0._wp) … … 356 351 IF(zflxu1(ji-1,jj) < 0._wp) zwdlmtu(ji-1,jj) = zcoef 357 352 IF(zflxv1(ji, jj) > 0._wp) zwdlmtv(ji ,jj) = zcoef 358 IF(zflxv1(ji,jj-1) < 0._wp) zwdlmtv(ji -1,jj) = zcoef353 IF(zflxv1(ji,jj-1) < 0._wp) zwdlmtv(ji,jj-1) = zcoef 359 354 END IF 360 355 END DO ! ji loop … … 379 374 IF(zflag == 1 .AND. lwp) WRITE(numout,*) 'Need more iterations in wad_lmt_bt!!!' 380 375 381 !IF( ln_rnf ) CALL sbc_rnf_div( hdivn ) ! runoffs (update hdivn field)382 !IF( nn_cla == 1 ) CALL cla_div ( kt ) ! Cross Land Advection (update hdivn field)383 376 ! 384 377 ! … … 390 383 IF( nn_timing == 1 ) CALL timing_stop('wad_lmt') 391 384 END SUBROUTINE wad_lmt_bt 385 386 SUBROUTINE wad_istate 387 !!---------------------------------------------------------------------- 388 !! *** ROUTINE wad_istate *** 389 !! 390 !! ** Purpose : Initialization of the dynamics and tracers for WAD test 391 !! configurations (channels or bowls with initial ssh gradients) 392 !! 393 !! ** Method : - set temperature field 394 !! - set salinity field 395 !! - set ssh slope (needs to be repeated in domvvl_rst_init to 396 !! set vertical metrics ) 397 !!---------------------------------------------------------------------- 398 ! 399 INTEGER :: ji, jj ! dummy loop indices 400 REAL(wp) :: zi, zj 401 !!---------------------------------------------------------------------- 402 ! 403 ! Uniform T & S in all test cases 404 tsn(:,:,:,jp_tem) = 10._wp 405 tsb(:,:,:,jp_tem) = 10._wp 406 tsn(:,:,:,jp_sal) = 35._wp 407 tsb(:,:,:,jp_sal) = 35._wp 408 SELECT CASE ( jp_cfg ) 409 ! ! ==================== 410 CASE ( 1 ) ! WAD 1 configuration 411 ! ! ==================== 412 ! 413 IF(lwp) WRITE(numout,*) 414 IF(lwp) WRITE(numout,*) 'istate_wad : Closed box with EW linear bottom slope' 415 IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 416 ! 417 do ji = 1,jpi 418 sshn(ji,:) = ( -5.5_wp + 5.5_wp*FLOAT(mig(ji))/FLOAT(jpidta-1))*tmask(ji,:,1) 419 end do 420 ! ! ==================== 421 CASE ( 2 ) ! WAD 2 configuration 422 ! ! ==================== 423 ! 424 IF(lwp) WRITE(numout,*) 425 IF(lwp) WRITE(numout,*) 'istate_wad : Parobolic EW channel, mid-range initial ssh slope' 426 IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 427 ! 428 do ji = 1,jpi 429 sshn(ji,:) = ( -5.5_wp + 3.9_wp*FLOAT(jpidta - mig(ji))/FLOAT(jpidta-1))*tmask(ji,:,1) 430 end do 431 ! ! ==================== 432 CASE ( 3 ) ! WAD 3 configuration 433 ! ! ==================== 434 ! 435 IF(lwp) WRITE(numout,*) 436 IF(lwp) WRITE(numout,*) 'istate_wad : Parobolic EW channel, extreme initial ssh slope' 437 IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 438 ! 439 do ji = 1,jpi 440 sshn(ji,:) = ( -7.5_wp + 6.9_wp*FLOAT(jpidta - mig(ji))/FLOAT(jpidta-1))*tmask(ji,:,1) 441 end do 442 443 ! 444 ! ! ==================== 445 CASE ( 4 ) ! WAD 4 configuration 446 ! ! ==================== 447 ! 448 IF(lwp) WRITE(numout,*) 449 IF(lwp) WRITE(numout,*) 'istate_wad : Parobolic bowl, mid-range initial ssh slope' 450 IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 451 ! 452 DO ji = 1, jpi 453 zi = MAX(1.0-FLOAT((mig(ji)-25)**2)/400.0, 0.0 ) 454 DO jj = 1, jpj 455 zj = MAX(1.0-FLOAT((mjg(jj)-17)**2)/144.0, 0.0 ) 456 sshn(ji,jj) = -8.5_wp + 8.5_wp*zi*zj 457 END DO 458 END DO 459 460 ! 461 ! ! =========================== 462 CASE ( 5 ) ! WAD 5 configuration 463 ! ! ==================== 464 ! 465 IF(lwp) WRITE(numout,*) 466 IF(lwp) WRITE(numout,*) 'istate_wad : Double slope with shelf' 467 IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 468 ! 469 ! Needed rn_wdmin2 increased to 0.01 for this case? 470 do ji = 1,jpi 471 sshn(ji,:) = ( -5.5_wp + 9.0_wp*FLOAT(mig(ji))/FLOAT(jpidta-1))*tmask(ji,:,1) 472 end do 473 474 ! 475 ! ! =========================== 476 CASE ( 6 ) ! WAD 6 configuration 477 ! ! ==================== 478 ! 479 IF(lwp) WRITE(numout,*) 480 IF(lwp) WRITE(numout,*) 'istate_wad : Parobolic EW channel with gaussian ridge' 481 IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 482 ! 483 do ji = 1,jpi 484 !6a 485 sshn(ji,:) = ( -5.5_wp + 9.0_wp*FLOAT(jpidta - mig(ji))/FLOAT(jpidta-1))*tmask(ji,:,1) 486 !Some variations in initial slope that have been tested 487 !6b 488 !sshn(ji,:) = ( -5.5_wp + 6.5_wp*FLOAT(jpidta - mig(ji))/FLOAT(jpidta-1))*tmask(ji,:,1) 489 !6c 490 !sshn(ji,:) = ( -5.5_wp + 7.5_wp*FLOAT(jpidta - mig(ji))/FLOAT(jpidta-1))*tmask(ji,:,1) 491 !6d 492 !sshn(ji,:) = ( -4.5_wp + 8.0_wp*FLOAT(jpidta - mig(ji))/FLOAT(jpidta-1))*tmask(ji,:,1) 493 end do 494 495 ! 496 ! ! =========================== 497 CASE DEFAULT ! NONE existing configuration 498 ! ! =========================== 499 WRITE(ctmp1,*) 'WAD test with a ', jp_cfg,' option is not coded' 500 ! 501 CALL ctl_stop( ctmp1 ) 502 ! 503 END SELECT 504 ! 505 ! Apply minimum wetdepth criterion 506 ! 507 do jj = 1,jpj 508 do ji = 1,jpi 509 IF( bathy(ji,jj) + sshn(ji,jj) < rn_wdmin1 ) THEN 510 sshn(ji,jj) = tmask(ji,jj,1)*( rn_wdmin1 - bathy(ji,jj) ) 511 ENDIF 512 end do 513 end do 514 sshb = sshn 515 ssha = sshn 516 ! 517 END SUBROUTINE wad_istate 518 519 !!===================================================================== 392 520 END MODULE wet_dry -
branches/2016/dev_NOC_CMCC_merge_2016/NEMOGCM/NEMO/OPA_SRC/IOM/iom.F90
r6519 r7397 9 9 !! 3.4 ! 2012-12 (R. Bourdalle-Badie and G. Reffray) add C1D case 10 10 !! 3.6 ! 2014-15 DIMG format removed 11 !! 3.6 ! 2015-15 (J. Harle) Added procedure to read REAL attributes 11 12 !!-------------------------------------------------------------------- 12 13 … … 67 68 END INTERFACE 68 69 INTERFACE iom_getatt 69 MODULE PROCEDURE iom_g0d_intatt 70 MODULE PROCEDURE iom_g0d_intatt, iom_g0d_ratt 70 71 END INTERFACE 71 72 INTERFACE iom_rstput … … 979 980 IF( iom_file(kiomid)%nfid > 0 ) THEN 980 981 SELECT CASE (iom_file(kiomid)%iolib) 981 CASE (jpnf90 ) ; CALL iom_nf90_getatt( kiomid, cdatt, pv ar )982 CASE (jpnf90 ) ; CALL iom_nf90_getatt( kiomid, cdatt, pv_i0d=pvar ) 982 983 CASE DEFAULT 983 984 CALL ctl_stop( 'iom_g0d_att: accepted IO library is only jpnf90' ) … … 987 988 END SUBROUTINE iom_g0d_intatt 988 989 990 SUBROUTINE iom_g0d_ratt( kiomid, cdatt, pvar, cdvar ) 991 INTEGER , INTENT(in ) :: kiomid ! Identifier of the file 992 CHARACTER(len=*), INTENT(in ) :: cdatt ! Name of the attribute 993 REAL(wp) , INTENT( out) :: pvar ! written field 994 CHARACTER(len=*), INTENT(in ), OPTIONAL :: cdvar ! Name of the variable 995 ! 996 IF( kiomid > 0 ) THEN 997 IF( iom_file(kiomid)%nfid > 0 ) THEN 998 SELECT CASE (iom_file(kiomid)%iolib) 999 CASE (jpnf90 ) ; IF( PRESENT(cdvar) ) THEN 1000 CALL iom_nf90_getatt( kiomid, cdatt, pv_r0d=pvar, cdvar=cdvar ) 1001 ELSE 1002 CALL iom_nf90_getatt( kiomid, cdatt, pv_r0d=pvar ) 1003 ENDIF 1004 CASE DEFAULT 1005 CALL ctl_stop( 'iom_g0d_att: accepted IO library is only jpnf90' ) 1006 END SELECT 1007 ENDIF 1008 ENDIF 1009 END SUBROUTINE iom_g0d_ratt 989 1010 990 1011 !!---------------------------------------------------------------------- -
branches/2016/dev_NOC_CMCC_merge_2016/NEMOGCM/NEMO/OPA_SRC/IOM/iom_nf90.F90
r6140 r7397 7 7 !! 9.0 ! 06 02 (S. Masson) Adaptation to NEMO 8 8 !! " ! 07 07 (D. Storkey) Changes to iom_nf90_gettime 9 !! 3.6 ! 2015-15 (J. Harle) Added procedure to read REAL attributes 9 10 !!-------------------------------------------------------------------- 10 11 !!gm caution add !DIR nec: improved performance to be checked as well as no result changes … … 35 36 END INTERFACE 36 37 INTERFACE iom_nf90_getatt 37 MODULE PROCEDURE iom_nf90_ intatt38 MODULE PROCEDURE iom_nf90_att 38 39 END INTERFACE 39 40 INTERFACE iom_nf90_rstput … … 313 314 314 315 315 SUBROUTINE iom_nf90_ intatt( kiomid, cdatt, pvar)316 !!----------------------------------------------------------------------- 317 !! *** ROUTINE iom_nf90_ intatt ***316 SUBROUTINE iom_nf90_att( kiomid, cdatt, pv_i0d, pv_r0d, cdvar) 317 !!----------------------------------------------------------------------- 318 !! *** ROUTINE iom_nf90_att *** 318 319 !! 319 320 !! ** Purpose : read an integer attribute with NF90 … … 321 322 INTEGER , INTENT(in ) :: kiomid ! Identifier of the file 322 323 CHARACTER(len=*), INTENT(in ) :: cdatt ! attribute name 323 INTEGER , INTENT( out) :: pvar ! read field 324 INTEGER , INTENT( out), OPTIONAL :: pv_i0d ! read field 325 REAL(wp), INTENT( out), OPTIONAL :: pv_r0d ! read field 326 CHARACTER(len=*), INTENT(in ), OPTIONAL :: cdvar ! name of the variable 324 327 ! 325 328 INTEGER :: if90id ! temporary integer 329 INTEGER :: ivarid ! NetCDF variable Id 326 330 LOGICAL :: llok ! temporary logical 327 331 CHARACTER(LEN=100) :: clinfo ! info character … … 329 333 ! 330 334 if90id = iom_file(kiomid)%nfid 331 llok = NF90_Inquire_attribute(if90id, NF90_GLOBAL, cdatt) == nf90_noerr 335 IF( PRESENT(cdvar) ) THEN 336 llok = NF90_INQ_VARID( if90id, TRIM(cdvar), ivarid ) == nf90_noerr ! does the variable exist in the file 337 IF( llok ) THEN 338 llok = NF90_Inquire_attribute(if90id, ivarid, cdatt) == nf90_noerr 339 ELSE 340 CALL ctl_warn('iom_nf90_getatt: no variable '//cdvar//' found') 341 ENDIF 342 ELSE 343 llok = NF90_Inquire_attribute(if90id, NF90_GLOBAL, cdatt) == nf90_noerr 344 ENDIF 345 ! 332 346 IF( llok) THEN 333 347 clinfo = 'iom_nf90_getatt, file: '//TRIM(iom_file(kiomid)%name)//', att: '//TRIM(cdatt) 334 CALL iom_nf90_check(NF90_GET_ATT(if90id, NF90_GLOBAL, cdatt, values=pvar), clinfo) 348 IF( PRESENT(pv_r0d) ) THEN 349 IF( PRESENT(cdvar) ) THEN 350 CALL iom_nf90_check(NF90_GET_ATT(if90id, ivarid, cdatt, values=pv_r0d), clinfo) 351 ELSE 352 CALL iom_nf90_check(NF90_GET_ATT(if90id, NF90_GLOBAL, cdatt, values=pv_r0d), clinfo) 353 ENDIF 354 ELSE 355 IF( PRESENT(cdvar) ) THEN 356 CALL iom_nf90_check(NF90_GET_ATT(if90id, ivarid, cdatt, values=pv_i0d), clinfo) 357 ELSE 358 CALL iom_nf90_check(NF90_GET_ATT(if90id, NF90_GLOBAL, cdatt, values=pv_i0d), clinfo) 359 ENDIF 360 ENDIF 335 361 ELSE 336 362 CALL ctl_warn('iom_nf90_getatt: no attribute '//cdatt//' found') 337 pvar = -999 363 IF( PRESENT(pv_r0d) ) THEN 364 pv_r0d = -999._wp 365 ELSE 366 pv_i0d = -999 367 ENDIF 338 368 ENDIF 339 369 ! 340 END SUBROUTINE iom_nf90_ intatt370 END SUBROUTINE iom_nf90_att 341 371 342 372 -
branches/2016/dev_NOC_CMCC_merge_2016/NEMOGCM/NEMO/OPA_SRC/LBC/mppini_2.h90
r6412 r7397 41 41 USE in_out_manager ! I/O Manager 42 42 USE iom 43 #if defined key_bdy 44 USE bdy_oce 45 #endif 43 46 !! 44 47 INTEGER :: ji, jj, jn, jproc, jarea ! dummy loop indices … … 73 76 ! read namelist for ln_zco 74 77 NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco, ln_isfcav, ln_linssh 75 78 #if defined key_bdy 79 NAMELIST/nambdy/ nb_bdy, ln_coords_file, cn_coords_file, & 80 & ln_mask_file, cn_mask_file, cn_dyn2d, nn_dyn2d_dta, & 81 & cn_dyn3d, nn_dyn3d_dta, cn_tra, nn_tra_dta, & 82 & ln_tra_dmp, ln_dyn3d_dmp, rn_time_dmp, rn_time_dmp_out, & 83 & cn_ice_lim, nn_ice_lim_dta, & 84 & rn_ice_tem, rn_ice_sal, rn_ice_age, & 85 & ln_vol, nn_volctl, nn_rimwidth, nb_jpk_bdy 86 #endif 76 87 !!---------------------------------------------------------------------- 77 88 !! OPA 9.0 , LOCEAN-IPSL (2005) … … 137 148 imask(:,:)=1 138 149 WHERE ( zdta(:,:) - zdtaisf(:,:) <= rn_isfhmin ) imask = 0 150 151 #if defined key_bdy 152 ! Adjust imask with bdy_msk if exists 153 154 REWIND( numnam_ref ) ! Namelist nambdy in reference namelist : BDY 155 READ ( numnam_ref, nambdy, IOSTAT = ios, ERR = 903) 156 903 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambdy in reference namelist (mppini_2)', lwp ) 157 158 REWIND( numnam_cfg ) ! Namelist nambdy in configuration namelist : BDY 159 READ ( numnam_cfg, nambdy, IOSTAT = ios, ERR = 904 ) 160 904 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambdy in configuration namelist (mppini_2)', lwp ) 161 IF(lwm) WRITE ( numond, namzgr ) 162 163 IF( ln_mask_file ) THEN 164 CALL iom_open( cn_mask_file, inum ) 165 CALL iom_get ( inum, jpdom_unknown, 'bdy_msk', zdta(:,:), kstart=(/jpizoom,jpjzoom/), kcount=(/jpiglo,jpjglo/) ) 166 CALL iom_close( inum ) 167 WHERE ( zdta(:,:) <= 0. ) imask = 0 168 ENDIF 169 #endif 139 170 140 171 ! 1. Dimension arrays for subdomains -
branches/2016/dev_NOC_CMCC_merge_2016/NEMOGCM/NEMO/OPA_SRC/LDF/ldfdyn.F90
r6140 r7397 42 42 REAL(wp), PUBLIC :: rn_ahm_b !: lateral laplacian background eddy viscosity [m2/s] 43 43 REAL(wp), PUBLIC :: rn_bhm_0 !: lateral bilaplacian eddy viscosity [m4/s] 44 !! If nn_ahm_ijk_t = 32 a time and space varying Smagorinsky viscosity 45 !! will be computed. 46 REAL(wp), PUBLIC :: rn_csmc !: Smagorinsky constant of proportionality 47 REAL(wp), PUBLIC :: rn_minfac !: Multiplicative factor of theorectical minimum Smagorinsky viscosity 48 REAL(wp), PUBLIC :: rn_maxfac !: Multiplicative factor of theorectical maximum Smagorinsky viscosity 44 49 45 50 LOGICAL , PUBLIC :: l_ldfdyn_time !: flag for time variation of the lateral eddy viscosity coef. 46 51 47 52 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ahmt, ahmf !: eddy diffusivity coef. at U- and V-points [m2/s or m4/s] 53 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: dtensq !: horizontal tension squared (Smagorinsky only) 54 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: dshesq !: horizontal shearing strain squared (Smagorinsky only) 55 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: esqt, esqf !: Square of the local gridscale (e1e2/(e1+e2))**2 48 56 49 57 REAL(wp) :: r1_12 = 1._wp / 12._wp ! =1/12 50 58 REAL(wp) :: r1_4 = 0.25_wp ! =1/4 59 REAL(wp) :: r1_8 = 0.125_wp ! =1/8 51 60 REAL(wp) :: r1_288 = 1._wp / 288._wp ! =1/( 12^2 * 2 ) 52 61 … … 79 88 !! = 31 = F(i,j,k,t) = F(local velocity) ( |u|e /12 laplacian operator 80 89 !! or |u|e^3/12 bilaplacian operator ) 81 !!---------------------------------------------------------------------- 82 INTEGER :: jk ! dummy loop indices 90 !! = 32 = F(i,j,k,t) = F(local deformation rate and gridscale) (D and L) (Smagorinsky) 91 !! ( L^2|D| laplacian operator 92 !! or L^4|D|/8 bilaplacian operator ) 93 !!---------------------------------------------------------------------- 94 INTEGER :: ji, jj, jk ! dummy loop indices 83 95 INTEGER :: ierr, inum, ios ! local integer 84 96 REAL(wp) :: zah0 ! local scalar … … 86 98 NAMELIST/namdyn_ldf/ ln_dynldf_lap, ln_dynldf_blp, & 87 99 & ln_dynldf_lev, ln_dynldf_hor, ln_dynldf_iso, & 88 & nn_ahm_ijk_t , rn_ahm_0, rn_ahm_b, rn_bhm_0 100 & nn_ahm_ijk_t , rn_ahm_0, rn_ahm_b, rn_bhm_0, & 101 & rn_csmc , rn_minfac, rn_maxfac 89 102 !!---------------------------------------------------------------------- 90 103 ! … … 115 128 WRITE(numout,*) ' coefficients :' 116 129 WRITE(numout,*) ' type of time-space variation nn_ahm_ijk_t = ', nn_ahm_ijk_t 117 WRITE(numout,*) ' lateral laplacian eddy viscosity rn_ahm_0 _lap= ', rn_ahm_0, ' m2/s'130 WRITE(numout,*) ' lateral laplacian eddy viscosity rn_ahm_0 = ', rn_ahm_0, ' m2/s' 118 131 WRITE(numout,*) ' background viscosity (iso case) rn_ahm_b = ', rn_ahm_b, ' m2/s' 119 WRITE(numout,*) ' lateral bilaplacian eddy viscosity rn_ahm_0_blp = ', rn_bhm_0, ' m4/s' 132 WRITE(numout,*) ' lateral bilaplacian eddy viscosity rn_bhm_0 = ', rn_bhm_0, ' m4/s' 133 WRITE(numout,*) ' smagorinsky settings (nn_ahm_ijk_t = 32) :' 134 WRITE(numout,*) ' Smagorinsky coefficient rn_csmc = ', rn_csmc 135 WRITE(numout,*) ' factor multiplier for theorectical lower limit for ' 136 WRITE(numout,*) ' Smagorinsky eddy visc (def. 1.0) rn_minfac = ', rn_minfac 137 WRITE(numout,*) ' factor multiplier for theorectical lower upper for ' 138 WRITE(numout,*) ' Smagorinsky eddy visc (def. 1.0) rn_maxfac = ', rn_maxfac 120 139 ENDIF 121 140 … … 208 227 l_ldfdyn_time = .TRUE. ! will be calculated by call to ldf_dyn routine in step.F90 209 228 ! 229 CASE( 32 ) !== time varying 3D field ==! 230 IF(lwp) WRITE(numout,*) ' momentum mixing coef. = F( latitude, longitude, depth , time )' 231 IF(lwp) WRITE(numout,*) ' proportional to the local deformation rate and gridscale (Smagorinsky)' 232 IF(lwp) WRITE(numout,*) ' : L^2|D| or L^4|D|/8' 233 ! 234 l_ldfdyn_time = .TRUE. ! will be calculated by call to ldf_dyn routine in step.F90 235 ! 236 ! allocate arrays used in ldf_dyn. 237 ALLOCATE( dtensq(jpi,jpj) , dshesq(jpi,jpj) , esqt(jpi,jpj) , esqf(jpi,jpj) , STAT=ierr ) 238 IF( ierr /= 0 ) CALL ctl_stop( 'STOP', 'ldf_dyn_init: failed to allocate Smagorinsky arrays') 239 ! 240 ! Set local gridscale values 241 DO jj = 2, jpjm1 242 DO ji = fs_2, fs_jpim1 243 esqt(ji,jj) = ( e1e2t(ji,jj) /( e1t(ji,jj) + e2t(ji,jj) ) )**2 244 esqf(ji,jj) = ( e1e2f(ji,jj) /( e1f(ji,jj) + e2f(ji,jj) ) )**2 245 END DO 246 END DO 247 ! 210 248 CASE DEFAULT 211 249 CALL ctl_stop('ldf_dyn_init: wrong choice for nn_ahm_ijk_t, the type of space-time variation of ahm') … … 232 270 !! nn_ahm_ijk_t = 31 ahmt, ahmf = F(i,j,k,t) = F(local velocity) 233 271 !! ( |u|e /12 or |u|e^3/12 for laplacian or bilaplacian operator ) 234 !! BLP case : sqrt of the eddy coef, since bilaplacian is en re-entrant laplacian235 272 !! 236 !! ** action : ahmt, ahmf update at each time step 273 !! nn_ahm_ijk_t = 32 ahmt, ahmf = F(i,j,k,t) = F(local deformation rate and gridscale) (D and L) (Smagorinsky) 274 !! ( L^2|D| or L^4|D|/8 for laplacian or bilaplacian operator ) 275 !! 276 !! ** note : in BLP cases the sqrt of the eddy coef is returned, since bilaplacian is en re-entrant laplacian 277 !! ** action : ahmt, ahmf updated at each time step 237 278 !!---------------------------------------------------------------------- 238 279 INTEGER, INTENT(in) :: kt ! time step index … … 240 281 INTEGER :: ji, jj, jk ! dummy loop indices 241 282 REAL(wp) :: zu2pv2_ij_p1, zu2pv2_ij, zu2pv2_ij_m1, zetmax, zefmax ! local scalar 283 REAL(wp) :: zcmsmag, zstabf_lo, zstabf_up, zdelta, zdb ! local scalar 242 284 !!---------------------------------------------------------------------- 243 285 ! … … 248 290 CASE( 31 ) !== time varying 3D field ==! = F( local velocity ) 249 291 ! 250 IF( ln_dynldf_lap ) THEN !laplacian operator : |u| e /12 = |u/144| e292 IF( ln_dynldf_lap ) THEN ! laplacian operator : |u| e /12 = |u/144| e 251 293 DO jk = 1, jpkm1 252 294 DO jj = 2, jpjm1 … … 280 322 CALL lbc_lnk( ahmt, 'T', 1. ) ; CALL lbc_lnk( ahmf, 'F', 1. ) 281 323 ! 324 ! 325 CASE( 32 ) !== time varying 3D field ==! = F( local deformation rate and gridscale ) (Smagorinsky) 326 ! 327 IF( ln_dynldf_lap .OR. ln_dynldf_blp ) THEN ! laplacian operator : (C_smag/pi)^2 L^2 |D| 328 ! 329 zcmsmag = (rn_csmc/rpi)**2 ! (C_smag/pi)^2 330 zstabf_lo = rn_minfac * rn_minfac / ( 2._wp * 4._wp * zcmsmag ) ! lower limit stability factor scaling 331 zstabf_up = rn_maxfac / ( 4._wp * zcmsmag * 2._wp * rdt ) ! upper limit stability factor scaling 332 IF( ln_dynldf_blp ) zstabf_lo = ( 16._wp / 9._wp ) * zstabf_lo ! provide |U|L^3/12 lower limit instead 333 ! ! of |U|L^3/16 in blp case 334 DO jk = 1, jpkm1 335 ! 336 DO jj = 2, jpj 337 DO ji = 2, jpi 338 zdb = ( ( ub(ji,jj,jk) * r1_e2u(ji,jj) - ub(ji-1,jj,jk) * r1_e2u(ji-1,jj) ) & 339 & * r1_e1t(ji,jj) * e2t(ji,jj) & 340 & - ( vb(ji,jj,jk) * r1_e1v(ji,jj) - vb(ji,jj-1,jk) * r1_e1v(ji,jj-1) ) & 341 & * r1_e2t(ji,jj) * e1t(ji,jj) ) * tmask(ji,jj,jk) 342 dtensq(ji,jj) = zdb*zdb 343 END DO 344 END DO 345 ! 346 DO jj = 1, jpjm1 347 DO ji = 1, jpim1 348 zdb = ( ( ub(ji,jj+1,jk) * r1_e1u(ji,jj+1) - ub(ji,jj,jk) * r1_e1u(ji,jj) ) & 349 & * r1_e2f(ji,jj) * e1f(ji,jj) & 350 & + ( vb(ji+1,jj,jk) * r1_e2v(ji+1,jj) - vb(ji,jj,jk) * r1_e2v(ji,jj) ) & 351 & * r1_e1f(ji,jj) * e2f(ji,jj) ) * fmask(ji,jj,jk) 352 dshesq(ji,jj) = zdb*zdb 353 END DO 354 END DO 355 ! 356 DO jj = 2, jpjm1 357 DO ji = fs_2, fs_jpim1 358 ! 359 zu2pv2_ij_p1 = ub(ji ,jj+1,jk) * ub(ji ,jj+1,jk) + vb(ji+1,jj ,jk) * vb(ji+1,jj ,jk) 360 zu2pv2_ij = ub(ji ,jj ,jk) * ub(ji ,jj ,jk) + vb(ji ,jj ,jk) * vb(ji ,jj ,jk) 361 zu2pv2_ij_m1 = ub(ji-1,jj ,jk) * ub(ji-1,jj ,jk) + vb(ji ,jj-1,jk) * vb(ji ,jj-1,jk) 362 ! T-point value 363 zdelta = zcmsmag * esqt(ji,jj) ! L^2 * (C_smag/pi)^2 364 ahmt(ji,jj,jk) = zdelta * sqrt( dtensq(ji,jj) + & 365 & r1_4 * ( dshesq(ji,jj) + dshesq(ji,jj-1) + & 366 & dshesq(ji-1,jj) + dshesq(ji-1,jj-1) ) ) 367 ahmt(ji,jj,jk) = MAX( ahmt(ji,jj,jk), & 368 & SQRT( (zu2pv2_ij + zu2pv2_ij_m1) * zdelta * zstabf_lo ) ) ! Impose lower limit == minfac * |U|L/2 369 ahmt(ji,jj,jk) = MIN( ahmt(ji,jj,jk), zdelta * zstabf_up ) ! Impose upper limit == maxfac * L^2/(4*2dt) 370 ! F-point value 371 zdelta = zcmsmag * esqf(ji,jj) ! L^2 * (C_smag/pi)^2 372 ahmf(ji,jj,jk) = zdelta * sqrt( dshesq(ji,jj) + & 373 & r1_4 * ( dtensq(ji,jj) + dtensq(ji,jj+1) + & 374 & dtensq(ji+1,jj) + dtensq(ji+1,jj+1) ) ) 375 ahmf(ji,jj,jk) = MAX( ahmf(ji,jj,jk), & 376 & SQRT( (zu2pv2_ij + zu2pv2_ij_p1) * zdelta * zstabf_lo ) ) ! Impose lower limit == minfac * |U|L/2 377 ahmf(ji,jj,jk) = MIN( ahmf(ji,jj,jk), zdelta * zstabf_up ) ! Impose upper limit == maxfac * L^2/(4*2dt) 378 ! 379 END DO 380 END DO 381 END DO 382 ! 383 ENDIF 384 ! 385 IF( ln_dynldf_blp ) THEN ! bilaplacian operator : sqrt( (C_smag/pi)^2 L^4 |D|/8) 386 ! = sqrt( A_lap_smag L^2/8 ) 387 ! stability limits already applied to laplacian values 388 ! effective default limits are |U|L^3/12 < B_hm < L^4/(32*2dt) 389 ! 390 DO jk = 1, jpkm1 391 ! 392 DO jj = 2, jpjm1 393 DO ji = fs_2, fs_jpim1 394 ! 395 ahmt(ji,jj,jk) = sqrt( r1_8 * esqt(ji,jj) * ahmt(ji,jj,jk) ) 396 ahmf(ji,jj,jk) = sqrt( r1_8 * esqf(ji,jj) * ahmf(ji,jj,jk) ) 397 ! 398 END DO 399 END DO 400 END DO 401 ! 402 ENDIF 403 ! 404 CALL lbc_lnk( ahmt, 'T', 1. ) ; CALL lbc_lnk( ahmf, 'F', 1. ) 405 ! 282 406 END SELECT 283 407 ! -
branches/2016/dev_NOC_CMCC_merge_2016/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90
r7299 r7397 7 7 !! ! 05-2008 (S. Alderson) Modified for Interpolation in memory from input grid to model grid 8 8 !! ! 10-2013 (D. Delrosso, P. Oddo) suppression of land point prior to interpolation 9 !! ! 12-2015 (J. Harle) Adding BDY on-the-fly interpolation 9 10 !!---------------------------------------------------------------------- 10 11 … … 67 68 INTEGER :: nreclast ! last record to be read in the current file 68 69 CHARACTER(len = 256) :: lsmname ! current name of the NetCDF mask file acting as a key 70 INTEGER :: igrd ! grid type for bdy data 71 INTEGER :: ibdy ! bdy set id number 69 72 END TYPE FLD 70 73 … … 114 117 CONTAINS 115 118 116 SUBROUTINE fld_read( kt, kn_fsbc, sd, map, kit, kt_offset )119 SUBROUTINE fld_read( kt, kn_fsbc, sd, map, kit, kt_offset, jpk_bdy, fvl ) 117 120 !!--------------------------------------------------------------------- 118 121 !! *** ROUTINE fld_read *** … … 135 138 ! ! kt_offset = +1 => fields at "after" time level 136 139 ! ! etc. 140 INTEGER , INTENT(in ), OPTIONAL :: jpk_bdy ! number of vertical levels in the BDY data 141 LOGICAL , INTENT(in ), OPTIONAL :: fvl ! number of vertical levels in the BDY data 142 !! 137 143 INTEGER :: itmp ! local variable 138 144 INTEGER :: imf ! size of the structure sd … … 171 177 DO jf = 1, imf 172 178 IF( PRESENT(map) ) imap = map(jf) 173 CALL fld_init( kn_fsbc, sd(jf), imap ) ! read each before field (put them in after as they will be swapped) 179 IF( PRESENT(jpk_bdy) ) THEN 180 CALL fld_init( kn_fsbc, sd(jf), imap, jpk_bdy, fvl ) ! read each before field (put them in after as they will be swapped) 181 ELSE 182 CALL fld_init( kn_fsbc, sd(jf), imap ) ! read each before field (put them in after as they will be swapped) 183 ENDIF 174 184 END DO 175 185 IF( lwp ) CALL wgt_print() ! control print … … 263 273 264 274 ! read after data 265 CALL fld_get( sd(jf), imap ) 266 275 IF( PRESENT(jpk_bdy) ) THEN 276 CALL fld_get( sd(jf), imap, jpk_bdy, fvl) 277 ELSE 278 CALL fld_get( sd(jf), imap ) 279 ENDIF 267 280 ENDIF ! read new data? 268 281 END DO ! --- end loop over field --- ! … … 302 315 303 316 304 SUBROUTINE fld_init( kn_fsbc, sdjf, map )317 SUBROUTINE fld_init( kn_fsbc, sdjf, map , jpk_bdy, fvl) 305 318 !!--------------------------------------------------------------------- 306 319 !! *** ROUTINE fld_init *** … … 309 322 !! - if time interpolation, read before data 310 323 !!---------------------------------------------------------------------- 311 INTEGER , INTENT(in ) :: kn_fsbc ! sbc computation period (in time step) 312 TYPE(FLD), INTENT(inout) :: sdjf ! input field related variables 313 TYPE(MAP_POINTER),INTENT(in) :: map ! global-to-local mapping indices 324 INTEGER , INTENT(in ) :: kn_fsbc ! sbc computation period (in time step) 325 TYPE(FLD), INTENT(inout) :: sdjf ! input field related variables 326 TYPE(MAP_POINTER),INTENT(in) :: map ! global-to-local mapping indices 327 INTEGER , INTENT(in), OPTIONAL :: jpk_bdy ! number of vertical levels in the BDY data 328 LOGICAL , INTENT(in), OPTIONAL :: fvl ! number of vertical levels in the BDY data 314 329 !! 315 330 LOGICAL :: llprevyr ! are we reading previous year file? … … 405 420 ENDIF 406 421 ! 407 CALL fld_get( sdjf, map ) ! read before data in after arrays(as we will swap it later) 422 ! read before data in after arrays(as we will swap it later) 423 IF( PRESENT(jpk_bdy) ) THEN 424 CALL fld_get( sdjf, map, jpk_bdy, fvl ) 425 ELSE 426 CALL fld_get( sdjf, map ) 427 ENDIF 408 428 ! 409 429 clfmt = "('fld_init : time-interpolation for ', a, ' read previous record = ', i6, ' at time = ', f7.2, ' days')" … … 581 601 582 602 583 SUBROUTINE fld_get( sdjf, map )603 SUBROUTINE fld_get( sdjf, map, jpk_bdy, fvl ) 584 604 !!--------------------------------------------------------------------- 585 605 !! *** ROUTINE fld_get *** … … 589 609 TYPE(FLD) , INTENT(inout) :: sdjf ! input field related variables 590 610 TYPE(MAP_POINTER), INTENT(in ) :: map ! global-to-local mapping indices 611 INTEGER , INTENT(in), OPTIONAL :: jpk_bdy ! number of vertical levels in the bdy data 612 LOGICAL , INTENT(in), OPTIONAL :: fvl ! number of vertical levels in the bdy data 591 613 ! 592 614 INTEGER :: ipk ! number of vertical levels of sdjf%fdta ( 2D: ipk=1 ; 3D: ipk=jpk ) … … 601 623 ! 602 624 IF( ASSOCIATED(map%ptr) ) THEN 603 IF( sdjf%ln_tint ) THEN ; CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fdta(:,:,:,2), sdjf%nrec_a(1), map ) 604 ELSE ; CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fnow(:,:,: ), sdjf%nrec_a(1), map ) 605 ENDIF 625 IF( PRESENT(jpk_bdy) ) THEN 626 IF( sdjf%ln_tint ) THEN ; CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fdta(:,:,:,2), & 627 sdjf%nrec_a(1), map, sdjf%igrd, sdjf%ibdy, jpk_bdy, fvl ) 628 ELSE ; CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fnow(:,:,: ), & 629 sdjf%nrec_a(1), map, sdjf%igrd, sdjf%ibdy, jpk_bdy, fvl ) 630 ENDIF 631 ELSE 632 IF( sdjf%ln_tint ) THEN ; CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fdta(:,:,:,2), sdjf%nrec_a(1), map ) 633 ELSE ; CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fnow(:,:,: ), sdjf%nrec_a(1), map ) 634 ENDIF 635 ENDIF 606 636 ELSE IF( LEN(TRIM(sdjf%wgtname)) > 0 ) THEN 607 637 CALL wgt_list( sdjf, iw ) … … 658 688 END SUBROUTINE fld_get 659 689 660 661 SUBROUTINE fld_map( num, clvar, dta, nrec, map ) 690 SUBROUTINE fld_map( num, clvar, dta, nrec, map, igrd, ibdy, jpk_bdy, fvl ) 662 691 !!--------------------------------------------------------------------- 663 692 !! *** ROUTINE fld_map *** … … 666 695 !! using a general mapping (for open boundaries) 667 696 !!---------------------------------------------------------------------- 668 USE bdy_oce, ONLY: ln_bdy, dta_global, dta_global2 ! workspace to read in global data arrays 697 698 USE bdy_oce, ONLY: ln_bdy, idx_bdy, dta_global, dta_global_z, dta_global_dz, dta_global2, dta_global2_z, dta_global2_dz ! workspace to read in global data arrays 669 699 670 700 INTEGER , INTENT(in ) :: num ! stream number 671 701 CHARACTER(LEN=*) , INTENT(in ) :: clvar ! variable name 672 REAL(wp), DIMENSION(:,:,:), INTENT(out) :: dta ! output field on model grid (2 dimensional)702 REAL(wp), DIMENSION(:,:,:), INTENT(out) :: dta ! output field on model grid (2 dimensional) 673 703 INTEGER , INTENT(in ) :: nrec ! record number to read (ie time slice) 674 704 TYPE(MAP_POINTER) , INTENT(in ) :: map ! global-to-local mapping indices 705 INTEGER , INTENT(in), OPTIONAL :: igrd, ibdy, jpk_bdy ! grid type, set number and number of vertical levels in the bdy data 706 LOGICAL , INTENT(in), OPTIONAL :: fvl ! grid type, set number and number of vertical levels in the bdy data 707 INTEGER :: jpkm1_bdy! number of vertical levels in the bdy data minus 1 675 708 !! 676 709 INTEGER :: ipi ! length of boundary data on local process … … 681 714 INTEGER :: ib, ik, ji, jj ! loop counters 682 715 INTEGER :: ierr 683 REAL(wp), POINTER, DIMENSION(:,:,:) :: dta_read ! work space for global data 716 REAL(wp) :: fv ! fillvalue 717 REAL(wp), POINTER, DIMENSION(:,:,:) :: dta_read ! work space for global data 718 REAL(wp), POINTER, DIMENSION(:,:,:) :: dta_read_z ! work space for global data 719 REAL(wp), POINTER, DIMENSION(:,:,:) :: dta_read_dz ! work space for global data 684 720 !!--------------------------------------------------------------------- 685 721 ! … … 695 731 IF( map%ll_unstruc) THEN ! unstructured open boundary data file 696 732 dta_read => dta_global 697 ELSE ! structured open boundary data file 733 IF( PRESENT(jpk_bdy) ) THEN 734 IF( jpk_bdy>0 ) THEN 735 dta_read_z => dta_global_z 736 dta_read_dz => dta_global_dz 737 jpkm1_bdy = jpk_bdy-1 738 ENDIF 739 ENDIF 740 ELSE ! structured open boundary file 698 741 dta_read => dta_global2 742 IF( PRESENT(jpk_bdy) ) THEN 743 IF( jpk_bdy>0 ) THEN 744 dta_read_z => dta_global2_z 745 dta_read_dz => dta_global2_dz 746 jpkm1_bdy = jpk_bdy-1 747 ENDIF 748 ENDIF 699 749 ENDIF 700 750 ENDIF … … 704 754 ! 705 755 SELECT CASE( ipk ) 706 CASE(1) ; CALL iom_get ( num, jpdom_unknown, clvar, dta_read(1:ilendta,1:ipj,1 ), nrec ) 707 CASE DEFAULT ; CALL iom_get ( num, jpdom_unknown, clvar, dta_read(1:ilendta,1:ipj,1:ipk), nrec ) 756 CASE(1) ; 757 CALL iom_get ( num, jpdom_unknown, clvar, dta_read(1:ilendta,1:ipj,1 ), nrec ) 758 IF ( map%ll_unstruc) THEN ! unstructured open boundary data file 759 DO ib = 1, ipi 760 DO ik = 1, ipk 761 dta(ib,1,ik) = dta_read(map%ptr(ib),1,ik) 762 END DO 763 END DO 764 ELSE ! we assume that this is a structured open boundary file 765 DO ib = 1, ipi 766 jj=1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta)) 767 ji=map%ptr(ib)-(jj-1)*ilendta 768 DO ik = 1, ipk 769 dta(ib,1,ik) = dta_read(ji,jj,ik) 770 END DO 771 END DO 772 ENDIF 773 774 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 775 ! Do we include something here to adjust barotropic velocities ! 776 ! in case of a depth difference between bdy files and ! 777 ! bathymetry in the case ln_full_vel = .false. and jpk_bdy>0? ! 778 ! [as the enveloping and parital cells could change H] ! 779 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 780 781 CASE DEFAULT ; 782 783 IF( PRESENT(jpk_bdy) .AND. jpk_bdy>0 ) THEN ! boundary data not on model grid: vertical interpolation 784 CALL iom_getatt(num, '_FillValue', fv, cdvar=clvar ) 785 dta_read(:,:,:) = -ABS(fv) 786 dta_read_z(:,:,:) = 0._wp 787 dta_read_dz(:,:,:) = 0._wp 788 CALL iom_get ( num, jpdom_unknown, clvar, dta_read(1:ilendta,1:ipj,1:jpk_bdy), nrec ) 789 SELECT CASE( igrd ) 790 CASE(1) 791 CALL iom_get ( num, jpdom_unknown, 'gdept', dta_read_z(1:ilendta,1:ipj,1:jpk_bdy) ) 792 CALL iom_get ( num, jpdom_unknown, 'e3t', dta_read_dz(1:ilendta,1:ipj,1:jpk_bdy) ) 793 CASE(2) 794 CALL iom_get ( num, jpdom_unknown, 'gdepu', dta_read_z(1:ilendta,1:ipj,1:jpk_bdy) ) 795 CALL iom_get ( num, jpdom_unknown, 'e3u', dta_read_dz(1:ilendta,1:ipj,1:jpk_bdy) ) 796 CASE(3) 797 CALL iom_get ( num, jpdom_unknown, 'gdepv', dta_read_z(1:ilendta,1:ipj,1:jpk_bdy) ) 798 CALL iom_get ( num, jpdom_unknown, 'e3v', dta_read_dz(1:ilendta,1:ipj,1:jpk_bdy) ) 799 END SELECT 800 801 IF ( ln_bdy ) & 802 CALL fld_bdy_interp(dta_read, dta_read_z, dta_read_dz, map, jpk_bdy, igrd, ibdy, fv, dta, fvl, ilendta) 803 804 ELSE ! boundary data assumed to be on model grid 805 806 CALL iom_get ( num, jpdom_unknown, clvar, dta_read(1:ilendta,1:ipj,1:ipk), nrec ) 807 IF ( map%ll_unstruc) THEN ! unstructured open boundary data file 808 DO ib = 1, ipi 809 DO ik = 1, ipk 810 dta(ib,1,ik) = dta_read(map%ptr(ib),1,ik) 811 END DO 812 END DO 813 ELSE ! we assume that this is a structured open boundary file 814 DO ib = 1, ipi 815 jj=1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta)) 816 ji=map%ptr(ib)-(jj-1)*ilendta 817 DO ik = 1, ipk 818 dta(ib,1,ik) = dta_read(ji,jj,ik) 819 END DO 820 END DO 821 ENDIF 822 ENDIF ! PRESENT(jpk_bdy) 708 823 END SELECT 709 ! 710 IF( map%ll_unstruc ) THEN ! unstructured open boundary data file 824 825 END SUBROUTINE fld_map 826 827 SUBROUTINE fld_bdy_interp(dta_read, dta_read_z, dta_read_dz, map, jpk_bdy, igrd, ibdy, fv, dta, fvl, ilendta) 828 829 !!--------------------------------------------------------------------- 830 !! *** ROUTINE fld_bdy_interp *** 831 !! 832 !! ** Purpose : on the fly vertical interpolation to allow the use of 833 !! boundary data from non-native vertical grid 834 !!---------------------------------------------------------------------- 835 USE bdy_oce, ONLY: idx_bdy ! indexing for map <-> ij transformation 836 837 REAL(wp), POINTER, DIMENSION(:,:,:), INTENT(in ) :: dta_read ! work space for global data 838 REAL(wp), POINTER, DIMENSION(:,:,:), INTENT(in ) :: dta_read_z ! work space for global data 839 REAL(wp), POINTER, DIMENSION(:,:,:), INTENT(in ) :: dta_read_dz ! work space for global data 840 REAL(wp) , INTENT(in) :: fv ! fillvalue and alternative -ABS(fv) 841 REAL(wp), DIMENSION(:,:,:), INTENT(out) :: dta ! output field on model grid (2 dimensional) 842 TYPE(MAP_POINTER) , INTENT(in ) :: map ! global-to-local mapping indices 843 LOGICAL , INTENT(in), OPTIONAL :: fvl ! grid type, set number and number of vertical levels in the bdy data 844 INTEGER , INTENT(in) :: igrd, ibdy, jpk_bdy ! number of levels in bdy data 845 INTEGER , INTENT(in) :: ilendta ! length of data in file 846 !! 847 INTEGER :: ipi ! length of boundary data on local process 848 INTEGER :: ipj ! length of dummy dimension ( = 1 ) 849 INTEGER :: ipk ! number of vertical levels of dta ( 2D: ipk=1 ; 3D: ipk=jpk ) 850 INTEGER :: jpkm1_bdy ! number of levels in bdy data minus 1 851 INTEGER :: ib, ik, ikk ! loop counters 852 INTEGER :: ji, jj, zij, zjj ! temporary indices 853 REAL(wp) :: zl, zi, zh ! tmp variable for current depth and interpolation factor 854 REAL(wp) :: fv_alt, ztrans, ztrans_new ! fillvalue and alternative -ABS(fv) 855 CHARACTER (LEN=10) :: ibstr 856 !!--------------------------------------------------------------------- 857 858 859 ipi = SIZE( dta, 1 ) 860 ipj = SIZE( dta_read, 2 ) 861 ipk = SIZE( dta, 3 ) 862 jpkm1_bdy = jpk_bdy-1 863 864 fv_alt = -ABS(fv) ! set _FillValue < 0 as we make use of MAXVAL and MAXLOC later 865 DO ib = 1, ipi 866 zij = idx_bdy(ibdy)%nbi(ib,igrd) 867 zjj = idx_bdy(ibdy)%nbj(ib,igrd) 868 IF(narea==2) WRITE(*,*) 'MAPI', ib, igrd, map%ptr(ib), narea-1, zij, zjj 869 ENDDO 870 ! 871 IF ( map%ll_unstruc ) THEN ! unstructured open boundary data file 872 711 873 DO ib = 1, ipi 712 DO ik = 1, ipk 713 dta(ib,1,ik) = dta_read(map%ptr(ib),1,ik) 874 DO ik = 1, jpk_bdy 875 IF( ( dta_read(map%ptr(ib),1,ik) == fv ) ) THEN 876 dta_read_z(map%ptr(ib),1,ik) = fv_alt ! safety: put fillvalue into external depth field so consistent with data 877 dta_read_dz(map%ptr(ib),1,ik) = 0._wp ! safety: put 0._wp into external thickness factors to ensure transport is correct 878 ENDIF 879 ENDDO 880 ENDDO 881 882 DO ib = 1, ipi 883 zij = idx_bdy(ibdy)%nbi(ib,igrd) 884 zjj = idx_bdy(ibdy)%nbj(ib,igrd) 885 zh = SUM(dta_read_dz(map%ptr(ib),1,:) ) 886 ! Warnings to flag differences in the input and model topgraphy - is this useful/necessary? 887 SELECT CASE( igrd ) 888 CASE(1) 889 IF( ABS( (zh - ht_n(zij,zjj)) / ht_n(zij,zjj)) * tmask(zij,zjj,1) > 0.01_wp ) THEN 890 WRITE(ibstr,"(I10.10)") map%ptr(ib) 891 CALL ctl_warn('fld_bdy_interp: T depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%') 892 ! IF(lwp) WRITE(*,*) 'DEPTHT', zh, sum(e3t_n(zij,zjj,:), mask=tmask(zij,zjj,:)==1), ht_n(zij,zjj), map%ptr(ib), ib, zij, zjj 893 ENDIF 894 CASE(2) 895 IF( ABS( (zh - hu_n(zij,zjj)) * r1_hu_n(zij,zjj)) * umask(zij,zjj,1) > 0.01_wp ) THEN 896 WRITE(ibstr,"(I10.10)") map%ptr(ib) 897 CALL ctl_warn('fld_bdy_interp: U depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%') 898 IF(lwp) WRITE(*,*) 'DEPTHU', zh, sum(e3u_n(zij,zjj,:), mask=umask(zij,zjj,:)==1), sum(umask(zij,zjj,:)), & 899 & hu_n(zij,zjj), map%ptr(ib), ib, zij, zjj, narea-1 , & 900 & dta_read(map%ptr(ib),1,:) 901 ENDIF 902 CASE(3) 903 IF( ABS( (zh - hv_n(zij,zjj)) * r1_hv_n(zij,zjj)) * vmask(zij,zjj,1) > 0.01_wp ) THEN 904 WRITE(ibstr,"(I10.10)") map%ptr(ib) 905 CALL ctl_warn('fld_bdy_interp: V depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%') 906 ENDIF 907 END SELECT 908 DO ik = 1, ipk 909 SELECT CASE( igrd ) 910 CASE(1) 911 zl = gdept_n(zij,zjj,ik) ! if using in step could use fsdept instead of gdept_n? 912 CASE(2) 913 IF(ln_sco) THEN 914 zl = ( gdept_n(zij,zjj,ik) + gdept_n(zij+1,zjj,ik) ) * 0.5_wp ! if using in step could use fsdept instead of gdept_n? 915 ELSE 916 zl = MIN( gdept_n(zij,zjj,ik), gdept_n(zij+1,zjj,ik) ) 917 ENDIF 918 CASE(3) 919 IF(ln_sco) THEN 920 zl = ( gdept_n(zij,zjj,ik) + gdept_n(zij,zjj+1,ik) ) * 0.5_wp ! if using in step could use fsdept instead of gdept_n? 921 ELSE 922 zl = MIN( gdept_n(zij,zjj,ik), gdept_n(zij,zjj+1,ik) ) 923 ENDIF 924 END SELECT 925 IF( zl < dta_read_z(map%ptr(ib),1,1) ) THEN ! above the first level of external data 926 dta(ib,1,ik) = dta_read(map%ptr(ib),1,1) 927 ELSEIF( zl > MAXVAL(dta_read_z(map%ptr(ib),1,:),1) ) THEN ! below the last level of external data 928 dta(ib,1,ik) = dta_read(map%ptr(ib),1,MAXLOC(dta_read_z(map%ptr(ib),1,:),1)) 929 ELSE ! inbetween : vertical interpolation between ikk & ikk+1 930 DO ikk = 1, jpkm1_bdy ! when gdept_n(ikk) < zl < gdept_n(ikk+1) 931 IF( ( (zl-dta_read_z(map%ptr(ib),1,ikk)) * (zl-dta_read_z(map%ptr(ib),1,ikk+1)) <= 0._wp) & 932 & .AND. (dta_read_z(map%ptr(ib),1,ikk+1) /= fv_alt)) THEN 933 zi = ( zl - dta_read_z(map%ptr(ib),1,ikk) ) / & 934 & ( dta_read_z(map%ptr(ib),1,ikk+1) - dta_read_z(map%ptr(ib),1,ikk) ) 935 dta(ib,1,ik) = dta_read(map%ptr(ib),1,ikk) + & 936 & ( dta_read(map%ptr(ib),1,ikk+1) - dta_read(map%ptr(ib),1,ikk) ) * zi 937 ENDIF 938 END DO 939 ENDIF 714 940 END DO 715 941 END DO 716 ELSE ! structured open boundary data file 942 943 IF(igrd == 2) THEN ! do we need to adjust the transport term? 944 DO ib = 1, ipi 945 zij = idx_bdy(ibdy)%nbi(ib,igrd) 946 zjj = idx_bdy(ibdy)%nbj(ib,igrd) 947 zh = SUM(dta_read_dz(map%ptr(ib),1,:) ) 948 ztrans = 0._wp 949 ztrans_new = 0._wp 950 DO ik = 1, jpk_bdy ! calculate transport on input grid 951 ztrans = ztrans + dta_read(map%ptr(ib),1,ik) * dta_read_dz(map%ptr(ib),1,ik) 952 ENDDO 953 DO ik = 1, ipk ! calculate transport on model grid 954 ztrans_new = ztrans_new + dta(ib,1,ik) * e3u_n(zij,zjj,ik) * umask(zij,zjj,ik) 955 ENDDO 956 DO ik = 1, ipk ! make transport correction 957 IF(fvl) THEN ! bdy data are total velocity so adjust bt transport term to match input data 958 dta(ib,1,ik) = ( dta(ib,1,ik) + ( ztrans - ztrans_new ) * r1_hu_n(zij,zjj) ) * umask(zij,zjj,ik) 959 ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero 960 IF( ABS(ztrans * r1_hu_n(zij,zjj)) > 0.01_wp ) & 961 & CALL ctl_warn('fld_bdy_interp: barotropic component of > 0.01 ms-1 found in baroclinic velocities at') 962 dta(ib,1,ik) = dta(ib,1,ik) + ( 0._wp - ztrans_new ) * r1_hu_n(zij,zjj) * umask(zij,zjj,ik) 963 ENDIF 964 ENDDO 965 ENDDO 966 ENDIF 967 968 IF(igrd == 3) THEN ! do we need to adjust the transport term? 969 DO ib = 1, ipi 970 zij = idx_bdy(ibdy)%nbi(ib,igrd) 971 zjj = idx_bdy(ibdy)%nbj(ib,igrd) 972 zh = SUM(dta_read_dz(map%ptr(ib),1,:) ) 973 ztrans = 0._wp 974 ztrans_new = 0._wp 975 DO ik = 1, jpk_bdy ! calculate transport on input grid 976 ztrans = ztrans + dta_read(map%ptr(ib),1,ik) * dta_read_dz(map%ptr(ib),1,ik) 977 ENDDO 978 DO ik = 1, ipk ! calculate transport on model grid 979 ztrans_new = ztrans_new + dta(ib,1,ik) * e3v_n(zij,zjj,ik) * vmask(zij,zjj,ik) 980 ENDDO 981 DO ik = 1, ipk ! make transport correction 982 IF(fvl) THEN ! bdy data are total velocity so adjust bt transport term to match input data 983 dta(ib,1,ik) = ( dta(ib,1,ik) + ( ztrans - ztrans_new ) * r1_hv_n(zij,zjj) ) * vmask(zij,zjj,ik) 984 ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero 985 dta(ib,1,ik) = dta(ib,1,ik) + ( 0._wp - ztrans_new ) * r1_hv_n(zij,zjj) * vmask(zij,zjj,ik) 986 ENDIF 987 ENDDO 988 ENDDO 989 ENDIF 990 991 ELSE ! structured open boundary file 992 717 993 DO ib = 1, ipi 718 994 jj=1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta)) 719 995 ji=map%ptr(ib)-(jj-1)*ilendta 720 DO ik = 1, ipk 721 dta(ib,1,ik) = dta_read(ji,jj,ik) 996 DO ik = 1, jpk_bdy 997 IF( ( dta_read(ji,jj,ik) == fv ) ) THEN 998 dta_read_z(ji,jj,ik) = fv_alt ! safety: put fillvalue into external depth field so consistent with data 999 dta_read_dz(ji,jj,ik) = 0._wp ! safety: put 0._wp into external thickness factors to ensure transport is correct 1000 ENDIF 1001 ENDDO 1002 ENDDO 1003 1004 1005 DO ib = 1, ipi 1006 jj=1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta)) 1007 ji=map%ptr(ib)-(jj-1)*ilendta 1008 zij = idx_bdy(ibdy)%nbi(ib,igrd) 1009 zjj = idx_bdy(ibdy)%nbj(ib,igrd) 1010 zh = SUM(dta_read_dz(ji,jj,:) ) 1011 ! Warnings to flag differences in the input and model topgraphy - is this useful/necessary? 1012 SELECT CASE( igrd ) 1013 CASE(1) 1014 IF( ABS( (zh - ht_n(zij,zjj)) / ht_n(zij,zjj)) * tmask(zij,zjj,1) > 0.01_wp ) THEN 1015 WRITE(ibstr,"(I10.10)") map%ptr(ib) 1016 CALL ctl_warn('fld_bdy_interp: T depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%') 1017 ! IF(lwp) WRITE(*,*) 'DEPTHT', zh, sum(e3t_n(zij,zjj,:), mask=tmask(zij,zjj,:)==1), ht_n(zij,zjj), map%ptr(ib), ib, zij, zjj 1018 ENDIF 1019 CASE(2) 1020 IF( ABS( (zh - hu_n(zij,zjj)) * r1_hu_n(zij,zjj)) * umask(zij,zjj,1) > 0.01_wp ) THEN 1021 WRITE(ibstr,"(I10.10)") map%ptr(ib) 1022 CALL ctl_warn('fld_bdy_interp: U depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%') 1023 ENDIF 1024 CASE(3) 1025 IF( ABS( (zh - hv_n(zij,zjj)) * r1_hv_n(zij,zjj)) * vmask(zij,zjj,1) > 0.01_wp ) THEN 1026 WRITE(ibstr,"(I10.10)") map%ptr(ib) 1027 CALL ctl_warn('fld_bdy_interp: V depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%') 1028 ENDIF 1029 END SELECT 1030 DO ik = 1, ipk 1031 SELECT CASE( igrd ) ! coded for sco - need zco and zps option using min 1032 CASE(1) 1033 zl = gdept_n(zij,zjj,ik) ! if using in step could use fsdept instead of gdept_n? 1034 CASE(2) 1035 IF(ln_sco) THEN 1036 zl = ( gdept_n(zij,zjj,ik) + gdept_n(zij+1,zjj,ik) ) * 0.5_wp ! if using in step could use fsdept instead of gdept_n? 1037 ELSE 1038 zl = MIN( gdept_n(zij,zjj,ik), gdept_n(zij+1,zjj,ik) ) 1039 ENDIF 1040 CASE(3) 1041 IF(ln_sco) THEN 1042 zl = ( gdept_n(zij,zjj,ik) + gdept_n(zij,zjj+1,ik) ) * 0.5_wp ! if using in step could use fsdept instead of gdept_n? 1043 ELSE 1044 zl = MIN( gdept_n(zij,zjj,ik), gdept_n(zij,zjj+1,ik) ) 1045 ENDIF 1046 END SELECT 1047 IF( zl < dta_read_z(ji,jj,1) ) THEN ! above the first level of external data 1048 dta(ib,1,ik) = dta_read(ji,jj,1) 1049 ELSEIF( zl > MAXVAL(dta_read_z(ji,jj,:),1) ) THEN ! below the last level of external data 1050 dta(ib,1,ik) = dta_read(ji,jj,MAXLOC(dta_read_z(ji,jj,:),1)) 1051 ELSE ! inbetween : vertical interpolation between ikk & ikk+1 1052 DO ikk = 1, jpkm1_bdy ! when gdept_n(ikk) < zl < gdept_n(ikk+1) 1053 IF( ( (zl-dta_read_z(ji,jj,ikk)) * (zl-dta_read_z(ji,jj,ikk+1)) <= 0._wp) & 1054 & .AND. (dta_read_z(ji,jj,ikk+1) /= fv_alt)) THEN 1055 zi = ( zl - dta_read_z(ji,jj,ikk) ) / & 1056 & ( dta_read_z(ji,jj,ikk+1) - dta_read_z(ji,jj,ikk) ) 1057 dta(ib,1,ik) = dta_read(ji,jj,ikk) + & 1058 & ( dta_read(ji,jj,ikk+1) - dta_read(ji,jj,ikk) ) * zi 1059 ENDIF 1060 END DO 1061 ENDIF 722 1062 END DO 723 1063 END DO 724 ENDIF 725 ! 726 END SUBROUTINE fld_map 1064 1065 IF(igrd == 2) THEN ! do we need to adjust the transport term? 1066 DO ib = 1, ipi 1067 jj=1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta)) 1068 ji=map%ptr(ib)-(jj-1)*ilendta 1069 zij = idx_bdy(ibdy)%nbi(ib,igrd) 1070 zjj = idx_bdy(ibdy)%nbj(ib,igrd) 1071 zh = SUM(dta_read_dz(ji,jj,:) ) 1072 ztrans = 0._wp 1073 ztrans_new = 0._wp 1074 DO ik = 1, jpk_bdy ! calculate transport on input grid 1075 ztrans = ztrans + dta_read(ji,jj,ik) * dta_read_dz(ji,jj,ik) 1076 ENDDO 1077 DO ik = 1, ipk ! calculate transport on model grid 1078 ztrans_new = ztrans_new + dta(ib,1,ik) * e3u_n(zij,zjj,ik) * umask(zij,zjj,ik) 1079 ENDDO 1080 DO ik = 1, ipk ! make transport correction 1081 IF(fvl) THEN ! bdy data are total velocity so adjust bt transport term to match input data 1082 dta(ib,1,ik) = ( dta(ib,1,ik) + ( ztrans - ztrans_new ) * r1_hu_n(zij,zjj) ) * umask(zij,zjj,ik) 1083 ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero 1084 dta(ib,1,ik) = ( dta(ib,1,ik) + ( 0._wp - ztrans_new ) * r1_hu_n(zij,zjj) ) * umask(zij,zjj,ik) 1085 ENDIF 1086 ENDDO 1087 ENDDO 1088 ENDIF 1089 1090 IF(igrd == 3) THEN ! do we need to adjust the transport term? 1091 DO ib = 1, ipi 1092 jj = 1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta)) 1093 ji = map%ptr(ib)-(jj-1)*ilendta 1094 zij = idx_bdy(ibdy)%nbi(ib,igrd) 1095 zjj = idx_bdy(ibdy)%nbj(ib,igrd) 1096 zh = SUM(dta_read_dz(ji,jj,:) ) 1097 ztrans = 0._wp 1098 ztrans_new = 0._wp 1099 DO ik = 1, jpk_bdy ! calculate transport on input grid 1100 ztrans = ztrans + dta_read(ji,jj,ik) * dta_read_dz(ji,jj,ik) 1101 ENDDO 1102 DO ik = 1, ipk ! calculate transport on model grid 1103 ztrans_new = ztrans_new + dta(ib,1,ik) * e3v_n(zij,zjj,ik) * vmask(zij,zjj,ik) 1104 ENDDO 1105 DO ik = 1, ipk ! make transport correction 1106 IF(fvl) THEN ! bdy data are total velocity so adjust bt transport term to match input data 1107 dta(ib,1,ik) = ( dta(ib,1,ik) + ( ztrans - ztrans_new ) * r1_hv_n(zij,zjj) ) * vmask(zij,zjj,ik) 1108 ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero 1109 dta(ib,1,ik) = ( dta(ib,1,ik) + ( 0._wp - ztrans_new ) * r1_hv_n(zij,zjj) ) * vmask(zij,zjj,ik) 1110 ENDIF 1111 ENDDO 1112 ENDDO 1113 ENDIF 1114 1115 ENDIF ! endif unstructured or structured 1116 1117 END SUBROUTINE fld_bdy_interp 727 1118 728 1119 -
branches/2016/dev_NOC_CMCC_merge_2016/NEMOGCM/NEMO/TOP_SRC/trcbdy.F90
r7299 r7397 46 46 INTEGER, INTENT( in ) :: kt ! Main time step counter 47 47 !! 48 INTEGER :: ib_bdy , jn! Loop indeces48 INTEGER :: ib_bdy ,jn ,igrd ! Loop indeces 49 49 REAL(wp), POINTER, DIMENSION(:,:) :: ztrc 50 50 REAL(wp), POINTER :: zfac … … 52 52 ! 53 53 IF( nn_timing == 1 ) CALL timing_start('trc_bdy') 54 ! 55 igrd = 1 54 56 ! 55 57 DO ib_bdy=1, nb_bdy … … 63 65 CASE('frs' ) ; CALL bdy_frs( idx_bdy(ib_bdy), tra(:,:,:,jn), ztrc*zfac ) 64 66 CASE('specified' ) ; CALL bdy_spe( idx_bdy(ib_bdy), tra(:,:,:,jn), ztrc*zfac ) 65 CASE('neumann' ) ; CALL bdy_nmn( idx_bdy(ib_bdy), 67 CASE('neumann' ) ; CALL bdy_nmn( idx_bdy(ib_bdy), igrd , tra(:,:,:,jn) ) 66 68 CASE('orlanski' ) ; CALL bdy_orl( idx_bdy(ib_bdy), trb(:,:,:,jn), tra(:,:,:,jn), ztrc*zfac, ll_npo=.false. ) 67 69 CASE('orlanski_npo') ; CALL bdy_orl( idx_bdy(ib_bdy), trb(:,:,:,jn), tra(:,:,:,jn), ztrc*zfac, ll_npo=.true. )
Note: See TracChangeset
for help on using the changeset viewer.
