Changeset 6667
- Timestamp:
- 2016-06-06T07:57:00+02:00 (9 years ago)
- Location:
- branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM
- Files:
-
- 1 added
- 23 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/CONFIG/AMM12/EXP00/namelist_cfg
r6624 r6667 20 20 &namcfg ! parameters of the configuration 21 21 !----------------------------------------------------------------------- 22 cp_cfg = "amm" ! name of the configuration 23 jp_cfg = 011 ! resolution of the configuration 24 / 25 !----------------------------------------------------------------------- 26 &namzgr ! vertical coordinate 27 !----------------------------------------------------------------------- 28 ln_sco = .true. ! s- or hybrid z-s-coordinate 29 / 30 !----------------------------------------------------------------------- 31 &namzgr_sco ! s-coordinate or hybrid z-s-coordinate 32 !----------------------------------------------------------------------- 33 ln_s_sh94 = .false. ! Song & Haidvogel 1994 hybrid S-sigma (T)| 34 ln_s_sf12 = .true. ! Siddorn & Furner 2012 hybrid S-z-sigma (T)| if both are false the NEMO tanh stretching is applied 35 ln_sigcrit = .true. ! use sigma coordinates below critical depth (T) or Z coordinates (F) for Siddorn & Furner stretch 36 ! stretching coefficients for all functions 37 rn_hc = 50.0 ! critical depth for transition to stretched coordinates 38 / 22 ln_read_cfg = .false. ! (=T) read the domain configuration in 'domain_cfg.nc" file 23 ! ! (=F) user defined configuration ==>>> see usrdef(_...) modules 24 ln_write_cfg= .false. ! (=T) create the domain configuration file 25 ! 26 cp_cfg = "amm" ! name of the configuration 27 jp_cfg = 011 ! resolution of the configuration 28 / 29 39 30 !----------------------------------------------------------------------- 40 31 &namdom ! space and time domain (bathymetry, mesh, timestep) 41 32 !----------------------------------------------------------------------- 42 33 rn_rdt = 600. ! time step for the dynamics (and tracer if nn_acc=0) 43 ppsur = 999999.0 ! ORCA r4, r2 and r05 coefficients44 ppa0 = 999999.0 ! (default coefficients)45 ppa1 = 999999.0 !46 ppkth = 23.563 !47 ppacr = 9.0 !48 ppdzmin = 6.0 ! Minimum vertical spacing49 pphmax = 5720. ! Maximum depth50 ldbletanh = .FALSE. ! Use/do not use double tanf function for vertical coordinates51 ppa2 = 999999. ! Double tanh function parameters52 ppkth2 = 999999. !53 ppacr2 = 999999.54 34 / 55 35 !----------------------------------------------------------------------- 56 36 &namcrs ! Grid coarsening for dynamics output and/or 57 37 ! ! passive tracer coarsened online simulations 58 38 !----------------------------------------------------------------------- 59 39 / -
branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/CONFIG/GYRE/EXP00/namelist_cfg
r6624 r6667 15 15 cn_exp = "GYRE" ! experience name 16 16 nn_it000 = 1 ! first time step 17 nn_itend = 4320 ! last time step17 nn_itend = 4320 !!gm 4320 ! last time step 18 18 nn_leapy = 30 ! Leap year calendar (1) or not (0) 19 19 nn_stock = 4320 ! frequency of creation of a restart file (modulo referenced to 1) … … 21 21 22 22 ln_clobber = .true. ! clobber (overwrite) an existing file 23 24 23 / 25 24 !----------------------------------------------------------------------- 26 25 &namcfg ! parameters of the configuration 27 26 !----------------------------------------------------------------------- 28 ln_read_cfg = .true. !!gm .false. ! flag to read (.true.) configuration definition files (coordinates, 29 cp_cfg = "gyre" ! name of the configuration 30 jp_cfg = 1 ! resolution of the configuration 31 / 32 !----------------------------------------------------------------------- 33 &namzgr ! vertical coordinate 34 !----------------------------------------------------------------------- 35 ln_zco = .true. ! z-coordinate - full steps 36 ln_linssh = .true. ! linear free surface 37 / 38 !----------------------------------------------------------------------- 39 &namzgr_sco ! s-coordinate or hybrid z-s-coordinate 40 !----------------------------------------------------------------------- 27 ln_read_cfg = .false. ! (=T) read the domain configuration in 'domain_cfg.nc" file 28 ! ! (=F) user defined configuration ==>>> see usrdef(_...) modules 29 ln_write_cfg= .true. ! (=T) create the domain configuration file 30 ! 31 cp_cfg = "default" ! name of the configuration 32 jp_cfg = 0 ! resolution of the configuration 33 ln_use_jattr = .false. ! use (T) the file attribute: open_ocean_jstart, if present 34 ! ! in netcdf input files, as the start j-row for reading 41 35 / 42 36 !----------------------------------------------------------------------- 43 37 &namdom ! space and time domain (bathymetry, mesh, timestep) 44 38 !----------------------------------------------------------------------- 45 nn_bathy = 0 ! compute (=0) or read (=1) the bathymetry file 46 rn_rdt = 7200. ! time step for the dynamics 47 nn_msh = -1 ! create (=1) a mesh file or not (=0) 48 ppsur = -2033.194295283385 ! ORCA r4, r2 and r05 coefficients 49 ppa0 = 155.8325369664153 ! (default coefficients) 50 ppa1 = 146.3615918601890 ! 51 ppkth = 17.28520372419791 ! 52 ppacr = 5.0 ! 53 ppdzmin = 999999.0 ! Minimum vertical spacing 54 pphmax = 999999.0 ! Maximum depth 55 ldbletanh = .FALSE. ! Use/do not use double tanf function for vertical coordinates 56 ppa2 = 999999.0 ! Double tanh function parameters 57 ppkth2 = 999999.0 ! 58 ppacr2 = 999999.0 ! 39 ln_linssh = .true. ! =T linear free surface ==>> model level are fixed in time 40 ! 41 nn_msh = 0 ! create (>0) a mesh file or not (=0) 42 ! 43 rn_rdt = 7200. ! time step for the dynamics (and tracer if nn_acc=0) 59 44 / 60 45 !----------------------------------------------------------------------- 61 46 &namcrs ! Grid coarsening for dynamics output and/or 62 47 ! ! passive tracer coarsened online simulations 63 48 !----------------------------------------------------------------------- 64 49 / -
branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/CONFIG/GYRE_BFM/EXP00/namelist_cfg
r6624 r6667 22 22 &namcfg ! parameters of the configuration 23 23 !----------------------------------------------------------------------- 24 cp_cfg = "gyre" ! name of the configuration 25 jp_cfg = 1 ! resolution of the configuration 26 / 27 &namzgr ! vertical coordinate 28 !----------------------------------------------------------------------- 29 ln_zco = .true. ! z-coordinate - full steps 30 ln_linssh = .true. ! linear free surface 31 / 32 !----------------------------------------------------------------------- 33 &namzgr_sco ! s-coordinate or hybrid z-s-coordinate 34 !----------------------------------------------------------------------- 24 ln_read_cfg = .false. ! (=T) read the domain configuration in 'domain_cfg.nc" file 25 ! ! (=F) user defined configuration ==>>> see usrdef(_...) modules 26 ln_write_cfg= .true. ! (=T) create the domain configuration file 27 ! 28 cp_cfg = "default" ! name of the configuration 29 jp_cfg = 0 ! resolution of the configuration 30 ln_use_jattr = .false. ! use (T) the file attribute: open_ocean_jstart, if present 31 ! ! in netcdf input files, as the start j-row for reading 35 32 / 36 33 !----------------------------------------------------------------------- 37 34 &namdom ! space and time domain (bathymetry, mesh, timestep) 38 35 !----------------------------------------------------------------------- 39 nn_bathy = 0 ! compute (=0) or read (=1) the bathymetry file 40 rn_rdt = 7200. ! time step for the dynamics 41 ppsur = -2033.194295283385 ! ORCA r4, r2 and r05 coefficients 42 ppa0 = 155.8325369664153 ! (default coefficients) 43 ppa1 = 146.3615918601890 ! 44 ppkth = 17.28520372419791 ! 45 ppacr = 5.0 ! 46 ppdzmin = 999999.0 ! Minimum vertical spacing 47 pphmax = 999999.0 ! Maximum depth 48 ldbletanh = .FALSE. ! Use/do not use double tanf function for vertical coordinates 49 ppa2 = 999999.0 ! Double tanh function parameters 50 ppkth2 = 999999.0 ! 51 ppacr2 = 999999.0 ! 36 ln_linssh = .true. ! =T linear free surface ==>> model level are fixed in time 37 ! 38 nn_msh = 0 ! create (>0) a mesh file or not (=0) 39 ! 40 rn_rdt = 7200. ! time step for the dynamics (and tracer if nn_acc=0) 52 41 / 53 42 !----------------------------------------------------------------------- -
branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/CONFIG/GYRE_PISCES/EXP00/namelist_cfg
r6624 r6667 15 15 &namcfg ! parameters of the configuration 16 16 !----------------------------------------------------------------------- 17 cp_cfg = "gyre" ! name of the configuration18 jp_cfg = 1 ! resolution of the configuration19 / 20 !----------------------------------------------------------------------- 21 &namzgr ! vertical coordinate 22 !----------------------------------------------------------------------- 23 ln_ zco = .true. ! z-coordinate - full steps24 ln_linssh = .true. ! linear free surface17 ln_read_cfg = .false. ! (=T) read the domain configuration in 'domain_cfg.nc" file 18 ! ! (=F) user defined configuration ==>>> see usrdef(_...) modules 19 ln_write_cfg= .true. ! (=T) create the domain configuration file 20 ! 21 cp_cfg = "default" ! name of the configuration 22 jp_cfg = 0 ! resolution of the configuration 23 ln_use_jattr = .false. ! use (T) the file attribute: open_ocean_jstart, if present 24 ! ! in netcdf input files, as the start j-row for reading 25 25 / 26 26 !----------------------------------------------------------------------- 27 27 &namdom ! space and time domain (bathymetry, mesh, timestep) 28 28 !----------------------------------------------------------------------- 29 nn_bathy = 0 ! compute (=0) or read (=1) the bathymetry file 30 rn_rdt = 7200. ! time step for the dynamics 31 ppsur = -2033.194295283385 ! ORCA r4, r2 and r05 coefficients 32 ppa0 = 155.8325369664153 ! (default coefficients) 33 ppa1 = 146.3615918601890 ! 34 ppkth = 17.28520372419791 ! 35 ppacr = 5.0 ! 36 ppdzmin = 999999.0 ! Minimum vertical spacing 37 pphmax = 999999.0 ! Maximum depth 38 ldbletanh = .FALSE. ! Use/do not use double tanf function for vertical coordinates 39 ppa2 = 999999.0 ! Double tanh function parameters 40 ppkth2 = 999999.0 ! 41 ppacr2 = 999999.0 ! 29 ln_linssh = .true. ! =T linear free surface ==>> model level are fixed in time 30 ! 31 nn_msh = 0 ! create (>0) a mesh file or not (=0) 32 ! 33 rn_rdt = 7200. ! time step for the dynamics (and tracer if nn_acc=0) 42 34 / 43 35 !----------------------------------------------------------------------- -
branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/CONFIG/GYRE_XIOS/EXP00/namelist_cfg
r6624 r6667 15 15 &namcfg ! parameters of the configuration 16 16 !----------------------------------------------------------------------- 17 cp_cfg = "gyre" ! name of the configuration 18 jp_cfg = 1 ! resolution of the configuration 19 / 20 !----------------------------------------------------------------------- 21 &namzgr ! vertical coordinate 22 !----------------------------------------------------------------------- 23 ln_zco = .true. ! z-coordinate - full steps 24 ln_linssh = .true. ! linear free surface 25 / 26 !----------------------------------------------------------------------- 27 &namzgr_sco ! s-coordinate or hybrid z-s-coordinate 28 !----------------------------------------------------------------------- 17 ln_read_cfg = .false. ! (=T) read the domain configuration in 'domain_cfg.nc" file 18 ! ! (=F) user defined configuration ==>>> see usrdef(_...) modules 19 ln_write_cfg= .true. ! (=T) create the domain configuration file 20 ! 21 cp_cfg = "default" ! name of the configuration 22 jp_cfg = 0 ! resolution of the configuration 23 ln_use_jattr = .false. ! use (T) the file attribute: open_ocean_jstart, if present 24 ! ! in netcdf input files, as the start j-row for reading 29 25 / 30 26 !----------------------------------------------------------------------- 31 27 &namdom ! space and time domain (bathymetry, mesh, timestep) 32 28 !----------------------------------------------------------------------- 33 nn_bathy = 0 ! compute (=0) or read (=1) the bathymetry file 34 rn_rdt = 7200. ! time step for the dynamics 35 ppsur = -2033.194295283385 ! ORCA r4, r2 and r05 coefficients 36 ppa0 = 155.8325369664153 ! (default coefficients) 37 ppa1 = 146.3615918601890 ! 38 ppkth = 17.28520372419791 ! 39 ppacr = 5.0 ! 40 ppdzmin = 999999.0 ! Minimum vertical spacing 41 pphmax = 999999.0 ! Maximum depth 42 ldbletanh = .FALSE. ! Use/do not use double tanf function for vertical coordinates 43 ppa2 = 999999.0 ! Double tanh function parameters 44 ppkth2 = 999999.0 ! 45 ppacr2 = 999999.0 ! 29 ln_linssh = .true. ! =T linear free surface ==>> model level are fixed in time 30 ! 31 nn_msh = 0 ! create (>0) a mesh file or not (=0) 32 ! 33 rn_rdt = 7200. ! time step for the dynamics (and tracer if nn_acc=0) 46 34 / 47 35 !----------------------------------------------------------------------- -
branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/CONFIG/SHARED/namelist_ref
r6624 r6667 3 3 !!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 4 4 !! NEMO/OPA : 1 - run manager (namrun) 5 !! namelists 2 - Domain (namcfg, namzgr, nam zgr_sco, namdom, namtsd)5 !! namelists 2 - Domain (namcfg, namzgr, namdom, namtsd) 6 6 !! 3 - Surface boundary (namsbc, namsbc_flx, namsbc_clio, namsbc_core, namsbc_sas 7 7 !! namsbc_cpl, namtra_qsr, namsbc_rnf, … … 60 60 !! namcfg parameters of the configuration 61 61 !! namzgr vertical coordinate 62 !! namzgr_sco s-coordinate or hybrid z-s-coordinate63 62 !! namdom space and time domain (bathymetry, mesh, timestep) 64 63 !! namcrs coarsened grid (for outputs and/or TOP) ("key_crs") … … 72 71 &namcfg ! parameters of the configuration 73 72 !----------------------------------------------------------------------- 74 ln_read_cfg = .false. ! flag to read (.true.) configuration definition files (coordinates, 75 ! bathymetry, boudary condition, initial state, sbc) or (.false.) to call user_defined.F90 module 73 ln_read_cfg = .false. ! (=T) read the domain configuration in 'domain_cfg.nc" file 74 ! ! (=F) user defined configuration ==>>> see usrdef(_...) modules 75 ln_write_cfg= .false. ! (=T) create the domain configuration file 76 ! 76 77 cp_cfg = "default" ! name of the configuration 77 78 jp_cfg = 0 ! resolution of the configuration 78 79 ln_use_jattr = .false. ! use (T) the file attribute: open_ocean_jstart, if present 79 80 ! ! in netcdf input files, as the start j-row for reading 80 81 / 81 82 !----------------------------------------------------------------------- … … 86 87 ln_sco = .false. ! s- or hybrid z-s-coordinate 87 88 ln_isfcav = .false. ! ice shelf cavity 88 ln_linssh = .false. ! linear free surface89 /90 !-----------------------------------------------------------------------91 &namzgr_sco ! s-coordinate or hybrid z-s-coordinate92 !-----------------------------------------------------------------------93 ln_s_sh94 = .false. ! Song & Haidvogel 1994 hybrid S-sigma (T)|94 ln_s_sf12 = .false. ! Siddorn & Furner 2012 hybrid S-z-sigma (T)| if both are false the NEMO tanh stretching is applied95 ln_sigcrit = .false. ! use sigma coordinates below critical depth (T) or Z coordinates (F) for Siddorn & Furner stretch96 ! stretching coefficients for all functions97 rn_sbot_min = 10.0 ! minimum depth of s-bottom surface (>0) (m)98 rn_sbot_max = 7000.0 ! maximum depth of s-bottom surface (= ocean depth) (>0) (m)99 rn_hc = 150.0 ! critical depth for transition to stretched coordinates100 !!!!!!! Envelop bathymetry101 rn_rmax = 0.3 ! maximum cut-off r-value allowed (0<r_max<1)102 !!!!!!! SH94 stretching coefficients (ln_s_sh94 = .true.)103 rn_theta = 6.0 ! surface control parameter (0<=theta<=20)104 rn_bb = 0.8 ! stretching with SH94 s-sigma105 !!!!!!! SF12 stretching coefficient (ln_s_sf12 = .true.)106 rn_alpha = 4.4 ! stretching with SF12 s-sigma107 rn_efold = 0.0 ! efold length scale for transition to stretched coord108 rn_zs = 1.0 ! depth of surface grid box109 ! bottom cell depth (Zb) is a linear function of water depth Zb = H*a + b110 rn_zb_a = 0.024 ! bathymetry scaling factor for calculating Zb111 rn_zb_b = -0.2 ! offset for calculating Zb112 !!!!!!!! Other stretching (not SH94 or SF12) [also uses rn_theta above]113 rn_thetb = 1.0 ! bottom control parameter (0<=thetb<= 1)114 89 / 115 90 !----------------------------------------------------------------------- 116 91 &namdom ! space and time domain (bathymetry, mesh, timestep) 117 92 !----------------------------------------------------------------------- 118 nn_bathy = 1 ! compute (=0) or read (=1) the bathymetry file 119 rn_bathy = 0. ! value of the bathymetry. if (=0) bottom flat at jpkm1 93 ln_linssh = .false. ! =T linear free surface ==>> model level are fixed in time 120 94 nn_closea = 0 ! remove (=0) or keep (=1) closed seas and lakes (ORCA) 121 nn_msh = 0 ! create (=1) a mesh file or not (=0)122 rn_hmin = -3. ! min depth of the ocean (>0) or min number of ocean level (<0)95 ! 96 nn_msh = 0 ! create (>0) a mesh file or not (=0) 123 97 rn_isfhmin = 1.00 ! treshold (m) to discriminate grounding ice to floating ice 124 rn_e3zps_min= 20. ! partial step thickness is set larger than the minimum of 125 rn_e3zps_rat= 0.1 ! rn_e3zps_min and rn_e3zps_rat*e3t, with 0<rn_e3zps_rat<1 126 ! 98 ! 127 99 rn_rdt = 5760. ! time step for the dynamics (and tracer if nn_acc=0) 128 100 rn_atfp = 0.1 ! asselin time filter parameter 101 ! 129 102 ln_crs = .false. ! Logical switch for coarsening module 130 !131 ppsur = -4762.96143546300 ! ORCA r4, r2 and r05 coefficients132 ppa0 = 255.58049070440 ! (default coefficients)133 ppa1 = 245.58132232490 !134 ppkth = 21.43336197938 !135 ppacr = 3.0 !136 ppdzmin = 10. ! Minimum vertical spacing137 pphmax = 5000. ! Maximum depth138 ldbletanh = .TRUE. ! Use/do not use double tanf function for vertical coordinates139 ppa2 = 100.760928500000 ! Double tanh function parameters140 ppkth2 = 48.029893720000 !141 ppacr2 = 13.000000000000 !142 103 / 143 104 !----------------------------------------------------------------------- -
branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/DIA/diatmb.F90
r6140 r6667 4 4 !! Harmonic analysis of tidal constituents 5 5 !!====================================================================== 6 !! History : 3.6 ! 2014 (E O'Dea) Original code 6 !! History : 3.6 ! 08-2014 (E O'Dea) Original code 7 !! 3.7 ! 05-2016 (G. Madec) use mbkt, mikt to account for ocean cavities 7 8 !!---------------------------------------------------------------------- 8 9 USE oce ! ocean dynamics and tracers variables 9 10 USE dom_oce ! ocean space and time domain 11 ! 10 12 USE in_out_manager ! I/O units 11 13 USE iom ! I/0 library … … 31 33 !! *** ROUTINE dia_tmb_init *** 32 34 !! 33 !! ** Purpose :Initialization of tmb namelist35 !! ** Purpose : Initialization of tmb namelist 34 36 !! 35 !! ** Method : Read namelist 36 !! History 37 !! 3.6 ! 08-14 (E. O'Dea) Routine to initialize dia_tmb 37 !! ** Method : Read namelist 38 38 !!--------------------------------------------------------------------------- 39 !!40 39 INTEGER :: ios ! Local integer output status for namelist read 41 40 ! … … 59 58 WRITE(numout,*) 'Switch for TMB diagnostics (T) or not (F) ln_diatmb = ', ln_diatmb 60 59 ENDIF 61 60 ! 62 61 END SUBROUTINE dia_tmb_init 63 62 64 SUBROUTINE dia_calctmb( pinfield,pouttmb ) 63 64 SUBROUTINE dia_calctmb( pfield, ptmb ) 65 65 !!--------------------------------------------------------------------- 66 66 !! *** ROUTINE dia_tmb *** … … 68 68 !! ** Purpose : Find the Top, Mid and Bottom fields of water Column 69 69 !! 70 !! ** Method : 71 !! use mbathy to find surface, mid and bottom of model levels70 !! ** Method : use mbkt, mikt to find surface, mid and bottom of 71 !! model levels due to potential existence of ocean cavities 72 72 !! 73 !! History :74 !! 3.6 ! 08-14 (E. O'Dea) Routine based on dia_wri_foam75 73 !!---------------------------------------------------------------------- 76 !! * Modules used 77 78 ! Routine to map 3d field to top, middle, bottom 79 IMPLICIT NONE 80 81 82 ! Routine arguments 83 REAL(wp), DIMENSION(jpi, jpj, jpk), INTENT(IN ) :: pinfield ! Input 3d field and mask 84 REAL(wp), DIMENSION(jpi, jpj, 3 ), INTENT( OUT) :: pouttmb ! Output top, middle, bottom 85 86 87 88 ! Local variables 89 INTEGER :: ji,jj,jk ! Dummy loop indices 90 91 ! Local Real 92 REAL(wp) :: zmdi ! set masked values 93 94 zmdi=1.e+20 !missing data indicator for masking 95 96 ! Calculate top 97 pouttmb(:,:,1) = pinfield(:,:,1)*tmask(:,:,1) + zmdi*(1.0-tmask(:,:,1)) 98 99 ! Calculate middle 100 DO jj = 1,jpj 101 DO ji = 1,jpi 102 jk = max(1,mbathy(ji,jj)/2) 103 pouttmb(ji,jj,2) = pinfield(ji,jj,jk)*tmask(ji,jj,jk) + zmdi*(1.0-tmask(ji,jj,jk)) 74 REAL(wp), DIMENSION(jpi, jpj, jpk), INTENT(in ) :: pfield ! Input 3d field and mask 75 REAL(wp), DIMENSION(jpi, jpj, 3 ), INTENT( out) :: ptmb ! top, middle, bottom extracted from pfield 76 ! 77 INTEGER :: ji, jj ! Dummy loop indices 78 INTEGER :: itop, imid, ibot ! local integers 79 REAL(wp) :: zmdi = 1.e+20_wp ! land value 80 !!--------------------------------------------------------------------- 81 ! 82 DO jj = 1, jpj 83 DO ji = 1, jpi 84 itop = mikt(ji,jj) ! top ocean 85 ibot = mbkt(ji,jj) ! bottom ocean 86 imid = itop + ( ibot - itop + 1 ) / 2 ! middle ocean 87 ! 88 ptmb(ji,jj,1) = pfield(ji,jj,itop)*tmask(ji,jj,itop) + zmdi*( 1._wp-tmask(ji,jj,itop) ) 89 ptmb(ji,jj,2) = pfield(ji,jj,imid)*tmask(ji,jj,imid) + zmdi*( 1._wp-tmask(ji,jj,imid) ) 90 ptmb(ji,jj,3) = pfield(ji,jj,ibot)*tmask(ji,jj,ibot) + zmdi*( 1._wp-tmask(ji,jj,ibot) ) 104 91 END DO 105 92 END DO 106 107 ! Calculate bottom 108 DO jj = 1,jpj 109 DO ji = 1,jpi 110 jk = max(1,mbathy(ji,jj) - 1) 111 pouttmb(ji,jj,3) = pinfield(ji,jj,jk)*tmask(ji,jj,jk) + zmdi*(1.0-tmask(ji,jj,jk)) 112 END DO 113 END DO 114 93 ! 115 94 END SUBROUTINE dia_calctmb 116 117 95 118 96 … … 122 100 !! ** Purpose : Write diagnostics for Top, Mid and Bottom of water Column 123 101 !! 124 !! ** Method : 125 !! use mbathy to find surface, mid and bottom of model levels 102 !! ** Method : use mikt,mbkt to find surface, mid and bottom of model levels 126 103 !! calls calctmb to retrieve TMB values before sending to iom_put 127 104 !! 128 !! History :129 !! 3.6 ! 08-14 (E. O'Dea)130 !!131 105 !!-------------------------------------------------------------------- 132 REAL(wp), POINTER, DIMENSION(:,:,:) :: zwtmb ! temporary workspace 133 REAL(wp) :: zmdi ! set masked values 134 135 zmdi=1.e+20 !missing data indicator for maskin 136 106 REAL(wp) :: zmdi =1.e+20 ! land value 107 REAL(wp), POINTER, DIMENSION(:,:,:) :: zwtmb ! workspace 108 !!-------------------------------------------------------------------- 109 ! 137 110 IF (ln_diatmb) THEN 138 CALL wrk_alloc( jpi , jpj, 3, zwtmb )111 CALL wrk_alloc( jpi,jpj,3 , zwtmb ) 139 112 CALL dia_calctmb( tsn(:,:,:,jp_tem),zwtmb ) 140 113 !ssh already output but here we output it masked 141 CALL iom_put( "sshnmasked" , sshn(:,:)*tmask(:,:,1) + zmdi*(1.0 - tmask(:,:,1)) ) ! tmb Temperature114 CALL iom_put( "sshnmasked" , sshn(:,:)*tmask(:,:,1) + zmdi*(1.0 - tmask(:,:,1)) ) 142 115 CALL iom_put( "top_temp" , zwtmb(:,:,1) ) ! tmb Temperature 143 116 CALL iom_put( "mid_temp" , zwtmb(:,:,2) ) ! tmb Temperature … … 161 134 CALL iom_put( "bot_v" , zwtmb(:,:,3) ) ! tmb V Velocity 162 135 !Called in dynspg_ts.F90 CALL iom_put( "baro_v" , vn_b ) ! Barotropic V Velocity 136 CALL wrk_dealloc( jpi,jpj,3 , zwtmb ) 163 137 ELSE 164 138 CALL ctl_warn('dia_tmb: tmb diagnostic is set to false you should not have seen this') 165 139 ENDIF 166 140 ! 167 141 END SUBROUTINE dia_tmb 168 142 !!====================================================================== -
branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90
r6624 r6667 30 30 !! ---------------------------- 31 31 ! !!* Namelist namdom : time & space domain * 32 INTEGER , PUBLIC :: nn_bathy !: = 0/1 ,compute/read the bathymetry file33 REAL(wp), PUBLIC :: rn_bathy !: depth of flat bottom (active if nn_bathy=0; if =0 depth=jpkm1)34 REAL(wp), PUBLIC :: rn_hmin !: minimum ocean depth (>0) or minimum number of ocean levels (<0)32 LOGICAL , PUBLIC :: ln_linssh !: =T linear free surface ==>> model level are fixed in time 33 INTEGER , PUBLIC :: nn_closea !: =0 suppress closed sea/lake from the ORCA domain or not (=1) 34 INTEGER , PUBLIC :: nn_msh !: >0 create a mesh-mask file (mesh_mask.nc) 35 35 REAL(wp), PUBLIC :: rn_isfhmin !: threshold to discriminate grounded ice to floating ice 36 REAL(wp), PUBLIC :: rn_e3zps_min !: miminum thickness for partial steps (meters) 37 REAL(wp), PUBLIC :: rn_e3zps_rat !: minimum thickness ration for partial steps 38 INTEGER , PUBLIC :: nn_msh !: = 1 create a mesh-mask file 36 REAL(wp), PUBLIC :: rn_rdt !: time step for the dynamics and tracer 39 37 REAL(wp), PUBLIC :: rn_atfp !: asselin time filter parameter 40 REAL(wp), PUBLIC :: rn_rdt !: time step for the dynamics and tracer41 INTEGER , PUBLIC :: nn_closea !: =0 suppress closed sea/lake from the ORCA domain or not (=1)42 38 INTEGER , PUBLIC :: nn_euler !: =0 start with forward time step or not (=1) 43 39 LOGICAL , PUBLIC :: ln_iscpl !: coupling with ice sheet … … 58 54 REAL(wp), PUBLIC :: rn_bt_cmax !: Maximum allowed courant number (used if ln_bt_auto=T) 59 55 60 !! Vertical grid parameter for domzgr61 !! ==================================62 REAL(wp) :: ppsur !: ORCA r4, r2 and r05 coefficients63 REAL(wp) :: ppa0 !: (default coefficients)64 REAL(wp) :: ppa1 !:65 REAL(wp) :: ppkth !:66 REAL(wp) :: ppacr !:67 !68 ! If both ppa0 ppa1 and ppsur are specified to 0, then69 ! they are computed from ppdzmin, pphmax , ppkth, ppacr in dom_zgr70 REAL(wp) :: ppdzmin !: Minimum vertical spacing71 REAL(wp) :: pphmax !: Maximum depth72 !73 LOGICAL :: ldbletanh !: Use/do not use double tanf function for vertical coordinates74 REAL(wp) :: ppa2 !: Double tanh function parameters75 REAL(wp) :: ppkth2 !:76 REAL(wp) :: ppacr2 !:77 56 78 57 ! !! old non-DOCTOR names still used in the model 79 INTEGER , PUBLIC :: ntopo !: = 0/1 ,compute/read the bathymetry file80 REAL(wp), PUBLIC :: e3zps_min !: miminum thickness for partial steps (meters)81 REAL(wp), PUBLIC :: e3zps_rat !: minimum thickness ration for partial steps82 INTEGER , PUBLIC :: nmsh !: = 1 create a mesh-mask file83 58 REAL(wp), PUBLIC :: atfp !: asselin time filter parameter 84 59 REAL(wp), PUBLIC :: rdt !: time step for the dynamics and tracer … … 123 98 INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: mig !: local ==> global domain i-index 124 99 INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: mjg !: local ==> global domain j-index 125 INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: mi0, mi1 !: global ==> local domain i-index !!bug ==> other solution?126 ! ! (mi0=1 and mi1=0 if the global indexis not in the local domain)127 INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: mj0, mj1 !: global ==> local domain j-index !!bug ==> other solution?128 ! ! (mi0=1 and mi1=0 if the global indexis not in the local domain)100 INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: mi0, mi1 !: global ==> local domain i-index (mi0=1 and mi1=0 if the global index 101 ! ! is not in the local domain) 102 INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: mj0, mj1 !: global ==> local domain j-index (mj0=1 and mj1=0 if the global index 103 ! ! is not in the local domain) 129 104 INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: nimppt, njmppt !: i-, j-indexes for each processor 130 105 INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: ibonit, ibonjt !: i-, j- processor neighbour existence … … 154 129 !! vertical coordinate and scale factors 155 130 !! --------------------------------------------------------------------- 156 ! !!* Namelist namzgr : vertical coordinate *157 131 LOGICAL, PUBLIC :: ln_zco !: z-coordinate - full step 158 132 LOGICAL, PUBLIC :: ln_zps !: z-coordinate - partial step 159 133 LOGICAL, PUBLIC :: ln_sco !: s-coordinate or hybrid z-s coordinate 160 134 LOGICAL, PUBLIC :: ln_isfcav !: presence of ISF 161 LOGICAL, PUBLIC :: ln_linssh !: variable grid flag162 163 135 ! ! ref. ! before ! now ! after ! 164 136 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3t_0 , e3t_b , e3t_n , e3t_a !: t- vert. scale factor [m] … … 192 164 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: e3tp , e3wp !: ocean bottom level thickness at T and W points 193 165 194 !!gm This should be removed from here.... ==>>> only used in domzgr at initialization phase195 !! s-coordinate and hybrid z-s-coordinate196 !! =----------------======---------------197 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: gsigt, gsigw !: model level depth coefficient at t-, w-levels (analytic)198 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: gsi3w !: model level depth coefficient at w-level (sum of gsigw)199 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: esigt, esigw !: vertical scale factor coef. at t-, w-levels200 201 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hbatv , hbatf !: ocean depth at the vertical of v--f202 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hbatt , hbatu !: t--u points (m)203 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: scosrf, scobot !: ocean surface and bottom topographies204 ! ! (if deviating from coordinate surfaces in HYBRID)205 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hifv , hiff !: interface depth between stretching at v--f206 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hift , hifu !: and quasi-uniform spacing t--u points (m)207 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: rx1 !: Maximum grid stiffness ratio208 !!gm end209 166 210 167 !!---------------------------------------------------------------------- 211 168 !! masks, bathymetry 212 169 !! --------------------------------------------------------------------- 213 INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: mbathy !: number of ocean level (=0, 1, ... , jpk-1)214 INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: mbkt, mbku, mbkv !: bottom last wet T-, U- and W-level215 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: bathy !: ocean depth (meters)170 !!gm INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: mbathy !: number of ocean level (=0, 1, ... , jpk-1) 171 INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: mbkt, mbku, mbkv !: bottom last wet T-, U- and V-level 172 !!gm REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: bathy !: ocean depth (meters) 216 173 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tmask_i !: interior domain T-point mask 217 174 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tmask_h !: internal domain T-point mask (Figure 8.5 NEMO book) … … 336 293 ! 337 294 ALLOCATE( gdept_1d(jpk) , gdepw_1d(jpk) , & 338 & e3t_1d (jpk) , e3w_1d (jpk) , e3tp (jpi,jpj), e3wp(jpi,jpj) , & 339 & gsigt (jpk) , gsigw (jpk) , gsi3w(jpk) , & 340 & esigt (jpk) , esigw (jpk) , STAT=ierr(7) ) 341 ! 342 ALLOCATE( hbatv (jpi,jpj) , hbatf (jpi,jpj) , & 343 & hbatt (jpi,jpj) , hbatu (jpi,jpj) , & 344 & scosrf(jpi,jpj) , scobot(jpi,jpj) , & 345 & hifv (jpi,jpj) , hiff (jpi,jpj) , & 346 & hift (jpi,jpj) , hifu (jpi,jpj) , rx1(jpi,jpj) , STAT=ierr(8) ) 347 348 ALLOCATE( mbathy (jpi,jpj) , bathy (jpi,jpj) , & 349 & tmask_i(jpi,jpj) , tmask_h(jpi,jpj) , & 295 & e3t_1d (jpk) , e3w_1d (jpk) , e3tp (jpi,jpj), e3wp(jpi,jpj) , STAT=ierr(7) ) 296 ! 297 ALLOCATE( tmask_i(jpi,jpj) , tmask_h(jpi,jpj) , & 350 298 & ssmask (jpi,jpj) , ssumask(jpi,jpj) , ssvmask(jpi,jpj) , ssfmask(jpi,jpj) , & 351 299 & mbkt (jpi,jpj) , mbku (jpi,jpj) , mbkv (jpi,jpj) , STAT=ierr(9) ) 352 300 ! 353 301 ALLOCATE( misfdep(jpi,jpj) , mikt(jpi,jpj) , miku(jpi,jpj) , & 354 302 & risfdep(jpi,jpj) , mikv(jpi,jpj) , mikf(jpi,jpj) , STAT=ierr(10) ) 355 303 ! 356 304 ALLOCATE( tmask(jpi,jpj,jpk) , umask(jpi,jpj,jpk) , & 357 305 & vmask(jpi,jpj,jpk) , fmask(jpi,jpj,jpk) , STAT=ierr(11) ) 358 306 ! 359 307 ALLOCATE( wmask(jpi,jpj,jpk) , wumask(jpi,jpj,jpk), wvmask(jpi,jpj,jpk) , STAT=ierr(12) ) 360 308 ! -
branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90
r6624 r6667 20 20 !! dom_nam : read and contral domain namelists 21 21 !! dom_ctl : control print for the ocean domain 22 !! dom_stiff : diagnose maximum grid stiffness/hydrostatic consistency (s-coordinate)23 22 !! cfg_wri : create the "domain_cfg.nc" file containing all required configuration information 24 23 !!---------------------------------------------------------------------- … … 67 66 !! and scale factors, and the coriolis factor 68 67 !! - dom_zgr: define the vertical coordinate and the bathymetry 69 !! - dom_wri: create the meshmask file if n msh=168 !! - dom_wri: create the meshmask file if nn_msh=1 70 69 !! - 1D configuration, move Coriolis, u and v at T-point 71 70 !!---------------------------------------------------------------------- 72 INTEGER :: jk ! dummy loop indices 73 INTEGER :: iconf = 0 ! local integers 74 CHARACTER (len=64) :: cform = "(A12, 3(A13, I7) )" 75 REAL(wp), POINTER, DIMENSION(:,:) :: z1_hu_0, z1_hv_0 71 INTEGER :: ji, jj, jk ! dummy loop indices 72 INTEGER :: iconf = 0 ! local integers 73 CHARACTER (len=64) :: cform = "(A12, 3(A13, I7))" 74 INTEGER, DIMENSION(jpi,jpj) :: ik_top, ik_bot ! top and bottom ocean level 75 REAL(wp), POINTER, DIMENSION(:,:) :: zht, z1_hu_0, z1_hv_0 76 76 !!---------------------------------------------------------------------- 77 77 ! 78 78 IF( nn_timing == 1 ) CALL timing_start('dom_init') 79 79 ! 80 IF(lwp) THEN 80 IF(lwp) THEN ! Ocean domain Parameters (control print) 81 81 WRITE(numout,*) 82 82 WRITE(numout,*) 'dom_init : domain initialization' 83 83 WRITE(numout,*) '~~~~~~~~' 84 ENDIF 85 ! 86 ! Ocean domain Parameters (control print) 87 ! ----------------------- 88 IF(lwp) THEN 84 ! 89 85 WRITE(numout,*) ' Domain info' 90 86 WRITE(numout,*) ' dimension of model' … … 100 96 WRITE(numout,*) ' lateral boundary of the Global domain : jperio = ', jperio 101 97 ENDIF 102 103 ! !== Reference coordinate system ==! 104 ! 105 CALL dom_nam ! read namelist ( namrun, namdom ) 106 CALL dom_clo ! Closed seas and lake 107 CALL dom_hgr ! Horizontal mesh 108 CALL dom_zgr ! Vertical mesh and bathymetry 109 CALL dom_msk ! Masks 110 ! 111 IF( ln_sco ) CALL dom_stiff ! Maximum stiffness ratio/hydrostatic consistency 98 ! 99 ! !== Reference coordinate system ==! 100 ! 101 CALL dom_nam ! read namelist ( namrun, namdom ) 102 CALL dom_clo ! Closed seas and lake 103 CALL dom_hgr ! Horizontal mesh 104 CALL dom_zgr( ik_top, ik_bot ) ! Vertical mesh and bathymetry 105 CALL dom_msk( ik_top, ik_bot ) ! Masks 106 ! 107 IF( ln_sco ) CALL dom_stiff ! Maximum stiffness ratio/hydrostatic consistency 108 ! 109 DO jj = 1, jpj ! depth of the iceshelves 110 DO ji = 1, jpj 111 risfdep(ji,jj) = gdepw_0(ji,jj,mikt(ji,jj)) 112 END DO 113 END DO 112 114 ! 113 115 ht_0(:,:) = e3t_0(:,:,1) * tmask(:,:,1) ! Reference ocean thickness … … 120 122 END DO 121 123 ! 122 ! 124 ! !== time varying part of coordinate system ==! 123 125 ! 124 126 IF( ln_linssh ) THEN ! Fix in time : set to the reference one for all … … 158 160 IF( lk_c1d ) CALL cor_c1d ! 1D configuration: Coriolis set at T-point 159 161 ! 160 IF( n msh > 0 .AND. .NOT. ln_iscpl ) CALL dom_wri ! Create a domain file161 IF( n msh > 0 .AND. ln_iscpl .AND. .NOT. ln_rstart ) CALL dom_wri ! Create a domain file162 IF( nn_msh > 0 .AND. .NOT. ln_iscpl ) CALL dom_wri ! Create a domain file 163 IF( nn_msh > 0 .AND. ln_iscpl .AND. .NOT. ln_rstart ) CALL dom_wri ! Create a domain file 162 164 IF( .NOT.ln_rstart ) CALL dom_ctl ! Domain control 163 165 ! … … 165 167 IF(lwp) THEN 166 168 WRITE(numout,*) 167 WRITE(numout,*) 'dom_init : end of domain initialization n msh=', nmsh169 WRITE(numout,*) 'dom_init : end of domain initialization nn_msh=', nn_msh 168 170 WRITE(numout,*) 169 171 ENDIF 170 172 ! 171 IF( nmsh == -1) CALL cfg_wri ! create the configuration file173 IF( ln_write_cfg ) CALL cfg_wri ! create the configuration file 172 174 ! 173 175 IF( nn_timing == 1 ) CALL timing_stop('dom_init') … … 192 194 & nn_stock, nn_write , ln_mskland , ln_clobber , nn_chunksz, nn_euler , & 193 195 & ln_cfmeta, ln_iscpl 194 NAMELIST/namdom/ nn_bathy, rn_bathy , rn_e3zps_min, rn_e3zps_rat, nn_msh, rn_hmin, rn_isfhmin, & 195 & rn_atfp , rn_rdt , nn_closea , ln_crs , & 196 & ppsur, ppa0, ppa1, ppkth, ppacr, ppdzmin, pphmax, ldbletanh, & 197 & ppa2, ppkth2, ppacr2 196 NAMELIST/namdom/ ln_linssh, nn_closea, nn_msh, rn_isfhmin, rn_rdt, rn_atfp, ln_crs 198 197 #if defined key_netcdf4 199 198 NAMELIST/namnc4/ nn_nchunks_i, nn_nchunks_j, nn_nchunks_k, ln_nc4zip … … 261 260 neuler = 0 262 261 ENDIF 263 264 262 ! ! control of output frequency 265 263 IF ( nstock == 0 .OR. nstock > nitend ) THEN … … 305 303 WRITE(numout,*) 306 304 WRITE(numout,*) ' Namelist namdom : space & time domain' 307 WRITE(numout,*) ' flag read/compute bathymetry nn_bathy = ', nn_bathy 308 WRITE(numout,*) ' Depth (if =0 bathy=jpkm1) rn_bathy = ', rn_bathy 309 WRITE(numout,*) ' min depth of the ocean (>0) or rn_hmin = ', rn_hmin 310 WRITE(numout,*) ' min number of ocean level (<0) ' 311 WRITE(numout,*) ' treshold to open the isf cavity rn_isfhmin = ', rn_isfhmin, ' (m)' 312 WRITE(numout,*) ' minimum thickness of partial rn_e3zps_min = ', rn_e3zps_min, ' (m)' 313 WRITE(numout,*) ' step level rn_e3zps_rat = ', rn_e3zps_rat 314 WRITE(numout,*) ' create mesh/mask file(s) nn_msh = ', nn_msh 305 WRITE(numout,*) ' linear free surface (=T) ln_linssh = ', ln_linssh 306 WRITE(numout,*) ' suppression of closed seas (=0) nn_closea = ', nn_closea 307 WRITE(numout,*) ' create mesh/mask file(s) nn_msh = ', nn_msh 315 308 WRITE(numout,*) ' = 0 no file created ' 316 309 WRITE(numout,*) ' = 1 mesh_mask ' 317 310 WRITE(numout,*) ' = 2 mesh and mask ' 318 311 WRITE(numout,*) ' = 3 mesh_hgr, msh_zgr and mask' 319 WRITE(numout,*) ' ocean time step rn_rdt = ', rn_rdt 320 WRITE(numout,*) ' asselin time filter parameter rn_atfp = ', rn_atfp 321 WRITE(numout,*) ' suppression of closed seas (=0) nn_closea = ', nn_closea 322 WRITE(numout,*) ' online coarsening of dynamical fields ln_crs = ', ln_crs 323 324 WRITE(numout,*) ' ORCA r4, r2 and r05 coefficients ppsur = ', ppsur 325 WRITE(numout,*) ' ppa0 = ', ppa0 326 WRITE(numout,*) ' ppa1 = ', ppa1 327 WRITE(numout,*) ' ppkth = ', ppkth 328 WRITE(numout,*) ' ppacr = ', ppacr 329 WRITE(numout,*) ' Minimum vertical spacing ppdzmin = ', ppdzmin 330 WRITE(numout,*) ' Maximum depth pphmax = ', pphmax 331 WRITE(numout,*) ' Use double tanf function for vertical coordinates ldbletanh = ', ldbletanh 332 WRITE(numout,*) ' Double tanh function parameters ppa2 = ', ppa2 333 WRITE(numout,*) ' ppkth2 = ', ppkth2 334 WRITE(numout,*) ' ppacr2 = ', ppacr2 312 WRITE(numout,*) ' treshold to open the isf cavity rn_isfhmin = ', rn_isfhmin, ' (m)' 313 WRITE(numout,*) ' ocean time step rn_rdt = ', rn_rdt 314 WRITE(numout,*) ' asselin time filter parameter rn_atfp = ', rn_atfp 315 WRITE(numout,*) ' online coarsening of dynamical fields ln_crs = ', ln_crs 335 316 ENDIF 336 317 337 318 call flush( numout ) 338 319 ! 339 ntopo = nn_bathy ! conversion DOCTOR names into model names (this should disappear soon) 340 e3zps_min = rn_e3zps_min 341 e3zps_rat = rn_e3zps_rat 342 nmsh = nn_msh 320 ! ! ! conversion DOCTOR names into model names (this should disappear soon) 343 321 atfp = rn_atfp 344 322 rdt = rn_rdt … … 373 351 snc4set%luse = .FALSE. ! No NetCDF 4 case 374 352 #endif 375 376 call flush( numout )377 378 353 ! 379 354 END SUBROUTINE dom_nam … … 403 378 ze1max = MAXVAL( e1t(:,:), mask = tmask_i(:,:) == 1._wp ) 404 379 ze2max = MAXVAL( e2t(:,:), mask = tmask_i(:,:) == 1._wp ) 405 380 ! 406 381 iloc = MINLOC( e1t(:,:), mask = tmask_i(:,:) == 1._wp ) 407 382 iimi1 = iloc(1) + nimpp - 1 … … 430 405 431 406 432 SUBROUTINE dom_stiff433 !!----------------------------------------------------------------------434 !! *** ROUTINE dom_stiff ***435 !!436 !! ** Purpose : Diagnose maximum grid stiffness/hydrostatic consistency437 !!438 !! ** Method : Compute Haney (1991) hydrostatic condition ratio439 !! Save the maximum in the vertical direction440 !! (this number is only relevant in s-coordinates)441 !!442 !! Haney, R. L., 1991: On the pressure gradient force443 !! over steep topography in sigma coordinate ocean models.444 !! J. Phys. Oceanogr., 21, 610-619.445 !!----------------------------------------------------------------------446 INTEGER :: ji, jj, jk447 REAL(wp) :: zrxmax448 REAL(wp), DIMENSION(4) :: zr1449 !!----------------------------------------------------------------------450 rx1(:,:) = 0._wp451 zrxmax = 0._wp452 zr1(:) = 0._wp453 !454 DO ji = 2, jpim1455 DO jj = 2, jpjm1456 DO jk = 1, jpkm1457 zr1(1) = ABS( ( gdepw_0(ji ,jj,jk )-gdepw_0(ji-1,jj,jk ) &458 & +gdepw_0(ji ,jj,jk+1)-gdepw_0(ji-1,jj,jk+1) ) &459 & / ( gdepw_0(ji ,jj,jk )+gdepw_0(ji-1,jj,jk ) &460 & -gdepw_0(ji ,jj,jk+1)-gdepw_0(ji-1,jj,jk+1) + rsmall ) ) * umask(ji-1,jj,jk)461 zr1(2) = ABS( ( gdepw_0(ji+1,jj,jk )-gdepw_0(ji ,jj,jk ) &462 & +gdepw_0(ji+1,jj,jk+1)-gdepw_0(ji ,jj,jk+1) ) &463 & / ( gdepw_0(ji+1,jj,jk )+gdepw_0(ji ,jj,jk ) &464 & -gdepw_0(ji+1,jj,jk+1)-gdepw_0(ji ,jj,jk+1) + rsmall ) ) * umask(ji ,jj,jk)465 zr1(3) = ABS( ( gdepw_0(ji,jj+1,jk )-gdepw_0(ji,jj ,jk ) &466 & +gdepw_0(ji,jj+1,jk+1)-gdepw_0(ji,jj ,jk+1) ) &467 & / ( gdepw_0(ji,jj+1,jk )+gdepw_0(ji,jj ,jk ) &468 & -gdepw_0(ji,jj+1,jk+1)-gdepw_0(ji,jj ,jk+1) + rsmall ) ) * vmask(ji,jj ,jk)469 zr1(4) = ABS( ( gdepw_0(ji,jj ,jk )-gdepw_0(ji,jj-1,jk ) &470 & +gdepw_0(ji,jj ,jk+1)-gdepw_0(ji,jj-1,jk+1) ) &471 & / ( gdepw_0(ji,jj ,jk )+gdepw_0(ji,jj-1,jk ) &472 & -gdepw_0(ji,jj ,jk+1)-gdepw_0(ji,jj-1,jk+1) + rsmall ) ) * vmask(ji,jj-1,jk)473 zrxmax = MAXVAL( zr1(1:4) )474 rx1(ji,jj) = MAX( rx1(ji,jj) , zrxmax )475 END DO476 END DO477 END DO478 CALL lbc_lnk( rx1, 'T', 1. )479 !480 zrxmax = MAXVAL( rx1 )481 !482 IF( lk_mpp ) CALL mpp_max( zrxmax ) ! max over the global domain483 !484 IF(lwp) THEN485 WRITE(numout,*)486 WRITE(numout,*) 'dom_stiff : maximum grid stiffness ratio: ', zrxmax487 WRITE(numout,*) '~~~~~~~~~'488 ENDIF489 !490 END SUBROUTINE dom_stiff491 492 493 407 SUBROUTINE cfg_wri 494 408 !!---------------------------------------------------------------------- … … 503 417 !! domhgr, domzgr, and dommsk. Note: the file contain depends on 504 418 !! the vertical coord. used (z-coord, partial steps, s-coord) 505 !! MOD(n msh, 3) = 1 : 'mesh_mask.nc' file419 !! MOD(nn_msh, 3) = 1 : 'mesh_mask.nc' file 506 420 !! = 2 : 'mesh.nc' and mask.nc' files 507 421 !! = 0 : 'mesh_hgr.nc', 'mesh_zgr.nc' and … … 510 424 !! vertical coordinate. 511 425 !! 512 !! if n msh <= 3: write full 3D arrays for e3[tuvw] and gdep[tuvw]513 !! if 3 < n msh <= 6: write full 3D arrays for e3[tuvw] and 2D arrays426 !! if nn_msh <= 3: write full 3D arrays for e3[tuvw] and gdep[tuvw] 427 !! if 3 < nn_msh <= 6: write full 3D arrays for e3[tuvw] and 2D arrays 514 428 !! corresponding to the depth of the bottom t- and w-points 515 !! if 6 < n msh <= 9: write 2D arrays corresponding to the depth and the429 !! if 6 < nn_msh <= 9: write 2D arrays corresponding to the depth and the 516 430 !! thickness (e3[tw]_ps) of the bottom points 517 431 !! … … 523 437 INTEGER :: inum ! temprary units for 'domain_cfg.nc' file 524 438 CHARACTER(len=21) :: clnam ! filename (mesh and mask informations) 439 REAL(wp), DIMENSION(jpi,jpj) :: z2d ! workspace 525 440 !!---------------------------------------------------------------------- 526 441 ! … … 537 452 538 453 ! !== global domain size ==! 454 ! 539 455 CALL iom_rstput( 0, 0, inum, 'jpiglo', REAL( jpiglo, wp), ktype = jp_i4 ) 540 456 CALL iom_rstput( 0, 0, inum, 'jpjglo', REAL( jpjglo, wp), ktype = jp_i4 ) 541 457 CALL iom_rstput( 0, 0, inum, 'jpkglo', REAL( jpk , wp), ktype = jp_i4 ) 542 458 ! 543 459 ! !== domain characteristics ==! 460 ! 544 461 ! ! lateral boundary of the global domain 545 462 CALL iom_rstput( 0, 0, inum, 'jperio', REAL( jperio, wp), ktype = jp_i4 ) 463 ! 546 464 ! ! type of vertical coordinate 547 465 IF( ln_zco ) THEN ; izco = 1 ; ELSE ; izco = 0 ; ENDIF … … 551 469 CALL iom_rstput( 0, 0, inum, 'ln_zps' , REAL( izps, wp), ktype = jp_i4 ) 552 470 CALL iom_rstput( 0, 0, inum, 'ln_sco' , REAL( isco, wp), ktype = jp_i4 ) 471 ! 553 472 ! ! ocean cavities under iceshelves 554 473 IF( ln_isfcav ) THEN ; icav = 1 ; ELSE ; icav = 0 ; ENDIF 555 474 CALL iom_rstput( 0, 0, inum, 'ln_isfcav', REAL( icav, wp), ktype = jp_i4 ) 556 475 ! 557 476 ! !== horizontal mesh ! 558 477 ! … … 579 498 CALL iom_rstput( 0, 0, inum, 'ff_f' , ff_f , ktype = jp_r8 ) ! coriolis factor 580 499 CALL iom_rstput( 0, 0, inum, 'ff_t' , ff_t , ktype = jp_r8 ) 581 582 500 ! 583 501 ! !== vertical mesh - 3D mask ==! 584 502 ! … … 594 512 CALL iom_rstput( 0, 0, inum, 'e3u_0' , e3u_0 , ktype = jp_r8 ) 595 513 CALL iom_rstput( 0, 0, inum, 'e3v_0' , e3v_0 , ktype = jp_r8 ) 514 CALL iom_rstput( 0, 0, inum, 'e3f_0' , e3f_0 , ktype = jp_r8 ) 596 515 CALL iom_rstput( 0, 0, inum, 'e3w_0' , e3w_0 , ktype = jp_r8 ) 597 ! 598 CALL iom_rstput( 0, 0, inum, 'tmask' , tmask , ktype = jp_i1 ) ! masks (in bytes) 599 CALL iom_rstput( 0, 0, inum, 'umask' , umask , ktype = jp_i1 ) 600 CALL iom_rstput( 0, 0, inum, 'vmask' , vmask , ktype = jp_i1 ) 601 CALL iom_rstput( 0, 0, inum, 'fmask' , fmask , ktype = jp_i1 ) 602 603 !!gm Probably not required fields : 604 CALL iom_rstput( 0, 0, inum, 'bathy' , bathy , ktype = jp_r8 ) ! depth of the ocean at T-points 605 CALL iom_rstput( 0, 0, inum, 'mbathy' , REAL( mbathy , wp ) , ktype = jp_i4 ) ! nb of ocean T-points 606 CALL iom_rstput( 0, 0, inum, 'mbkt' , REAL( mbkt , wp ) , ktype = jp_i4 ) ! nb of ocean T-points 607 608 ! 609 CALL iom_rstput( 0, 0, inum, 'bathy' , risfdep , ktype = jp_r8 ) ! depth of the iceshelves at T-points 610 CALL iom_rstput( 0, 0, inum, 'misfdep', REAL( misfdep, wp ) , ktype = jp_i4 ) ! nb of ocean T-points (ISF) 611 612 !!gm end 613 614 !!gm ? 615 ! CALL iom_rstput( 0, 0, inum, 'hbatt', hbatt ) 616 ! CALL iom_rstput( 0, 0, inum, 'hbatu', hbatu ) 617 ! CALL iom_rstput( 0, 0, inum, 'hbatv', hbatv ) 618 ! CALL iom_rstput( 0, 0, inum, 'hbatf', hbatf ) 619 !!gm ? 620 !!gm ? 621 ! CALL iom_rstput( 0, 0, inum, 'rx1' , rx1 ) ! Max. grid stiffness ratio 622 !!gm ? 623 ! 624 ! DO jk = 1,jpk 625 ! DO jj = 1, jpjm1 626 ! DO ji = 1, fs_jpim1 ! vector opt. 627 ! zdepu(ji,jj,jk) = MIN( gdept_0(ji,jj,jk) , gdept_0(ji+1,jj ,jk) ) 628 ! zdepv(ji,jj,jk) = MIN( gdept_0(ji,jj,jk) , gdept_0(ji ,jj+1,jk) ) 629 ! END DO 630 ! END DO 631 ! END DO 632 ! CALL lbc_lnk( zdepu, 'U', 1. ) ; CALL lbc_lnk( zdepv, 'V', 1. ) 633 ! CALL iom_rstput( 0, 0, inum, 'gdepu' , zdepu , ktype = jp_r8 ) 634 ! CALL iom_rstput( 0, 0, inum, 'gdepv' , zdepv , ktype = jp_r8 ) 635 516 CALL iom_rstput( 0, 0, inum, 'e3uw_0' , e3uw_0 , ktype = jp_r8 ) 517 CALL iom_rstput( 0, 0, inum, 'e3vw_0' , e3vw_0 , ktype = jp_r8 ) 518 ! 519 ! !== ocean top and bottom level ==! 520 ! 521 CALL iom_rstput( 0, 0, inum, 'bottom level' , REAL( mbkt, wp )*ssmask , ktype = jp_i4 ) ! nb of ocean T-points 522 CALL iom_rstput( 0, 0, inum, 'top level' , REAL( mikt, wp )*ssmask , ktype = jp_i4 ) ! nb of ocean T-points (ISF) 523 ! 524 IF( ln_sco ) THEN ! s-coordinate: store grid stiffness ratio (Not required anyway) 525 CALL dom_stiff( z2d ) 526 CALL iom_rstput( 0, 0, inum, 'stiffness', z2d ) ! ! Max. grid stiffness ratio 527 ENDIF 528 ! 636 529 ! ! ============================ 637 530 ! ! close the files -
branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90
r6624 r6667 91 91 IF( ln_read_cfg ) THEN !== read in mesh_mask.nc file ==! 92 92 IF(lwp) WRITE(numout,*) 93 IF(lwp) WRITE(numout,*) ' read horizontal mesh in " mesh_mask" file'93 IF(lwp) WRITE(numout,*) ' read horizontal mesh in "domain_cfg" file' 94 94 ! 95 95 CALL hgr_read ( glamt , glamu , glamv , glamf , & ! geographic position (required) -
branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90
r6140 r6667 9 9 !! - ! 1996-05 (G. Madec) mask computed from tmask 10 10 !! 8.0 ! 1997-02 (G. Madec) mesh information put in domhgr.F 11 !! 8.1 ! 1997-07 (G. Madec) modification of mbathyand fmask11 !! 8.1 ! 1997-07 (G. Madec) modification of kbat and fmask 12 12 !! - ! 1998-05 (G. Roullet) free surface 13 13 !! 8.2 ! 2000-03 (G. Madec) no slip accurate … … 50 50 CONTAINS 51 51 52 SUBROUTINE dom_msk 52 SUBROUTINE dom_msk( k_top, k_bot ) 53 53 !!--------------------------------------------------------------------- 54 54 !! *** ROUTINE dom_msk *** … … 57 57 !! zontal velocity points (u & v), vorticity points (f) points. 58 58 !! 59 !! ** Method : The ocean/land mask is computed from the basin bathy- 60 !! metry in level (mbathy) which is defined or read in dommba. 61 !! mbathy equals 0 over continental T-point 62 !! and the number of ocean level over the ocean. 63 !! 64 !! At a given position (ji,jj,jk) the ocean/land mask is given by: 65 !! t-point : 0. IF mbathy( ji ,jj) =< 0 66 !! 1. IF mbathy( ji ,jj) >= jk 67 !! u-point : 0. IF mbathy( ji ,jj) or mbathy(ji+1, jj ) =< 0 68 !! 1. IF mbathy( ji ,jj) and mbathy(ji+1, jj ) >= jk. 69 !! v-point : 0. IF mbathy( ji ,jj) or mbathy( ji ,jj+1) =< 0 70 !! 1. IF mbathy( ji ,jj) and mbathy( ji ,jj+1) >= jk. 71 !! f-point : 0. IF mbathy( ji ,jj) or mbathy( ji ,jj+1) 72 !! or mbathy(ji+1,jj) or mbathy(ji+1,jj+1) =< 0 73 !! 1. IF mbathy( ji ,jj) and mbathy( ji ,jj+1) 74 !! and mbathy(ji+1,jj) and mbathy(ji+1,jj+1) >= jk. 75 !! tmask_i : interior ocean mask at t-point, i.e. excluding duplicated 76 !! rows/lines due to cyclic or North Fold boundaries as well 77 !! as MPP halos. 78 !! 79 !! The lateral friction is set through the value of fmask along 80 !! the coast and topography. This value is defined by rn_shlat, a 81 !! namelist parameter: 59 !! ** Method : The ocean/land mask at t-point is deduced from ko_top 60 !! and ko_bot, the indices of the fist and last ocean t-levels which 61 !! are either defined in usrdef_zgr or read in zgr_read. 62 !! The velocity masks (umask, vmask, wmask, wumask, wvmask) 63 !! are deduced from a product of the two neighboring tmask. 64 !! The vorticity mask (fmask) is deduced from tmask taking 65 !! into account the choice of lateral boundary condition (rn_shlat) : 82 66 !! rn_shlat = 0, free slip (no shear along the coast) 83 67 !! rn_shlat = 2, no slip (specified zero velocity at the coast) … … 85 69 !! 2 < rn_shlat, strong slip | in the lateral boundary layer 86 70 !! 87 !! N.B. If nperio not equal to 0, the land/ocean mask arrays 88 !! are defined with the proper value at lateral domain boundaries. 71 !! tmask_i : interior ocean mask at t-point, i.e. excluding duplicated 72 !! rows/lines due to cyclic or North Fold boundaries as well 73 !! as MPP halos. 74 !! tmask_h : halo mask at t-point, i.e. excluding duplicated rows/lines 75 !! due to cyclic or North Fold boundaries as well 76 !! as MPP halos. 89 77 !! 90 78 !! In case of open boundaries (lk_bdy=T): 91 !! - tmask is set to 1 on the points to be computed b ay the open79 !! - tmask is set to 1 on the points to be computed by the open 92 80 !! boundaries routines. 93 81 !! 94 !! ** Action : tmask : land/ocean mask at t-point (=0. or 1.) 95 !! umask : land/ocean mask at u-point (=0. or 1.) 96 !! vmask : land/ocean mask at v-point (=0. or 1.) 97 !! fmask : land/ocean mask at f-point (=0. or 1.) 98 !! =rn_shlat along lateral boundaries 99 !! tmask_i : interior ocean mask 82 !! ** Action : tmask, umask, vmask, wmask, wumask, wvmask : land/ocean mask 83 !! at t-, u-, v- w, wu-, and wv-points (=0. or 1.) 84 !! fmask : land/ocean mask at f-point (=0., or =1., or 85 !! =rn_shlat along lateral boundaries) 86 !! tmask_i : interior ocean mask 87 !! tmask_h : halo mask 88 !! ssmask , ssumask, ssvmask, ssfmask : 2D ocean mask 100 89 !!---------------------------------------------------------------------- 90 INTEGER, DIMENSION(:,:), INTENT(in) :: k_top, k_bot ! first and last ocean level 91 ! 101 92 INTEGER :: ji, jj, jk ! dummy loop indices 102 93 INTEGER :: iif, iil, ii0, ii1, ii ! local integers 103 94 INTEGER :: ijf, ijl, ij0, ij1 ! - - 95 INTEGER :: iktop, ikbot ! - - 104 96 INTEGER :: ios 105 97 INTEGER :: isrow ! index for ORCA1 starting row 106 INTEGER , POINTER, DIMENSION(:,:) :: imsk107 98 REAL(wp), POINTER, DIMENSION(:,:) :: zwf 108 99 !! … … 111 102 ! 112 103 IF( nn_timing == 1 ) CALL timing_start('dom_msk') 113 !114 CALL wrk_alloc( jpi, jpj, imsk )115 CALL wrk_alloc( jpi, jpj, zwf )116 104 ! 117 105 REWIND( numnam_ref ) ! Namelist namlbc in reference namelist : Lateral momentum boundary condition … … 142 130 ENDIF 143 131 144 ! 1. Ocean/land mask at t-point (computed from mbathy) 145 ! ----------------------------- 146 ! N.B. tmask has already the right boundary conditions since mbathy is ok 147 ! 148 tmask(:,:,:) = 0._wp 149 DO jk = 1, jpk 150 DO jj = 1, jpj 151 DO ji = 1, jpi 152 IF( REAL( mbathy(ji,jj) - jk, wp ) + 0.1_wp >= 0._wp ) tmask(ji,jj,jk) = 1._wp 153 END DO 132 133 ! Ocean/land mask at t-point (computed from ko_top and ko_bot) 134 ! ---------------------------- 135 ! 136 DO jj = 1, jpj 137 DO ji = 1, jpi 138 iktop = k_top(ji,jj) 139 ikbot = k_bot(ji,jj) 140 tmask(ji,jj, 1:iktop-1) = 0._wp ! mask the iceshelves 141 tmask(ji,jj,iktop :ikbot ) = 1._wp 142 tmask(ji,jj,ikbot+1:jpkglo ) = 0._wp ! mask the ocean topography 154 143 END DO 155 144 END DO 145 146 ! 2D ocean mask (=1 if at least one level of the water column is ocean, =0 otherwise) 147 WHERE( k_bot(:,:) > 0 ) ; ssmask(:,:) = 1._wp 148 ELSEWHERE ; ssmask(:,:) = 0._wp 149 END WHERE 156 150 157 ! (ISF) define barotropic mask and mask the ice shelf point158 ssmask(:,:)=tmask(:,:,1) ! at this stage ice shelf is not masked159 151 160 DO jk = 1, jpk 161 DO jj = 1, jpj 162 DO ji = 1, jpi 163 IF( REAL( misfdep(ji,jj) - jk, wp ) - 0.1_wp >= 0._wp ) THEN 164 tmask(ji,jj,jk) = 0._wp 165 END IF 166 END DO 167 END DO 168 END DO 169 170 ! Interior domain mask (used for global sum) 152 ! Interior domain mask (used for global sum) 171 153 ! -------------------- 172 tmask_i(:,:) = ssmask(:,:) ! (ISH) tmask_i = 1 even on the ice shelf 173 174 tmask_h(:,:) = 1._wp ! 0 on the halo and 1 elsewhere 175 iif = jpreci ! ??? 176 iil = nlci - jpreci + 1 177 ijf = jprecj ! ??? 178 ijl = nlcj - jprecj + 1 179 154 ! 155 iif = jpreci ; iil = nlci - jpreci + 1 156 ijf = jprecj ; ijl = nlcj - jprecj + 1 157 ! 158 ! ! halo mask : 0 on the halo and 1 elsewhere 159 tmask_h(:,:) = 1._wp 180 160 tmask_h( 1 :iif, : ) = 0._wp ! first columns 181 161 tmask_h(iil:jpi, : ) = 0._wp ! last columns (including mpp extra columns) 182 162 tmask_h( : , 1 :ijf) = 0._wp ! first rows 183 163 tmask_h( : ,ijl:jpj) = 0._wp ! last rows (including mpp extra rows) 184 185 ! north fold mask 186 ! --------------- 164 ! 165 ! ! north fold mask 187 166 tpol(1:jpiglo) = 1._wp 188 167 fpol(1:jpiglo) = 1._wp … … 190 169 tpol(jpiglo/2+1:jpiglo) = 0._wp 191 170 fpol( 1 :jpiglo) = 0._wp 192 IF( mjg(nlej) == jpjglo ) THEN ! only half of the nlcj-1 row 171 IF( mjg(nlej) == jpjglo ) THEN ! only half of the nlcj-1 row for tmask_h 193 172 DO ji = iif+1, iil-1 194 173 tmask_h(ji,nlej-1) = tmask_h(ji,nlej-1) * tpol(mig(ji)) … … 196 175 ENDIF 197 176 ENDIF 198 199 tmask_i(:,:) = tmask_i(:,:) * tmask_h(:,:) 200 177 ! 201 178 IF( jperio == 5 .OR. jperio == 6 ) THEN ! F-point pivot 202 179 tpol( 1 :jpiglo) = 0._wp 203 180 fpol(jpiglo/2+1:jpiglo) = 0._wp 204 181 ENDIF 205 206 ! 2. Ocean/land mask at u-, v-, and z-points (computed from tmask) 207 ! ------------------------------------------- 182 ! 183 ! ! interior mask : 2D ocean mask x halo mask 184 tmask_i(:,:) = ssmask(:,:) * tmask_h(:,:) 185 186 187 ! Ocean/land mask at u-, v-, and z-points (computed from tmask) 188 ! ---------------------------------------- 208 189 DO jk = 1, jpk 209 190 DO jj = 1, jpjm1 … … 221 202 DO jj = 1, jpjm1 222 203 DO ji = 1, fs_jpim1 ! vector loop 204 !!gm simpler : 205 ! ssumask(ji,jj) = MIN( 1._wp , SUM( umask(ji,jj,:) ) ) 206 ! ssvmask(ji,jj) = MIN( 1._wp , SUM( vmask(ji,jj,:) ) ) 207 !!gm 208 !!gm faster : 209 ! ssumask(ji,jj) = ssmask(ji,jj) * tmask(ji+1,jj ) 210 ! ssvmask(ji,jj) = ssmask(ji,jj) * tmask(ji ,jj+1) 211 !!gm 223 212 ssumask(ji,jj) = ssmask(ji,jj) * ssmask(ji+1,jj ) * MIN(1._wp,SUM(umask(ji,jj,:))) 224 213 ssvmask(ji,jj) = ssmask(ji,jj) * ssmask(ji ,jj+1) * MIN(1._wp,SUM(vmask(ji,jj,:))) 214 !!end 225 215 END DO 226 216 DO ji = 1, jpim1 ! NO vector opt. 217 !!gm faster 218 ! ssfmask(ji,jj) = ssmask(ji,jj ) * ssmask(ji+1,jj ) & 219 ! & * ssmask(ji,jj+1) * ssmask(ji+1,jj+1) 220 !!gm 227 221 ssfmask(ji,jj) = ssmask(ji,jj ) * ssmask(ji+1,jj ) & 228 222 & * ssmask(ji,jj+1) * ssmask(ji+1,jj+1) * MIN(1._wp,SUM(fmask(ji,jj,:))) 223 !!gm 229 224 END DO 230 225 END DO 231 226 CALL lbc_lnk( umask , 'U', 1._wp ) ! Lateral boundary conditions 232 227 CALL lbc_lnk( vmask , 'V', 1._wp ) 233 CALL lbc_lnk( fmask , 'F', 1._wp ) 234 CALL lbc_lnk( ssumask, 'U', 1._wp ) ! Lateral boundary conditions228 ! CALL lbc_lnk( fmask , 'F', 1._wp ) ! applied after the specification of lateral b.c. 229 CALL lbc_lnk( ssumask, 'U', 1._wp ) 235 230 CALL lbc_lnk( ssvmask, 'V', 1._wp ) 236 231 CALL lbc_lnk( ssfmask, 'F', 1._wp ) 237 232 238 ! 3. Ocean/land mask at wu-, wv- and w points 233 234 ! Ocean/land mask at wu-, wv- and w points 239 235 !---------------------------------------------- 240 236 wmask (:,:,1) = tmask(:,:,1) ! surface … … 247 243 END DO 248 244 245 249 246 ! Lateral boundary conditions on velocity (modify fmask) 250 247 ! --------------------------------------- 248 CALL wrk_alloc( jpi,jpj, zwf ) 249 ! 251 250 DO jk = 1, jpk 252 251 zwf(:,:) = fmask(:,:,jk) … … 277 276 END DO 278 277 ! 278 CALL wrk_dealloc( jpi,jpj, zwf ) 279 ! 279 280 IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! ORCA_R2 configuration 280 281 ! ! Increased lateral friction near of some straits … … 349 350 ! CAUTION : The fmask may be further modified in dyn_vor_init ( dynvor.F90 ) 350 351 ! 351 CALL wrk_dealloc( jpi, jpj, imsk )352 CALL wrk_dealloc( jpi, jpj, zwf )353 !354 352 IF( nn_timing == 1 ) CALL timing_stop('dom_msk') 355 353 ! -
branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90
r6624 r6667 885 885 e3t_n(ji,jj,:) = 0.5_wp * rn_wdmin1 886 886 e3t_a(ji,jj,:) = 0.5_wp * rn_wdmin1 887 sshb(ji,jj) = rn_wdmin1 - bathy(ji,jj)888 sshn(ji,jj) = rn_wdmin1 - bathy(ji,jj)889 ssha(ji,jj) = rn_wdmin1 - bathy(ji,jj)887 sshb(ji,jj) = rn_wdmin1 - ht_0(ji,jj) !!gm I don't understand that ! 888 sshn(ji,jj) = rn_wdmin1 - ht_0(ji,jj) 889 ssha(ji,jj) = rn_wdmin1 - ht_0(ji,jj) 890 890 ENDIF 891 891 ENDDO … … 894 894 895 895 IF( ln_vvl_ztilde .OR. ln_vvl_layer) THEN 896 tilde_e3t_b(:,:,:) = 0. 0_wp897 tilde_e3t_n(:,:,:) = 0. 0_wp896 tilde_e3t_b(:,:,:) = 0._wp 897 tilde_e3t_n(:,:,:) = 0._wp 898 898 IF( ln_vvl_ztilde ) hdiv_lf(:,:,:) = 0._wp 899 899 END IF -
branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/DOM/domwri.F90
r6624 r6667 13 13 !! dom_wri : create and write mesh and mask file(s) 14 14 !! dom_uniq : identify unique point of a grid (TUVF) 15 !! dom_stiff : diagnose maximum grid stiffness/hydrostatic consistency (s-coordinate) 15 16 !!---------------------------------------------------------------------- 16 17 USE dom_oce ! ocean space and time domain … … 27 28 PUBLIC dom_wri ! routine called by inidom.F90 28 29 PUBLIC dom_wri_coordinate ! routine called by domhgr.F90 29 30 PUBLIC dom_stiff ! routine called by inidom.F90 31 30 32 !! * Substitutions 31 33 # include "vectopt_loop_substitute.h90" … … 114 116 !! domhgr, domzgr, and dommsk. Note: the file contain depends on 115 117 !! the vertical coord. used (z-coord, partial steps, s-coord) 116 !! MOD(n msh, 3) = 1 : 'mesh_mask.nc' file118 !! MOD(nn_msh, 3) = 1 : 'mesh_mask.nc' file 117 119 !! = 2 : 'mesh.nc' and mask.nc' files 118 120 !! = 0 : 'mesh_hgr.nc', 'mesh_zgr.nc' and … … 121 123 !! vertical coordinate. 122 124 !! 123 !! if n msh <= 3: write full 3D arrays for e3[tuvw] and gdep[tuvw]124 !! if 3 < n msh <= 6: write full 3D arrays for e3[tuvw] and 2D arrays125 !! if nn_msh <= 3: write full 3D arrays for e3[tuvw] and gdep[tuvw] 126 !! if 3 < nn_msh <= 6: write full 3D arrays for e3[tuvw] and 2D arrays 125 127 !! corresponding to the depth of the bottom t- and w-points 126 !! if 6 < n msh <= 9: write 2D arrays corresponding to the depth and the128 !! if 6 < nn_msh <= 9: write 2D arrays corresponding to the depth and the 127 129 !! thickness (e3[tw]_ps) of the bottom points 128 130 !! … … 149 151 IF( nn_timing == 1 ) CALL timing_start('dom_wri') 150 152 ! 151 CALL wrk_alloc( jpi, jpj, zprt, zprw)152 CALL wrk_alloc( jpi, jpj, jpk,zdepu, zdepv )153 CALL wrk_alloc( jpi,jpj, zprt , zprw ) 154 CALL wrk_alloc( jpi,jpj,jpk, zdepu, zdepv ) 153 155 ! 154 156 IF(lwp) WRITE(numout,*) … … 162 164 clnam4 = 'mesh_zgr' ! filename (vertical mesh informations) 163 165 164 SELECT CASE ( MOD(n msh, 3) )166 SELECT CASE ( MOD(nn_msh, 3) ) 165 167 ! ! ============================ 166 168 CASE ( 1 ) ! create 'mesh_mask.nc' file … … 201 203 IF( ln_zps ) THEN ; izps = 1 ; ELSE ; izps = 0 ; ENDIF 202 204 IF( ln_sco ) THEN ; isco = 1 ; ELSE ; isco = 0 ; ENDIF 203 CALL iom_rstput( 0, 0, inum , 'ln_zco' , REAL( izco, wp), ktype = jp_i4 )204 CALL iom_rstput( 0, 0, inum , 'ln_zps' , REAL( izps, wp), ktype = jp_i4 )205 CALL iom_rstput( 0, 0, inum , 'ln_sco' , REAL( isco, wp), ktype = jp_i4 )205 CALL iom_rstput( 0, 0, inum2, 'ln_zco' , REAL( izco, wp), ktype = jp_i4 ) 206 CALL iom_rstput( 0, 0, inum2, 'ln_zps' , REAL( izps, wp), ktype = jp_i4 ) 207 CALL iom_rstput( 0, 0, inum2, 'ln_sco' , REAL( isco, wp), ktype = jp_i4 ) 206 208 ! ! ocean cavities under iceshelves 207 209 IF( ln_isfcav ) THEN ; icav = 1 ; ELSE ; icav = 0 ; ENDIF 208 CALL iom_rstput( 0, 0, inum , 'ln_isfcav', REAL( icav, wp), ktype = jp_i4 )210 CALL iom_rstput( 0, 0, inum2, 'ln_isfcav', REAL( icav, wp), ktype = jp_i4 ) 209 211 210 212 ! ! masks (inum2) … … 280 282 281 283 IF( ln_sco ) THEN ! s-coordinate 282 CALL iom_rstput( 0, 0, inum4, 'hbatt', hbatt )283 CALL iom_rstput( 0, 0, inum4, 'hbatu', hbatu )284 CALL iom_rstput( 0, 0, inum4, 'hbatv', hbatv )285 CALL iom_rstput( 0, 0, inum4, 'hbatf', hbatf )286 !287 CALL iom_rstput( 0, 0, inum4, 'gsigt', gsigt ) ! ! scaling coef.288 CALL iom_rstput( 0, 0, inum4, 'gsigw', gsigw )289 CALL iom_rstput( 0, 0, inum4, 'gsi3w', gsi3w )290 CALL iom_rstput( 0, 0, inum4, 'esigt', esigt )291 CALL iom_rstput( 0, 0, inum4, 'esigw', esigw )292 !293 284 CALL iom_rstput( 0, 0, inum4, 'e3t_0', e3t_0 ) ! ! scale factors 294 285 CALL iom_rstput( 0, 0, inum4, 'e3u_0', e3u_0 ) 295 286 CALL iom_rstput( 0, 0, inum4, 'e3v_0', e3v_0 ) 296 287 CALL iom_rstput( 0, 0, inum4, 'e3w_0', e3w_0 ) 297 CALL iom_rstput( 0, 0, inum4, 'rx1', rx1 ) ! ! Max. grid stiffness ratio 288 ! 289 CALL dom_stiff( zprt ) 290 CALL iom_rstput( 0, 0, inum4, 'stiffness', zprt ) ! ! Max. grid stiffness ratio 298 291 ! 299 292 CALL iom_rstput( 0, 0, inum4, 'gdept_1d' , gdept_1d ) ! ! stretched system … … 305 298 IF( ln_zps ) THEN ! z-coordinate - partial steps 306 299 ! 307 IF( n msh <= 6 ) THEN ! ! 3D vertical scale factors300 IF( nn_msh <= 6 ) THEN ! ! 3D vertical scale factors 308 301 CALL iom_rstput( 0, 0, inum4, 'e3t_0', e3t_0 ) 309 302 CALL iom_rstput( 0, 0, inum4, 'e3u_0', e3u_0 ) … … 321 314 END IF 322 315 ! 323 IF( n msh <= 3 ) THEN ! ! 3D depth316 IF( nn_msh <= 3 ) THEN ! ! 3D depth 324 317 CALL iom_rstput( 0, 0, inum4, 'gdept_0', gdept_0, ktype = jp_r8 ) 325 318 DO jk = 1,jpk … … 362 355 ! ! close the files 363 356 ! ! ============================ 364 SELECT CASE ( MOD(n msh, 3) )357 SELECT CASE ( MOD(nn_msh, 3) ) 365 358 CASE ( 1 ) 366 359 CALL iom_close( inum0 ) … … 425 418 END SUBROUTINE dom_uniq 426 419 420 421 SUBROUTINE dom_stiff( px1 ) 422 !!---------------------------------------------------------------------- 423 !! *** ROUTINE dom_stiff *** 424 !! 425 !! ** Purpose : Diagnose maximum grid stiffness/hydrostatic consistency 426 !! 427 !! ** Method : Compute Haney (1991) hydrostatic condition ratio 428 !! Save the maximum in the vertical direction 429 !! (this number is only relevant in s-coordinates) 430 !! 431 !! Haney, 1991, J. Phys. Oceanogr., 21, 610-619. 432 !!---------------------------------------------------------------------- 433 REAL(wp), DIMENSION(:,:), INTENT(out), OPTIONAL :: px1 ! stiffness 434 ! 435 INTEGER :: ji, jj, jk 436 REAL(wp) :: zrxmax 437 REAL(wp), DIMENSION(4) :: zr1 438 REAL(wp), DIMENSION(jpi,jpj) :: zx1 439 !!---------------------------------------------------------------------- 440 zx1(:,:) = 0._wp 441 zrxmax = 0._wp 442 zr1(:) = 0._wp 443 ! 444 DO ji = 2, jpim1 445 DO jj = 2, jpjm1 446 DO jk = 1, jpkm1 447 !!gm remark: dk(gdepw) = e3t ===>>> possible simplification of the following calculation.... 448 !! especially since it is gde3w which is used to compute the pressure gradient 449 !! furthermore, I think gdept_0 should be used below instead of w point in the numerator 450 !! so that the ratio is computed at the same point (i.e. uw and vw) .... 451 zr1(1) = ABS( ( gdepw_0(ji ,jj,jk )-gdepw_0(ji-1,jj,jk ) & 452 & +gdepw_0(ji ,jj,jk+1)-gdepw_0(ji-1,jj,jk+1) ) & 453 & / ( gdepw_0(ji ,jj,jk )+gdepw_0(ji-1,jj,jk ) & 454 & -gdepw_0(ji ,jj,jk+1)-gdepw_0(ji-1,jj,jk+1) + rsmall ) ) * umask(ji-1,jj,jk) 455 zr1(2) = ABS( ( gdepw_0(ji+1,jj,jk )-gdepw_0(ji ,jj,jk ) & 456 & +gdepw_0(ji+1,jj,jk+1)-gdepw_0(ji ,jj,jk+1) ) & 457 & / ( gdepw_0(ji+1,jj,jk )+gdepw_0(ji ,jj,jk ) & 458 & -gdepw_0(ji+1,jj,jk+1)-gdepw_0(ji ,jj,jk+1) + rsmall ) ) * umask(ji ,jj,jk) 459 zr1(3) = ABS( ( gdepw_0(ji,jj+1,jk )-gdepw_0(ji,jj ,jk ) & 460 & +gdepw_0(ji,jj+1,jk+1)-gdepw_0(ji,jj ,jk+1) ) & 461 & / ( gdepw_0(ji,jj+1,jk )+gdepw_0(ji,jj ,jk ) & 462 & -gdepw_0(ji,jj+1,jk+1)-gdepw_0(ji,jj ,jk+1) + rsmall ) ) * vmask(ji,jj ,jk) 463 zr1(4) = ABS( ( gdepw_0(ji,jj ,jk )-gdepw_0(ji,jj-1,jk ) & 464 & +gdepw_0(ji,jj ,jk+1)-gdepw_0(ji,jj-1,jk+1) ) & 465 & / ( gdepw_0(ji,jj ,jk )+gdepw_0(ji,jj-1,jk ) & 466 & -gdepw_0(ji,jj ,jk+1)-gdepw_0(ji,jj-1,jk+1) + rsmall ) ) * vmask(ji,jj-1,jk) 467 zrxmax = MAXVAL( zr1(1:4) ) 468 zx1(ji,jj) = MAX( zx1(ji,jj) , zrxmax ) 469 END DO 470 END DO 471 END DO 472 CALL lbc_lnk( zx1, 'T', 1. ) 473 ! 474 IF( PRESENT( px1 ) ) px1 = zx1 475 ! 476 zrxmax = MAXVAL( zx1 ) 477 ! 478 IF( lk_mpp ) CALL mpp_max( zrxmax ) ! max over the global domain 479 ! 480 IF(lwp) THEN 481 WRITE(numout,*) 482 WRITE(numout,*) 'dom_stiff : maximum grid stiffness ratio: ', zrxmax 483 WRITE(numout,*) '~~~~~~~~~' 484 ENDIF 485 ! 486 END SUBROUTINE dom_stiff 487 427 488 !!====================================================================== 428 489 END MODULE domwri -
branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
r6624 r6667 22 22 23 23 !!---------------------------------------------------------------------- 24 !! dom_zgr : defined the ocean vertical coordinate system 25 !! zgr_read : read the vertical domain coordinate and mask in domain_cfg file 26 !! zgr_bat : bathymetry fields (levels and meters) 27 !! zgr_bat_ctl : check the bathymetry files 28 !! zgr_bot_level: deepest ocean level for t-, u, and v-points 29 !! zgr_z : reference z-coordinate 30 !! zgr_zco : z-coordinate 31 !! zgr_zps : z-coordinate with partial steps 32 !! zgr_sco : s-coordinate 33 !! fssig : tanh stretch function 34 !! fssig1 : Song and Haidvogel 1994 stretch function 35 !! fgamma : Siddorn and Furner 2012 stretching function 24 !! dom_zgr : read or set the ocean vertical coordinate system 25 !! zgr_read : read the vertical domain coordinate and mask in domain_cfg file 26 !! zgr_tb_level : ocean top and bottom level for t-, u, and v-points with 1 as minimum value 36 27 !!--------------------------------------------------------------------- 37 USE oce ! ocean variables 38 USE dom_oce ! ocean domain 39 USE wet_dry ! wetting and drying 40 USE closea ! closed seas 41 USE c1d ! 1D vertical configuration 28 USE oce ! ocean variables 29 USE dom_oce ! ocean domain 30 USE usrdef_zgr ! user defined vertical coordinate system 42 31 ! 43 USE in_out_manager 44 USE iom 45 USE lbclnk 46 USE lib_mpp 47 USE wrk_nemo 48 USE timing 32 USE in_out_manager ! I/O manager 33 USE iom ! I/O library 34 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 35 USE lib_mpp ! distributed memory computing library 36 USE wrk_nemo ! Memory allocation 37 USE timing ! Timing 49 38 50 39 IMPLICIT NONE … … 52 41 53 42 PUBLIC dom_zgr ! called by dom_init.F90 54 55 ! !!* Namelist namzgr_sco *56 LOGICAL :: ln_s_sh94 ! use hybrid s-sig Song and Haidvogel 1994 stretching function fssig1 (ln_sco=T)57 LOGICAL :: ln_s_sf12 ! use hybrid s-z-sig Siddorn and Furner 2012 stretching function fgamma (ln_sco=T)58 !59 REAL(wp) :: rn_sbot_min ! minimum depth of s-bottom surface (>0) (m)60 REAL(wp) :: rn_sbot_max ! maximum depth of s-bottom surface (= ocean depth) (>0) (m)61 REAL(wp) :: rn_rmax ! maximum cut-off r-value allowed (0<rn_rmax<1)62 REAL(wp) :: rn_hc ! Critical depth for transition from sigma to stretched coordinates63 ! Song and Haidvogel 1994 stretching parameters64 REAL(wp) :: rn_theta ! surface control parameter (0<=rn_theta<=20)65 REAL(wp) :: rn_thetb ! bottom control parameter (0<=rn_thetb<= 1)66 REAL(wp) :: rn_bb ! stretching parameter67 ! ! ( rn_bb=0; top only, rn_bb =1; top and bottom)68 ! Siddorn and Furner stretching parameters69 LOGICAL :: ln_sigcrit ! use sigma coordinates below critical depth (T) or Z coordinates (F) for Siddorn & Furner stretch70 REAL(wp) :: rn_alpha ! control parameter ( > 1 stretch towards surface, < 1 towards seabed)71 REAL(wp) :: rn_efold ! efold length scale for transition to stretched coord72 REAL(wp) :: rn_zs ! depth of surface grid box73 ! bottom cell depth (Zb) is a linear function of water depth Zb = H*a + b74 REAL(wp) :: rn_zb_a ! bathymetry scaling factor for calculating Zb75 REAL(wp) :: rn_zb_b ! offset for calculating Zb76 43 77 44 !! * Substitutions … … 84 51 CONTAINS 85 52 86 SUBROUTINE dom_zgr 53 SUBROUTINE dom_zgr( k_top, k_bot ) 87 54 !!---------------------------------------------------------------------- 88 55 !! *** ROUTINE dom_zgr *** … … 101 68 !! ** Action : define gdep., e3., mbathy and bathy 102 69 !!---------------------------------------------------------------------- 103 INTEGER :: ioptio, ibat ! local integer 104 INTEGER :: ios 105 ! 106 NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco, ln_isfcav, ln_linssh 70 INTEGER, DIMENSION(:,:), INTENT(out) :: k_top, k_bot ! ocean first and last level indices 71 ! 72 INTEGER :: jk ! dummy loop index 73 INTEGER :: ioptio, ibat, ios ! local integer 74 REAL(wp) :: zrefdep ! depth of the reference level (~10m) 75 !! 76 ! NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco, ln_isfcav, ln_linssh 107 77 !!---------------------------------------------------------------------- 108 78 ! 109 79 IF( nn_timing == 1 ) CALL timing_start('dom_zgr') 110 80 ! 111 REWIND( numnam_ref ) ! Namelist namzgr in reference namelist : Vertical coordinate112 READ ( numnam_ref, namzgr, IOSTAT = ios, ERR = 901 )113 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr in reference namelist', lwp )114 115 REWIND( numnam_cfg ) ! Namelist namzgr in configuration namelist : Vertical coordinate116 READ ( numnam_cfg, namzgr, IOSTAT = ios, ERR = 902 )117 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr in configuration namelist', lwp )118 IF(lwm) WRITE ( numond, namzgr )81 ! REWIND( numnam_ref ) ! Namelist namzgr in reference namelist : Vertical coordinate 82 ! READ ( numnam_ref, namzgr, IOSTAT = ios, ERR = 901 ) 83 !901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr in reference namelist', lwp ) 84 ! 85 ! REWIND( numnam_cfg ) ! Namelist namzgr in configuration namelist : Vertical coordinate 86 ! READ ( numnam_cfg, namzgr, IOSTAT = ios, ERR = 902 ) 87 !902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr in configuration namelist', lwp ) 88 ! IF(lwm) WRITE ( numond, namzgr ) 119 89 120 90 IF(lwp) THEN ! Control print … … 122 92 WRITE(numout,*) 'dom_zgr : vertical coordinate' 123 93 WRITE(numout,*) '~~~~~~~' 124 WRITE(numout,*) ' Namelist namzgr : set vertical coordinate' 94 ENDIF 95 96 IF( ln_linssh .AND. lwp) WRITE(numout,*) ' linear free surface: the vertical mesh does not change in time' 97 98 99 IF( ln_read_cfg ) THEN !== read in mesh_mask.nc file ==! 100 IF(lwp) WRITE(numout,*) 101 IF(lwp) WRITE(numout,*) ' Read vertical mesh in "domain_cfg" file' 102 ! 103 CALL zgr_read ( ln_zco , ln_zps , ln_sco, ln_isfcav, & 104 & gdept_1d, gdepw_1d, e3t_1d, e3w_1d , & ! 1D gridpoints depth 105 & gdept_0 , gdepw_0 , & ! gridpoints depth 106 & e3t_0 , e3u_0 , e3v_0 , e3f_0 , & ! vertical scale factors 107 & e3w_0 , e3uw_0 , e3vw_0 , & ! vertical scale factors 108 & k_top , k_bot ) ! 1st & last ocean level 109 ! 110 ELSE !== User defined configuration ==! 111 IF(lwp) WRITE(numout,*) 112 IF(lwp) WRITE(numout,*) ' User defined horizontal mesh (usr_def_hgr)' 113 ! 114 CALL usr_def_zgr( ln_zco , ln_zps , ln_sco, ln_isfcav, & 115 & gdept_1d, gdepw_1d, e3t_1d, e3w_1d , & ! 1D gridpoints depth 116 & gdept_0 , gdepw_0 , & ! gridpoints depth 117 & e3t_0 , e3u_0 , e3v_0 , e3f_0 , & ! vertical scale factors 118 & e3w_0 , e3uw_0 , e3vw_0 , & ! vertical scale factors 119 & k_top , k_bot ) ! 1st & last ocean level 120 ! 121 ENDIF 122 ! 123 !!gm to be remove when changing the definition of e3 scale factors so that gde3w disappears 124 ! Compute gde3w_0 (vertical sum of e3w) 125 gde3w_0(:,:,1) = 0.5_wp * e3w_0(:,:,1) 126 DO jk = 2, jpk 127 gde3w_0(:,:,jk) = gde3w_0(:,:,jk-1) + e3w_0(:,:,jk) 128 END DO 129 ! 130 IF(lwp) THEN ! Control print 131 WRITE(numout,*) 132 WRITE(numout,*) ' Read in domain_cfg.nc or user defined type of vertical coordinate:' 125 133 WRITE(numout,*) ' z-coordinate - full steps ln_zco = ', ln_zco 126 134 WRITE(numout,*) ' z-coordinate - partial steps ln_zps = ', ln_zps 127 135 WRITE(numout,*) ' s- or hybrid z-s-coordinate ln_sco = ', ln_sco 128 136 WRITE(numout,*) ' ice shelf cavities ln_isfcav = ', ln_isfcav 129 WRITE(numout,*) ' linear free surface ln_linssh = ', ln_linssh 130 ENDIF 131 132 IF( ln_linssh .AND. lwp) WRITE(numout,*) ' linear free surface: the vertical mesh does not change in time' 137 ENDIF 133 138 134 139 ioptio = 0 ! Check Vertical coordinate options … … 137 142 IF( ln_sco ) ioptio = ioptio + 1 138 143 IF( ioptio /= 1 ) CALL ctl_stop( ' none or several vertical coordinate options used' ) 139 ! 140 ! Build the vertical coordinate system 141 ! ------------------------------------ 142 CALL zgr_z ! Reference z-coordinate system (always called) 143 CALL zgr_bat ! Bathymetry fields (levels and meters) 144 IF( lk_c1d ) CALL lbc_lnk( bathy , 'T', 1._wp ) ! 1D config.: same bathy value over the 3x3 domain 145 IF( ln_zco ) CALL zgr_zco ! z-coordinate 146 IF( ln_zps ) CALL zgr_zps ! Partial step z-coordinate 147 IF( ln_sco ) CALL zgr_sco ! s-coordinate or hybrid z-s coordinate 148 ! 149 ! final adjustment of mbathy & check 150 ! ----------------------------------- 151 IF( .NOT.lk_c1d ) CALL zgr_bat_ctl ! check bathymetry (mbathy) and suppress isolated ocean points 152 CALL zgr_bot_level ! deepest ocean level for t-, u- and v-points 153 CALL zgr_top_level ! shallowest ocean level for T-, U-, V- points 154 ! 155 IF( lk_c1d ) THEN ! 1D config.: same mbathy value over the 3x3 domain 156 ibat = mbathy(2,2) 157 mbathy(:,:) = ibat 158 END IF 144 145 146 ! ! top/bottom ocean level indices for t-, u- and v-points (f-point also for top) 147 CALL zgr_tb_level( k_top, k_bot ) ! with a minimum value set to 1 148 149 150 ! ! deepest/shallowest W level Above/Below ~10m 151 !!gm BUG in s-coordinate this does not work! 152 zrefdep = 10._wp - 0.1_wp * MINVAL( e3w_1d ) ! ref. depth with tolerance (10% of minimum layer thickness) 153 nlb10 = MINLOC( gdepw_1d, mask = gdepw_1d > zrefdep, dim = 1 ) ! shallowest W level Below ~10m 154 nla10 = nlb10 - 1 ! deepest W level Above ~10m 155 !!gm end bug 156 ! 159 157 ! 160 158 IF( nprint == 1 .AND. lwp ) THEN 161 WRITE(numout,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) ) 159 WRITE(numout,*) ' MIN val k_top ', MINVAL( k_top(:,:) ), ' MAX ', MAXVAL( k_top(:,:) ) 160 WRITE(numout,*) ' MIN val k_bot ', MINVAL( k_bot(:,:) ), ' MAX ', MAXVAL( k_bot(:,:) ) 162 161 WRITE(numout,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ), & 163 162 & ' w ', MINVAL( gdepw_0(:,:,:) ), '3w ', MINVAL( gde3w_0(:,:,:) ) … … 180 179 181 180 182 SUBROUTINE zgr_read( pdept_1d, pdepw_1d, pe3t_1d , pe3w_1d , & 183 & pdept , pdepw , & ! gridpoints depth (required) 184 & pe3t , pe3u , pe3v , pe3w , & ! scale factors (required) 185 & ptmask , pumask , pvmask , pfmask, & 186 & pbathy , prisfdep, & 187 & kbathy, kbkt, kisfdep ) 188 ! & pe3f , pe3uw , pe3vw ) ! u- & v-surfaces (if gridsize reduction in some straits) 181 SUBROUTINE zgr_read( ld_zco , ld_zps , ld_sco , ld_isfcav, & ! type of vertical coordinate 182 & pdept_1d, pdepw_1d, pe3t_1d , pe3w_1d , & ! 1D reference vertical coordinate 183 & pdept , pdepw , & ! 3D t & w-points depth 184 & pe3t , pe3u , pe3v , pe3f , & ! vertical scale factors 185 & pe3w , pe3uw , pe3vw , & ! - - - 186 & k_top , k_bot ) ! top & bottom ocean level 189 187 !!--------------------------------------------------------------------- 190 188 !! *** ROUTINE zgr_read *** … … 193 191 !! 194 192 !!---------------------------------------------------------------------- 195 REAL(wp), DIMENSION(:,:,:), INTENT(out) :: pdept_1d, pdepw_1d ! 1D grid-point depth [m]196 REAL(wp), DIMENSION(:,:,:), INTENT(out) :: pe3t_1d , pe3w_1d ! 1D grid-point depth [m]197 REAL(wp), DIMENSION(: ,:,:), INTENT(out) :: pdept, pdepw !, pde3w ! grid-point depth[m]198 REAL(wp), DIMENSION(: ,:,:), INTENT(out) :: pe3t , pe3u , pe3v , pe3w ! vertical scale factors[m]199 ! REAL(wp), DIMENSION(:,:,:), INTENT(out) :: pe3f , pe3uw, pe3vw ! i-scale factors 200 REAL(wp), DIMENSION(:,:,:), INTENT(out) :: p tmask , pumask , pvmask , pfmask ! masks [-]201 REAL(wp), DIMENSION(:,: ) , INTENT(out) :: pbathy , prisfdep ! bathymetry, iceshelf depth [m]202 INTEGER , DIMENSION(:,:) , INTENT(out) :: k bathy, kbkt, kisfdep ! bathymetry, iceshelf depth [m]193 LOGICAL , INTENT(out) :: ld_zco, ld_zps, ld_sco ! vertical coordinate flags 194 LOGICAL , INTENT(out) :: ld_isfcav ! under iceshelf cavity flag 195 REAL(wp), DIMENSION(:) , INTENT(out) :: pdept_1d, pdepw_1d ! 1D grid-point depth [m] 196 REAL(wp), DIMENSION(:) , INTENT(out) :: pe3t_1d , pe3w_1d ! 1D vertical scale factors [m] 197 REAL(wp), DIMENSION(:,:,:), INTENT(out) :: pdept, pdepw ! grid-point depth [m] 198 REAL(wp), DIMENSION(:,:,:), INTENT(out) :: pe3t , pe3u , pe3v , pe3f ! vertical scale factors [m] 199 REAL(wp), DIMENSION(:,:,:), INTENT(out) :: pe3w , pe3uw, pe3vw ! - - - 200 INTEGER , DIMENSION(:,:) , INTENT(out) :: k_top , k_bot ! first & last ocean level 203 201 ! 204 202 INTEGER :: inum ! local logical unit 203 REAL(WP) :: z_zco, z_zps, z_sco, z_cav 205 204 REAL(wp), DIMENSION(jpi,jpj) :: z2d ! 2D workspace 206 205 !!---------------------------------------------------------------------- … … 208 207 IF(lwp) THEN 209 208 WRITE(numout,*) 210 WRITE(numout,*) 'hgr_read : read the vertical coordinates in mesh_mask' 211 WRITE(numout,*) '~~~~~~~~ jpiglo = ', jpiglo, ' jpjglo = ', jpjglo, ' jpk = ', jpk 212 ENDIF 213 ! 214 CALL iom_open( 'mesh_mask', inum ) 215 ! 216 CALL iom_get( inum, jpdom_unknown, 'gdept_1d', pdept_1d ) ! reference 1D-coordinate 209 WRITE(numout,*) 'hgr_read : read the vertical coordinates in "domain_cfg.nc" file' 210 WRITE(numout,*) '~~~~~~~~ jpiglo = ', jpiglo, ' jpjglo = ', jpjglo, ' jpkglo = ', jpkglo 211 ENDIF 212 ! 213 CALL iom_open( 'domain_cfg', inum ) 214 ! 215 ! ! type of vertical coordinate 216 CALL iom_get( inum, 'ln_zco' , z_zco ) 217 CALL iom_get( inum, 'ln_zps' , z_zps ) 218 CALL iom_get( inum, 'ln_sco' , z_sco ) 219 IF( z_zco == 0._wp ) THEN ; ld_zco = .false. ; ELSE ; ld_zco = .true. ; ENDIF 220 IF( z_zps == 0._wp ) THEN ; ld_zps = .false. ; ELSE ; ld_zps = .true. ; ENDIF 221 IF( z_sco == 0._wp ) THEN ; ld_sco = .false. ; ELSE ; ld_sco = .true. ; ENDIF 222 ! 223 ! ! ocean cavities under iceshelves 224 CALL iom_get( inum, 'ln_isfcav', z_cav ) 225 IF( z_cav == 0._wp ) THEN ; ld_isfcav = .false. ; ELSE ; ld_isfcav = .true. ; ENDIF 226 ! 227 ! ! reference 1D-coordinate 228 CALL iom_get( inum, jpdom_unknown, 'gdept_1d', pdept_1d ) 217 229 CALL iom_get( inum, jpdom_unknown, 'gdepw_1d', pdepw_1d ) 218 230 CALL iom_get( inum, jpdom_unknown, 'e3t_1d' , pe3t_1d ) 219 231 CALL iom_get( inum, jpdom_unknown, 'e3w_1d' , pe3w_1d ) 220 232 ! 221 CALL iom_get( inum, jpdom_data, 'gdept_0', pdept, lrowattr=ln_use_jattr ) ! depth 222 CALL iom_get( inum, jpdom_data, 'gdepw_0', pdepw, lrowattr=ln_use_jattr ) 223 ! CALL iom_get( inum, jpdom_data, 'gde3w_0', pde3w, lrowattr=ln_use_jattr ) 224 ! 225 CALL iom_get( inum, jpdom_data, 'e3t_0' , pe3t , lrowattr=ln_use_jattr ) ! vertical scale factors 226 CALL iom_get( inum, jpdom_data, 'e3u_0' , pe3u , lrowattr=ln_use_jattr ) 227 CALL iom_get( inum, jpdom_data, 'e3v_0' , pe3v , lrowattr=ln_use_jattr ) 228 ! CALL iom_get( inum, jpdom_data, 'e3f_0' , pe3f , lrowattr=ln_use_jattr ) 229 CALL iom_get( inum, jpdom_data, 'e3w_0' , pe3w , lrowattr=ln_use_jattr ) 230 ! CALL iom_get( inum, jpdom_data, 'e3uw_0' , pe3uw, lrowattr=ln_use_jattr ) 231 ! CALL iom_get( inum, jpdom_data, 'e3vw_0' , pe3vw, lrowattr=ln_use_jattr ) 232 233 234 235 236 237 CALL iom_get( inum, jpdom_data, 'tmask' , ptmask ) ! masks 238 CALL iom_get( inum, jpdom_data, 'umask' , pumask ) 239 CALL iom_get( inum, jpdom_data, 'vmask' , pvmask ) 240 CALL iom_get( inum, jpdom_data, 'fmask' , pfmask ) 241 242 243 !!gm Probably not required fields : 244 CALL iom_get( inum, jpdom_data, 'bathy' , bathy ) ! depth of the ocean at T-points 245 CALL iom_get( inum, jpdom_data, 'mbathy' , z2d ) ! nb of ocean T-points 246 kbathy(:,:) = INT( z2d(:,:) ) 247 CALL iom_get( inum, jpdom_data, 'mbkt' , z2d ) ! nb of ocean T-points 248 kbkt(:,:) = INT( z2d(:,:) ) 249 250 ! 251 CALL iom_get( inum, jpdom_data, 'bathy' , risfdep ) ! depth of the iceshelves at T-points 252 CALL iom_get( inum, jpdom_data, 'misfdep', z2d ) ! nb of ocean T-points (ISF) 253 kisfdep(:,:) = INT( z2d(:,:) ) 254 255 233 ! ! depth 234 CALL iom_get( inum, jpdom_data, 'gdept_0', pdept ) 235 CALL iom_get( inum, jpdom_data, 'gdepw_0', pdepw ) 236 ! CALL iom_get( inum, jpdom_data, 'gde3w_0', pde3w ) 237 ! 238 ! ! vertical scale factors 239 CALL iom_get( inum, jpdom_data, 'e3t_0' , pe3t ) 240 CALL iom_get( inum, jpdom_data, 'e3u_0' , pe3u ) 241 CALL iom_get( inum, jpdom_data, 'e3v_0' , pe3v ) 242 CALL iom_get( inum, jpdom_data, 'e3f_0' , pe3f ) 243 CALL iom_get( inum, jpdom_data, 'e3w_0' , pe3w ) 244 CALL iom_get( inum, jpdom_data, 'e3uw_0' , pe3uw ) 245 CALL iom_get( inum, jpdom_data, 'e3vw_0' , pe3vw ) 246 247 ! ! ocean top and bottom level 248 CALL iom_get( inum, jpdom_data, 'bottom level' , z2d ) ! nb of ocean T-points 249 k_bot(:,:) = INT( z2d(:,:) ) 250 CALL iom_get( inum, jpdom_data, 'top level' , z2d ) ! nb of ocean T-points (ISF) 251 k_top(:,:) = INT( z2d(:,:) ) 256 252 ! 257 253 CALL iom_close( inum ) … … 260 256 261 257 262 SUBROUTINE zgr_z 263 !!---------------------------------------------------------------------- 264 !! *** ROUTINE zgr_z *** 265 !! 266 !! ** Purpose : set the depth of model levels and the resulting 267 !! vertical scale factors. 268 !! 269 !! ** Method : z-coordinate system (use in all type of coordinate) 270 !! The depth of model levels is defined from an analytical 271 !! function the derivative of which gives the scale factors. 272 !! both depth and scale factors only depend on k (1d arrays). 273 !! w-level: gdepw_1d = gdep(k) 274 !! e3w_1d(k) = dk(gdep)(k) = e3(k) 275 !! t-level: gdept_1d = gdep(k+0.5) 276 !! e3t_1d(k) = dk(gdep)(k+0.5) = e3(k+0.5) 277 !! 278 !! ** Action : - gdept_1d, gdepw_1d : depth of T- and W-point (m) 279 !! - e3t_1d , e3w_1d : scale factors at T- and W-levels (m) 280 !! 281 !! Reference : Marti, Madec & Delecluse, 1992, JGR, 97, No8, 12,763-12,766. 282 !!---------------------------------------------------------------------- 283 INTEGER :: jk ! dummy loop indices 284 REAL(wp) :: zt, zw ! temporary scalars 285 REAL(wp) :: zsur, za0, za1, zkth ! Values set from parameters in 286 REAL(wp) :: zacr, zdzmin, zhmax ! par_CONFIG_Rxx.h90 287 REAL(wp) :: zrefdep ! depth of the reference level (~10m) 288 REAL(wp) :: za2, zkth2, zacr2 ! Values for optional double tanh function set from parameters 289 !!---------------------------------------------------------------------- 290 ! 291 IF( nn_timing == 1 ) CALL timing_start('zgr_z') 292 ! 293 ! Set variables from parameters 294 ! ------------------------------ 295 zkth = ppkth ; zacr = ppacr 296 zdzmin = ppdzmin ; zhmax = pphmax 297 zkth2 = ppkth2 ; zacr2 = ppacr2 ! optional (ldbletanh=T) double tanh parameters 298 299 ! If ppa1 and ppa0 and ppsur are et to pp_to_be_computed 300 ! za0, za1, zsur are computed from ppdzmin , pphmax, ppkth, ppacr 301 IF( ppa1 == pp_to_be_computed .AND. & 302 & ppa0 == pp_to_be_computed .AND. & 303 & ppsur == pp_to_be_computed ) THEN 304 ! 305 #if defined key_agrif 306 za1 = ( ppdzmin - pphmax / REAL( jpkdta-1 , wp ) ) & 307 & / ( TANH((1-ppkth)/ppacr) - ppacr/REAL( jpkdta-1 , wp ) * ( LOG( COSH( (jpkdta - ppkth) / ppacr) )& 308 & - LOG( COSH( ( 1 - ppkth) / ppacr) ) ) ) 309 #else 310 za1 = ( ppdzmin - pphmax / FLOAT(jpkm1) ) & 311 & / ( TANH((1-ppkth)/ppacr) - ppacr/FLOAT(jpk-1) * ( LOG( COSH( (jpk - ppkth) / ppacr) ) & 312 & - LOG( COSH( ( 1 - ppkth) / ppacr) ) ) ) 313 #endif 314 za0 = ppdzmin - za1 * TANH( (1-ppkth) / ppacr ) 315 zsur = - za0 - za1 * ppacr * LOG( COSH( (1-ppkth) / ppacr ) ) 316 ELSE 317 za1 = ppa1 ; za0 = ppa0 ; zsur = ppsur 318 za2 = ppa2 ! optional (ldbletanh=T) double tanh parameter 319 ENDIF 320 321 IF(lwp) THEN ! Parameter print 322 WRITE(numout,*) 323 WRITE(numout,*) ' zgr_z : Reference vertical z-coordinates' 324 WRITE(numout,*) ' ~~~~~~~' 325 IF( ppkth == 0._wp ) THEN 326 WRITE(numout,*) ' Uniform grid with ',jpk-1,' layers' 327 WRITE(numout,*) ' Total depth :', zhmax 328 #if defined key_agrif 329 WRITE(numout,*) ' Layer thickness:', zhmax/(jpkdta-1) 330 #else 331 WRITE(numout,*) ' Layer thickness:', zhmax/(jpk-1) 332 #endif 333 ELSE 334 IF( ppa1 == 0._wp .AND. ppa0 == 0._wp .AND. ppsur == 0._wp ) THEN 335 WRITE(numout,*) ' zsur, za0, za1 computed from ' 336 WRITE(numout,*) ' zdzmin = ', zdzmin 337 WRITE(numout,*) ' zhmax = ', zhmax 338 ENDIF 339 WRITE(numout,*) ' Value of coefficients for vertical mesh:' 340 WRITE(numout,*) ' zsur = ', zsur 341 WRITE(numout,*) ' za0 = ', za0 342 WRITE(numout,*) ' za1 = ', za1 343 WRITE(numout,*) ' zkth = ', zkth 344 WRITE(numout,*) ' zacr = ', zacr 345 IF( ldbletanh ) THEN 346 WRITE(numout,*) ' (Double tanh za2 = ', za2 347 WRITE(numout,*) ' parameters) zkth2= ', zkth2 348 WRITE(numout,*) ' zacr2= ', zacr2 349 ENDIF 350 ENDIF 351 ENDIF 352 353 354 ! Reference z-coordinate (depth - scale factor at T- and W-points) 355 ! ====================== 356 IF( ppkth == 0._wp ) THEN ! uniform vertical grid 357 #if defined key_agrif 358 za1 = zhmax / FLOAT(jpkdta-1) 359 #else 360 za1 = zhmax / FLOAT(jpk-1) 361 #endif 362 DO jk = 1, jpk 363 zw = FLOAT( jk ) 364 zt = FLOAT( jk ) + 0.5_wp 365 gdepw_1d(jk) = ( zw - 1 ) * za1 366 gdept_1d(jk) = ( zt - 1 ) * za1 367 e3w_1d (jk) = za1 368 e3t_1d (jk) = za1 369 END DO 370 ELSE ! Madec & Imbard 1996 function 371 IF( .NOT. ldbletanh ) THEN 372 DO jk = 1, jpk 373 zw = REAL( jk , wp ) 374 zt = REAL( jk , wp ) + 0.5_wp 375 gdepw_1d(jk) = ( zsur + za0 * zw + za1 * zacr * LOG ( COSH( (zw-zkth) / zacr ) ) ) 376 gdept_1d(jk) = ( zsur + za0 * zt + za1 * zacr * LOG ( COSH( (zt-zkth) / zacr ) ) ) 377 e3w_1d (jk) = za0 + za1 * TANH( (zw-zkth) / zacr ) 378 e3t_1d (jk) = za0 + za1 * TANH( (zt-zkth) / zacr ) 379 END DO 380 ELSE 381 DO jk = 1, jpk 382 zw = FLOAT( jk ) 383 zt = FLOAT( jk ) + 0.5_wp 384 ! Double tanh function 385 gdepw_1d(jk) = ( zsur + za0 * zw + za1 * zacr * LOG ( COSH( (zw-zkth ) / zacr ) ) & 386 & + za2 * zacr2* LOG ( COSH( (zw-zkth2) / zacr2 ) ) ) 387 gdept_1d(jk) = ( zsur + za0 * zt + za1 * zacr * LOG ( COSH( (zt-zkth ) / zacr ) ) & 388 & + za2 * zacr2* LOG ( COSH( (zt-zkth2) / zacr2 ) ) ) 389 e3w_1d (jk) = za0 + za1 * TANH( (zw-zkth ) / zacr ) & 390 & + za2 * TANH( (zw-zkth2) / zacr2 ) 391 e3t_1d (jk) = za0 + za1 * TANH( (zt-zkth ) / zacr ) & 392 & + za2 * TANH( (zt-zkth2) / zacr2 ) 393 END DO 394 ENDIF 395 gdepw_1d(1) = 0._wp ! force first w-level to be exactly at zero 396 ENDIF 397 398 IF ( ln_isfcav ) THEN 399 ! need to be like this to compute the pressure gradient with ISF. If not, level beneath the ISF are not aligned (sum(e3t) /= depth) 400 ! define e3t_0 and e3w_0 as the differences between gdept and gdepw respectively 401 DO jk = 1, jpkm1 402 e3t_1d(jk) = gdepw_1d(jk+1)-gdepw_1d(jk) 403 END DO 404 e3t_1d(jpk) = e3t_1d(jpk-1) ! we don't care because this level is masked in NEMO 405 406 DO jk = 2, jpk 407 e3w_1d(jk) = gdept_1d(jk) - gdept_1d(jk-1) 408 END DO 409 e3w_1d(1 ) = 2._wp * (gdept_1d(1) - gdepw_1d(1)) 410 END IF 411 412 !!gm BUG in s-coordinate this does not work! 413 ! deepest/shallowest W level Above/Below ~10m 414 zrefdep = 10._wp - 0.1_wp * MINVAL( e3w_1d ) ! ref. depth with tolerance (10% of minimum layer thickness) 415 nlb10 = MINLOC( gdepw_1d, mask = gdepw_1d > zrefdep, dim = 1 ) ! shallowest W level Below ~10m 416 nla10 = nlb10 - 1 ! deepest W level Above ~10m 417 !!gm end bug 418 419 IF(lwp) THEN ! control print 420 WRITE(numout,*) 421 WRITE(numout,*) ' Reference z-coordinate depth and scale factors:' 422 WRITE(numout, "(9x,' level gdept_1d gdepw_1d e3t_1d e3w_1d ')" ) 423 WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, gdept_1d(jk), gdepw_1d(jk), e3t_1d(jk), e3w_1d(jk), jk = 1, jpk ) 424 ENDIF 425 DO jk = 1, jpk ! control positivity 426 IF( e3w_1d (jk) <= 0._wp .OR. e3t_1d (jk) <= 0._wp ) CALL ctl_stop( 'dom:zgr_z: e3w_1d or e3t_1d =< 0 ' ) 427 IF( gdepw_1d(jk) < 0._wp .OR. gdept_1d(jk) < 0._wp ) CALL ctl_stop( 'dom:zgr_z: gdepw_1d or gdept_1d < 0 ' ) 428 END DO 429 ! 430 IF( nn_timing == 1 ) CALL timing_stop('zgr_z') 431 ! 432 END SUBROUTINE zgr_z 433 434 435 SUBROUTINE zgr_bat 436 !!---------------------------------------------------------------------- 437 !! *** ROUTINE zgr_bat *** 438 !! 439 !! ** Purpose : set bathymetry both in levels and meters 440 !! 441 !! ** Method : read or define mbathy and bathy arrays 442 !! * level bathymetry: 443 !! The ocean basin geometry is given by a two-dimensional array, 444 !! mbathy, which is defined as follow : 445 !! mbathy(ji,jj) = 1, ..., jpk-1, the number of ocean level 446 !! at t-point (ji,jj). 447 !! = 0 over the continental t-point. 448 !! The array mbathy is checked to verified its consistency with 449 !! model option. in particular: 450 !! mbathy must have at least 1 land grid-points (mbathy<=0) 451 !! along closed boundary. 452 !! mbathy must be cyclic IF jperio=1. 453 !! mbathy must be lower or equal to jpk-1. 454 !! isolated ocean grid points are suppressed from mbathy 455 !! since they are only connected to remaining 456 !! ocean through vertical diffusion. 457 !! ntopo=-1 : rectangular channel or bassin with a bump 458 !! ntopo= 0 : flat rectangular channel or basin 459 !! ntopo= 1 : mbathy is read in 'bathy_level.nc' NetCDF file 460 !! bathy is read in 'bathy_meter.nc' NetCDF file 461 !! 462 !! ** Action : - mbathy: level bathymetry (in level index) 463 !! - bathy : meter bathymetry (in meters) 464 !!---------------------------------------------------------------------- 465 INTEGER :: ji, jj, jk ! dummy loop indices 466 INTEGER :: inum ! temporary logical unit 467 INTEGER :: ierror ! error flag 468 INTEGER :: ii_bump, ij_bump, ih ! bump center position 469 INTEGER :: ii0, ii1, ij0, ij1, ik ! local indices 470 REAL(wp) :: r_bump , h_bump , h_oce ! bump characteristics 471 REAL(wp) :: zi, zj, zh, zhmin ! local scalars 472 INTEGER , ALLOCATABLE, DIMENSION(:,:) :: idta ! global domain integer data 473 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdta ! global domain scalar data 474 !!---------------------------------------------------------------------- 475 ! 476 IF( nn_timing == 1 ) CALL timing_start('zgr_bat') 477 ! 478 IF(lwp) WRITE(numout,*) 479 IF(lwp) WRITE(numout,*) ' zgr_bat : defines level and meter bathymetry' 480 IF(lwp) WRITE(numout,*) ' ~~~~~~~' 481 ! ! ================== ! 482 IF( ntopo == 0 .OR. ntopo == -1 ) THEN ! defined by hand ! 483 ! ! ================== ! 484 ! ! global domain level and meter bathymetry (idta,zdta) 485 ! 486 ALLOCATE( idta(jpidta,jpjdta), STAT=ierror ) 487 IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'zgr_bat: unable to allocate idta array' ) 488 ALLOCATE( zdta(jpidta,jpjdta), STAT=ierror ) 489 IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'zgr_bat: unable to allocate zdta array' ) 490 ! 491 IF( ntopo == 0 ) THEN ! flat basin 492 IF(lwp) WRITE(numout,*) 493 IF(lwp) WRITE(numout,*) ' bathymetry field: flat basin' 494 IF( rn_bathy > 0.01 ) THEN 495 IF(lwp) WRITE(numout,*) ' Depth = rn_bathy read in namelist' 496 zdta(:,:) = rn_bathy 497 IF( ln_sco ) THEN ! s-coordinate (zsc ): idta()=jpk 498 idta(:,:) = jpkm1 499 ELSE ! z-coordinate (zco or zps): step-like topography 500 idta(:,:) = jpkm1 501 DO jk = 1, jpkm1 502 WHERE( gdept_1d(jk) < zdta(:,:) .AND. zdta(:,:) <= gdept_1d(jk+1) ) idta(:,:) = jk 503 END DO 504 ENDIF 505 ELSE 506 IF(lwp) WRITE(numout,*) ' Depth = depthw(jpkm1)' 507 idta(:,:) = jpkm1 ! before last level 508 zdta(:,:) = gdepw_1d(jpk) ! last w-point depth 509 h_oce = gdepw_1d(jpk) 510 ENDIF 511 ENDIF 512 ! ! set GLOBAL boundary conditions 513 ! ! Caution : idta on the global domain: use of jperio, not nperio 514 IF( jperio == 1 .OR. jperio == 4 .OR. jperio == 6 ) THEN 515 idta( : , 1 ) = -1 ; zdta( : , 1 ) = -1._wp 516 idta( : ,jpjdta) = 0 ; zdta( : ,jpjdta) = 0._wp 517 ELSEIF( jperio == 2 ) THEN 518 idta( : , 1 ) = idta( : , 3 ) ; zdta( : , 1 ) = zdta( : , 3 ) 519 idta( : ,jpjdta) = 0 ; zdta( : ,jpjdta) = 0._wp 520 idta( 1 , : ) = 0 ; zdta( 1 , : ) = 0._wp 521 idta(jpidta, : ) = 0 ; zdta(jpidta, : ) = 0._wp 522 ELSE 523 ih = 0 ; zh = 0._wp 524 IF( ln_sco ) ih = jpkm1 ; IF( ln_sco ) zh = h_oce 525 idta( : , 1 ) = ih ; zdta( : , 1 ) = zh 526 idta( : ,jpjdta) = ih ; zdta( : ,jpjdta) = zh 527 idta( 1 , : ) = ih ; zdta( 1 , : ) = zh 528 idta(jpidta, : ) = ih ; zdta(jpidta, : ) = zh 529 ENDIF 530 531 ! ! local domain level and meter bathymetries (mbathy,bathy) 532 mbathy(:,:) = 0 ! set to zero extra halo points 533 bathy (:,:) = 0._wp ! (require for mpp case) 534 DO jj = 1, nlcj ! interior values 535 DO ji = 1, nlci 536 mbathy(ji,jj) = idta( mig(ji), mjg(jj) ) 537 bathy (ji,jj) = zdta( mig(ji), mjg(jj) ) 538 END DO 539 END DO 540 risfdep(:,:) = 0._wp 541 misfdep(:,:) = 1 542 ! 543 DEALLOCATE( idta, zdta ) 544 ! 545 ! ! ================ ! 546 ELSEIF( ntopo == 1 ) THEN ! read in file ! (over the local domain) 547 ! ! ================ ! 548 ! 549 IF( ln_zco ) THEN ! zco : read level bathymetry 550 CALL iom_open ( 'bathy_level.nc', inum ) 551 CALL iom_get ( inum, jpdom_data, 'Bathy_level', bathy ) 552 CALL iom_close( inum ) 553 mbathy(:,:) = INT( bathy(:,:) ) 554 ! ! ===================== 555 IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! ORCA R2 configuration 556 ! ! ===================== 557 ! 558 ii0 = 140 ; ii1 = 140 ! Gibraltar Strait open 559 ij0 = 102 ; ij1 = 102 ! (Thomson, Ocean Modelling, 1995) 560 DO ji = mi0(ii0), mi1(ii1) 561 DO jj = mj0(ij0), mj1(ij1) 562 mbathy(ji,jj) = 15 563 END DO 564 END DO 565 IF(lwp) WRITE(numout,*) 566 IF(lwp) WRITE(numout,*) ' orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0 567 ! 568 ii0 = 160 ; ii1 = 160 ! Bab el mandeb Strait open 569 ij0 = 88 ; ij1 = 88 ! (Thomson, Ocean Modelling, 1995) 570 DO ji = mi0(ii0), mi1(ii1) 571 DO jj = mj0(ij0), mj1(ij1) 572 mbathy(ji,jj) = 12 573 END DO 574 END DO 575 IF(lwp) WRITE(numout,*) 576 IF(lwp) WRITE(numout,*) ' orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0 577 ! 578 ENDIF 579 ! 580 ENDIF 581 IF( ln_zps .OR. ln_sco ) THEN ! zps or sco : read meter bathymetry 582 CALL iom_open ( 'bathy_meter.nc', inum ) 583 IF ( ln_isfcav ) THEN 584 CALL iom_get ( inum, jpdom_data, 'Bathymetry_isf', bathy, lrowattr=.false. ) 585 ELSE 586 CALL iom_get ( inum, jpdom_data, 'Bathymetry' , bathy, lrowattr=ln_use_jattr ) 587 END IF 588 CALL iom_close( inum ) 589 ! 590 risfdep(:,:) = 0._wp 591 misfdep(:,:) = 1 592 IF ( ln_isfcav ) THEN 593 CALL iom_open ( 'isf_draft_meter.nc', inum ) 594 CALL iom_get ( inum, jpdom_data, 'isf_draft', risfdep ) 595 CALL iom_close( inum ) 596 WHERE( bathy(:,:) <= 0._wp ) risfdep(:,:) = 0._wp 597 598 ! set grounded point to 0 599 ! (a treshold could be set here if needed, or set it offline based on the grounded fraction) 600 WHERE ( bathy(:,:) <= risfdep(:,:) + rn_isfhmin ) 601 misfdep(:,:) = 0 ; risfdep(:,:) = 0._wp 602 mbathy (:,:) = 0 ; bathy (:,:) = 0._wp 603 END WHERE 604 END IF 605 ! 606 IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! ORCA R2 configuration 607 ! 608 ii0 = 140 ; ii1 = 140 ! Gibraltar Strait open 609 ij0 = 102 ; ij1 = 102 ! (Thomson, Ocean Modelling, 1995) 610 DO ji = mi0(ii0), mi1(ii1) 611 DO jj = mj0(ij0), mj1(ij1) 612 bathy(ji,jj) = 284._wp 613 END DO 614 END DO 615 IF(lwp) WRITE(numout,*) 616 IF(lwp) WRITE(numout,*) ' orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0 617 ! 618 ii0 = 160 ; ii1 = 160 ! Bab el mandeb Strait open 619 ij0 = 88 ; ij1 = 88 ! (Thomson, Ocean Modelling, 1995) 620 DO ji = mi0(ii0), mi1(ii1) 621 DO jj = mj0(ij0), mj1(ij1) 622 bathy(ji,jj) = 137._wp 623 END DO 624 END DO 625 IF(lwp) WRITE(numout,*) 626 IF(lwp) WRITE(numout,*) ' orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0 627 ! 628 ENDIF 629 ! 630 ENDIF 631 ! ! =============== ! 632 ELSE ! error ! 633 ! ! =============== ! 634 WRITE(ctmp1,*) 'parameter , ntopo = ', ntopo 635 CALL ctl_stop( ' zgr_bat : '//trim(ctmp1) ) 636 ENDIF 637 ! 638 IF( nn_closea == 0 ) CALL clo_bat( bathy, mbathy ) !== NO closed seas or lakes ==! 639 ! 640 IF ( .not. ln_sco ) THEN !== set a minimum depth ==! 641 IF( rn_hmin < 0._wp ) THEN ; ik = - INT( rn_hmin ) ! from a nb of level 642 ELSE ; ik = MINLOC( gdepw_1d, mask = gdepw_1d > rn_hmin, dim = 1 ) ! from a depth 643 ENDIF 644 zhmin = gdepw_1d(ik+1) ! minimum depth = ik+1 w-levels 645 WHERE( bathy(:,:) <= 0._wp ) ; bathy(:,:) = 0._wp ! min=0 over the lands 646 ELSE WHERE ; bathy(:,:) = MAX( zhmin , bathy(:,:) ) ! min=zhmin over the oceans 647 END WHERE 648 IF(lwp) write(numout,*) 'Minimum ocean depth: ', zhmin, ' minimum number of ocean levels : ', ik 649 ENDIF 650 ! 651 IF( nn_timing == 1 ) CALL timing_stop('zgr_bat') 652 ! 653 END SUBROUTINE zgr_bat 654 655 656 657 658 SUBROUTINE zgr_bat_ctl 659 !!---------------------------------------------------------------------- 660 !! *** ROUTINE zgr_bat_ctl *** 661 !! 662 !! ** Purpose : check the bathymetry in levels 663 !! 664 !! ** Method : The array mbathy is checked to verified its consistency 665 !! with the model options. in particular: 666 !! mbathy must have at least 1 land grid-points (mbathy<=0) 667 !! along closed boundary. 668 !! mbathy must be cyclic IF jperio=1. 669 !! mbathy must be lower or equal to jpk-1. 670 !! isolated ocean grid points are suppressed from mbathy 671 !! since they are only connected to remaining 672 !! ocean through vertical diffusion. 673 !! C A U T I O N : mbathy will be modified during the initializa- 674 !! tion phase to become the number of non-zero w-levels of a water 675 !! column, with a minimum value of 1. 676 !! 677 !! ** Action : - update mbathy: level bathymetry (in level index) 678 !! - update bathy : meter bathymetry (in meters) 679 !!---------------------------------------------------------------------- 680 INTEGER :: ji, jj, jl ! dummy loop indices 681 INTEGER :: icompt, ibtest, ikmax ! temporary integers 682 REAL(wp), POINTER, DIMENSION(:,:) :: zbathy 683 !!---------------------------------------------------------------------- 684 ! 685 IF( nn_timing == 1 ) CALL timing_start('zgr_bat_ctl') 686 ! 687 CALL wrk_alloc( jpi, jpj, zbathy ) 688 ! 689 IF(lwp) WRITE(numout,*) 690 IF(lwp) WRITE(numout,*) ' zgr_bat_ctl : check the bathymetry' 691 IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~' 692 ! ! Suppress isolated ocean grid points 693 IF(lwp) WRITE(numout,*) 694 IF(lwp) WRITE(numout,*)' suppress isolated ocean grid points' 695 IF(lwp) WRITE(numout,*)' -----------------------------------' 696 icompt = 0 697 DO jl = 1, 2 698 IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN 699 mbathy( 1 ,:) = mbathy(jpim1,:) ! local domain is cyclic east-west 700 mbathy(jpi,:) = mbathy( 2 ,:) 701 ENDIF 702 DO jj = 2, jpjm1 703 DO ji = 2, jpim1 704 ibtest = MAX( mbathy(ji-1,jj), mbathy(ji+1,jj), & 705 & mbathy(ji,jj-1), mbathy(ji,jj+1) ) 706 IF( ibtest < mbathy(ji,jj) ) THEN 707 IF(lwp) WRITE(numout,*) ' the number of ocean level at ', & 708 & 'grid-point (i,j) = ',ji,jj,' is changed from ', mbathy(ji,jj),' to ', ibtest 709 mbathy(ji,jj) = ibtest 710 icompt = icompt + 1 711 ENDIF 712 END DO 713 END DO 714 END DO 715 IF( lk_mpp ) CALL mpp_sum( icompt ) 716 IF( icompt == 0 ) THEN 717 IF(lwp) WRITE(numout,*)' no isolated ocean grid points' 718 ELSE 719 IF(lwp) WRITE(numout,*)' ',icompt,' ocean grid points suppressed' 720 ENDIF 721 IF( lk_mpp ) THEN 722 zbathy(:,:) = FLOAT( mbathy(:,:) ) 723 CALL lbc_lnk( zbathy, 'T', 1._wp ) 724 mbathy(:,:) = INT( zbathy(:,:) ) 725 ENDIF 726 ! ! East-west cyclic boundary conditions 727 IF( nperio == 0 ) THEN 728 IF(lwp) WRITE(numout,*) ' mbathy set to 0 along east and west boundary: nperio = ', nperio 729 IF( lk_mpp ) THEN 730 IF( nbondi == -1 .OR. nbondi == 2 ) THEN 731 IF( jperio /= 1 ) mbathy(1,:) = 0 732 ENDIF 733 IF( nbondi == 1 .OR. nbondi == 2 ) THEN 734 IF( jperio /= 1 ) mbathy(nlci,:) = 0 735 ENDIF 736 ELSE 737 IF( ln_zco .OR. ln_zps ) THEN 738 mbathy( 1 ,:) = 0 739 mbathy(jpi,:) = 0 740 ELSE 741 mbathy( 1 ,:) = jpkm1 742 mbathy(jpi,:) = jpkm1 743 ENDIF 744 ENDIF 745 ELSEIF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN 746 IF(lwp) WRITE(numout,*)' east-west cyclic boundary conditions on mbathy: nperio = ', nperio 747 mbathy( 1 ,:) = mbathy(jpim1,:) 748 mbathy(jpi,:) = mbathy( 2 ,:) 749 ELSEIF( nperio == 2 ) THEN 750 IF(lwp) WRITE(numout,*) ' equatorial boundary conditions on mbathy: nperio = ', nperio 751 ELSE 752 IF(lwp) WRITE(numout,*) ' e r r o r' 753 IF(lwp) WRITE(numout,*) ' parameter , nperio = ', nperio 754 ! STOP 'dom_mba' 755 ENDIF 756 ! Boundary condition on mbathy 757 IF( .NOT.lk_mpp ) THEN 758 !!gm !!bug ??? think about it ! 759 ! ... mono- or macro-tasking: T-point, >0, 2D array, no slab 760 zbathy(:,:) = FLOAT( mbathy(:,:) ) 761 CALL lbc_lnk( zbathy, 'T', 1._wp ) 762 mbathy(:,:) = INT( zbathy(:,:) ) 763 ENDIF 764 ! Number of ocean level inferior or equal to jpkm1 765 ikmax = 0 766 DO jj = 1, jpj 767 DO ji = 1, jpi 768 ikmax = MAX( ikmax, mbathy(ji,jj) ) 769 END DO 770 END DO 771 !!gm !!! test to do: ikmax = MAX( mbathy(:,:) ) ??? 772 IF( ikmax > jpkm1 ) THEN 773 IF(lwp) WRITE(numout,*) ' maximum number of ocean level = ', ikmax,' > jpk-1' 774 IF(lwp) WRITE(numout,*) ' change jpk to ',ikmax+1,' to use the exact ead bathymetry' 775 ELSE IF( ikmax < jpkm1 ) THEN 776 IF(lwp) WRITE(numout,*) ' maximum number of ocean level = ', ikmax,' < jpk-1' 777 IF(lwp) WRITE(numout,*) ' you can decrease jpk to ', ikmax+1 778 ENDIF 779 ! 780 CALL wrk_dealloc( jpi, jpj, zbathy ) 781 ! 782 IF( nn_timing == 1 ) CALL timing_stop('zgr_bat_ctl') 783 ! 784 END SUBROUTINE zgr_bat_ctl 785 786 787 SUBROUTINE zgr_bot_level 788 !!---------------------------------------------------------------------- 789 !! *** ROUTINE zgr_bot_level *** 258 SUBROUTINE zgr_tb_level( k_top, k_bot ) 259 !!---------------------------------------------------------------------- 260 !! *** ROUTINE zgr_tb_level *** 790 261 !! 791 262 !! ** Purpose : defines the vertical index of ocean bottom (mbk. arrays) 792 263 !! 793 !! ** Method : computes from mbathy with a minimum value of 1 over land 794 !! 264 !! ** Method : computes from k_top and k_bot with a minimum value of 1 over land 265 !! 266 !! ** Action : mikt, miku, mikv : vertical indices of the shallowest 267 !! ocean level at t-, u- & v-points 268 !! (min value = 1) 795 269 !! ** Action : mbkt, mbku, mbkv : vertical indices of the deeptest 796 270 !! ocean level at t-, u- & v-points 797 271 !! (min value = 1 over land) 798 272 !!---------------------------------------------------------------------- 273 INTEGER , DIMENSION(:,:) , INTENT(out) :: k_top, k_bot ! top & bottom ocean level indices 274 ! 799 275 INTEGER :: ji, jj ! dummy loop indices 800 REAL(wp), POINTER, DIMENSION(:,:) :: z mbk276 REAL(wp), POINTER, DIMENSION(:,:) :: zk 801 277 !!---------------------------------------------------------------------- 802 278 ! 803 279 IF( nn_timing == 1 ) CALL timing_start('zgr_bot_level') 804 280 ! 805 CALL wrk_alloc( jpi, jpj, zmbk )281 CALL wrk_alloc( jpi,jpj, zk ) 806 282 ! 807 283 IF(lwp) WRITE(numout,*) 808 IF(lwp) WRITE(numout,*) ' zgr_bot_level : ocean bottom k-index of T-, U-, V- and W-levels ' 809 IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~' 810 ! 811 mbkt(:,:) = MAX( mbathy(:,:) , 1 ) ! bottom k-index of T-level (=1 over land) 284 IF(lwp) WRITE(numout,*) ' zgr_tb_level : ocean top and bottom k-index of T-, U-, V- and W-levels ' 285 IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~' 286 ! 287 mikt(:,:) = MAX( k_top(:,:) , 1 ) ! top ocean k-index of T-level (=1 over land) 288 ! 289 mbkt(:,:) = MAX( k_bot(:,:) , 1 ) ! bottom ocean k-index of T-level (=1 over land) 812 290 813 ! ! bottom k-index of W-level = mbkt+1 814 DO jj = 1, jpjm1 ! bottom k-index of u- (v-) level 291 ! ! N.B. top k-index of W-level = mikt 292 ! ! bottom k-index of W-level = mbkt+1 293 DO jj = 1, jpjm1 815 294 DO ji = 1, jpim1 295 miku(ji,jj) = MAX( mikt(ji+1,jj ) , mikt(ji,jj) ) 296 mikv(ji,jj) = MAX( mikt(ji ,jj+1) , mikt(ji,jj) ) 297 mikf(ji,jj) = MAX( mikt(ji ,jj+1) , mikt(ji,jj), mikt(ji+1,jj ), mikt(ji+1,jj+1) ) 298 ! 816 299 mbku(ji,jj) = MIN( mbkt(ji+1,jj ) , mbkt(ji,jj) ) 817 300 mbkv(ji,jj) = MIN( mbkt(ji ,jj+1) , mbkt(ji,jj) ) … … 819 302 END DO 820 303 ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk 821 zmbk(:,:) = REAL( mbku(:,:), wp ) ; CALL lbc_lnk(zmbk,'U',1.) ; mbku (:,:) = MAX( INT( zmbk(:,:) ), 1 ) 822 zmbk(:,:) = REAL( mbkv(:,:), wp ) ; CALL lbc_lnk(zmbk,'V',1.) ; mbkv (:,:) = MAX( INT( zmbk(:,:) ), 1 ) 823 ! 824 CALL wrk_dealloc( jpi, jpj, zmbk ) 825 ! 826 IF( nn_timing == 1 ) CALL timing_stop('zgr_bot_level') 827 ! 828 END SUBROUTINE zgr_bot_level 829 830 831 SUBROUTINE zgr_top_level 832 !!---------------------------------------------------------------------- 833 !! *** ROUTINE zgr_top_level *** 834 !! 835 !! ** Purpose : defines the vertical index of ocean top (mik. arrays) 836 !! 837 !! ** Method : computes from misfdep with a minimum value of 1 838 !! 839 !! ** Action : mikt, miku, mikv : vertical indices of the shallowest 840 !! ocean level at t-, u- & v-points 841 !! (min value = 1) 842 !!---------------------------------------------------------------------- 843 INTEGER :: ji, jj ! dummy loop indices 844 REAL(wp), POINTER, DIMENSION(:,:) :: zmik 845 !!---------------------------------------------------------------------- 846 ! 847 IF( nn_timing == 1 ) CALL timing_start('zgr_top_level') 848 ! 849 CALL wrk_alloc( jpi, jpj, zmik ) 850 ! 851 IF(lwp) WRITE(numout,*) 852 IF(lwp) WRITE(numout,*) ' zgr_top_level : ocean top k-index of T-, U-, V- and W-levels ' 853 IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~' 854 ! 855 mikt(:,:) = MAX( misfdep(:,:) , 1 ) ! top k-index of T-level (=1) 856 ! ! top k-index of W-level (=mikt) 857 DO jj = 1, jpjm1 ! top k-index of U- (U-) level 858 DO ji = 1, jpim1 859 miku(ji,jj) = MAX( mikt(ji+1,jj ) , mikt(ji,jj) ) 860 mikv(ji,jj) = MAX( mikt(ji ,jj+1) , mikt(ji,jj) ) 861 mikf(ji,jj) = MAX( mikt(ji ,jj+1) , mikt(ji,jj), mikt(ji+1,jj ), mikt(ji+1,jj+1) ) 862 END DO 863 END DO 864 865 ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk 866 zmik(:,:) = REAL( miku(:,:), wp ) ; CALL lbc_lnk(zmik,'U',1.) ; miku (:,:) = MAX( INT( zmik(:,:) ), 1 ) 867 zmik(:,:) = REAL( mikv(:,:), wp ) ; CALL lbc_lnk(zmik,'V',1.) ; mikv (:,:) = MAX( INT( zmik(:,:) ), 1 ) 868 zmik(:,:) = REAL( mikf(:,:), wp ) ; CALL lbc_lnk(zmik,'F',1.) ; mikf (:,:) = MAX( INT( zmik(:,:) ), 1 ) 869 ! 870 CALL wrk_dealloc( jpi, jpj, zmik ) 304 zk(:,:) = REAL( miku(:,:), wp ) ; CALL lbc_lnk( zk, 'U', 1. ) ; miku(:,:) = MAX( INT( zk(:,:) ), 1 ) 305 zk(:,:) = REAL( mikv(:,:), wp ) ; CALL lbc_lnk( zk, 'V', 1. ) ; mikv(:,:) = MAX( INT( zk(:,:) ), 1 ) 306 zk(:,:) = REAL( mikf(:,:), wp ) ; CALL lbc_lnk( zk, 'F', 1. ) ; mikf(:,:) = MAX( INT( zk(:,:) ), 1 ) 307 ! 308 zk(:,:) = REAL( mbku(:,:), wp ) ; CALL lbc_lnk( zk, 'U', 1. ) ; mbku(:,:) = MAX( INT( zk(:,:) ), 1 ) 309 zk(:,:) = REAL( mbkv(:,:), wp ) ; CALL lbc_lnk( zk, 'V', 1. ) ; mbkv(:,:) = MAX( INT( zk(:,:) ), 1 ) 310 ! 311 CALL wrk_dealloc( jpi,jpj, zk ) 871 312 ! 872 313 IF( nn_timing == 1 ) CALL timing_stop('zgr_top_level') 873 314 ! 874 END SUBROUTINE zgr_top_level 875 876 877 SUBROUTINE zgr_zco 878 !!---------------------------------------------------------------------- 879 !! *** ROUTINE zgr_zco *** 880 !! 881 !! ** Purpose : define the reference z-coordinate system 882 !! 883 !! ** Method : set 3D coord. arrays to reference 1D array 884 !!---------------------------------------------------------------------- 885 INTEGER :: jk 886 !!---------------------------------------------------------------------- 887 ! 888 IF( nn_timing == 1 ) CALL timing_start('zgr_zco') 889 ! 890 DO jk = 1, jpk 891 gdept_0(:,:,jk) = gdept_1d(jk) 892 gdepw_0(:,:,jk) = gdepw_1d(jk) 893 gde3w_0(:,:,jk) = gdepw_1d(jk) 894 e3t_0 (:,:,jk) = e3t_1d (jk) 895 e3u_0 (:,:,jk) = e3t_1d (jk) 896 e3v_0 (:,:,jk) = e3t_1d (jk) 897 e3f_0 (:,:,jk) = e3t_1d (jk) 898 e3w_0 (:,:,jk) = e3w_1d (jk) 899 e3uw_0 (:,:,jk) = e3w_1d (jk) 900 e3vw_0 (:,:,jk) = e3w_1d (jk) 901 END DO 902 ! 903 IF( nn_timing == 1 ) CALL timing_stop('zgr_zco') 904 ! 905 END SUBROUTINE zgr_zco 906 907 908 SUBROUTINE zgr_zps 909 !!---------------------------------------------------------------------- 910 !! *** ROUTINE zgr_zps *** 911 !! 912 !! ** Purpose : the depth and vertical scale factor in partial step 913 !! reference z-coordinate case 914 !! 915 !! ** Method : Partial steps : computes the 3D vertical scale factors 916 !! of T-, U-, V-, W-, UW-, VW and F-points that are associated with 917 !! a partial step representation of bottom topography. 918 !! 919 !! The reference depth of model levels is defined from an analytical 920 !! function the derivative of which gives the reference vertical 921 !! scale factors. 922 !! From depth and scale factors reference, we compute there new value 923 !! with partial steps on 3d arrays ( i, j, k ). 924 !! 925 !! w-level: gdepw_0(i,j,k) = gdep(k) 926 !! e3w_0(i,j,k) = dk(gdep)(k) = e3(i,j,k) 927 !! t-level: gdept_0(i,j,k) = gdep(k+0.5) 928 !! e3t_0(i,j,k) = dk(gdep)(k+0.5) = e3(i,j,k+0.5) 929 !! 930 !! With the help of the bathymetric file ( bathymetry_depth_ORCA_R2.nc), 931 !! we find the mbathy index of the depth at each grid point. 932 !! This leads us to three cases: 933 !! 934 !! - bathy = 0 => mbathy = 0 935 !! - 1 < mbathy < jpkm1 936 !! - bathy > gdepw_0(jpk) => mbathy = jpkm1 937 !! 938 !! Then, for each case, we find the new depth at t- and w- levels 939 !! and the new vertical scale factors at t-, u-, v-, w-, uw-, vw- 940 !! and f-points. 941 !! 942 !! This routine is given as an example, it must be modified 943 !! following the user s desiderata. nevertheless, the output as 944 !! well as the way to compute the model levels and scale factors 945 !! must be respected in order to insure second order accuracy 946 !! schemes. 947 !! 948 !! c a u t i o n : gdept_1d, gdepw_1d and e3._1d are positives 949 !! - - - - - - - gdept_0, gdepw_0 and e3. are positives 950 !! 951 !! Reference : Pacanowsky & Gnanadesikan 1997, Mon. Wea. Rev., 126, 3248-3270. 952 !!---------------------------------------------------------------------- 953 INTEGER :: ji, jj, jk ! dummy loop indices 954 INTEGER :: ik, it, ikb, ikt ! temporary integers 955 REAL(wp) :: ze3tp , ze3wp ! Last ocean level thickness at T- and W-points 956 REAL(wp) :: zdepwp, zdepth ! Ajusted ocean depth to avoid too small e3t 957 REAL(wp) :: zdiff ! temporary scalar 958 REAL(wp) :: zmax ! temporary scalar 959 REAL(wp), POINTER, DIMENSION(:,:,:) :: zprt 960 !!--------------------------------------------------------------------- 961 ! 962 IF( nn_timing == 1 ) CALL timing_start('zgr_zps') 963 ! 964 CALL wrk_alloc( jpi,jpj,jpk, zprt ) 965 ! 966 IF(lwp) WRITE(numout,*) 967 IF(lwp) WRITE(numout,*) ' zgr_zps : z-coordinate with partial steps' 968 IF(lwp) WRITE(numout,*) ' ~~~~~~~ ' 969 IF(lwp) WRITE(numout,*) ' mbathy is recomputed : bathy_level file is NOT used' 970 971 ! bathymetry in level (from bathy_meter) 972 ! =================== 973 zmax = gdepw_1d(jpk) + e3t_1d(jpk) ! maximum depth (i.e. the last ocean level thickness <= 2*e3t_1d(jpkm1) ) 974 bathy(:,:) = MIN( zmax , bathy(:,:) ) ! bounded value of bathy (min already set at the end of zgr_bat) 975 WHERE( bathy(:,:) == 0._wp ) ; mbathy(:,:) = 0 ! land : set mbathy to 0 976 ELSE WHERE ; mbathy(:,:) = jpkm1 ! ocean : initialize mbathy to the max ocean level 977 END WHERE 978 979 ! Compute mbathy for ocean points (i.e. the number of ocean levels) 980 ! find the number of ocean levels such that the last level thickness 981 ! is larger than the minimum of e3zps_min and e3zps_rat * e3t_1d (where 982 ! e3t_1d is the reference level thickness 983 DO jk = jpkm1, 1, -1 984 zdepth = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 985 WHERE( 0._wp < bathy(:,:) .AND. bathy(:,:) <= zdepth ) mbathy(:,:) = jk-1 986 END DO 987 988 ! Scale factors and depth at T- and W-points 989 DO jk = 1, jpk ! intitialization to the reference z-coordinate 990 gdept_0(:,:,jk) = gdept_1d(jk) 991 gdepw_0(:,:,jk) = gdepw_1d(jk) 992 e3t_0 (:,:,jk) = e3t_1d (jk) 993 e3w_0 (:,:,jk) = e3w_1d (jk) 994 END DO 995 996 ! Bathy, iceshelf draft, scale factor and depth at T- and W- points in case of isf 997 IF ( ln_isfcav ) CALL zgr_isf 998 999 ! Scale factors and depth at T- and W-points 1000 IF ( .NOT. ln_isfcav ) THEN 1001 DO jj = 1, jpj 1002 DO ji = 1, jpi 1003 ik = mbathy(ji,jj) 1004 IF( ik > 0 ) THEN ! ocean point only 1005 ! max ocean level case 1006 IF( ik == jpkm1 ) THEN 1007 zdepwp = bathy(ji,jj) 1008 ze3tp = bathy(ji,jj) - gdepw_1d(ik) 1009 ze3wp = 0.5_wp * e3w_1d(ik) * ( 1._wp + ( ze3tp/e3t_1d(ik) ) ) 1010 e3t_0(ji,jj,ik ) = ze3tp 1011 e3t_0(ji,jj,ik+1) = ze3tp 1012 e3w_0(ji,jj,ik ) = ze3wp 1013 e3w_0(ji,jj,ik+1) = ze3tp 1014 gdepw_0(ji,jj,ik+1) = zdepwp 1015 gdept_0(ji,jj,ik ) = gdept_1d(ik-1) + ze3wp 1016 gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + ze3tp 1017 ! 1018 ELSE ! standard case 1019 IF( bathy(ji,jj) <= gdepw_1d(ik+1) ) THEN ; gdepw_0(ji,jj,ik+1) = bathy(ji,jj) 1020 ELSE ; gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1) 1021 ENDIF 1022 !gm Bug? check the gdepw_1d 1023 ! ... on ik 1024 gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) ) & 1025 & * ((gdept_1d( ik ) - gdepw_1d(ik) ) & 1026 & / ( gdepw_1d( ik+1) - gdepw_1d(ik) )) 1027 e3t_0 (ji,jj,ik) = e3t_1d (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) ) & 1028 & / ( gdepw_1d( ik+1) - gdepw_1d(ik) ) 1029 e3w_0(ji,jj,ik) = 0.5_wp * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2._wp * gdepw_1d(ik) ) & 1030 & * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) ) 1031 ! ... on ik+1 1032 e3w_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik) 1033 e3t_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik) 1034 gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik) 1035 ENDIF 1036 ENDIF 1037 END DO 1038 END DO 1039 ! 1040 it = 0 1041 DO jj = 1, jpj 1042 DO ji = 1, jpi 1043 ik = mbathy(ji,jj) 1044 IF( ik > 0 ) THEN ! ocean point only 1045 e3tp (ji,jj) = e3t_0(ji,jj,ik) 1046 e3wp (ji,jj) = e3w_0(ji,jj,ik) 1047 ! test 1048 zdiff= gdepw_0(ji,jj,ik+1) - gdept_0(ji,jj,ik ) 1049 IF( zdiff <= 0._wp .AND. lwp ) THEN 1050 it = it + 1 1051 WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj 1052 WRITE(numout,*) ' bathy = ', bathy(ji,jj) 1053 WRITE(numout,*) ' gdept_0 = ', gdept_0(ji,jj,ik), ' gdepw_0 = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff 1054 WRITE(numout,*) ' e3tp = ', e3t_0 (ji,jj,ik), ' e3wp = ', e3w_0 (ji,jj,ik ) 1055 ENDIF 1056 ENDIF 1057 END DO 1058 END DO 1059 END IF 1060 ! 1061 ! Scale factors and depth at U-, V-, UW and VW-points 1062 DO jk = 1, jpk ! initialisation to z-scale factors 1063 e3u_0 (:,:,jk) = e3t_1d(jk) 1064 e3v_0 (:,:,jk) = e3t_1d(jk) 1065 e3uw_0(:,:,jk) = e3w_1d(jk) 1066 e3vw_0(:,:,jk) = e3w_1d(jk) 1067 END DO 1068 1069 DO jk = 1,jpk ! Computed as the minimum of neighbooring scale factors 1070 DO jj = 1, jpjm1 1071 DO ji = 1, fs_jpim1 ! vector opt. 1072 e3u_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji+1,jj,jk) ) 1073 e3v_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji,jj+1,jk) ) 1074 e3uw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji+1,jj,jk) ) 1075 e3vw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji,jj+1,jk) ) 1076 END DO 1077 END DO 1078 END DO 1079 IF ( ln_isfcav ) THEN 1080 ! (ISF) define e3uw (adapted for 2 cells in the water column) 1081 DO jj = 2, jpjm1 1082 DO ji = 2, fs_jpim1 ! vector opt. 1083 ikb = MAX(mbathy (ji,jj),mbathy (ji+1,jj)) 1084 ikt = MAX(misfdep(ji,jj),misfdep(ji+1,jj)) 1085 IF (ikb == ikt+1) e3uw_0(ji,jj,ikb) = MIN( gdept_0(ji,jj,ikb ), gdept_0(ji+1,jj ,ikb ) ) & 1086 & - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji+1,jj ,ikb-1) ) 1087 ikb = MAX(mbathy (ji,jj),mbathy (ji,jj+1)) 1088 ikt = MAX(misfdep(ji,jj),misfdep(ji,jj+1)) 1089 IF (ikb == ikt+1) e3vw_0(ji,jj,ikb) = MIN( gdept_0(ji,jj,ikb ), gdept_0(ji ,jj+1,ikb ) ) & 1090 & - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji ,jj+1,ikb-1) ) 1091 END DO 1092 END DO 1093 END IF 1094 1095 CALL lbc_lnk( e3u_0 , 'U', 1._wp ) ; CALL lbc_lnk( e3uw_0, 'U', 1._wp ) ! lateral boundary conditions 1096 CALL lbc_lnk( e3v_0 , 'V', 1._wp ) ; CALL lbc_lnk( e3vw_0, 'V', 1._wp ) 1097 ! 1098 1099 DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries) 1100 WHERE( e3u_0 (:,:,jk) == 0._wp ) e3u_0 (:,:,jk) = e3t_1d(jk) 1101 WHERE( e3v_0 (:,:,jk) == 0._wp ) e3v_0 (:,:,jk) = e3t_1d(jk) 1102 WHERE( e3uw_0(:,:,jk) == 0._wp ) e3uw_0(:,:,jk) = e3w_1d(jk) 1103 WHERE( e3vw_0(:,:,jk) == 0._wp ) e3vw_0(:,:,jk) = e3w_1d(jk) 1104 END DO 1105 1106 ! Scale factor at F-point 1107 DO jk = 1, jpk ! initialisation to z-scale factors 1108 e3f_0(:,:,jk) = e3t_1d(jk) 1109 END DO 1110 DO jk = 1, jpk ! Computed as the minimum of neighbooring V-scale factors 1111 DO jj = 1, jpjm1 1112 DO ji = 1, fs_jpim1 ! vector opt. 1113 e3f_0(ji,jj,jk) = MIN( e3v_0(ji,jj,jk), e3v_0(ji+1,jj,jk) ) 1114 END DO 1115 END DO 1116 END DO 1117 CALL lbc_lnk( e3f_0, 'F', 1._wp ) ! Lateral boundary conditions 1118 ! 1119 DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries) 1120 WHERE( e3f_0(:,:,jk) == 0._wp ) e3f_0(:,:,jk) = e3t_1d(jk) 1121 END DO 1122 !!gm bug ? : must be a do loop with mj0,mj1 1123 ! 1124 e3t_0(:,mj0(1),:) = e3t_0(:,mj0(2),:) ! we duplicate factor scales for jj = 1 and jj = 2 1125 e3w_0(:,mj0(1),:) = e3w_0(:,mj0(2),:) 1126 e3u_0(:,mj0(1),:) = e3u_0(:,mj0(2),:) 1127 e3v_0(:,mj0(1),:) = e3v_0(:,mj0(2),:) 1128 e3f_0(:,mj0(1),:) = e3f_0(:,mj0(2),:) 1129 1130 ! Control of the sign 1131 IF( MINVAL( e3t_0 (:,:,:) ) <= 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r e3t_0 <= 0' ) 1132 IF( MINVAL( e3w_0 (:,:,:) ) <= 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r e3w_0 <= 0' ) 1133 IF( MINVAL( gdept_0(:,:,:) ) < 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r gdept_0 < 0' ) 1134 IF( MINVAL( gdepw_0(:,:,:) ) < 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r gdepw_0 < 0' ) 1135 1136 ! Compute gde3w_0 (vertical sum of e3w) 1137 IF ( ln_isfcav ) THEN ! if cavity 1138 WHERE( misfdep == 0 ) misfdep = 1 1139 DO jj = 1,jpj 1140 DO ji = 1,jpi 1141 gde3w_0(ji,jj,1) = 0.5_wp * e3w_0(ji,jj,1) 1142 DO jk = 2, misfdep(ji,jj) 1143 gde3w_0(ji,jj,jk) = gde3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk) 1144 END DO 1145 IF( misfdep(ji,jj) >= 2 ) gde3w_0(ji,jj,misfdep(ji,jj)) = risfdep(ji,jj) + 0.5_wp * e3w_0(ji,jj,misfdep(ji,jj)) 1146 DO jk = misfdep(ji,jj) + 1, jpk 1147 gde3w_0(ji,jj,jk) = gde3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk) 1148 END DO 1149 END DO 1150 END DO 1151 ELSE ! no cavity 1152 gde3w_0(:,:,1) = 0.5_wp * e3w_0(:,:,1) 1153 DO jk = 2, jpk 1154 gde3w_0(:,:,jk) = gde3w_0(:,:,jk-1) + e3w_0(:,:,jk) 1155 END DO 1156 END IF 1157 ! 1158 CALL wrk_dealloc( jpi,jpj,jpk, zprt ) 1159 ! 1160 IF( nn_timing == 1 ) CALL timing_stop('zgr_zps') 1161 ! 1162 END SUBROUTINE zgr_zps 1163 1164 1165 SUBROUTINE zgr_isf 1166 !!---------------------------------------------------------------------- 1167 !! *** ROUTINE zgr_isf *** 1168 !! 1169 !! ** Purpose : check the bathymetry in levels 1170 !! 1171 !! ** Method : THe water column have to contained at least 2 cells 1172 !! Bathymetry and isfdraft are modified (dig/close) to respect 1173 !! this criterion. 1174 !! 1175 !! ** Action : - test compatibility between isfdraft and bathy 1176 !! - bathy and isfdraft are modified 1177 !!---------------------------------------------------------------------- 1178 INTEGER :: ji, jj, jl, jk ! dummy loop indices 1179 INTEGER :: ik, it ! temporary integers 1180 INTEGER :: icompt, ibtest ! (ISF) 1181 INTEGER :: ibtestim1, ibtestip1 ! (ISF) 1182 INTEGER :: ibtestjm1, ibtestjp1 ! (ISF) 1183 REAL(wp) :: zdepth ! Ajusted ocean depth to avoid too small e3t 1184 REAL(wp) :: zmax ! Maximum and minimum depth 1185 REAL(wp) :: zbathydiff ! isf temporary scalar 1186 REAL(wp) :: zrisfdepdiff ! isf temporary scalar 1187 REAL(wp) :: ze3tp , ze3wp ! Last ocean level thickness at T- and W-points 1188 REAL(wp) :: zdepwp ! Ajusted ocean depth to avoid too small e3t 1189 REAL(wp) :: zdiff ! temporary scalar 1190 REAL(wp), POINTER, DIMENSION(:,:) :: zrisfdep, zbathy, zmask ! 2D workspace (ISH) 1191 INTEGER , POINTER, DIMENSION(:,:) :: zmbathy, zmisfdep ! 2D workspace (ISH) 1192 !!--------------------------------------------------------------------- 1193 ! 1194 IF( nn_timing == 1 ) CALL timing_start('zgr_isf') 1195 ! 1196 CALL wrk_alloc( jpi,jpj, zbathy, zmask, zrisfdep) 1197 CALL wrk_alloc( jpi,jpj, zmisfdep, zmbathy ) 1198 1199 ! (ISF) compute misfdep 1200 WHERE( risfdep(:,:) == 0._wp .AND. bathy(:,:) /= 0 ) ; misfdep(:,:) = 1 ! open water : set misfdep to 1 1201 ELSEWHERE ; misfdep(:,:) = 2 ! iceshelf : initialize misfdep to second level 1202 END WHERE 1203 1204 ! Compute misfdep for ocean points (i.e. first wet level) 1205 ! find the first ocean level such that the first level thickness 1206 ! is larger than the bot_level of e3zps_min and e3zps_rat * e3t_0 (where 1207 ! e3t_0 is the reference level thickness 1208 DO jk = 2, jpkm1 1209 zdepth = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 1210 WHERE( 0._wp < risfdep(:,:) .AND. risfdep(:,:) >= zdepth ) misfdep(:,:) = jk+1 1211 END DO 1212 WHERE ( 0._wp < risfdep(:,:) .AND. risfdep(:,:) <= e3t_1d(1) ) 1213 risfdep(:,:) = 0._wp ; misfdep(:,:) = 1 1214 END WHERE 1215 1216 ! remove very shallow ice shelf (less than ~ 10m if 75L) 1217 WHERE (risfdep(:,:) <= 10._wp .AND. misfdep(:,:) > 1) 1218 misfdep = 0; risfdep = 0.0_wp; 1219 mbathy = 0; bathy = 0.0_wp; 1220 END WHERE 1221 WHERE (bathy(:,:) <= 30.0_wp .AND. gphit < -60._wp) 1222 misfdep = 0; risfdep = 0.0_wp; 1223 mbathy = 0; bathy = 0.0_wp; 1224 END WHERE 1225 1226 ! basic check for the compatibility of bathy and risfdep. I think it should be offline because it is not perfect and cannot solved all the situation 1227 icompt = 0 1228 ! run the bathy check 10 times to be sure all the modif in the bathy or iceshelf draft are compatible together 1229 DO jl = 1, 10 1230 ! check at each iteration if isf is grounded or not (1cm treshold have to be update after first coupling experiments) 1231 WHERE (bathy(:,:) <= risfdep(:,:) + rn_isfhmin) 1232 misfdep(:,:) = 0 ; risfdep(:,:) = 0._wp 1233 mbathy (:,:) = 0 ; bathy (:,:) = 0._wp 1234 END WHERE 1235 WHERE (mbathy(:,:) <= 0) 1236 misfdep(:,:) = 0 ; risfdep(:,:) = 0._wp 1237 mbathy (:,:) = 0 ; bathy (:,:) = 0._wp 1238 END WHERE 1239 IF( lk_mpp ) THEN 1240 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1241 CALL lbc_lnk( zbathy, 'T', 1. ) 1242 misfdep(:,:) = INT( zbathy(:,:) ) 1243 1244 CALL lbc_lnk( risfdep,'T', 1. ) 1245 CALL lbc_lnk( bathy, 'T', 1. ) 1246 1247 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1248 CALL lbc_lnk( zbathy, 'T', 1. ) 1249 mbathy(:,:) = INT( zbathy(:,:) ) 1250 ENDIF 1251 IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN 1252 misfdep( 1 ,:) = misfdep(jpim1,:) ! local domain is cyclic east-west 1253 misfdep(jpi,:) = misfdep( 2 ,:) 1254 mbathy( 1 ,:) = mbathy(jpim1,:) ! local domain is cyclic east-west 1255 mbathy(jpi,:) = mbathy( 2 ,:) 1256 ENDIF 1257 1258 ! split last cell if possible (only where water column is 2 cell or less) 1259 ! if coupled to ice sheet, we do not modify the bathymetry (can be discuss). 1260 IF ( .NOT. ln_iscpl) THEN 1261 DO jk = jpkm1, 1, -1 1262 zmax = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 1263 WHERE( gdepw_1d(jk) < bathy(:,:) .AND. bathy(:,:) <= zmax .AND. misfdep + 1 >= mbathy) 1264 mbathy(:,:) = jk 1265 bathy(:,:) = zmax 1266 END WHERE 1267 END DO 1268 END IF 1269 1270 ! split top cell if possible (only where water column is 2 cell or less) 1271 DO jk = 2, jpkm1 1272 zmax = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 1273 WHERE( gdepw_1d(jk+1) > risfdep(:,:) .AND. risfdep(:,:) >= zmax .AND. misfdep + 1 >= mbathy) 1274 misfdep(:,:) = jk 1275 risfdep(:,:) = zmax 1276 END WHERE 1277 END DO 1278 1279 1280 ! Case where bathy and risfdep compatible but not the level variable mbathy/misfdep because of partial cell condition 1281 DO jj = 1, jpj 1282 DO ji = 1, jpi 1283 ! find the minimum change option: 1284 ! test bathy 1285 IF (risfdep(ji,jj) > 1) THEN 1286 IF ( .NOT. ln_iscpl ) THEN 1287 zbathydiff =ABS(bathy(ji,jj) - (gdepw_1d(mbathy (ji,jj)+1) & 1288 & + MIN( e3zps_min, e3t_1d(mbathy (ji,jj)+1)*e3zps_rat ))) 1289 zrisfdepdiff=ABS(risfdep(ji,jj) - (gdepw_1d(misfdep(ji,jj) ) & 1290 & - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ))) 1291 IF (bathy(ji,jj) > risfdep(ji,jj) .AND. mbathy(ji,jj) < misfdep(ji,jj)) THEN 1292 IF (zbathydiff <= zrisfdepdiff) THEN 1293 bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj)+1)*e3zps_rat ) 1294 mbathy(ji,jj)= mbathy(ji,jj) + 1 1295 ELSE 1296 risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ) 1297 misfdep(ji,jj) = misfdep(ji,jj) - 1 1298 END IF 1299 ENDIF 1300 ELSE 1301 IF (bathy(ji,jj) > risfdep(ji,jj) .AND. mbathy(ji,jj) < misfdep(ji,jj)) THEN 1302 risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ) 1303 misfdep(ji,jj) = misfdep(ji,jj) - 1 1304 END IF 1305 END IF 1306 END IF 1307 END DO 1308 END DO 1309 1310 ! At least 2 levels for water thickness at T, U, and V point. 1311 DO jj = 1, jpj 1312 DO ji = 1, jpi 1313 ! find the minimum change option: 1314 ! test bathy 1315 IF( misfdep(ji,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) > 1) THEN 1316 IF ( .NOT. ln_iscpl ) THEN 1317 zbathydiff =ABS(bathy(ji,jj) - ( gdepw_1d(mbathy (ji,jj)+1) & 1318 & + MIN( e3zps_min,e3t_1d(mbathy (ji,jj)+1)*e3zps_rat ))) 1319 zrisfdepdiff=ABS(risfdep(ji,jj) - ( gdepw_1d(misfdep(ji,jj) ) & 1320 & - MIN( e3zps_min,e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ))) 1321 IF (zbathydiff <= zrisfdepdiff) THEN 1322 mbathy(ji,jj) = mbathy(ji,jj) + 1 1323 bathy(ji,jj) = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat ) 1324 ELSE 1325 misfdep(ji,jj)= misfdep(ji,jj) - 1 1326 risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj))*e3zps_rat ) 1327 END IF 1328 ELSE 1329 misfdep(ji,jj)= misfdep(ji,jj) - 1 1330 risfdep(ji,jj)= gdepw_1d(misfdep(ji,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj))*e3zps_rat ) 1331 END IF 1332 ENDIF 1333 END DO 1334 END DO 1335 1336 ! point V mbathy(ji,jj) == misfdep(ji,jj+1) 1337 DO jj = 1, jpjm1 1338 DO ji = 1, jpim1 1339 IF( misfdep(ji,jj+1) == mbathy(ji,jj) .AND. mbathy(ji,jj) > 1) THEN 1340 IF ( .NOT. ln_iscpl ) THEN 1341 zbathydiff =ABS(bathy(ji,jj ) - ( gdepw_1d(mbathy (ji,jj)+1) & 1342 & + MIN( e3zps_min, e3t_1d(mbathy (ji,jj )+1)*e3zps_rat ))) 1343 zrisfdepdiff=ABS(risfdep(ji,jj+1) - ( gdepw_1d(misfdep(ji,jj+1)) & 1344 & - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1)-1)*e3zps_rat ))) 1345 IF (zbathydiff <= zrisfdepdiff) THEN 1346 mbathy(ji,jj) = mbathy(ji,jj) + 1 1347 bathy(ji,jj) = gdepw_1d(mbathy (ji,jj )) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj )+1)*e3zps_rat ) 1348 ELSE 1349 misfdep(ji,jj+1) = misfdep(ji,jj+1) - 1 1350 risfdep (ji,jj+1) = gdepw_1d(misfdep(ji,jj+1)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1))*e3zps_rat ) 1351 END IF 1352 ELSE 1353 misfdep(ji,jj+1) = misfdep(ji,jj+1) - 1 1354 risfdep (ji,jj+1) = gdepw_1d(misfdep(ji,jj+1)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1))*e3zps_rat ) 1355 END IF 1356 ENDIF 1357 END DO 1358 END DO 1359 1360 IF( lk_mpp ) THEN 1361 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1362 CALL lbc_lnk( zbathy, 'T', 1. ) 1363 misfdep(:,:) = INT( zbathy(:,:) ) 1364 1365 CALL lbc_lnk( risfdep,'T', 1. ) 1366 CALL lbc_lnk( bathy, 'T', 1. ) 1367 1368 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1369 CALL lbc_lnk( zbathy, 'T', 1. ) 1370 mbathy(:,:) = INT( zbathy(:,:) ) 1371 ENDIF 1372 ! point V misdep(ji,jj) == mbathy(ji,jj+1) 1373 DO jj = 1, jpjm1 1374 DO ji = 1, jpim1 1375 IF( misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) > 1) THEN 1376 IF ( .NOT. ln_iscpl ) THEN 1377 zbathydiff =ABS( bathy(ji,jj+1) - ( gdepw_1d(mbathy (ji,jj+1)+1) & 1378 & + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat ))) 1379 zrisfdepdiff=ABS(risfdep(ji,jj ) - ( gdepw_1d(misfdep(ji,jj ) ) & 1380 & - MIN( e3zps_min, e3t_1d(misfdep(ji,jj )-1)*e3zps_rat ))) 1381 IF (zbathydiff <= zrisfdepdiff) THEN 1382 mbathy (ji,jj+1) = mbathy(ji,jj+1) + 1 1383 bathy (ji,jj+1) = gdepw_1d(mbathy (ji,jj+1) ) + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat ) 1384 ELSE 1385 misfdep(ji,jj) = misfdep(ji,jj) - 1 1386 risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj )+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj ) )*e3zps_rat ) 1387 END IF 1388 ELSE 1389 misfdep(ji,jj) = misfdep(ji,jj) - 1 1390 risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj )+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj ) )*e3zps_rat ) 1391 END IF 1392 ENDIF 1393 END DO 1394 END DO 1395 1396 1397 IF( lk_mpp ) THEN 1398 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1399 CALL lbc_lnk( zbathy, 'T', 1. ) 1400 misfdep(:,:) = INT( zbathy(:,:) ) 1401 1402 CALL lbc_lnk( risfdep,'T', 1. ) 1403 CALL lbc_lnk( bathy, 'T', 1. ) 1404 1405 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1406 CALL lbc_lnk( zbathy, 'T', 1. ) 1407 mbathy(:,:) = INT( zbathy(:,:) ) 1408 ENDIF 1409 1410 ! point U mbathy(ji,jj) == misfdep(ji,jj+1) 1411 DO jj = 1, jpjm1 1412 DO ji = 1, jpim1 1413 IF( misfdep(ji+1,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) > 1) THEN 1414 IF ( .NOT. ln_iscpl ) THEN 1415 zbathydiff =ABS( bathy(ji ,jj) - ( gdepw_1d(mbathy (ji,jj)+1) & 1416 & + MIN( e3zps_min, e3t_1d(mbathy (ji ,jj)+1)*e3zps_rat ))) 1417 zrisfdepdiff=ABS(risfdep(ji+1,jj) - ( gdepw_1d(misfdep(ji+1,jj)) & 1418 & - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj)-1)*e3zps_rat ))) 1419 IF (zbathydiff <= zrisfdepdiff) THEN 1420 mbathy(ji,jj) = mbathy(ji,jj) + 1 1421 bathy(ji,jj) = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat ) 1422 ELSE 1423 misfdep(ji+1,jj)= misfdep(ji+1,jj) - 1 1424 risfdep(ji+1,jj) = gdepw_1d(misfdep(ji+1,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj))*e3zps_rat ) 1425 END IF 1426 ELSE 1427 misfdep(ji+1,jj)= misfdep(ji+1,jj) - 1 1428 risfdep(ji+1,jj) = gdepw_1d(misfdep(ji+1,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj))*e3zps_rat ) 1429 ENDIF 1430 ENDIF 1431 ENDDO 1432 ENDDO 1433 1434 IF( lk_mpp ) THEN 1435 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1436 CALL lbc_lnk( zbathy, 'T', 1. ) 1437 misfdep(:,:) = INT( zbathy(:,:) ) 1438 1439 CALL lbc_lnk( risfdep,'T', 1. ) 1440 CALL lbc_lnk( bathy, 'T', 1. ) 1441 1442 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1443 CALL lbc_lnk( zbathy, 'T', 1. ) 1444 mbathy(:,:) = INT( zbathy(:,:) ) 1445 ENDIF 1446 1447 ! point U misfdep(ji,jj) == bathy(ji,jj+1) 1448 DO jj = 1, jpjm1 1449 DO ji = 1, jpim1 1450 IF( misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) > 1) THEN 1451 IF ( .NOT. ln_iscpl ) THEN 1452 zbathydiff =ABS( bathy(ji+1,jj) - ( gdepw_1d(mbathy (ji+1,jj)+1) & 1453 & + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj)+1)*e3zps_rat ))) 1454 zrisfdepdiff=ABS(risfdep(ji ,jj) - ( gdepw_1d(misfdep(ji ,jj) ) & 1455 & - MIN( e3zps_min, e3t_1d(misfdep(ji ,jj)-1)*e3zps_rat ))) 1456 IF (zbathydiff <= zrisfdepdiff) THEN 1457 mbathy(ji+1,jj) = mbathy (ji+1,jj) + 1 1458 bathy (ji+1,jj) = gdepw_1d(mbathy (ji+1,jj) ) + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj) +1)*e3zps_rat ) 1459 ELSE 1460 misfdep(ji,jj) = misfdep(ji ,jj) - 1 1461 risfdep(ji,jj) = gdepw_1d(misfdep(ji ,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji ,jj) )*e3zps_rat ) 1462 END IF 1463 ELSE 1464 misfdep(ji,jj) = misfdep(ji ,jj) - 1 1465 risfdep(ji,jj) = gdepw_1d(misfdep(ji ,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji ,jj) )*e3zps_rat ) 1466 ENDIF 1467 ENDIF 1468 ENDDO 1469 ENDDO 1470 1471 IF( lk_mpp ) THEN 1472 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1473 CALL lbc_lnk( zbathy, 'T', 1. ) 1474 misfdep(:,:) = INT( zbathy(:,:) ) 1475 1476 CALL lbc_lnk( risfdep,'T', 1. ) 1477 CALL lbc_lnk( bathy, 'T', 1. ) 1478 1479 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1480 CALL lbc_lnk( zbathy, 'T', 1. ) 1481 mbathy(:,:) = INT( zbathy(:,:) ) 1482 ENDIF 1483 END DO 1484 ! end dig bathy/ice shelf to be compatible 1485 ! now fill single point in "coastline" of ice shelf, bathy, hole, and test again one cell tickness 1486 DO jl = 1,20 1487 1488 ! remove single point "bay" on isf coast line in the ice shelf draft' 1489 DO jk = 2, jpk 1490 WHERE (misfdep==0) misfdep=jpk 1491 zmask=0._wp 1492 WHERE (misfdep <= jk) zmask=1 1493 DO jj = 2, jpjm1 1494 DO ji = 2, jpim1 1495 IF (misfdep(ji,jj) == jk) THEN 1496 ibtest = zmask(ji-1,jj) + zmask(ji+1,jj) + zmask(ji,jj-1) + zmask(ji,jj+1) 1497 IF (ibtest <= 1) THEN 1498 risfdep(ji,jj)=gdepw_1d(jk+1) ; misfdep(ji,jj)=jk+1 1499 IF (misfdep(ji,jj) > mbathy(ji,jj)) misfdep(ji,jj) = jpk 1500 END IF 1501 END IF 1502 END DO 1503 END DO 1504 END DO 1505 WHERE (misfdep==jpk) 1506 misfdep=0 ; risfdep=0._wp ; mbathy=0 ; bathy=0._wp 1507 END WHERE 1508 IF( lk_mpp ) THEN 1509 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1510 CALL lbc_lnk( zbathy, 'T', 1. ) 1511 misfdep(:,:) = INT( zbathy(:,:) ) 1512 1513 CALL lbc_lnk( risfdep,'T', 1. ) 1514 CALL lbc_lnk( bathy, 'T', 1. ) 1515 1516 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1517 CALL lbc_lnk( zbathy, 'T', 1. ) 1518 mbathy(:,:) = INT( zbathy(:,:) ) 1519 ENDIF 1520 1521 ! remove single point "bay" on bathy coast line beneath an ice shelf' 1522 DO jk = jpk,1,-1 1523 zmask=0._wp 1524 WHERE (mbathy >= jk ) zmask=1 1525 DO jj = 2, jpjm1 1526 DO ji = 2, jpim1 1527 IF (mbathy(ji,jj) == jk .AND. misfdep(ji,jj) >= 2) THEN 1528 ibtest = zmask(ji-1,jj) + zmask(ji+1,jj) + zmask(ji,jj-1) + zmask(ji,jj+1) 1529 IF (ibtest <= 1) THEN 1530 bathy(ji,jj)=gdepw_1d(jk) ; mbathy(ji,jj)=jk-1 1531 IF (misfdep(ji,jj) > mbathy(ji,jj)) mbathy(ji,jj) = 0 1532 END IF 1533 END IF 1534 END DO 1535 END DO 1536 END DO 1537 WHERE (mbathy==0) 1538 misfdep=0 ; risfdep=0._wp ; mbathy=0 ; bathy=0._wp 1539 END WHERE 1540 IF( lk_mpp ) THEN 1541 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1542 CALL lbc_lnk( zbathy, 'T', 1. ) 1543 misfdep(:,:) = INT( zbathy(:,:) ) 1544 1545 CALL lbc_lnk( risfdep,'T', 1. ) 1546 CALL lbc_lnk( bathy, 'T', 1. ) 1547 1548 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1549 CALL lbc_lnk( zbathy, 'T', 1. ) 1550 mbathy(:,:) = INT( zbathy(:,:) ) 1551 ENDIF 1552 1553 ! fill hole in ice shelf 1554 zmisfdep = misfdep 1555 zrisfdep = risfdep 1556 WHERE (zmisfdep <= 1._wp) zmisfdep=jpk 1557 DO jj = 2, jpjm1 1558 DO ji = 2, jpim1 1559 ibtestim1 = zmisfdep(ji-1,jj ) ; ibtestip1 = zmisfdep(ji+1,jj ) 1560 ibtestjm1 = zmisfdep(ji ,jj-1) ; ibtestjp1 = zmisfdep(ji ,jj+1) 1561 IF( zmisfdep(ji,jj) >= mbathy(ji-1,jj ) ) ibtestim1 = jpk 1562 IF( zmisfdep(ji,jj) >= mbathy(ji+1,jj ) ) ibtestip1 = jpk 1563 IF( zmisfdep(ji,jj) >= mbathy(ji ,jj-1) ) ibtestjm1 = jpk 1564 IF( zmisfdep(ji,jj) >= mbathy(ji ,jj+1) ) ibtestjp1 = jpk 1565 ibtest=MIN(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1) 1566 IF( ibtest == jpk .AND. misfdep(ji,jj) >= 2) THEN 1567 mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp 1568 END IF 1569 IF( zmisfdep(ji,jj) < ibtest .AND. misfdep(ji,jj) >= 2) THEN 1570 misfdep(ji,jj) = ibtest 1571 risfdep(ji,jj) = gdepw_1d(ibtest) 1572 ENDIF 1573 ENDDO 1574 ENDDO 1575 1576 IF( lk_mpp ) THEN 1577 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1578 CALL lbc_lnk( zbathy, 'T', 1. ) 1579 misfdep(:,:) = INT( zbathy(:,:) ) 1580 1581 CALL lbc_lnk( risfdep, 'T', 1. ) 1582 CALL lbc_lnk( bathy, 'T', 1. ) 1583 1584 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1585 CALL lbc_lnk( zbathy, 'T', 1. ) 1586 mbathy(:,:) = INT( zbathy(:,:) ) 1587 ENDIF 1588 ! 1589 !! fill hole in bathymetry 1590 zmbathy (:,:)=mbathy (:,:) 1591 DO jj = 2, jpjm1 1592 DO ji = 2, jpim1 1593 ibtestim1 = zmbathy(ji-1,jj ) ; ibtestip1 = zmbathy(ji+1,jj ) 1594 ibtestjm1 = zmbathy(ji ,jj-1) ; ibtestjp1 = zmbathy(ji ,jj+1) 1595 IF( zmbathy(ji,jj) < misfdep(ji-1,jj ) ) ibtestim1 = 0 1596 IF( zmbathy(ji,jj) < misfdep(ji+1,jj ) ) ibtestip1 = 0 1597 IF( zmbathy(ji,jj) < misfdep(ji ,jj-1) ) ibtestjm1 = 0 1598 IF( zmbathy(ji,jj) < misfdep(ji ,jj+1) ) ibtestjp1 = 0 1599 ibtest=MAX(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1) 1600 IF( ibtest == 0 .AND. misfdep(ji,jj) >= 2) THEN 1601 mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp ; 1602 END IF 1603 IF( ibtest < zmbathy(ji,jj) .AND. misfdep(ji,jj) >= 2) THEN 1604 mbathy(ji,jj) = ibtest 1605 bathy(ji,jj) = gdepw_1d(ibtest+1) 1606 ENDIF 1607 END DO 1608 END DO 1609 IF( lk_mpp ) THEN 1610 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1611 CALL lbc_lnk( zbathy, 'T', 1. ) 1612 misfdep(:,:) = INT( zbathy(:,:) ) 1613 1614 CALL lbc_lnk( risfdep, 'T', 1. ) 1615 CALL lbc_lnk( bathy, 'T', 1. ) 1616 1617 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1618 CALL lbc_lnk( zbathy, 'T', 1. ) 1619 mbathy(:,:) = INT( zbathy(:,:) ) 1620 ENDIF 1621 ! if not compatible after all check (ie U point water column less than 2 cells), mask U 1622 DO jj = 1, jpjm1 1623 DO ji = 1, jpim1 1624 IF (mbathy(ji,jj) == misfdep(ji+1,jj) .AND. mbathy(ji,jj) >= 1 .AND. mbathy(ji+1,jj) >= 1) THEN 1625 mbathy(ji,jj) = mbathy(ji,jj) - 1 ; bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)+1) ; 1626 END IF 1627 END DO 1628 END DO 1629 IF( lk_mpp ) THEN 1630 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1631 CALL lbc_lnk( zbathy, 'T', 1. ) 1632 misfdep(:,:) = INT( zbathy(:,:) ) 1633 1634 CALL lbc_lnk( risfdep, 'T', 1. ) 1635 CALL lbc_lnk( bathy, 'T', 1. ) 1636 1637 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1638 CALL lbc_lnk( zbathy, 'T', 1. ) 1639 mbathy(:,:) = INT( zbathy(:,:) ) 1640 ENDIF 1641 ! if not compatible after all check (ie U point water column less than 2 cells), mask U 1642 DO jj = 1, jpjm1 1643 DO ji = 1, jpim1 1644 IF (misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) >= 1 .AND. mbathy(ji+1,jj) >= 1) THEN 1645 mbathy(ji+1,jj) = mbathy(ji+1,jj) - 1; bathy(ji+1,jj) = gdepw_1d(mbathy(ji+1,jj)+1) ; 1646 END IF 1647 END DO 1648 END DO 1649 IF( lk_mpp ) THEN 1650 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1651 CALL lbc_lnk( zbathy, 'T', 1. ) 1652 misfdep(:,:) = INT( zbathy(:,:) ) 1653 1654 CALL lbc_lnk( risfdep,'T', 1. ) 1655 CALL lbc_lnk( bathy, 'T', 1. ) 1656 1657 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1658 CALL lbc_lnk( zbathy, 'T', 1. ) 1659 mbathy(:,:) = INT( zbathy(:,:) ) 1660 ENDIF 1661 ! if not compatible after all check (ie V point water column less than 2 cells), mask V 1662 DO jj = 1, jpjm1 1663 DO ji = 1, jpi 1664 IF (mbathy(ji,jj) == misfdep(ji,jj+1) .AND. mbathy(ji,jj) >= 1 .AND. mbathy(ji,jj+1) >= 1) THEN 1665 mbathy(ji,jj) = mbathy(ji,jj) - 1 ; bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)+1) ; 1666 END IF 1667 END DO 1668 END DO 1669 IF( lk_mpp ) THEN 1670 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1671 CALL lbc_lnk( zbathy, 'T', 1. ) 1672 misfdep(:,:) = INT( zbathy(:,:) ) 1673 1674 CALL lbc_lnk( risfdep,'T', 1. ) 1675 CALL lbc_lnk( bathy, 'T', 1. ) 1676 1677 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1678 CALL lbc_lnk( zbathy, 'T', 1. ) 1679 mbathy(:,:) = INT( zbathy(:,:) ) 1680 ENDIF 1681 ! if not compatible after all check (ie V point water column less than 2 cells), mask V 1682 DO jj = 1, jpjm1 1683 DO ji = 1, jpi 1684 IF (misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) >= 1 .AND. mbathy(ji,jj+1) >= 1) THEN 1685 mbathy(ji,jj+1) = mbathy(ji,jj+1) - 1 ; bathy(ji,jj+1) = gdepw_1d(mbathy(ji,jj+1)+1) ; 1686 END IF 1687 END DO 1688 END DO 1689 IF( lk_mpp ) THEN 1690 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1691 CALL lbc_lnk( zbathy, 'T', 1. ) 1692 misfdep(:,:) = INT( zbathy(:,:) ) 1693 1694 CALL lbc_lnk( risfdep,'T', 1. ) 1695 CALL lbc_lnk( bathy, 'T', 1. ) 1696 1697 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1698 CALL lbc_lnk( zbathy, 'T', 1. ) 1699 mbathy(:,:) = INT( zbathy(:,:) ) 1700 ENDIF 1701 ! if not compatible after all check, mask T 1702 DO jj = 1, jpj 1703 DO ji = 1, jpi 1704 IF (mbathy(ji,jj) <= misfdep(ji,jj)) THEN 1705 misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0._wp ; mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0._wp ; 1706 END IF 1707 END DO 1708 END DO 1709 1710 WHERE (mbathy(:,:) == 1) 1711 mbathy = 0; bathy = 0.0_wp ; misfdep = 0 ; risfdep = 0.0_wp 1712 END WHERE 1713 END DO 1714 ! end check compatibility ice shelf/bathy 1715 ! remove very shallow ice shelf (less than ~ 10m if 75L) 1716 WHERE (risfdep(:,:) <= 10._wp) 1717 misfdep = 1; risfdep = 0.0_wp; 1718 END WHERE 1719 1720 IF( icompt == 0 ) THEN 1721 IF(lwp) WRITE(numout,*)' no points with ice shelf too close to bathymetry' 1722 ELSE 1723 IF(lwp) WRITE(numout,*)' ',icompt,' ocean grid points with ice shelf thickness reduced to avoid bathymetry' 1724 ENDIF 1725 1726 ! compute scale factor and depth at T- and W- points 1727 DO jj = 1, jpj 1728 DO ji = 1, jpi 1729 ik = mbathy(ji,jj) 1730 IF( ik > 0 ) THEN ! ocean point only 1731 ! max ocean level case 1732 IF( ik == jpkm1 ) THEN 1733 zdepwp = bathy(ji,jj) 1734 ze3tp = bathy(ji,jj) - gdepw_1d(ik) 1735 ze3wp = 0.5_wp * e3w_1d(ik) * ( 1._wp + ( ze3tp/e3t_1d(ik) ) ) 1736 e3t_0(ji,jj,ik ) = ze3tp 1737 e3t_0(ji,jj,ik+1) = ze3tp 1738 e3w_0(ji,jj,ik ) = ze3wp 1739 e3w_0(ji,jj,ik+1) = ze3tp 1740 gdepw_0(ji,jj,ik+1) = zdepwp 1741 gdept_0(ji,jj,ik ) = gdept_1d(ik-1) + ze3wp 1742 gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + ze3tp 1743 ! 1744 ELSE ! standard case 1745 IF( bathy(ji,jj) <= gdepw_1d(ik+1) ) THEN ; gdepw_0(ji,jj,ik+1) = bathy(ji,jj) 1746 ELSE ; gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1) 1747 ENDIF 1748 ! gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1) 1749 !gm Bug? check the gdepw_1d 1750 ! ... on ik 1751 gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) ) & 1752 & * ((gdept_1d( ik ) - gdepw_1d(ik) ) & 1753 & / ( gdepw_1d( ik+1) - gdepw_1d(ik) )) 1754 e3t_0 (ji,jj,ik ) = gdepw_0(ji,jj,ik+1) - gdepw_1d(ik ) 1755 e3w_0 (ji,jj,ik ) = gdept_0(ji,jj,ik ) - gdept_1d(ik-1) 1756 ! ... on ik+1 1757 e3w_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik) 1758 e3t_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik) 1759 ENDIF 1760 ENDIF 1761 END DO 1762 END DO 1763 ! 1764 it = 0 1765 DO jj = 1, jpj 1766 DO ji = 1, jpi 1767 ik = mbathy(ji,jj) 1768 IF( ik > 0 ) THEN ! ocean point only 1769 e3tp (ji,jj) = e3t_0(ji,jj,ik) 1770 e3wp (ji,jj) = e3w_0(ji,jj,ik) 1771 ! test 1772 zdiff= gdepw_0(ji,jj,ik+1) - gdept_0(ji,jj,ik ) 1773 IF( zdiff <= 0._wp .AND. lwp ) THEN 1774 it = it + 1 1775 WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj 1776 WRITE(numout,*) ' bathy = ', bathy(ji,jj) 1777 WRITE(numout,*) ' gdept_0 = ', gdept_0(ji,jj,ik), ' gdepw_0 = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff 1778 WRITE(numout,*) ' e3tp = ', e3t_0 (ji,jj,ik), ' e3wp = ', e3w_0 (ji,jj,ik ) 1779 ENDIF 1780 ENDIF 1781 END DO 1782 END DO 1783 ! 1784 ! (ISF) Definition of e3t, u, v, w for ISF case 1785 DO jj = 1, jpj 1786 DO ji = 1, jpi 1787 ik = misfdep(ji,jj) 1788 IF( ik > 1 ) THEN ! ice shelf point only 1789 IF( risfdep(ji,jj) < gdepw_1d(ik) ) risfdep(ji,jj)= gdepw_1d(ik) 1790 gdepw_0(ji,jj,ik) = risfdep(ji,jj) 1791 !gm Bug? check the gdepw_0 1792 ! ... on ik 1793 gdept_0(ji,jj,ik) = gdepw_1d(ik+1) - ( gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) ) & 1794 & * ( gdepw_1d(ik+1) - gdept_1d(ik) ) & 1795 & / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) 1796 e3t_0 (ji,jj,ik ) = gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) 1797 e3w_0 (ji,jj,ik+1) = gdept_1d(ik+1) - gdept_0(ji,jj,ik) 1798 1799 IF( ik + 1 == mbathy(ji,jj) ) THEN ! ice shelf point only (2 cell water column) 1800 e3w_0 (ji,jj,ik+1) = gdept_0(ji,jj,ik+1) - gdept_0(ji,jj,ik) 1801 ENDIF 1802 ! ... on ik / ik-1 1803 e3w_0 (ji,jj,ik ) = e3t_0 (ji,jj,ik) !2._wp * (gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik)) 1804 e3t_0 (ji,jj,ik-1) = gdepw_0(ji,jj,ik) - gdepw_1d(ik-1) 1805 ! The next line isn't required and doesn't affect results - included for consistency with bathymetry code 1806 gdept_0(ji,jj,ik-1) = gdept_1d(ik-1) 1807 ENDIF 1808 END DO 1809 END DO 1810 1811 it = 0 1812 DO jj = 1, jpj 1813 DO ji = 1, jpi 1814 ik = misfdep(ji,jj) 1815 IF( ik > 1 ) THEN ! ice shelf point only 1816 e3tp (ji,jj) = e3t_0(ji,jj,ik ) 1817 e3wp (ji,jj) = e3w_0(ji,jj,ik+1 ) 1818 ! test 1819 zdiff= gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik ) 1820 IF( zdiff <= 0. .AND. lwp ) THEN 1821 it = it + 1 1822 WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj 1823 WRITE(numout,*) ' risfdep = ', risfdep(ji,jj) 1824 WRITE(numout,*) ' gdept = ', gdept_0(ji,jj,ik), ' gdepw = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff 1825 WRITE(numout,*) ' e3tp = ', e3tp(ji,jj), ' e3wp = ', e3wp(ji,jj) 1826 ENDIF 1827 ENDIF 1828 END DO 1829 END DO 1830 1831 CALL wrk_dealloc( jpi, jpj, zmask, zbathy, zrisfdep ) 1832 CALL wrk_dealloc( jpi, jpj, zmisfdep, zmbathy ) 1833 ! 1834 IF( nn_timing == 1 ) CALL timing_stop('zgr_isf') 1835 ! 1836 END SUBROUTINE zgr_isf 1837 1838 1839 SUBROUTINE zgr_sco 1840 !!---------------------------------------------------------------------- 1841 !! *** ROUTINE zgr_sco *** 1842 !! 1843 !! ** Purpose : define the s-coordinate system 1844 !! 1845 !! ** Method : s-coordinate 1846 !! The depth of model levels is defined as the product of an 1847 !! analytical function by the local bathymetry, while the vertical 1848 !! scale factors are defined as the product of the first derivative 1849 !! of the analytical function by the bathymetry. 1850 !! (this solution save memory as depth and scale factors are not 1851 !! 3d fields) 1852 !! - Read bathymetry (in meters) at t-point and compute the 1853 !! bathymetry at u-, v-, and f-points. 1854 !! hbatu = mi( hbatt ) 1855 !! hbatv = mj( hbatt ) 1856 !! hbatf = mi( mj( hbatt ) ) 1857 !! - Compute z_gsigt, z_gsigw, z_esigt, z_esigw from an analytical 1858 !! function and its derivative given as function. 1859 !! z_gsigt(k) = fssig (k ) 1860 !! z_gsigw(k) = fssig (k-0.5) 1861 !! z_esigt(k) = fsdsig(k ) 1862 !! z_esigw(k) = fsdsig(k-0.5) 1863 !! Three options for stretching are give, and they can be modified 1864 !! following the users requirements. Nevertheless, the output as 1865 !! well as the way to compute the model levels and scale factors 1866 !! must be respected in order to insure second order accuracy 1867 !! schemes. 1868 !! 1869 !! The three methods for stretching available are: 1870 !! 1871 !! s_sh94 (Song and Haidvogel 1994) 1872 !! a sinh/tanh function that allows sigma and stretched sigma 1873 !! 1874 !! s_sf12 (Siddorn and Furner 2012?) 1875 !! allows the maintenance of fixed surface and or 1876 !! bottom cell resolutions (cf. geopotential coordinates) 1877 !! within an analytically derived stretched S-coordinate framework. 1878 !! 1879 !! s_tanh (Madec et al 1996) 1880 !! a cosh/tanh function that gives stretched coordinates 1881 !! 1882 !!---------------------------------------------------------------------- 1883 INTEGER :: ji, jj, jk, jl ! dummy loop argument 1884 INTEGER :: iip1, ijp1, iim1, ijm1 ! temporary integers 1885 INTEGER :: ios ! Local integer output status for namelist read 1886 REAL(wp) :: zrmax, ztaper ! temporary scalars 1887 REAL(wp) :: zrfact 1888 ! 1889 REAL(wp), POINTER, DIMENSION(:,: ) :: ztmpi1, ztmpi2, ztmpj1, ztmpj2 1890 REAL(wp), POINTER, DIMENSION(:,: ) :: zenv, ztmp, zmsk, zri, zrj, zhbat 1891 !! 1892 NAMELIST/namzgr_sco/ln_s_sh94, ln_s_sf12, ln_sigcrit, rn_sbot_min, rn_sbot_max, rn_hc, rn_rmax,rn_theta, & 1893 & rn_thetb, rn_bb, rn_alpha, rn_efold, rn_zs, rn_zb_a, rn_zb_b 1894 !!---------------------------------------------------------------------- 1895 ! 1896 IF( nn_timing == 1 ) CALL timing_start('zgr_sco') 1897 ! 1898 CALL wrk_alloc( jpi,jpj, zenv, ztmp, zmsk, zri, zrj, zhbat , ztmpi1, ztmpi2, ztmpj1, ztmpj2 ) 1899 ! 1900 REWIND( numnam_ref ) ! Namelist namzgr_sco in reference namelist : Sigma-stretching parameters 1901 READ ( numnam_ref, namzgr_sco, IOSTAT = ios, ERR = 901) 1902 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr_sco in reference namelist', lwp ) 1903 1904 REWIND( numnam_cfg ) ! Namelist namzgr_sco in configuration namelist : Sigma-stretching parameters 1905 READ ( numnam_cfg, namzgr_sco, IOSTAT = ios, ERR = 902 ) 1906 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr_sco in configuration namelist', lwp ) 1907 IF(lwm) WRITE ( numond, namzgr_sco ) 1908 1909 IF(lwp) THEN ! control print 1910 WRITE(numout,*) 1911 WRITE(numout,*) 'domzgr_sco : s-coordinate or hybrid z-s-coordinate' 1912 WRITE(numout,*) '~~~~~~~~~~~' 1913 WRITE(numout,*) ' Namelist namzgr_sco' 1914 WRITE(numout,*) ' stretching coeffs ' 1915 WRITE(numout,*) ' maximum depth of s-bottom surface (>0) rn_sbot_max = ',rn_sbot_max 1916 WRITE(numout,*) ' minimum depth of s-bottom surface (>0) rn_sbot_min = ',rn_sbot_min 1917 WRITE(numout,*) ' Critical depth rn_hc = ',rn_hc 1918 WRITE(numout,*) ' maximum cut-off r-value allowed rn_rmax = ',rn_rmax 1919 WRITE(numout,*) ' Song and Haidvogel 1994 stretching ln_s_sh94 = ',ln_s_sh94 1920 WRITE(numout,*) ' Song and Haidvogel 1994 stretching coefficients' 1921 WRITE(numout,*) ' surface control parameter (0<=rn_theta<=20) rn_theta = ',rn_theta 1922 WRITE(numout,*) ' bottom control parameter (0<=rn_thetb<= 1) rn_thetb = ',rn_thetb 1923 WRITE(numout,*) ' stretching parameter (song and haidvogel) rn_bb = ',rn_bb 1924 WRITE(numout,*) ' Siddorn and Furner 2012 stretching ln_s_sf12 = ',ln_s_sf12 1925 WRITE(numout,*) ' switching to sigma (T) or Z (F) at H<Hc ln_sigcrit = ',ln_sigcrit 1926 WRITE(numout,*) ' Siddorn and Furner 2012 stretching coefficients' 1927 WRITE(numout,*) ' stretchin parameter ( >1 surface; <1 bottom) rn_alpha = ',rn_alpha 1928 WRITE(numout,*) ' e-fold length scale for transition region rn_efold = ',rn_efold 1929 WRITE(numout,*) ' Surface cell depth (Zs) (m) rn_zs = ',rn_zs 1930 WRITE(numout,*) ' Bathymetry multiplier for Zb rn_zb_a = ',rn_zb_a 1931 WRITE(numout,*) ' Offset for Zb rn_zb_b = ',rn_zb_b 1932 WRITE(numout,*) ' Bottom cell (Zb) (m) = H*rn_zb_a + rn_zb_b' 1933 ENDIF 1934 1935 hift(:,:) = rn_sbot_min ! set the minimum depth for the s-coordinate 1936 hifu(:,:) = rn_sbot_min 1937 hifv(:,:) = rn_sbot_min 1938 hiff(:,:) = rn_sbot_min 1939 1940 ! ! set maximum ocean depth 1941 bathy(:,:) = MIN( rn_sbot_max, bathy(:,:) ) 1942 1943 IF( .NOT.ln_wd ) THEN 1944 DO jj = 1, jpj 1945 DO ji = 1, jpi 1946 IF( bathy(ji,jj) > 0._wp ) bathy(ji,jj) = MAX( rn_sbot_min, bathy(ji,jj) ) 1947 END DO 1948 END DO 1949 END IF 1950 ! ! ============================= 1951 ! ! Define the envelop bathymetry (hbatt) 1952 ! ! ============================= 1953 ! use r-value to create hybrid coordinates 1954 zenv(:,:) = bathy(:,:) 1955 ! 1956 IF( .NOT.ln_wd ) THEN 1957 ! set first land point adjacent to a wet cell to sbot_min as this needs to be included in smoothing 1958 DO jj = 1, jpj 1959 DO ji = 1, jpi 1960 IF( bathy(ji,jj) == 0._wp ) THEN 1961 iip1 = MIN( ji+1, jpi ) 1962 ijp1 = MIN( jj+1, jpj ) 1963 iim1 = MAX( ji-1, 1 ) 1964 ijm1 = MAX( jj-1, 1 ) 1965 !!gm BUG fix see ticket #1617 1966 IF( ( + bathy(iim1,ijm1) + bathy(ji,ijp1) + bathy(iip1,ijp1) & 1967 & + bathy(iim1,jj ) + bathy(iip1,jj ) & 1968 & + bathy(iim1,ijm1) + bathy(ji,ijm1) + bathy(iip1,ijp1) ) > 0._wp ) & 1969 & zenv(ji,jj) = rn_sbot_min 1970 !!gm 1971 !!gm IF( ( bathy(iip1,jj ) + bathy(iim1,jj ) + bathy(ji,ijp1 ) + bathy(ji,ijm1) + & 1972 !!gm & bathy(iip1,ijp1) + bathy(iim1,ijm1) + bathy(iip1,ijp1) + bathy(iim1,ijm1)) > 0._wp ) THEN 1973 !!gm zenv(ji,jj) = rn_sbot_min 1974 !!gm ENDIF 1975 !!gm end 1976 ENDIF 1977 END DO 1978 END DO 1979 END IF 1980 1981 ! apply lateral boundary condition CAUTION: keep the value when the lbc field is zero 1982 CALL lbc_lnk( zenv, 'T', 1._wp, 'no0' ) 1983 ! 1984 ! smooth the bathymetry (if required) 1985 scosrf(:,:) = 0._wp ! ocean surface depth (here zero: no under ice-shelf sea) 1986 scobot(:,:) = bathy(:,:) ! ocean bottom depth 1987 ! 1988 jl = 0 1989 zrmax = 1._wp 1990 ! 1991 ! 1992 ! set scaling factor used in reducing vertical gradients 1993 zrfact = ( 1._wp - rn_rmax ) / ( 1._wp + rn_rmax ) 1994 ! 1995 ! initialise temporary evelope depth arrays 1996 ztmpi1(:,:) = zenv(:,:) 1997 ztmpi2(:,:) = zenv(:,:) 1998 ztmpj1(:,:) = zenv(:,:) 1999 ztmpj2(:,:) = zenv(:,:) 2000 ! 2001 ! initialise temporary r-value arrays 2002 zri(:,:) = 1._wp 2003 zrj(:,:) = 1._wp 2004 ! ! ================ ! 2005 DO WHILE( jl <= 10000 .AND. ( zrmax - rn_rmax ) > 1.e-8_wp ) ! Iterative loop ! 2006 ! ! ================ ! 2007 jl = jl + 1 2008 zrmax = 0._wp 2009 ! we set zrmax from previous r-values (zri and zrj) first 2010 ! if set after current r-value calculation (as previously) 2011 ! we could exit DO WHILE prematurely before checking r-value 2012 ! of current zenv 2013 DO jj = 1, nlcj 2014 DO ji = 1, nlci 2015 zrmax = MAX( zrmax, ABS(zri(ji,jj)), ABS(zrj(ji,jj)) ) 2016 END DO 2017 END DO 2018 zri(:,:) = 0._wp 2019 zrj(:,:) = 0._wp 2020 DO jj = 1, nlcj 2021 DO ji = 1, nlci 2022 iip1 = MIN( ji+1, nlci ) ! force zri = 0 on last line (ji=ncli+1 to jpi) 2023 ijp1 = MIN( jj+1, nlcj ) ! force zrj = 0 on last raw (jj=nclj+1 to jpj) 2024 IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(iip1,jj) > 0._wp)) THEN 2025 zri(ji,jj) = ( zenv(iip1,jj ) - zenv(ji,jj) ) / ( zenv(iip1,jj ) + zenv(ji,jj) ) 2026 END IF 2027 IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(ji,ijp1) > 0._wp)) THEN 2028 zrj(ji,jj) = ( zenv(ji ,ijp1) - zenv(ji,jj) ) / ( zenv(ji ,ijp1) + zenv(ji,jj) ) 2029 END IF 2030 IF( zri(ji,jj) > rn_rmax ) ztmpi1(ji ,jj ) = zenv(iip1,jj ) * zrfact 2031 IF( zri(ji,jj) < -rn_rmax ) ztmpi2(iip1,jj ) = zenv(ji ,jj ) * zrfact 2032 IF( zrj(ji,jj) > rn_rmax ) ztmpj1(ji ,jj ) = zenv(ji ,ijp1) * zrfact 2033 IF( zrj(ji,jj) < -rn_rmax ) ztmpj2(ji ,ijp1) = zenv(ji ,jj ) * zrfact 2034 END DO 2035 END DO 2036 IF( lk_mpp ) CALL mpp_max( zrmax ) ! max over the global domain 2037 ! 2038 IF(lwp)WRITE(numout,*) 'zgr_sco : iter= ',jl, ' rmax= ', zrmax 2039 ! 2040 DO jj = 1, nlcj 2041 DO ji = 1, nlci 2042 zenv(ji,jj) = MAX(zenv(ji,jj), ztmpi1(ji,jj), ztmpi2(ji,jj), ztmpj1(ji,jj), ztmpj2(ji,jj) ) 2043 END DO 2044 END DO 2045 ! apply lateral boundary condition CAUTION: keep the value when the lbc field is zero 2046 CALL lbc_lnk( zenv, 'T', 1._wp, 'no0' ) 2047 ! ! ================ ! 2048 END DO ! End loop ! 2049 ! ! ================ ! 2050 DO jj = 1, jpj 2051 DO ji = 1, jpi 2052 zenv(ji,jj) = MAX( zenv(ji,jj), rn_sbot_min ) ! set all points to avoid undefined scale value warnings 2053 END DO 2054 END DO 2055 ! 2056 ! Envelope bathymetry saved in hbatt 2057 hbatt(:,:) = zenv(:,:) 2058 IF( MINVAL( gphit(:,:) ) * MAXVAL( gphit(:,:) ) <= 0._wp ) THEN 2059 CALL ctl_warn( ' s-coordinates are tapered in vicinity of the Equator' ) 2060 DO jj = 1, jpj 2061 DO ji = 1, jpi 2062 ztaper = EXP( -(gphit(ji,jj)/8._wp)**2._wp ) 2063 hbatt(ji,jj) = rn_sbot_max * ztaper + hbatt(ji,jj) * ( 1._wp - ztaper ) 2064 END DO 2065 END DO 2066 ENDIF 2067 ! 2068 ! ! ============================== 2069 ! ! hbatu, hbatv, hbatf fields 2070 ! ! ============================== 2071 IF(lwp) THEN 2072 WRITE(numout,*) 2073 IF( .NOT.ln_wd ) THEN 2074 WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', rn_sbot_min 2075 ELSE 2076 WRITE(numout,*) ' zgr_sco: minimum positive depth of the envelop topography set to : ', rn_sbot_min 2077 WRITE(numout,*) ' zgr_sco: minimum negative depth of the envelop topography set to : ', -rn_wdld 2078 ENDIF 2079 ENDIF 2080 hbatu(:,:) = rn_sbot_min 2081 hbatv(:,:) = rn_sbot_min 2082 hbatf(:,:) = rn_sbot_min 2083 DO jj = 1, jpjm1 2084 DO ji = 1, jpim1 ! NO vector opt. 2085 hbatu(ji,jj) = 0.50_wp * ( hbatt(ji ,jj) + hbatt(ji+1,jj ) ) 2086 hbatv(ji,jj) = 0.50_wp * ( hbatt(ji ,jj) + hbatt(ji ,jj+1) ) 2087 hbatf(ji,jj) = 0.25_wp * ( hbatt(ji ,jj) + hbatt(ji ,jj+1) & 2088 & + hbatt(ji+1,jj) + hbatt(ji+1,jj+1) ) 2089 END DO 2090 END DO 2091 2092 IF( ln_wd ) THEN !avoid the zero depth on T- (U-,V-,F-) points 2093 DO jj = 1, jpj 2094 DO ji = 1, jpi 2095 IF(ABS(hbatt(ji,jj)) < rn_wdmin1) & 2096 & hbatt(ji,jj) = SIGN(1._wp, hbatt(ji,jj)) * rn_wdmin1 2097 IF(ABS(hbatu(ji,jj)) < rn_wdmin1) & 2098 & hbatu(ji,jj) = SIGN(1._wp, hbatu(ji,jj)) * rn_wdmin1 2099 IF(ABS(hbatv(ji,jj)) < rn_wdmin1) & 2100 & hbatv(ji,jj) = SIGN(1._wp, hbatv(ji,jj)) * rn_wdmin1 2101 IF(ABS(hbatf(ji,jj)) < rn_wdmin1) & 2102 & hbatf(ji,jj) = SIGN(1._wp, hbatf(ji,jj)) * rn_wdmin1 2103 END DO 2104 END DO 2105 END IF 2106 ! 2107 ! Apply lateral boundary condition 2108 !!gm ! CAUTION: retain non zero value in the initial file this should be OK for orca cfg, not for EEL 2109 zhbat(:,:) = hbatu(:,:) ; CALL lbc_lnk( hbatu, 'U', 1._wp ) 2110 DO jj = 1, jpj 2111 DO ji = 1, jpi 2112 IF( hbatu(ji,jj) == 0._wp ) THEN 2113 !No worries about the following line when ln_wd == .true. 2114 IF( zhbat(ji,jj) == 0._wp ) hbatu(ji,jj) = rn_sbot_min 2115 IF( zhbat(ji,jj) /= 0._wp ) hbatu(ji,jj) = zhbat(ji,jj) 2116 ENDIF 2117 END DO 2118 END DO 2119 zhbat(:,:) = hbatv(:,:) ; CALL lbc_lnk( hbatv, 'V', 1._wp ) 2120 DO jj = 1, jpj 2121 DO ji = 1, jpi 2122 IF( hbatv(ji,jj) == 0._wp ) THEN 2123 IF( zhbat(ji,jj) == 0._wp ) hbatv(ji,jj) = rn_sbot_min 2124 IF( zhbat(ji,jj) /= 0._wp ) hbatv(ji,jj) = zhbat(ji,jj) 2125 ENDIF 2126 END DO 2127 END DO 2128 zhbat(:,:) = hbatf(:,:) ; CALL lbc_lnk( hbatf, 'F', 1._wp ) 2129 DO jj = 1, jpj 2130 DO ji = 1, jpi 2131 IF( hbatf(ji,jj) == 0._wp ) THEN 2132 IF( zhbat(ji,jj) == 0._wp ) hbatf(ji,jj) = rn_sbot_min 2133 IF( zhbat(ji,jj) /= 0._wp ) hbatf(ji,jj) = zhbat(ji,jj) 2134 ENDIF 2135 END DO 2136 END DO 2137 2138 !!bug: key_helsinki a verifer 2139 IF( .NOT.ln_wd ) THEN 2140 hift(:,:) = MIN( hift(:,:), hbatt(:,:) ) 2141 hifu(:,:) = MIN( hifu(:,:), hbatu(:,:) ) 2142 hifv(:,:) = MIN( hifv(:,:), hbatv(:,:) ) 2143 hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) ) 2144 END IF 2145 2146 IF( nprint == 1 .AND. lwp ) THEN 2147 WRITE(numout,*) ' MAX val hif t ', MAXVAL( hift (:,:) ), ' f ', MAXVAL( hiff (:,:) ), & 2148 & ' u ', MAXVAL( hifu (:,:) ), ' v ', MAXVAL( hifv (:,:) ) 2149 WRITE(numout,*) ' MIN val hif t ', MINVAL( hift (:,:) ), ' f ', MINVAL( hiff (:,:) ), & 2150 & ' u ', MINVAL( hifu (:,:) ), ' v ', MINVAL( hifv (:,:) ) 2151 WRITE(numout,*) ' MAX val hbat t ', MAXVAL( hbatt(:,:) ), ' f ', MAXVAL( hbatf(:,:) ), & 2152 & ' u ', MAXVAL( hbatu(:,:) ), ' v ', MAXVAL( hbatv(:,:) ) 2153 WRITE(numout,*) ' MIN val hbat t ', MINVAL( hbatt(:,:) ), ' f ', MINVAL( hbatf(:,:) ), & 2154 & ' u ', MINVAL( hbatu(:,:) ), ' v ', MINVAL( hbatv(:,:) ) 2155 ENDIF 2156 !! helsinki 2157 2158 ! ! ======================= 2159 ! ! s-ccordinate fields (gdep., e3.) 2160 ! ! ======================= 2161 ! 2162 ! non-dimensional "sigma" for model level depth at w- and t-levels 2163 2164 2165 !======================================================================== 2166 ! Song and Haidvogel 1994 (ln_s_sh94=T) 2167 ! Siddorn and Furner 2012 (ln_sf12=T) 2168 ! or tanh function (both false) 2169 !======================================================================== 2170 IF ( ln_s_sh94 ) THEN 2171 CALL s_sh94() 2172 ELSE IF ( ln_s_sf12 ) THEN 2173 CALL s_sf12() 2174 ELSE 2175 CALL s_tanh() 2176 ENDIF 2177 2178 CALL lbc_lnk( e3t_0 , 'T', 1._wp ) 2179 CALL lbc_lnk( e3u_0 , 'U', 1._wp ) 2180 CALL lbc_lnk( e3v_0 , 'V', 1._wp ) 2181 CALL lbc_lnk( e3f_0 , 'F', 1._wp ) 2182 CALL lbc_lnk( e3w_0 , 'W', 1._wp ) 2183 CALL lbc_lnk( e3uw_0, 'U', 1._wp ) 2184 CALL lbc_lnk( e3vw_0, 'V', 1._wp ) 2185 ! 2186 IF( .NOT.ln_wd ) THEN 2187 WHERE( e3t_0 (:,:,:) == 0._wp ) e3t_0 (:,:,:) = 1._wp 2188 WHERE( e3u_0 (:,:,:) == 0._wp ) e3u_0 (:,:,:) = 1._wp 2189 WHERE( e3v_0 (:,:,:) == 0._wp ) e3v_0 (:,:,:) = 1._wp 2190 WHERE( e3f_0 (:,:,:) == 0._wp ) e3f_0 (:,:,:) = 1._wp 2191 WHERE( e3w_0 (:,:,:) == 0._wp ) e3w_0 (:,:,:) = 1._wp 2192 WHERE( e3uw_0(:,:,:) == 0._wp ) e3uw_0(:,:,:) = 1._wp 2193 WHERE( e3vw_0(:,:,:) == 0._wp ) e3vw_0(:,:,:) = 1._wp 2194 END IF 2195 2196 #if defined key_agrif 2197 IF( .NOT. Agrif_Root() ) THEN ! Ensure meaningful vertical scale factors in ghost lines/columns 2198 IF( nbondi == -1 .OR. nbondi == 2 ) e3u_0( 1 , : ,:) = e3u_0( 2 , : ,:) 2199 IF( nbondi == 1 .OR. nbondi == 2 ) e3u_0(nlci-1, : ,:) = e3u_0(nlci-2, : ,:) 2200 IF( nbondj == -1 .OR. nbondj == 2 ) e3v_0( : , 1 ,:) = e3v_0( : , 2 ,:) 2201 IF( nbondj == 1 .OR. nbondj == 2 ) e3v_0( : ,nlcj-1,:) = e3v_0( : ,nlcj-2,:) 2202 ENDIF 2203 #endif 2204 2205 !!gm I don't like that HERE we are supposed to set the reference coordinate (i.e. _0 arrays) 2206 !!gm and only that !!!!! 2207 !!gm THIS should be removed from here ! 2208 gdept_n(:,:,:) = gdept_0(:,:,:) 2209 gdepw_n(:,:,:) = gdepw_0(:,:,:) 2210 gde3w_n(:,:,:) = gde3w_0(:,:,:) 2211 e3t_n (:,:,:) = e3t_0 (:,:,:) 2212 e3u_n (:,:,:) = e3u_0 (:,:,:) 2213 e3v_n (:,:,:) = e3v_0 (:,:,:) 2214 e3f_n (:,:,:) = e3f_0 (:,:,:) 2215 e3w_n (:,:,:) = e3w_0 (:,:,:) 2216 e3uw_n (:,:,:) = e3uw_0 (:,:,:) 2217 e3vw_n (:,:,:) = e3vw_0 (:,:,:) 2218 !!gm and obviously in the following, use the _0 arrays until the end of this subroutine 2219 !! gm end 2220 !! 2221 ! HYBRID : 2222 DO jj = 1, jpj 2223 DO ji = 1, jpi 2224 DO jk = 1, jpkm1 2225 IF( scobot(ji,jj) >= gdept_n(ji,jj,jk) ) mbathy(ji,jj) = MAX( 2, jk ) 2226 END DO 2227 IF( ln_wd ) THEN 2228 IF( scobot(ji,jj) <= -(rn_wdld - rn_wdmin2)) THEN 2229 mbathy(ji,jj) = 0 2230 ELSEIF(scobot(ji,jj) <= rn_wdmin1) THEN 2231 mbathy(ji,jj) = 2 2232 ENDIF 2233 ELSE 2234 IF( scobot(ji,jj) == 0._wp ) mbathy(ji,jj) = 0 2235 ENDIF 2236 END DO 2237 END DO 2238 IF( nprint == 1 .AND. lwp ) WRITE(numout,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ), & 2239 & ' MAX ', MAXVAL( mbathy(:,:) ) 2240 2241 IF( nprint == 1 .AND. lwp ) THEN ! min max values over the local domain 2242 WRITE(numout,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) ) 2243 WRITE(numout,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ), & 2244 & ' w ', MINVAL( gdepw_0(:,:,:) ), '3w ' , MINVAL( gde3w_0(:,:,:) ) 2245 WRITE(numout,*) ' MIN val e3 t ', MINVAL( e3t_0 (:,:,:) ), ' f ' , MINVAL( e3f_0 (:,:,:) ), & 2246 & ' u ', MINVAL( e3u_0 (:,:,:) ), ' u ' , MINVAL( e3v_0 (:,:,:) ), & 2247 & ' uw', MINVAL( e3uw_0 (:,:,:) ), ' vw' , MINVAL( e3vw_0 (:,:,:) ), & 2248 & ' w ', MINVAL( e3w_0 (:,:,:) ) 2249 2250 WRITE(numout,*) ' MAX val depth t ', MAXVAL( gdept_0(:,:,:) ), & 2251 & ' w ', MAXVAL( gdepw_0(:,:,:) ), '3w ' , MAXVAL( gde3w_0(:,:,:) ) 2252 WRITE(numout,*) ' MAX val e3 t ', MAXVAL( e3t_0 (:,:,:) ), ' f ' , MAXVAL( e3f_0 (:,:,:) ), & 2253 & ' u ', MAXVAL( e3u_0 (:,:,:) ), ' u ' , MAXVAL( e3v_0 (:,:,:) ), & 2254 & ' uw', MAXVAL( e3uw_0 (:,:,:) ), ' vw' , MAXVAL( e3vw_0 (:,:,:) ), & 2255 & ' w ', MAXVAL( e3w_0 (:,:,:) ) 2256 ENDIF 2257 ! END DO 2258 IF(lwp) THEN ! selected vertical profiles 2259 WRITE(numout,*) 2260 WRITE(numout,*) ' domzgr: vertical coordinates : point (1,1,k) bathy = ', bathy(1,1), hbatt(1,1) 2261 WRITE(numout,*) ' ~~~~~~ --------------------' 2262 WRITE(numout,"(9x,' level gdept_0 gdepw_0 e3t_0 e3w_0')") 2263 WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(1,1,jk), gdepw_0(1,1,jk), & 2264 & e3t_0 (1,1,jk) , e3w_0 (1,1,jk) , jk=1,jpk ) 2265 DO jj = mj0(20), mj1(20) 2266 DO ji = mi0(20), mi1(20) 2267 WRITE(numout,*) 2268 WRITE(numout,*) ' domzgr: vertical coordinates : point (20,20,k) bathy = ', bathy(ji,jj), hbatt(ji,jj) 2269 WRITE(numout,*) ' ~~~~~~ --------------------' 2270 WRITE(numout,"(9x,' level gdept_0 gdepw_0 e3t_0 e3w_0')") 2271 WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(ji,jj,jk), gdepw_0(ji,jj,jk), & 2272 & e3t_0 (ji,jj,jk) , e3w_0 (ji,jj,jk) , jk=1,jpk ) 2273 END DO 2274 END DO 2275 DO jj = mj0(74), mj1(74) 2276 DO ji = mi0(100), mi1(100) 2277 WRITE(numout,*) 2278 WRITE(numout,*) ' domzgr: vertical coordinates : point (100,74,k) bathy = ', bathy(ji,jj), hbatt(ji,jj) 2279 WRITE(numout,*) ' ~~~~~~ --------------------' 2280 WRITE(numout,"(9x,' level gdept_0 gdepw_0 e3t_0 e3w_0')") 2281 WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(ji,jj,jk), gdepw_0(ji,jj,jk), & 2282 & e3t_0 (ji,jj,jk) , e3w_0 (ji,jj,jk) , jk=1,jpk ) 2283 END DO 2284 END DO 2285 ENDIF 2286 ! 2287 !================================================================================ 2288 ! check the coordinate makes sense 2289 !================================================================================ 2290 DO ji = 1, jpi 2291 DO jj = 1, jpj 2292 ! 2293 IF( hbatt(ji,jj) > 0._wp) THEN 2294 DO jk = 1, mbathy(ji,jj) 2295 ! check coordinate is monotonically increasing 2296 IF (e3w_n(ji,jj,jk) <= 0._wp .OR. e3t_n(ji,jj,jk) <= 0._wp ) THEN 2297 WRITE(ctmp1,*) 'ERROR zgr_sco : e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk 2298 WRITE(numout,*) 'ERROR zgr_sco : e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk 2299 WRITE(numout,*) 'e3w',e3w_n(ji,jj,:) 2300 WRITE(numout,*) 'e3t',e3t_n(ji,jj,:) 2301 CALL ctl_stop( ctmp1 ) 2302 ENDIF 2303 ! and check it has never gone negative 2304 IF( gdepw_n(ji,jj,jk) < 0._wp .OR. gdept_n(ji,jj,jk) < 0._wp ) THEN 2305 WRITE(ctmp1,*) 'ERROR zgr_sco : gdepw or gdept =< 0 at point (i,j,k)= ', ji, jj, jk 2306 WRITE(numout,*) 'ERROR zgr_sco : gdepw or gdept =< 0 at point (i,j,k)= ', ji, jj, jk 2307 WRITE(numout,*) 'gdepw',gdepw_n(ji,jj,:) 2308 WRITE(numout,*) 'gdept',gdept_n(ji,jj,:) 2309 CALL ctl_stop( ctmp1 ) 2310 ENDIF 2311 ! and check it never exceeds the total depth 2312 IF( gdepw_n(ji,jj,jk) > hbatt(ji,jj) ) THEN 2313 WRITE(ctmp1,*) 'ERROR zgr_sco : gdepw > hbatt at point (i,j,k)= ', ji, jj, jk 2314 WRITE(numout,*) 'ERROR zgr_sco : gdepw > hbatt at point (i,j,k)= ', ji, jj, jk 2315 WRITE(numout,*) 'gdepw',gdepw_n(ji,jj,:) 2316 CALL ctl_stop( ctmp1 ) 2317 ENDIF 2318 END DO 2319 ! 2320 DO jk = 1, mbathy(ji,jj)-1 2321 ! and check it never exceeds the total depth 2322 IF( gdept_n(ji,jj,jk) > hbatt(ji,jj) ) THEN 2323 WRITE(ctmp1,*) 'ERROR zgr_sco : gdept > hbatt at point (i,j,k)= ', ji, jj, jk 2324 WRITE(numout,*) 'ERROR zgr_sco : gdept > hbatt at point (i,j,k)= ', ji, jj, jk 2325 WRITE(numout,*) 'gdept',gdept_n(ji,jj,:) 2326 CALL ctl_stop( ctmp1 ) 2327 ENDIF 2328 END DO 2329 ENDIF 2330 END DO 2331 END DO 2332 ! 2333 CALL wrk_dealloc( jpi, jpj, zenv, ztmp, zmsk, zri, zrj, zhbat , ztmpi1, ztmpi2, ztmpj1, ztmpj2 ) 2334 ! 2335 IF( nn_timing == 1 ) CALL timing_stop('zgr_sco') 2336 ! 2337 END SUBROUTINE zgr_sco 2338 2339 2340 SUBROUTINE s_sh94() 2341 !!---------------------------------------------------------------------- 2342 !! *** ROUTINE s_sh94 *** 2343 !! 2344 !! ** Purpose : stretch the s-coordinate system 2345 !! 2346 !! ** Method : s-coordinate stretch using the Song and Haidvogel 1994 2347 !! mixed S/sigma coordinate 2348 !! 2349 !! Reference : Song and Haidvogel 1994. 2350 !!---------------------------------------------------------------------- 2351 INTEGER :: ji, jj, jk ! dummy loop argument 2352 REAL(wp) :: zcoeft, zcoefw ! temporary scalars 2353 REAL(wp) :: ztmpu, ztmpv, ztmpf 2354 REAL(wp) :: ztmpu1, ztmpv1, ztmpf1 2355 ! 2356 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3 2357 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 2358 !!---------------------------------------------------------------------- 2359 2360 CALL wrk_alloc( jpi,jpj,jpk, z_gsigw3, z_gsigt3, z_gsi3w3 ) 2361 CALL wrk_alloc( jpi,jpj,jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 ) 2362 2363 z_gsigw3 = 0._wp ; z_gsigt3 = 0._wp ; z_gsi3w3 = 0._wp 2364 z_esigt3 = 0._wp ; z_esigw3 = 0._wp 2365 z_esigtu3 = 0._wp ; z_esigtv3 = 0._wp ; z_esigtf3 = 0._wp 2366 z_esigwu3 = 0._wp ; z_esigwv3 = 0._wp 2367 ! 2368 DO ji = 1, jpi 2369 DO jj = 1, jpj 2370 ! 2371 IF( hbatt(ji,jj) > rn_hc ) THEN !deep water, stretched sigma 2372 DO jk = 1, jpk 2373 z_gsigw3(ji,jj,jk) = -fssig1( REAL(jk,wp)-0.5_wp, rn_bb ) 2374 z_gsigt3(ji,jj,jk) = -fssig1( REAL(jk,wp) , rn_bb ) 2375 END DO 2376 ELSE ! shallow water, uniform sigma 2377 DO jk = 1, jpk 2378 z_gsigw3(ji,jj,jk) = REAL(jk-1,wp) / REAL(jpk-1,wp) 2379 z_gsigt3(ji,jj,jk) = ( REAL(jk-1,wp) + 0.5_wp ) / REAL(jpk-1,wp) 2380 END DO 2381 ENDIF 2382 ! 2383 DO jk = 1, jpkm1 2384 z_esigt3(ji,jj,jk ) = z_gsigw3(ji,jj,jk+1) - z_gsigw3(ji,jj,jk) 2385 z_esigw3(ji,jj,jk+1) = z_gsigt3(ji,jj,jk+1) - z_gsigt3(ji,jj,jk) 2386 END DO 2387 z_esigw3(ji,jj,1 ) = 2._wp * ( z_gsigt3(ji,jj,1 ) - z_gsigw3(ji,jj,1 ) ) 2388 z_esigt3(ji,jj,jpk) = 2._wp * ( z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk) ) 2389 ! 2390 ! Coefficients for vertical depth as the sum of e3w scale factors 2391 z_gsi3w3(ji,jj,1) = 0.5_wp * z_esigw3(ji,jj,1) 2392 DO jk = 2, jpk 2393 z_gsi3w3(ji,jj,jk) = z_gsi3w3(ji,jj,jk-1) + z_esigw3(ji,jj,jk) 2394 END DO 2395 ! 2396 DO jk = 1, jpk 2397 zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp) 2398 zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp) 2399 gdept_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigt3(ji,jj,jk)+rn_hc*zcoeft ) 2400 gdepw_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigw3(ji,jj,jk)+rn_hc*zcoefw ) 2401 gde3w_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsi3w3(ji,jj,jk)+rn_hc*zcoeft ) 2402 END DO 2403 ! 2404 END DO ! for all jj's 2405 END DO ! for all ji's 2406 2407 DO ji = 1, jpim1 2408 DO jj = 1, jpjm1 2409 ! extended for Wetting/Drying case 2410 ztmpu = hbatt(ji,jj)+hbatt(ji+1,jj) 2411 ztmpv = hbatt(ji,jj)+hbatt(ji,jj+1) 2412 ztmpf = hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) 2413 ztmpu1 = hbatt(ji,jj)*hbatt(ji+1,jj) 2414 ztmpv1 = hbatt(ji,jj)*hbatt(ji,jj+1) 2415 ztmpf1 = MIN(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1)) * & 2416 & MAX(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1)) 2417 DO jk = 1, jpk 2418 IF( ln_wd .AND. (ztmpu1 < 0._wp.OR.ABS(ztmpu) < rn_wdmin1) ) THEN 2419 z_esigtu3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji+1,jj,jk) ) 2420 z_esigwu3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji+1,jj,jk) ) 2421 ELSE 2422 z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) & 2423 & / ztmpu 2424 z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) & 2425 & / ztmpu 2426 END IF 2427 2428 IF( ln_wd .AND. (ztmpv1 < 0._wp.OR.ABS(ztmpv) < rn_wdmin1) ) THEN 2429 z_esigtv3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji,jj+1,jk) ) 2430 z_esigwv3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji,jj+1,jk) ) 2431 ELSE 2432 z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) & 2433 & / ztmpv 2434 z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) & 2435 & / ztmpv 2436 END IF 2437 2438 IF( ln_wd .AND. (ztmpf1 < 0._wp.OR.ABS(ztmpf) < rn_wdmin1) ) THEN 2439 z_esigtf3(ji,jj,jk) = 0.25_wp * ( z_esigt3(ji,jj ,jk) + z_esigt3(ji+1,jj ,jk) & 2440 & + z_esigt3(ji,jj+1,jk) + z_esigt3(ji+1,jj+1,jk) ) 2441 ELSE 2442 z_esigtf3(ji,jj,jk) = ( hbatt(ji ,jj )*z_esigt3(ji ,jj ,jk) & 2443 & + hbatt(ji+1,jj )*z_esigt3(ji+1,jj ,jk) & 2444 & + hbatt(ji ,jj+1)*z_esigt3(ji ,jj+1,jk) & 2445 & + hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) / ztmpf 2446 END IF 2447 2448 ! 2449 e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigt3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 2450 e3u_0(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigtu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 2451 e3v_0(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigtv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 2452 e3f_0(ji,jj,jk) = ( (hbatf(ji,jj)-rn_hc)*z_esigtf3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 2453 ! 2454 e3w_0 (ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigw3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 2455 e3uw_0(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigwu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 2456 e3vw_0(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigwv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 2457 END DO 2458 END DO 2459 END DO 2460 ! 2461 CALL wrk_dealloc( jpi,jpj,jpk, z_gsigw3, z_gsigt3, z_gsi3w3 ) 2462 CALL wrk_dealloc( jpi,jpj,jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 ) 2463 ! 2464 END SUBROUTINE s_sh94 2465 2466 2467 SUBROUTINE s_sf12 2468 !!---------------------------------------------------------------------- 2469 !! *** ROUTINE s_sf12 *** 2470 !! 2471 !! ** Purpose : stretch the s-coordinate system 2472 !! 2473 !! ** Method : s-coordinate stretch using the Siddorn and Furner 2012? 2474 !! mixed S/sigma/Z coordinate 2475 !! 2476 !! This method allows the maintenance of fixed surface and or 2477 !! bottom cell resolutions (cf. geopotential coordinates) 2478 !! within an analytically derived stretched S-coordinate framework. 2479 !! 2480 !! 2481 !! Reference : Siddorn and Furner 2012 (submitted Ocean modelling). 2482 !!---------------------------------------------------------------------- 2483 INTEGER :: ji, jj, jk ! dummy loop argument 2484 REAL(wp) :: zsmth ! smoothing around critical depth 2485 REAL(wp) :: zzs, zzb ! Surface and bottom cell thickness in sigma space 2486 REAL(wp) :: ztmpu, ztmpv, ztmpf 2487 REAL(wp) :: ztmpu1, ztmpv1, ztmpf1 2488 ! 2489 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3 2490 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 2491 !!---------------------------------------------------------------------- 2492 ! 2493 CALL wrk_alloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3 ) 2494 CALL wrk_alloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 ) 2495 2496 z_gsigw3 = 0._wp ; z_gsigt3 = 0._wp ; z_gsi3w3 = 0._wp 2497 z_esigt3 = 0._wp ; z_esigw3 = 0._wp 2498 z_esigtu3 = 0._wp ; z_esigtv3 = 0._wp ; z_esigtf3 = 0._wp 2499 z_esigwu3 = 0._wp ; z_esigwv3 = 0._wp 2500 2501 DO ji = 1, jpi 2502 DO jj = 1, jpj 2503 2504 IF (hbatt(ji,jj)>rn_hc) THEN !deep water, stretched sigma 2505 2506 zzb = hbatt(ji,jj)*rn_zb_a + rn_zb_b ! this forces a linear bottom cell depth relationship with H,. 2507 ! could be changed by users but care must be taken to do so carefully 2508 zzb = 1.0_wp-(zzb/hbatt(ji,jj)) 2509 2510 zzs = rn_zs / hbatt(ji,jj) 2511 2512 IF (rn_efold /= 0.0_wp) THEN 2513 zsmth = tanh( (hbatt(ji,jj)- rn_hc ) / rn_efold ) 2514 ELSE 2515 zsmth = 1.0_wp 2516 ENDIF 2517 2518 DO jk = 1, jpk 2519 z_gsigw3(ji,jj,jk) = REAL(jk-1,wp) /REAL(jpk-1,wp) 2520 z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp) 2521 ENDDO 2522 z_gsigw3(ji,jj,:) = fgamma( z_gsigw3(ji,jj,:), zzb, zzs, zsmth ) 2523 z_gsigt3(ji,jj,:) = fgamma( z_gsigt3(ji,jj,:), zzb, zzs, zsmth ) 2524 2525 ELSE IF (ln_sigcrit) THEN ! shallow water, uniform sigma 2526 2527 DO jk = 1, jpk 2528 z_gsigw3(ji,jj,jk) = REAL(jk-1,wp) /REAL(jpk-1,wp) 2529 z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp) 2530 END DO 2531 2532 ELSE ! shallow water, z coordinates 2533 2534 DO jk = 1, jpk 2535 z_gsigw3(ji,jj,jk) = REAL(jk-1,wp) /REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj)) 2536 z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj)) 2537 END DO 2538 2539 ENDIF 2540 2541 DO jk = 1, jpkm1 2542 z_esigt3(ji,jj,jk) = z_gsigw3(ji,jj,jk+1) - z_gsigw3(ji,jj,jk) 2543 z_esigw3(ji,jj,jk+1) = z_gsigt3(ji,jj,jk+1) - z_gsigt3(ji,jj,jk) 2544 END DO 2545 z_esigw3(ji,jj,1 ) = 2.0_wp * (z_gsigt3(ji,jj,1 ) - z_gsigw3(ji,jj,1 )) 2546 z_esigt3(ji,jj,jpk) = 2.0_wp * (z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk)) 2547 2548 ! Coefficients for vertical depth as the sum of e3w scale factors 2549 z_gsi3w3(ji,jj,1) = 0.5 * z_esigw3(ji,jj,1) 2550 DO jk = 2, jpk 2551 z_gsi3w3(ji,jj,jk) = z_gsi3w3(ji,jj,jk-1) + z_esigw3(ji,jj,jk) 2552 END DO 2553 2554 DO jk = 1, jpk 2555 gdept_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigt3(ji,jj,jk) 2556 gdepw_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigw3(ji,jj,jk) 2557 gde3w_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsi3w3(ji,jj,jk) 2558 END DO 2559 2560 ENDDO ! for all jj's 2561 ENDDO ! for all ji's 2562 2563 DO ji=1,jpi-1 2564 DO jj=1,jpj-1 2565 2566 ! extend to suit for Wetting/Drying case 2567 ztmpu = hbatt(ji,jj)+hbatt(ji+1,jj) 2568 ztmpv = hbatt(ji,jj)+hbatt(ji,jj+1) 2569 ztmpf = hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) 2570 ztmpu1 = hbatt(ji,jj)*hbatt(ji+1,jj) 2571 ztmpv1 = hbatt(ji,jj)*hbatt(ji,jj+1) 2572 ztmpf1 = MIN(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1)) * & 2573 & MAX(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1)) 2574 DO jk = 1, jpk 2575 IF( ln_wd .AND. (ztmpu1 < 0._wp.OR.ABS(ztmpu) < rn_wdmin1) ) THEN 2576 z_esigtu3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji+1,jj,jk) ) 2577 z_esigwu3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji+1,jj,jk) ) 2578 ELSE 2579 z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) & 2580 & / ztmpu 2581 z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) & 2582 & / ztmpu 2583 END IF 2584 2585 IF( ln_wd .AND. (ztmpv1 < 0._wp.OR.ABS(ztmpv) < rn_wdmin1) ) THEN 2586 z_esigtv3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji,jj+1,jk) ) 2587 z_esigwv3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji,jj+1,jk) ) 2588 ELSE 2589 z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) & 2590 & / ztmpv 2591 z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) & 2592 & / ztmpv 2593 END IF 2594 2595 IF( ln_wd .AND. (ztmpf1 < 0._wp.OR.ABS(ztmpf) < rn_wdmin1) ) THEN 2596 z_esigtf3(ji,jj,jk) = 0.25_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji+1,jj,jk) & 2597 & + z_esigt3(ji,jj+1,jk) + z_esigt3(ji+1,jj+1,jk) ) 2598 ELSE 2599 z_esigtf3(ji,jj,jk) = ( hbatt(ji ,jj )*z_esigt3(ji ,jj ,jk) & 2600 & + hbatt(ji+1,jj )*z_esigt3(ji+1,jj ,jk) & 2601 & + hbatt(ji ,jj+1)*z_esigt3(ji ,jj+1,jk) & 2602 & + hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) / ztmpf 2603 END IF 2604 2605 ! Code prior to wetting and drying option (for reference) 2606 !z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) & 2607 ! /( hbatt(ji,jj)+hbatt(ji+1,jj) ) 2608 ! 2609 !z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) & 2610 ! /( hbatt(ji,jj)+hbatt(ji+1,jj) ) 2611 ! 2612 !z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) & 2613 ! /( hbatt(ji,jj)+hbatt(ji,jj+1) ) 2614 ! 2615 !z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) & 2616 ! /( hbatt(ji,jj)+hbatt(ji,jj+1) ) 2617 ! 2618 !z_esigtf3(ji,jj,jk) = ( hbatt(ji ,jj )*z_esigt3(ji ,jj ,jk) & 2619 ! & +hbatt(ji+1,jj )*z_esigt3(ji+1,jj ,jk) & 2620 ! +hbatt(ji ,jj+1)*z_esigt3(ji ,jj+1,jk) & 2621 ! & +hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) & 2622 ! /( hbatt(ji ,jj )+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) ) 2623 2624 e3t_0(ji,jj,jk)=(scosrf(ji,jj)+hbatt(ji,jj))*z_esigt3(ji,jj,jk) 2625 e3u_0(ji,jj,jk)=(scosrf(ji,jj)+hbatu(ji,jj))*z_esigtu3(ji,jj,jk) 2626 e3v_0(ji,jj,jk)=(scosrf(ji,jj)+hbatv(ji,jj))*z_esigtv3(ji,jj,jk) 2627 e3f_0(ji,jj,jk)=(scosrf(ji,jj)+hbatf(ji,jj))*z_esigtf3(ji,jj,jk) 2628 ! 2629 e3w_0 (ji,jj,jk)=hbatt(ji,jj)*z_esigw3(ji,jj,jk) 2630 e3uw_0(ji,jj,jk)=hbatu(ji,jj)*z_esigwu3(ji,jj,jk) 2631 e3vw_0(ji,jj,jk)=hbatv(ji,jj)*z_esigwv3(ji,jj,jk) 2632 END DO 2633 2634 ENDDO 2635 ENDDO 2636 ! 2637 CALL lbc_lnk(e3t_0 ,'T',1.) ; CALL lbc_lnk(e3u_0 ,'T',1.) 2638 CALL lbc_lnk(e3v_0 ,'T',1.) ; CALL lbc_lnk(e3f_0 ,'T',1.) 2639 CALL lbc_lnk(e3w_0 ,'T',1.) 2640 CALL lbc_lnk(e3uw_0,'T',1.) ; CALL lbc_lnk(e3vw_0,'T',1.) 2641 ! 2642 CALL wrk_dealloc( jpi,jpj,jpk, z_gsigw3, z_gsigt3, z_gsi3w3 ) 2643 CALL wrk_dealloc( jpi,jpj,jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 ) 2644 ! 2645 END SUBROUTINE s_sf12 2646 2647 2648 SUBROUTINE s_tanh() 2649 !!---------------------------------------------------------------------- 2650 !! *** ROUTINE s_tanh*** 2651 !! 2652 !! ** Purpose : stretch the s-coordinate system 2653 !! 2654 !! ** Method : s-coordinate stretch 2655 !! 2656 !! Reference : Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408. 2657 !!---------------------------------------------------------------------- 2658 INTEGER :: ji, jj, jk ! dummy loop argument 2659 REAL(wp) :: zcoeft, zcoefw ! temporary scalars 2660 REAL(wp), POINTER, DIMENSION(:) :: z_gsigw, z_gsigt, z_gsi3w 2661 REAL(wp), POINTER, DIMENSION(:) :: z_esigt, z_esigw 2662 !!---------------------------------------------------------------------- 2663 2664 CALL wrk_alloc( jpk, z_gsigw, z_gsigt, z_gsi3w ) 2665 CALL wrk_alloc( jpk, z_esigt, z_esigw ) 2666 2667 z_gsigw = 0._wp ; z_gsigt = 0._wp ; z_gsi3w = 0._wp 2668 z_esigt = 0._wp ; z_esigw = 0._wp 2669 2670 DO jk = 1, jpk 2671 z_gsigw(jk) = -fssig( REAL(jk,wp)-0.5_wp ) 2672 z_gsigt(jk) = -fssig( REAL(jk,wp) ) 2673 END DO 2674 IF( nprint == 1 .AND. lwp ) WRITE(numout,*) 'z_gsigw 1 jpk ', z_gsigw(1), z_gsigw(jpk) 2675 ! 2676 ! Coefficients for vertical scale factors at w-, t- levels 2677 !!gm bug : define it from analytical function, not like juste bellow.... 2678 !!gm or betteroffer the 2 possibilities.... 2679 DO jk = 1, jpkm1 2680 z_esigt(jk ) = z_gsigw(jk+1) - z_gsigw(jk) 2681 z_esigw(jk+1) = z_gsigt(jk+1) - z_gsigt(jk) 2682 END DO 2683 z_esigw( 1 ) = 2._wp * ( z_gsigt(1 ) - z_gsigw(1 ) ) 2684 z_esigt(jpk) = 2._wp * ( z_gsigt(jpk) - z_gsigw(jpk) ) 2685 ! 2686 ! Coefficients for vertical depth as the sum of e3w scale factors 2687 z_gsi3w(1) = 0.5_wp * z_esigw(1) 2688 DO jk = 2, jpk 2689 z_gsi3w(jk) = z_gsi3w(jk-1) + z_esigw(jk) 2690 END DO 2691 !!gm: depuw, depvw can be suppressed (modif in ldfslp) and depw=dep3w can be set (save 3 3D arrays) 2692 DO jk = 1, jpk 2693 zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp) 2694 zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp) 2695 gdept_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigt(jk) + hift(:,:)*zcoeft ) 2696 gdepw_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigw(jk) + hift(:,:)*zcoefw ) 2697 gde3w_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsi3w(jk) + hift(:,:)*zcoeft ) 2698 END DO 2699 !!gm: e3uw, e3vw can be suppressed (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays) 2700 DO jj = 1, jpj 2701 DO ji = 1, jpi 2702 DO jk = 1, jpk 2703 e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigt(jk) + hift(ji,jj)/REAL(jpkm1,wp) ) 2704 e3u_0(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigt(jk) + hifu(ji,jj)/REAL(jpkm1,wp) ) 2705 e3v_0(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigt(jk) + hifv(ji,jj)/REAL(jpkm1,wp) ) 2706 e3f_0(ji,jj,jk) = ( (hbatf(ji,jj)-hiff(ji,jj))*z_esigt(jk) + hiff(ji,jj)/REAL(jpkm1,wp) ) 2707 ! 2708 e3w_0 (ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigw(jk) + hift(ji,jj)/REAL(jpkm1,wp) ) 2709 e3uw_0(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigw(jk) + hifu(ji,jj)/REAL(jpkm1,wp) ) 2710 e3vw_0(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigw(jk) + hifv(ji,jj)/REAL(jpkm1,wp) ) 2711 END DO 2712 END DO 2713 END DO 2714 ! 2715 CALL wrk_dealloc( jpk, z_gsigw, z_gsigt, z_gsi3w ) 2716 CALL wrk_dealloc( jpk, z_esigt, z_esigw ) 2717 ! 2718 END SUBROUTINE s_tanh 2719 2720 2721 FUNCTION fssig( pk ) RESULT( pf ) 2722 !!---------------------------------------------------------------------- 2723 !! *** ROUTINE fssig *** 2724 !! 2725 !! ** Purpose : provide the analytical function in s-coordinate 2726 !! 2727 !! ** Method : the function provide the non-dimensional position of 2728 !! T and W (i.e. between 0 and 1) 2729 !! T-points at integer values (between 1 and jpk) 2730 !! W-points at integer values - 1/2 (between 0.5 and jpk-0.5) 2731 !!---------------------------------------------------------------------- 2732 REAL(wp), INTENT(in) :: pk ! continuous "k" coordinate 2733 REAL(wp) :: pf ! sigma value 2734 !!---------------------------------------------------------------------- 2735 ! 2736 pf = ( TANH( rn_theta * ( -(pk-0.5_wp) / REAL(jpkm1) + rn_thetb ) ) & 2737 & - TANH( rn_thetb * rn_theta ) ) & 2738 & * ( COSH( rn_theta ) & 2739 & + COSH( rn_theta * ( 2._wp * rn_thetb - 1._wp ) ) ) & 2740 & / ( 2._wp * SINH( rn_theta ) ) 2741 ! 2742 END FUNCTION fssig 2743 2744 2745 FUNCTION fssig1( pk1, pbb ) RESULT( pf1 ) 2746 !!---------------------------------------------------------------------- 2747 !! *** ROUTINE fssig1 *** 2748 !! 2749 !! ** Purpose : provide the Song and Haidvogel version of the analytical function in s-coordinate 2750 !! 2751 !! ** Method : the function provides the non-dimensional position of 2752 !! T and W (i.e. between 0 and 1) 2753 !! T-points at integer values (between 1 and jpk) 2754 !! W-points at integer values - 1/2 (between 0.5 and jpk-0.5) 2755 !!---------------------------------------------------------------------- 2756 REAL(wp), INTENT(in) :: pk1 ! continuous "k" coordinate 2757 REAL(wp), INTENT(in) :: pbb ! Stretching coefficient 2758 REAL(wp) :: pf1 ! sigma value 2759 !!---------------------------------------------------------------------- 2760 ! 2761 IF ( rn_theta == 0 ) then ! uniform sigma 2762 pf1 = - ( pk1 - 0.5_wp ) / REAL( jpkm1 ) 2763 ELSE ! stretched sigma 2764 pf1 = ( 1._wp - pbb ) * ( SINH( rn_theta*(-(pk1-0.5_wp)/REAL(jpkm1)) ) ) / SINH( rn_theta ) & 2765 & + pbb * ( (TANH( rn_theta*( (-(pk1-0.5_wp)/REAL(jpkm1)) + 0.5_wp) ) - TANH( 0.5_wp * rn_theta ) ) & 2766 & / ( 2._wp * TANH( 0.5_wp * rn_theta ) ) ) 2767 ENDIF 2768 ! 2769 END FUNCTION fssig1 2770 2771 2772 FUNCTION fgamma( pk1, pzb, pzs, psmth) RESULT( p_gamma ) 2773 !!---------------------------------------------------------------------- 2774 !! *** ROUTINE fgamma *** 2775 !! 2776 !! ** Purpose : provide analytical function for the s-coordinate 2777 !! 2778 !! ** Method : the function provides the non-dimensional position of 2779 !! T and W (i.e. between 0 and 1) 2780 !! T-points at integer values (between 1 and jpk) 2781 !! W-points at integer values - 1/2 (between 0.5 and jpk-0.5) 2782 !! 2783 !! This method allows the maintenance of fixed surface and or 2784 !! bottom cell resolutions (cf. geopotential coordinates) 2785 !! within an analytically derived stretched S-coordinate framework. 2786 !! 2787 !! Reference : Siddorn and Furner, in prep 2788 !!---------------------------------------------------------------------- 2789 REAL(wp), INTENT(in ) :: pk1(jpk) ! continuous "k" coordinate 2790 REAL(wp) :: p_gamma(jpk) ! stretched coordinate 2791 REAL(wp), INTENT(in ) :: pzb ! Bottom box depth 2792 REAL(wp), INTENT(in ) :: pzs ! surface box depth 2793 REAL(wp), INTENT(in ) :: psmth ! Smoothing parameter 2794 ! 2795 INTEGER :: jk ! dummy loop index 2796 REAL(wp) :: za1,za2,za3 ! local scalar 2797 REAL(wp) :: zn1,zn2 ! - - 2798 REAL(wp) :: za,zb,zx ! - - 2799 !!---------------------------------------------------------------------- 2800 ! 2801 zn1 = 1._wp / REAL( jpkm1, wp ) 2802 zn2 = 1._wp - zn1 2803 ! 2804 za1 = (rn_alpha+2.0_wp)*zn1**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*zn1**(rn_alpha+2.0_wp) 2805 za2 = (rn_alpha+2.0_wp)*zn2**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*zn2**(rn_alpha+2.0_wp) 2806 za3 = (zn2**3.0_wp - za2)/( zn1**3.0_wp - za1) 2807 ! 2808 za = pzb - za3*(pzs-za1)-za2 2809 za = za/( zn2-0.5_wp*(za2+zn2**2.0_wp) - za3*(zn1-0.5_wp*(za1+zn1**2.0_wp) ) ) 2810 zb = (pzs - za1 - za*( zn1-0.5_wp*(za1+zn1**2.0_wp ) ) ) / (zn1**3.0_wp - za1) 2811 zx = 1.0_wp-za/2.0_wp-zb 2812 ! 2813 DO jk = 1, jpk 2814 p_gamma(jk) = za*(pk1(jk)*(1.0_wp-pk1(jk)/2.0_wp))+zb*pk1(jk)**3.0_wp + & 2815 & zx*( (rn_alpha+2.0_wp)*pk1(jk)**(rn_alpha+1.0_wp)- & 2816 & (rn_alpha+1.0_wp)*pk1(jk)**(rn_alpha+2.0_wp) ) 2817 p_gamma(jk) = p_gamma(jk)*psmth+pk1(jk)*(1.0_wp-psmth) 2818 END DO 2819 ! 2820 END FUNCTION fgamma 315 END SUBROUTINE zgr_tb_level 2821 316 2822 317 !!====================================================================== -
branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/DOM/iscplrst.F90
r6140 r6667 50 50 !!---------------------------------------------------------------------- 51 51 INTEGER :: inum0 52 REAL(wp), DIMENSION(:,: ), POINTER :: zsmask_b53 REAL(wp), DIMENSION(:,:,:), POINTER :: ztmask_b, zumask_b, zvmask_b54 REAL(wp), DIMENSION(:,:,:), POINTER :: ze3t_b , ze3u_b , ze3v_b55 REAL(wp), DIMENSION(:,:,:), POINTER :: zdepw_b52 REAL(wp), DIMENSION(:,: ), POINTER :: zsmask_b 53 REAL(wp), DIMENSION(:,:,:), POINTER :: ztmask_b, zumask_b, zvmask_b 54 REAL(wp), DIMENSION(:,:,:), POINTER :: ze3t_b , ze3u_b , ze3v_b 55 REAL(wp), DIMENSION(:,:,:), POINTER :: zdepw_b 56 56 CHARACTER(20) :: cfile 57 57 !!---------------------------------------------------------------------- 58 58 59 CALL wrk_alloc(jpi,jpj,jpk, ztmask_b, zumask_b, zvmask_b) ! mask before60 CALL wrk_alloc(jpi,jpj,jpk, ze3t_b , ze3u_b , ze3v_b ) ! e3 before61 CALL wrk_alloc(jpi,jpj,jpk, zdepw_b )62 CALL wrk_alloc(jpi,jpj, zsmask_b )59 CALL wrk_alloc(jpi,jpj,jpk, ztmask_b, zumask_b, zvmask_b) ! mask before 60 CALL wrk_alloc(jpi,jpj,jpk, ze3t_b , ze3u_b , ze3v_b ) ! e3 before 61 CALL wrk_alloc(jpi,jpj,jpk, zdepw_b ) 62 CALL wrk_alloc(jpi,jpj, zsmask_b ) 63 63 64 64 … … 86 86 87 87 !! print mesh/mask 88 IF( n msh /= 0 .AND. ln_iscpl ) CALL dom_wri ! Create a domain file88 IF( nn_msh /= 0 .AND. ln_iscpl ) CALL dom_wri ! Create a domain file 89 89 90 90 IF ( ln_hsb ) THEN … … 98 98 END IF 99 99 100 CALL wrk_dealloc(jpi,jpj,jpk, ztmask_b,zumask_b,zvmask_b )101 CALL wrk_dealloc(jpi,jpj,jpk, ze3t_b ,ze3u_b ,ze3v_b )102 CALL wrk_dealloc(jpi,jpj,jpk, zdepw_b )103 CALL wrk_dealloc(jpi,jpj, zsmask_b )100 CALL wrk_dealloc(jpi,jpj,jpk, ztmask_b,zumask_b,zvmask_b ) 101 CALL wrk_dealloc(jpi,jpj,jpk, ze3t_b ,ze3u_b ,ze3v_b ) 102 CALL wrk_dealloc(jpi,jpj,jpk, zdepw_b ) 103 CALL wrk_dealloc(jpi,jpj, zsmask_b ) 104 104 105 105 !! next step is an euler time step … … 108 108 !! set _b and _n variables equal 109 109 tsb (:,:,:,:) = tsn (:,:,:,:) 110 ub (:,:,: ) = un (:,:,:)111 vb (:,:,: ) = vn (:,:,:)112 sshb(:,: )= sshn(:,:)110 ub (:,:,:) = un (:,:,:) 111 vb (:,:,:) = vn (:,:,:) 112 sshb(:,:) = sshn(:,:) 113 113 114 114 !! set _b and _n vertical scale factor equal … … 117 117 e3v_b (:,:,:) = e3v_n (:,:,:) 118 118 119 e3uw_b (:,:,:) = e3uw_n(:,:,:)120 e3vw_b (:,:,:) = e3vw_n(:,:,:)121 gdept_b(:,:,:) 119 e3uw_b (:,:,:) = e3uw_n (:,:,:) 120 e3vw_b (:,:,:) = e3vw_n (:,:,:) 121 gdept_b(:,:,:) = gdept_n(:,:,:) 122 122 gdepw_b(:,:,:) = gdepw_n(:,:,:) 123 hu_b (:,:) = hu_n(:,:)124 hv_b (:,:) = hv_n(:,:)125 r1_hu_b(:,:) = r1_hu_n(:,:)126 r1_hv_b(:,:) = r1_hv_n(:,:)123 hu_b (:,:) = hu_n (:,:) 124 hv_b (:,:) = hv_n (:,:) 125 r1_hu_b(:,:) = r1_hu_n(:,:) 126 r1_hv_b(:,:) = r1_hv_n(:,:) 127 127 ! 128 128 END SUBROUTINE iscpl_stp 129 129 130 130 131 SUBROUTINE iscpl_rst_interpol (ptmask_b, pumask_b, pvmask_b, psmask_b, pe3t_b, pe3u_b, pe3v_b, pdepw_b) 131 132 !!---------------------------------------------------------------------- … … 436 437 CALL wrk_dealloc(jpi,jpj, zbub , zbvb , zbun , zbvn ) 437 438 CALL wrk_dealloc(jpi,jpj, zssh0 , zssh1 , zhu1 , zhv1 ) 438 439 ! 439 440 END SUBROUTINE iscpl_rst_interpol 440 441 442 !!====================================================================== 441 443 END MODULE iscplrst -
branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
r6152 r6667 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_wdmin2458 ll_tmp3 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(- bathy(ji,jj), -bathy(ji+1,jj)) + &456 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-ht_0(ji,jj), -ht_0(ji+1,jj)) 457 ll_tmp2 = MAX(sshn(ji,jj) + ht_0(ji,jj), sshn(ji+1,jj) + ht_0(ji+1,jj)) > rn_wdmin1 + rn_wdmin2 458 ll_tmp3 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-ht_0(ji,jj), -ht_0(ji+1,jj)) + & 459 459 & rn_wdmin1 + rn_wdmin2 460 460 … … 464 464 ELSE IF(ll_tmp3) THEN 465 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)) / &466 zcpx(ji,jj) = ABS((sshn(ji+1,jj) + ht_0(ji+1,jj) - sshn(ji,jj) - ht_0(ji,jj)) / & 467 467 & (sshn(ji+1,jj) - sshn(ji,jj))) 468 468 wduflt(ji,jj) = 1.0_wp … … 472 472 END IF 473 473 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_wdmin2476 ll_tmp3 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(- bathy(ji,jj), -bathy(ji,jj+1)) + &474 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-ht_0(ji,jj), -ht_0(ji,jj+1)) 475 ll_tmp2 = MAX(sshn(ji,jj) + ht_0(ji,jj), sshn(ji,jj+1) + ht_0(ji,jj+1)) > rn_wdmin1 + rn_wdmin2 476 ll_tmp3 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-ht_0(ji,jj), -ht_0(ji,jj+1)) + & 477 477 & rn_wdmin1 + rn_wdmin2 478 478 … … 482 482 ELSE IF(ll_tmp3) THEN 483 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)) / &484 zcpy(ji,jj) = ABS((sshn(ji,jj+1) + ht_0(ji,jj+1) - sshn(ji,jj) - ht_0(ji,jj)) / & 485 485 & (sshn(ji,jj+1) - sshn(ji,jj))) 486 486 wdvflt(ji,jj) = 1.0_wp … … 707 707 DO jj = 2, jpjm1 708 708 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)) &709 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-ht_0(ji,jj), -ht_0(ji+1,jj)) & 710 & .and. MAX(sshn(ji,jj) + ht_0(ji,jj), sshn(ji+1,jj) + ht_0(ji+1,jj)) & 711 711 & > rn_wdmin1 + rn_wdmin2 712 ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(- bathy(ji,jj), -bathy(ji+1,jj)) +&712 ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-ht_0(ji,jj), -ht_0(ji+1,jj)) +& 713 713 & rn_wdmin1 + rn_wdmin2 714 714 … … 717 717 ELSE IF(ll_tmp2) THEN 718 718 ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen ! here 719 zcpx(ji,jj) = ABS((sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) /&719 zcpx(ji,jj) = ABS((sshn(ji+1,jj) + ht_0(ji+1,jj) - sshn(ji,jj) - ht_0(ji,jj)) /& 720 720 & (sshn(ji+1,jj) - sshn(ji,jj))) 721 721 ELSE … … 723 723 END IF 724 724 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)) &725 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-ht_0(ji,jj), -ht_0(ji,jj+1)) & 726 & .and. MAX(sshn(ji,jj) + ht_0(ji,jj), sshn(ji,jj+1) + ht_0(ji,jj+1)) & 727 727 & > rn_wdmin1 + rn_wdmin2 728 ll_tmp2 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(- bathy(ji,jj), -bathy(ji,jj+1)) +&728 ll_tmp2 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-ht_0(ji,jj), -ht_0(ji,jj+1)) +& 729 729 & rn_wdmin1 + rn_wdmin2 730 730 … … 733 733 ELSE IF(ll_tmp2) THEN 734 734 ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen ! here 735 zcpy(ji,jj) = ABS((sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) /&735 zcpy(ji,jj) = ABS((sshn(ji,jj+1) + ht_0(ji,jj+1) - sshn(ji,jj) - ht_0(ji,jj)) /& 736 736 & (sshn(ji,jj+1) - sshn(ji,jj))) 737 737 ELSE … … 1003 1003 DO jj = 2, jpjm1 1004 1004 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)) &1005 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-ht_0(ji,jj), -ht_0(ji+1,jj)) & 1006 & .and. MAX(sshn(ji,jj) + ht_0(ji,jj), sshn(ji+1,jj) + ht_0(ji+1,jj)) & 1007 1007 & > rn_wdmin1 + rn_wdmin2 1008 ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(- bathy(ji,jj), -bathy(ji+1,jj)) +&1008 ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-ht_0(ji,jj), -ht_0(ji+1,jj)) +& 1009 1009 & rn_wdmin1 + rn_wdmin2 1010 1010 … … 1013 1013 ELSE IF(ll_tmp2) THEN 1014 1014 ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen ! here 1015 zcpx(ji,jj) = ABS((sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) /&1015 zcpx(ji,jj) = ABS((sshn(ji+1,jj) + ht_0(ji+1,jj) - sshn(ji,jj) - ht_0(ji,jj)) /& 1016 1016 & (sshn(ji+1,jj) - sshn(ji,jj))) 1017 1017 ELSE … … 1019 1019 END IF 1020 1020 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)) &1021 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-ht_0(ji,jj), -ht_0(ji,jj+1)) & 1022 & .and. MAX(sshn(ji,jj) + ht_0(ji,jj), sshn(ji,jj+1) + ht_0(ji,jj+1)) & 1023 1023 & > rn_wdmin1 + rn_wdmin2 1024 ll_tmp2 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(- bathy(ji,jj), -bathy(ji,jj+1)) +&1024 ll_tmp2 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-ht_0(ji,jj), -ht_0(ji,jj+1)) +& 1025 1025 & rn_wdmin1 + rn_wdmin2 1026 1026 … … 1029 1029 ELSE IF(ll_tmp2) THEN 1030 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)) /&1031 zcpy(ji,jj) = ABS((sshn(ji,jj+1) + ht_0(ji,jj+1) - sshn(ji,jj) - ht_0(ji,jj)) /& 1032 1032 & (sshn(ji,jj+1) - sshn(ji,jj))) 1033 1033 ELSE … … 1046 1046 DO jj = 1, jpj 1047 1047 DO ji = 1, jpi 1048 jk = mb athy(ji,jj)1048 jk = mbkt(ji,jj)+1 1049 1049 IF( jk <= 0 ) THEN ; zrhh(ji,jj, : ) = 0._wp 1050 1050 ELSEIF( jk == 1 ) THEN ; zrhh(ji,jj,jk+1:jpk) = rhd(ji,jj,jk) -
branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
r6596 r6667 255 255 zwz(:,:) = 0._wp 256 256 zhf(:,:) = 0._wp 257 IF ( .not. ln_sco ) THEN 258 259 !!gm agree the JC comment : this should be done in a much clear way 260 261 ! JC: It not clear yet what should be the depth at f-points over land in z-coordinate case 262 ! Set it to zero for the time being 263 ! IF( rn_hmin < 0._wp ) THEN ; jk = - INT( rn_hmin ) ! from a nb of level 264 ! ELSE ; jk = MINLOC( gdepw_0, mask = gdepw_0 > rn_hmin, dim = 1 ) ! from a depth 265 ! ENDIF 266 ! zhf(:,:) = gdepw_0(:,:,jk+1) 267 ELSE 268 zhf(:,:) = hbatf(:,:) 269 END IF 270 271 DO jj = 1, jpjm1 272 zhf(:,jj) = zhf(:,jj) * (1._wp- umask(:,jj,1) * umask(:,jj+1,1)) 273 END DO 257 258 !!gm assume 0 in both cases (xhich is almost surely WRONG ! ) as hvatf has been removed 259 !!gm A priori a better value should be something like : 260 !!gm zhf(i,j) = masked sum of ht(i,j) , ht(i+1,j) , ht(i,j+1) , (i+1,j+1) 261 !!gm divided by the sum of the corresponding mask 262 !!gm 263 !! 264 !! IF ( .not. ln_sco ) THEN 265 !! 266 !! !!gm agree the JC comment : this should be done in a much clear way 267 !! 268 !! ! JC: It not clear yet what should be the depth at f-points over land in z-coordinate case 269 !! ! Set it to zero for the time being 270 !! ! IF( rn_hmin < 0._wp ) THEN ; jk = - INT( rn_hmin ) ! from a nb of level 271 !! ! ELSE ; jk = MINLOC( gdepw_0, mask = gdepw_0 > rn_hmin, dim = 1 ) ! from a depth 272 !! ! ENDIF 273 !! ! zhf(:,:) = gdepw_0(:,:,jk+1) 274 !!