| 1 | MODULE ldfeke |
|---|
| 2 | !!====================================================================== |
|---|
| 3 | !! *** MODULE ldfeke *** |
|---|
| 4 | !! Ocean physics: Eddy induced velocity coefficient according to the |
|---|
| 5 | !! GEOMETRIC parameterization scheme |
|---|
| 6 | !!===================================================================== |
|---|
| 7 | !! History : 4.0 ! 2017-11 (J. Mak, G. Madec) original code |
|---|
| 8 | !!---------------------------------------------------------------------- |
|---|
| 9 | |
|---|
| 10 | !!---------------------------------------------------------------------- |
|---|
| 11 | !! ldf_eke : time step depth-integrated EKE and update aeiw (and |
|---|
| 12 | !! from that aeiu and aeiv) according to the GEOMETRIC |
|---|
| 13 | !! parameterization scheme |
|---|
| 14 | !! ldf_eke_init : initialization, namelist read, and parameters control |
|---|
| 15 | !! eke_rst : read/write eke restart in ocean restart file |
|---|
| 16 | !!---------------------------------------------------------------------- |
|---|
| 17 | USE oce ! ocean: dynamics and active tracers variables |
|---|
| 18 | USE phycst ! physical constants |
|---|
| 19 | USE dom_oce ! domain: ocean |
|---|
| 20 | USE ldfslp ! lateral physics: slope of iso-neutral surfaces |
|---|
| 21 | USE ldftra ! lateral physics: eddy coefficients |
|---|
| 22 | USE dynldf_lap_blp, ONLY : eke_keS ! source term of eke due to KE dissipation |
|---|
| 23 | USE zdfmxl ! vertical physics: mixed layer |
|---|
| 24 | |
|---|
| 25 | USE dynspg_ts , ONLY : un_adv, vn_adv |
|---|
| 26 | |
|---|
| 27 | ! |
|---|
| 28 | USE in_out_manager ! I/O manager |
|---|
| 29 | USE iom ! I/O manager library |
|---|
| 30 | USE lib_mpp ! MPP library |
|---|
| 31 | USE lbclnk ! ocean lateral boundary conditions (or mpp link) |
|---|
| 32 | USE prtctl ! Print control |
|---|
| 33 | USE timing ! Timing |
|---|
| 34 | USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) |
|---|
| 35 | |
|---|
| 36 | IMPLICIT NONE |
|---|
| 37 | PRIVATE |
|---|
| 38 | |
|---|
| 39 | PUBLIC ldf_eke ! routine called in step module |
|---|
| 40 | PUBLIC ldf_eke_init ! routine called in opa module |
|---|
| 41 | |
|---|
| 42 | |
|---|
| 43 | ! !!** Namelist namldf_eke ** |
|---|
| 44 | REAL(wp) :: rn_ekedis ! dissipation time scale of EKE [days] |
|---|
| 45 | REAL(wp) :: rn_geom ! geometric parameterization master coefficient [-] |
|---|
| 46 | REAL(wp) :: rn_eke_lap ! diffusion of EKE [m2/s] |
|---|
| 47 | REAL(wp) :: rn_eke_init ! initial value of total EKE [m2/s2] |
|---|
| 48 | REAL(wp) :: rn_eke_min ! background value of total EKE [m2/s2] |
|---|
| 49 | REAL(wp) :: rn_ross_min ! tapering based of minimum Rossby radius [m] |
|---|
| 50 | REAL(wp) :: rn_aeiv_min, rn_aeiv_max ! min and max bounds of aeiv coefficient [m2/s] |
|---|
| 51 | REAL(wp) :: rn_SFmin, rn_SFmax ! min and max bounds of Structure Function [-] |
|---|
| 52 | REAL(wp) :: zf0, zbeta ! f0 and beta for computing Rossby speed |
|---|
| 53 | ! |
|---|
| 54 | INTEGER :: nn_eke_opt ! option for what term to include in eddy energy budget |
|---|
| 55 | INTEGER :: nn_eke_dis ! option for taking constant or spatially varying linear dissipation |
|---|
| 56 | ! |
|---|
| 57 | LOGICAL :: ln_adv_wav ! option for having advection by Rossby wave or not |
|---|
| 58 | LOGICAL :: ln_beta_plane ! option for computing long Rossby phase speed |
|---|
| 59 | ! |
|---|
| 60 | REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: eke_b, eke_n, eke_a ! vertical sum of total Eddy Kinetic Energy [m3/s2] |
|---|
| 61 | REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: r1_ekedis ! linear dissipation rate (= 1/rn_ekedis) [ /s ] |
|---|
| 62 | REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zc1, zc_ros, zadv_wav ! 1st baroclinic mode and long Rossby speed [m /s ] |
|---|
| 63 | |
|---|
| 64 | !! * Substitutions |
|---|
| 65 | # include "vectopt_loop_substitute.h90" |
|---|
| 66 | !!---------------------------------------------------------------------- |
|---|
| 67 | !! NEMO/OPA 4.0 , NEMO Consortium (2017) |
|---|
| 68 | !! $Id: ldfeke.F90 7813 2017-03-20 16:17:45Z clem $ |
|---|
| 69 | !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) |
|---|
| 70 | !!---------------------------------------------------------------------- |
|---|
| 71 | CONTAINS |
|---|
| 72 | |
|---|
| 73 | SUBROUTINE ldf_eke( kt ) |
|---|
| 74 | !!---------------------------------------------------------------------- |
|---|
| 75 | !! *** ROUTINE ldf_eke *** |
|---|
| 76 | !! |
|---|
| 77 | !! ** Purpose : implements the GEOMETRIC parameterization |
|---|
| 78 | !! |
|---|
| 79 | !! ** Notes : If nn_aei_ijk_t = 32 then eke and aeiw are BOTH updated |
|---|
| 80 | !! If ln_eke_equ = .true. in namtra_ldfeiv but nn_aei_ijk_t |
|---|
| 81 | !! is something else, then ONLY eddy equation is updated |
|---|
| 82 | !! (but the eddy energy is passive and doesn't do anything) |
|---|
| 83 | !! |
|---|
| 84 | !! ** Method : GEOMETRIC calculates the Gent-McWilliams / eddy induced |
|---|
| 85 | !! velocity coefficient according to |
|---|
| 86 | !! |
|---|
| 87 | !! aeiw = alpha * (\int E dz) / (\int S M^2 / N dz), |
|---|
| 88 | !! |
|---|
| 89 | !! where (\int E dz) is the depth-integrated eddy energy |
|---|
| 90 | !! (at the previous time level), informed by a parameterized |
|---|
| 91 | !! depth-integrated eddy energy, where |
|---|
| 92 | !! |
|---|
| 93 | !! nn_eke_opt = 0 => default: just PE->EKE growth and linear dissipation |
|---|
| 94 | !! ! |
|---|
| 95 | !! = 1 => default + advection by depth-averaged flow |
|---|
| 96 | !! ! |
|---|
| 97 | !! = 2 => default + advection + contribution to EKE from |
|---|
| 98 | !! momentum dissipation |
|---|
| 99 | !! ! |
|---|
| 100 | !! = 88 => ONLY advection |
|---|
| 101 | !! ! |
|---|
| 102 | !! = 99 => ONLY Laplacian diffusion |
|---|
| 103 | !! |
|---|
| 104 | !! S is a structure function, and M and N are horizontal and |
|---|
| 105 | !! vertical buoyancy frequencies |
|---|
| 106 | !! |
|---|
| 107 | !! linear dissipation may be specified by |
|---|
| 108 | !! |
|---|
| 109 | !! nn_eke_dus = 0 => constant |
|---|
| 110 | !! ! |
|---|
| 111 | !! =-20 => read in a geom_diss_2D.nc field (give it in days) |
|---|
| 112 | !! |
|---|
| 113 | !! ** Action : * time step depth-integrated eddy energy budget |
|---|
| 114 | !! * calculate aeiw |
|---|
| 115 | !! |
|---|
| 116 | !! References : Marshall, Maddison, Berloff JPO 2012 ; Mak et al JPO 2018 |
|---|
| 117 | !!---------------------------------------------------------------------- |
|---|
| 118 | INTEGER, INTENT(in) :: kt ! ocean time step |
|---|
| 119 | ! |
|---|
| 120 | INTEGER :: ji, jj, jk ! dummy loop arguments |
|---|
| 121 | INTEGER :: ik ! local integer |
|---|
| 122 | REAL(wp) :: zn2_min = 1.e-8 ! minimum value of N2 used to compute structure function |
|---|
| 123 | REAL(wp) :: z1_f20, zfw, ze3w ! local scalar |
|---|
| 124 | REAL(wp) :: zfp_ui, zfp_vj ! - - |
|---|
| 125 | REAL(wp) :: zfm_ui, zfm_vj ! - - |
|---|
| 126 | REAL(wp) :: zn_slp2, zn2 ! - - |
|---|
| 127 | REAL(wp) :: zmsku, zaeiu_w ! - - |
|---|
| 128 | REAL(wp) :: zmskv, zaeiv_w ! - - |
|---|
| 129 | REAL(wp) :: zen, zed ! - - |
|---|
| 130 | REAL(wp) :: zeke_rhs ! - - |
|---|
| 131 | REAL(wp) :: zck, zwslpi, zwslpj ! - - tapering near coasts |
|---|
| 132 | REAL(wp) :: zc_rosu ! - - |
|---|
| 133 | REAL(wp), DIMENSION(jpi,jpj) :: zeke_peS, zn_slp ! 2D workspace, PE-KE conversion |
|---|
| 134 | REAL(wp), DIMENSION(jpi,jpj) :: zadv_ubt, zwx, zwy ! - - barotropic advection |
|---|
| 135 | REAL(wp), DIMENSION(jpi,jpj) :: zlap, zaheeu, zaheev, ztu, ztv ! - - diffusion |
|---|
| 136 | REAL(wp), DIMENSION(jpi,jpj) :: zdis ! - - linear dissipation |
|---|
| 137 | REAL(wp), DIMENSION(jpi,jpj) :: zn, zross ! - - tapering near coasts |
|---|
| 138 | REAL(wp), DIMENSION(jpi,jpj,jpk) :: zaeiw, zzSF ! 3D workspace |
|---|
| 139 | !!-------------------------------------------------------------------- |
|---|
| 140 | ! |
|---|
| 141 | IF( ln_timing ) CALL timing_start('ldf_eke') |
|---|
| 142 | ! |
|---|
| 143 | IF( kt == nit000 .AND. lwp) THEN !* Control print |
|---|
| 144 | WRITE(numout,*) |
|---|
| 145 | WRITE(numout,*) 'ldf_eke : GEOMETRIC parameterization (total EKE time evolution equation)' |
|---|
| 146 | WRITE(numout,*) '~~~~~~~' |
|---|
| 147 | ENDIF |
|---|
| 148 | |
|---|
| 149 | ! !== EIV mean conversion to EKE (ah N^2 slp2) & N slp ==! |
|---|
| 150 | ! |
|---|
| 151 | ! work out the 3D structure function here |
|---|
| 152 | ! |
|---|
| 153 | ! current: Ferreria et al, zzSF = N^2 / N^2_ref, W-points |
|---|
| 154 | ! capped between rn_SFmax and rn_SFmin |
|---|
| 155 | zzSF(:,:,:) = 0._wp |
|---|
| 156 | DO jk = 2, jpkm1 |
|---|
| 157 | DO jj = 1, jpj |
|---|
| 158 | DO ji = 1, jpi |
|---|
| 159 | IF( jk <= nmln(ji,jj) ) THEN ! above and at mixed layer |
|---|
| 160 | zzSF(ji,jj,jk) = rn_SFmax |
|---|
| 161 | ELSE |
|---|
| 162 | ik = MIN( nmln(ji,jj), mbkt(ji,jj) ) + 1 ! one level below the mixed layer (MIN in case ML depth is the ocean depth) |
|---|
| 163 | zzSF(ji,jj,jk) = MAX( 0._wp , rn2b(ji,jj,jk) ) / MAX( zn2_min , rn2b(ji,jj,ik) ) ! Structure Function : N^2 / N^2_ref |
|---|
| 164 | zzSF(ji,jj,jk) = MAX( rn_SFmin , MIN( zzSF(ji,jj,jk) , rn_SFmax ) ) |
|---|
| 165 | ENDIF |
|---|
| 166 | END DO |
|---|
| 167 | END DO |
|---|
| 168 | END DO |
|---|
| 169 | ! |
|---|
| 170 | ! !* parametrized PE_EKE conversion due to eddy induced velocity |
|---|
| 171 | zeke_peS(:,:) = 0._wp |
|---|
| 172 | zn_slp (:,:) = 0._wp |
|---|
| 173 | zn (:,:) = 0._wp |
|---|
| 174 | DO jk = 2, jpkm1 ! query: index? |
|---|
| 175 | DO jj = 2, jpjm1 |
|---|
| 176 | DO ji = fs_2, fs_jpim1 ! vector opt. ! NB: ah_slp2 is w-masked |
|---|
| 177 | zmsku = wmask(ji,jj,jk) / MAX( umask(ji ,jj,jk-1) + umask(ji-1,jj,jk) & |
|---|
| 178 | & + umask(ji-1,jj,jk-1) + umask(ji ,jj,jk) , 1._wp ) |
|---|
| 179 | zmskv = wmask(ji,jj,jk) / MAX( vmask(ji,jj ,jk-1) + vmask(ji,jj-1,jk) & |
|---|
| 180 | & + vmask(ji,jj-1,jk-1) + vmask(ji,jj ,jk) , 1._wp ) |
|---|
| 181 | ! |
|---|
| 182 | zaeiu_w = ( aeiu(ji ,jj,jk-1) + aeiu(ji-1,jj,jk) & |
|---|
| 183 | & + aeiu(ji-1,jj,jk-1) + aeiu(ji ,jj,jk) ) * zmsku |
|---|
| 184 | zaeiv_w = ( aeiv(ji,jj ,jk-1) + aeiv(ji,jj-1,jk) & |
|---|
| 185 | & + aeiv(ji,jj-1,jk-1) + aeiv(ji,jj ,jk) ) * zmskv |
|---|
| 186 | ! |
|---|
| 187 | zn_slp2 = ( zaeiu_w * wslpi(ji,jj,jk) * wslpi(ji,jj,jk) & ! (slope ** 2) * aeiv |
|---|
| 188 | & + zaeiv_w * wslpj(ji,jj,jk) * wslpj(ji,jj,jk) ) ! JM 28 Jun: undo slope reduction here too? |
|---|
| 189 | zn2 = MAX( 0._wp , rn2b(ji,jj,jk) ) |
|---|
| 190 | ! |
|---|
| 191 | ze3w = e3w_b(ji,jj,jk) * tmask(ji,jj,jk) |
|---|
| 192 | zn(ji,jj) = zn(ji,jj) + SQRT( zn2 ) * ze3w ! for working out taper at small rossby radius regions later |
|---|
| 193 | ! |
|---|
| 194 | zeke_peS(ji,jj) = zeke_peS(ji,jj) + ze3w * zn2 * zn_slp2 ! note this is >=0 |
|---|
| 195 | ! |
|---|
| 196 | zck = ( umask(ji,jj,jk) + umask(ji-1,jj,jk) ) & ! taken from ldfslp, undo the slope reduction |
|---|
| 197 | & * ( vmask(ji,jj,jk) + vmask(ji,jj-1,jk) ) * 0.25 ! near topographic features |
|---|
| 198 | zwslpi = wslpi(ji,jj,jk) / MAX( zck, 0.1_wp) ! (just to avoid dividing by zeros) |
|---|
| 199 | zwslpj = wslpj(ji,jj,jk) / MAX( zck, 0.1_wp) |
|---|
| 200 | ! |
|---|
| 201 | zn_slp(ji,jj) = zn_slp(ji,jj) + ze3w * zzSF(ji,jj,jk) & ! note this >=0 and structure function weighted |
|---|
| 202 | & * SQRT( zn2 * ( zwslpi * zwslpi + zwslpj * zwslpj ) ) |
|---|
| 203 | END DO |
|---|
| 204 | END DO |
|---|
| 205 | END DO |
|---|
| 206 | |
|---|
| 207 | ! !* upstream advection with initial mass fluxes & intermediate update |
|---|
| 208 | ! !* upstream tracer flux in the i and j direction |
|---|
| 209 | DO jj = 1, jpjm1 |
|---|
| 210 | DO ji = 1, fs_jpim1 ! vector opt. |
|---|
| 211 | ! upstream scheme |
|---|
| 212 | zfp_ui = un_adv(ji,jj) + ABS( un_adv(ji,jj) ) |
|---|
| 213 | zfm_ui = un_adv(ji,jj) - ABS( un_adv(ji,jj) ) |
|---|
| 214 | zfp_vj = vn_adv(ji,jj) + ABS( vn_adv(ji,jj) ) |
|---|
| 215 | zfm_vj = vn_adv(ji,jj) - ABS( vn_adv(ji,jj) ) |
|---|
| 216 | zwx(ji,jj) = 0.5 * ( zfp_ui * eke_b(ji,jj) + zfm_ui * eke_b(ji+1,jj ) ) |
|---|
| 217 | zwy(ji,jj) = 0.5 * ( zfp_vj * eke_b(ji,jj) + zfm_vj * eke_b(ji ,jj+1) ) |
|---|
| 218 | END DO |
|---|
| 219 | END DO |
|---|
| 220 | ! !* divergence of ubt advective fluxes |
|---|
| 221 | DO jj = 2, jpjm1 |
|---|
| 222 | DO ji = fs_2, fs_jpim1 ! vector opt. |
|---|
| 223 | zadv_ubt(ji,jj) = - ( zwx(ji,jj) - zwx(ji-1,jj ) & |
|---|
| 224 | & + zwy(ji,jj) - zwy(ji ,jj-1) ) * r1_e1e2t(ji,jj) |
|---|
| 225 | END DO |
|---|
| 226 | END DO |
|---|
| 227 | IF( ln_linssh ) THEN !* top value (linear free surf. only as zwz is multiplied by wmask) |
|---|
| 228 | DO jj = 2, jpjm1 |
|---|
| 229 | DO ji = fs_2, fs_jpim1 ! vector opt. |
|---|
| 230 | zadv_ubt(ji,jj) = - wn(ji,jj,1) * eke_b(ji,jj) / ( ht_0(ji,jj) + 1._wp-tmask(ji,jj,1) ) ! jm (03 Mar 18): was eke_n |
|---|
| 231 | END DO |
|---|
| 232 | END DO |
|---|
| 233 | ENDIF |
|---|
| 234 | ! |
|---|
| 235 | ! !* same as above but for advection by Rossby waves |
|---|
| 236 | zadv_wav(:,:) = 0._wp |
|---|
| 237 | IF ( ln_adv_wav ) THEN |
|---|
| 238 | ! |
|---|
| 239 | DO jj = 2, jpjm1 ! first work out the Rossby phase speeds |
|---|
| 240 | DO ji = 2, fs_jpim1 |
|---|
| 241 | ! compute only for deep enough places |
|---|
| 242 | IF ( ht_0(ji,jj) > 300._wp ) THEN ! jm: should use something like ht_b really... |
|---|
| 243 | ! compute vertical mode phase speed on T point |
|---|
| 244 | zc1(ji,jj) = MIN( 10._wp, zn(ji,jj) / rpi ) |
|---|
| 245 | ! compute long Rossby phase speed on T point (minus sign later) |
|---|
| 246 | IF ( ln_beta_plane ) THEN |
|---|
| 247 | zc_ros(ji,jj) = zc1(ji,jj) * zc1(ji,jj) * zbeta / (zf0 * zf0) |
|---|
| 248 | ELSE |
|---|
| 249 | zc_ros(ji,jj) = zc1(ji,jj) * zc1(ji,jj) * COS( rad * gphit(ji,jj) ) & |
|---|
| 250 | / ( ra * ff_t(ji,jj) * SIN( rad * gphit(ji,jj) ) & |
|---|
| 251 | + rsmall ) |
|---|
| 252 | END IF |
|---|
| 253 | ! cap the Rossby phase speeds by the fastest equatorial Rossby wave speed |
|---|
| 254 | ! the minus sign for westward propagation goes here |
|---|
| 255 | zc_ros(ji,jj) = -MIN( zc1(ji,jj) / 3._wp, zc_ros(ji,jj) ) |
|---|
| 256 | END IF |
|---|
| 257 | END DO |
|---|
| 258 | END DO |
|---|
| 259 | ! |
|---|
| 260 | zwx(:,:) = 0._wp ! wipe the advective contribution from above |
|---|
| 261 | ! |
|---|
| 262 | DO jj = 2, jpjm1 ! average Rossby phase speeds onto U points |
|---|
| 263 | DO ji = 1, fs_jpim1 |
|---|
| 264 | zc_rosu = 0.5 * ( zc_ros(ji,jj) + zc_ros(ji+1,jj) ) * umask(ji,jj,1) |
|---|
| 265 | zfp_ui = zc_rosu + ABS( zc_rosu ) |
|---|
| 266 | zfm_ui = zc_rosu - ABS( zc_rosu ) |
|---|
| 267 | zwx(ji,jj) = 0.5 * ( zfp_ui * eke_b(ji,jj) + zfm_ui * eke_b(ji+1,jj) ) |
|---|
| 268 | END DO |
|---|
| 269 | END DO |
|---|
| 270 | ! !* divergence of wav advective fluxes (only e1 here) |
|---|
| 271 | z1_f20 = 1._wp / ( 2._wp * omega * sin( rad * 20._wp ) ) |
|---|
| 272 | DO jj = 2, jpjm1 |
|---|
| 273 | DO ji = fs_2, fs_jpim1 ! vector opt. |
|---|
| 274 | zadv_wav(ji,jj) = - ( zwx(ji,jj) - zwx(ji-1,jj) ) * r1_e1t(ji,jj) |
|---|
| 275 | zadv_wav(ji,jj) = zadv_wav(ji,jj) * MIN( 1._wp, ABS( ff_t(ji,jj) * z1_f20 ) ) ! tropical decrease |
|---|
| 276 | END DO |
|---|
| 277 | END DO |
|---|
| 278 | END IF |
|---|
| 279 | ! |
|---|
| 280 | !* divergence of diffusive fluxes |
|---|
| 281 | ! code adapted from traldf_lap_blp.F90 |
|---|
| 282 | IF ( rn_eke_lap >= 0._wp ) THEN |
|---|
| 283 | DO jj = 1, jpjm1 |
|---|
| 284 | DO ji = 1, fs_jpim1 ! vector opt. |
|---|
| 285 | zaheeu(ji,jj) = rn_eke_lap * e2_e1u(ji,jj) * umask(ji,jj,1) ! rn_eke_lap is constant (for now) and NOT masked |
|---|
| 286 | zaheev(ji,jj) = rn_eke_lap * e1_e2v(ji,jj) * vmask(ji,jj,1) ! before it is pahu and pahv which IS maked |
|---|
| 287 | END DO |
|---|
| 288 | END DO |
|---|
| 289 | DO jj = 1, jpjm1 !== First derivative (gradient) ==! |
|---|
| 290 | DO ji = 1, fs_jpim1 |
|---|
| 291 | ztu(ji,jj) = zaheeu(ji,jj) * ( eke_b(ji+1,jj ) - eke_b(ji,jj) ) |
|---|
| 292 | ztv(ji,jj) = zaheev(ji,jj) * ( eke_b(ji ,jj+1) - eke_b(ji,jj) ) |
|---|
| 293 | END DO |
|---|
| 294 | END DO |
|---|
| 295 | DO jj = 2, jpjm1 !== Second derivative (divergence), form trend ==! |
|---|
| 296 | DO ji = fs_2, fs_jpim1 |
|---|
| 297 | zlap(ji,jj) = ( ztu(ji,jj) - ztu(ji-1,jj ) & |
|---|
| 298 | & + ztv(ji,jj) - ztv(ji ,jj-1) ) & |
|---|
| 299 | & / e1e2t(ji,jj) |
|---|
| 300 | END DO |
|---|
| 301 | END DO |
|---|
| 302 | ELSE |
|---|
| 303 | zlap(:,:) = 0._wp |
|---|
| 304 | END IF |
|---|
| 305 | !* form the trend for linear dissipation |
|---|
| 306 | DO jj = 2, jpjm1 |
|---|
| 307 | DO ji = fs_2, fs_jpim1 |
|---|
| 308 | zdis(ji,jj) = - r1_ekedis(ji,jj) * (eke_b(ji,jj) - rn_eke_min) * tmask(ji,jj,1) |
|---|
| 309 | END DO |
|---|
| 310 | END DO |
|---|
| 311 | ! !== time stepping of EKE Eq. ==! |
|---|
| 312 | ! |
|---|
| 313 | ! note: the rn_eke_min term is a forcing in eddy equation, thus damping in mean equation |
|---|
| 314 | ! added to prevent overshoots and oscillations of energy from exponential growth/decay |
|---|
| 315 | ! maintains a background (depth-integrated) eddy energy level |
|---|
| 316 | ! |
|---|
| 317 | zeke_rhs = 0._wp |
|---|
| 318 | DO jj = 2, jpjm1 |
|---|
| 319 | DO ji = fs_2, fs_jpim1 ! vector opt. |
|---|
| 320 | ! |
|---|
| 321 | SELECT CASE( nn_eke_opt ) ! Specification of what to include in EKE budget |
|---|
| 322 | ! |
|---|
| 323 | CASE( 0 ) ! default: just PE->EKE growth and linear dissipation |
|---|
| 324 | zeke_rhs = zeke_peS(ji,jj) + zdis(ji,jj) + zlap(ji,jj) |
|---|
| 325 | CASE( 1 ) ! as default but with full advection |
|---|
| 326 | zeke_rhs = - zadv_ubt(ji,jj) + zadv_wav(ji,jj) + zeke_peS(ji,jj) + zdis(ji,jj) + zlap(ji,jj) |
|---|
| 327 | CASE( 2 ) ! full thing with additional KE->EKE growth |
|---|
| 328 | zeke_rhs = - zadv_ubt(ji,jj) + zadv_wav(ji,jj) + zeke_peS(ji,jj) + zdis(ji,jj) + zlap(ji,jj) + eke_keS(ji,jj) |
|---|
| 329 | CASE( 88 ) ! ONLY advection by mean flow |
|---|
| 330 | zeke_rhs = - zadv_ubt(ji,jj) |
|---|
| 331 | CASE( 99 ) ! ONLY diffusion |
|---|
| 332 | zeke_rhs = zlap(ji,jj) |
|---|
| 333 | CASE DEFAULT |
|---|
| 334 | CALL ctl_stop('ldf_eke: wrong choice nn_eke_opt, set at least to 0 (default)') |
|---|
| 335 | END SELECT |
|---|
| 336 | ! |
|---|
| 337 | ! leap-frog |
|---|
| 338 | eke_a(ji,jj) = eke_b(ji,jj) + r2dt * zeke_rhs * ssmask(ji,jj) |
|---|
| 339 | ! |
|---|
| 340 | END DO |
|---|
| 341 | END DO |
|---|
| 342 | CALL lbc_lnk( 'ldfeke', eke_a, 'T', 1. ) ! Lateral boundary conditions on zwi (unchanged sign) |
|---|
| 343 | |
|---|
| 344 | IF( neuler == 0 .AND. kt == nit000 ) THEN ! Euler time-stepping at first time-step (only swap) |
|---|
| 345 | eke_n(:,:) = eke_a(:,:) |
|---|
| 346 | ! |
|---|
| 347 | ELSE ! Leap-Frog + Asselin filter time stepping |
|---|
| 348 | DO jj = 1, jpj |
|---|
| 349 | DO ji = 1, jpi |
|---|
| 350 | zen = eke_n(ji,jj) |
|---|
| 351 | zed = eke_a(ji,jj) - 2._wp * zen + eke_b(ji,jj) ! time laplacian on tracers |
|---|
| 352 | ! |
|---|
| 353 | eke_b(ji,jj) = zen + atfp * zed ! eke_b <-- filtered eke_n |
|---|
| 354 | eke_n(ji,jj) = eke_a(ji,jj) ! ele_n <-- eke_a |
|---|
| 355 | END DO |
|---|
| 356 | END DO |
|---|
| 357 | ENDIF |
|---|
| 358 | |
|---|
| 359 | ! initialise it here so XIOS stops complaining... |
|---|
| 360 | zross(:,:) = 0._wp |
|---|
| 361 | zaeiw(:,:,:) = 0._wp |
|---|
| 362 | IF( l_eke_eiv ) THEN |
|---|
| 363 | ! !== resulting EIV coefficient ==! |
|---|
| 364 | ! |
|---|
| 365 | ! ! surface value |
|---|
| 366 | z1_f20 = 1._wp / ( 2._wp * omega * sin( rad * 20._wp ) ) |
|---|
| 367 | DO jj = 2, jpjm1 |
|---|
| 368 | DO ji = fs_2, fs_jpim1 ! vector opt. |
|---|
| 369 | ! Rossby radius at w-point, Ro = .5 * sum_jpk(N) / f |
|---|
| 370 | zfw = MAX( ABS( ff_t(ji,jj) ) , 1.e-10 ) |
|---|
| 371 | zross(ji,jj) = 0.5 * zn(ji,jj) / zfw |
|---|
| 372 | ! |
|---|
| 373 | zaeiw(ji,jj,1) = rn_geom * eke_b(ji,jj) / MAX( 1.e-10 , zn_slp(ji,jj) ) * tmask(ji,jj,1) ! zn_slp has zzSF multiplied to it |
|---|
| 374 | zaeiw(ji,jj,1) = MIN( rn_aeiv_max, zaeiw(ji,jj,1) ) ! bound aeiv from above |
|---|
| 375 | zaeiw(ji,jj,1) = zaeiw(ji,jj,1) & |
|---|
| 376 | ! tanh taper to deal with some some large values near coast |
|---|
| 377 | & * 0.5 * ( 1._wp - TANH( ( -ht_0(ji,jj) + 800._wp ) / 300._wp ) ) & |
|---|
| 378 | ! tanh taper of aeiv on internal Rossby radius |
|---|
| 379 | & * 0.5 * ( 1._wp + TANH( ( zross(ji,jj) - rn_ross_min ) / 2._wp ) ) |
|---|
| 380 | zaeiw(ji,jj,1) = zaeiw(ji,jj,1) * MIN( 1._wp, ABS( ff_t(ji,jj) * z1_f20 ) ) ! tropical decrease |
|---|
| 381 | zaeiw(ji,jj,1) = MAX( rn_aeiv_min, zaeiw(ji,jj,1) ) ! bound aeiv from below |
|---|
| 382 | END DO |
|---|
| 383 | END DO |
|---|
| 384 | CALL lbc_lnk( 'ldfeke', zaeiw(:,:,1), 'W', 1. ) |
|---|
| 385 | ! |
|---|
| 386 | ! ! inner value |
|---|
| 387 | ! set bottom to zero and use the un-masked zaeiw(ji,jj,1) first... |
|---|
| 388 | zaeiw(:,:,jpk) = 0._wp |
|---|
| 389 | DO jk = 2, jpkm1 |
|---|
| 390 | DO jj = 1, jpj |
|---|
| 391 | DO ji = 1, jpi |
|---|
| 392 | zaeiw(ji,jj,jk) = zaeiw(ji,jj,1) * zzSF(ji,jj,jk) * wmask(ji,jj,jk) ! zzSF has already been capped |
|---|
| 393 | END DO |
|---|
| 394 | END DO |
|---|
| 395 | END DO |
|---|
| 396 | ! ! aei at u- and v-points |
|---|
| 397 | aeiu(:,:,jpk) = 0._wp |
|---|
| 398 | aeiv(:,:,jpk) = 0._wp |
|---|
| 399 | DO jk = 1, jpkm1 |
|---|
| 400 | DO jj = 2, jpjm1 |
|---|
| 401 | DO ji = fs_2, fs_jpim1 ! vector opt. |
|---|
| 402 | aeiu(ji,jj,jk) = ( zaeiw(ji,jj,jk ) + zaeiw(ji+1,jj,jk ) & |
|---|
| 403 | & + zaeiw(ji,jj,jk+1) + zaeiw(ji+1,jj,jk+1) ) & |
|---|
| 404 | & / MAX( wmask(ji,jj,jk ) + wmask(ji+1,jj,jk ) & |
|---|
| 405 | & + wmask(ji,jj,jk+1) + wmask(ji+1,jj,jk+1) , 1._wp ) * umask(ji,jj,jk) |
|---|
| 406 | aeiv(ji,jj,jk) = ( zaeiw(ji,jj,jk ) + zaeiw(ji,jj+1,jk ) & |
|---|
| 407 | & + zaeiw(ji,jj,jk+1) + zaeiw(ji,jj+1,jk+1) ) & |
|---|
| 408 | & / MAX( wmask(ji,jj,jk ) + wmask(ji,jj+1,jk ) & |
|---|
| 409 | & + wmask(ji,jj,jk+1) + wmask(ji,jj+1,jk+1) , 1._wp ) * vmask(ji,jj,jk) |
|---|
| 410 | END DO |
|---|
| 411 | END DO |
|---|
| 412 | END DO |
|---|
| 413 | ! !-- diagnostics --! |
|---|
| 414 | CALL lbc_lnk_multi( 'ldfeke', aeiu , 'U', 1. , aeiv , 'V', 1. ) ! lateral boundary condition |
|---|
| 415 | CALL lbc_lnk( 'ldfeke', zross, 'W', 1. ) |
|---|
| 416 | END IF |
|---|
| 417 | ! |
|---|
| 418 | IF( lrst_oce ) CALL eke_rst( kt, 'WRITE' ) |
|---|
| 419 | |
|---|
| 420 | ! !== output EKE related variables ==! |
|---|
| 421 | CALL lbc_lnk_multi( 'ldfeke', & ! lateral boundary condition |
|---|
| 422 | & zadv_ubt, 'T', 1. , zadv_wav, 'T', 1. , zdis, 'T', 1. , & |
|---|
| 423 | & zlap , 'T', 1. , zeke_peS, 'T', 1. , & |
|---|
| 424 | & zc1 , 'W', 1. , zc_ros , 'W', 1. ) |
|---|
| 425 | ! |
|---|
| 426 | CALL iom_put( "eke" , eke_n ) ! parameterized total EKE (EPE+ EKE) |
|---|
| 427 | CALL iom_put( "trd_eke_adv_ubt", zadv_ubt ) ! ubt advective trend of EKE(LHS) |
|---|
| 428 | CALL iom_put( "trd_eke_adv_wav", zadv_wav ) ! rossby advective trend of EKE(LHS) |
|---|
| 429 | CALL iom_put( "trd_eke_dis" , zdis ) ! dissipative trend of EKE (RHS) |
|---|
| 430 | CALL iom_put( "trd_eke_lap" , zlap ) ! diffusive trend of EKE (RHS) |
|---|
| 431 | CALL iom_put( "trd_eke_peS" , zeke_peS ) ! PE to EKE source trend (RHS) |
|---|
| 432 | CALL iom_put( "trd_eke_keS" , eke_keS ) ! KE to EKE source trend (RHS) |
|---|
| 433 | CALL iom_put( "aeiv_geom" , zaeiw ) ! eddy induced coefficient from GEOMETRIC param |
|---|
| 434 | CALL iom_put( "rossby_rad" , zross ) ! internal Rossby deformation radius |
|---|
| 435 | CALL iom_put( "c1_vert" , zc1 ) ! 1st baroclinic mode phase speed |
|---|
| 436 | CALL iom_put( "c_ros" , zc_ros ) ! long Rossby phase speed |
|---|
| 437 | !! |
|---|
| 438 | !! jm: modified ldf_tra in ldftra.F90 to include the ln_eke_equ flag |
|---|
| 439 | !! into consideration, otherwise iom_put is called twice aeiu and aeiv |
|---|
| 440 | !! |
|---|
| 441 | CALL iom_put( "aeiu_2d", aeiu(:,:,1) ) ! surface u-EIV coeff. |
|---|
| 442 | CALL iom_put( "aeiv_2d", aeiv(:,:,1) ) ! surface v-EIV coeff. |
|---|
| 443 | CALL iom_put( "aeiu_3d", aeiu(:,:,:) ) ! 3D u-EIV coeff. |
|---|
| 444 | CALL iom_put( "aeiv_3d", aeiv(:,:,:) ) ! 3D v-EIV coeff. |
|---|
| 445 | ! |
|---|
| 446 | IF( ln_timing ) CALL timing_stop('ldf_eke') |
|---|
| 447 | ! |
|---|
| 448 | END SUBROUTINE ldf_eke |
|---|
| 449 | |
|---|
| 450 | SUBROUTINE ldf_eke_init |
|---|
| 451 | !!---------------------------------------------------------------------- |
|---|
| 452 | !! *** ROUTINE ldf_eke_init *** |
|---|
| 453 | !! |
|---|
| 454 | !! ** Purpose : Initialization of the depth-integrated eke, vertical |
|---|
| 455 | !! structure function and max/min caps of aeiw when using the |
|---|
| 456 | !! GEOMETRIC parameterization |
|---|
| 457 | !! |
|---|
| 458 | !! ** Method : Read the namldf_eke namelist and check the parameters |
|---|
| 459 | !! called at the first timestep (nit000) |
|---|
| 460 | !! |
|---|
| 461 | !! ** input : Namlist namldf_eke |
|---|
| 462 | !!---------------------------------------------------------------------- |
|---|
| 463 | INTEGER :: ios, inum |
|---|
| 464 | INTEGER :: ji, jj |
|---|
| 465 | !! |
|---|
| 466 | NAMELIST/namldf_eke/ rn_ekedis , nn_eke_dis , & ! GEOM master params (lambda and option, alpha, e0) |
|---|
| 467 | & rn_geom , rn_eke_init, & |
|---|
| 468 | & rn_eke_lap , & ! coeff of Laplacian diffusion of EKE |
|---|
| 469 | & rn_eke_min , & ! background value of (depth-integrated) EKE |
|---|
| 470 | & rn_ross_min, & ! taper aeiv based on Rossby internal radius |
|---|
| 471 | & rn_aeiv_max, rn_aeiv_min, & ! caps and options on result aeiv |
|---|
| 472 | & rn_SFmin, rn_SFmax, & ! caps on vertical structure function |
|---|
| 473 | & nn_eke_opt , & ! option for EKE budget terms |
|---|
| 474 | & ln_adv_wav , & ! option for advection by Rossby wave |
|---|
| 475 | & ln_beta_plane ! beta plane option for computing Rossby wave speed |
|---|
| 476 | !!---------------------------------------------------------------------- |
|---|
| 477 | ! |
|---|
| 478 | ! |
|---|
| 479 | ! note: GEOMETRIC is *** off *** by default (so no ref list) |
|---|
| 480 | ! REWIND( numnam_ref ) ! Namelist namldf_eke in reference namelist : Turbulent Kinetic Energy |
|---|
| 481 | ! READ ( numnam_ref, namldf_eke, IOSTAT = ios, ERR = 901) |
|---|
| 482 | !901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namldf_eke in reference namelist' ) |
|---|
| 483 | |
|---|
| 484 | REWIND( numnam_cfg ) ! Namelist namldf_eke in configuration namelist : GEOMETRIC parameterization |
|---|
| 485 | READ ( numnam_cfg, namldf_eke, IOSTAT = ios, ERR = 902 ) |
|---|
| 486 | 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namldf_eke in configuration namelist' ) |
|---|
| 487 | IF(lwm) WRITE ( numond, namldf_eke ) |
|---|
| 488 | ! |
|---|
| 489 | IF(lwp) THEN !* Control print |
|---|
| 490 | WRITE(numout,*) |
|---|
| 491 | WRITE(numout,*) 'ldf_eke_init : GEOMETRIC parameterization (total EKE time evolution equation)' |
|---|
| 492 | WRITE(numout,*) '~~~~~~~~~~~~' |
|---|
| 493 | WRITE(numout,*) ' Namelist namldf_eke : set the GEOMETRIC parameters' |
|---|
| 494 | WRITE(numout,*) ' aeiw updated according to GEOMETRIC ln_eke_eiv = ', l_eke_eiv |
|---|
| 495 | WRITE(numout,*) ' dissipation time scale of EKE in days rn_ekedis = ', rn_ekedis, ' days' |
|---|
| 496 | WRITE(numout,*) ' option for dissipation field nn_eke_dis = ', nn_eke_dis |
|---|
| 497 | SELECT CASE( nn_eke_dis ) |
|---|
| 498 | ! |
|---|
| 499 | CASE( 0 ) |
|---|
| 500 | WRITE(numout,*) ' spatially constant ' |
|---|
| 501 | CASE( 2 ) |
|---|
| 502 | WRITE(numout,*) ' read in a 2d varying file geom_diss_2D.nc ' |
|---|
| 503 | END SELECT |
|---|
| 504 | ! |
|---|
| 505 | WRITE(numout,*) ' geometric parameterization master coefficient rn_geom = ', rn_geom |
|---|
| 506 | WRITE(numout,*) ' coeff of Lapican diffusion of EKE rn_eke_lap = ', rn_eke_lap |
|---|
| 507 | WRITE(numout,*) ' initial total EKE value rn_eke_init = ', rn_eke_init |
|---|
| 508 | WRITE(numout,*) ' background values of total EKE rn_eke_min = ', rn_eke_min |
|---|
| 509 | WRITE(numout,*) ' taper aeiv subject to min Rossby radius rn_ross_min = ', rn_ross_min |
|---|
| 510 | WRITE(numout,*) ' maximum bound of aeiv coeff rn_aeiv_max = ', rn_aeiv_max |
|---|
| 511 | WRITE(numout,*) ' minimum bound of aeiv coeff rn_aeiv_min = ', rn_aeiv_min |
|---|
| 512 | WRITE(numout,*) ' minimum bound of Structure Function rn_SFmin = ', rn_SFmin |
|---|
| 513 | WRITE(numout,*) ' maximum bound of Structure Function rn_SFmax = ', rn_SFmax |
|---|
| 514 | WRITE(numout,*) ' [set rnSFmin = rnSFmax = 1 to use aeiv = aeiv(x,y,t)] ' |
|---|
| 515 | WRITE(numout,*) ' option for terms in eddy energy budget nn_eke_opt = ', nn_eke_opt |
|---|
| 516 | WRITE(numout,*) ' default: just PE->EKE growth and linear dissipation ' |
|---|
| 517 | SELECT CASE( nn_eke_opt ) |
|---|
| 518 | ! |
|---|
| 519 | CASE( 0 ) |
|---|
| 520 | WRITE(numout,*) ' default ' |
|---|
| 521 | CASE( 1 ) |
|---|
| 522 | WRITE(numout,*) ' default + advection ' |
|---|
| 523 | CASE( 2 ) |
|---|
| 524 | WRITE(numout,*) ' default + advection + KE->EKE ' |
|---|
| 525 | CASE( 88 ) |
|---|
| 526 | WRITE(numout,*) ' ONLY advection by z-avg mean flow (no growth or dissipation)' |
|---|
| 527 | CASE( 99 ) |
|---|
| 528 | WRITE(numout,*) ' ONLY Laplacian diffusion (no growth or dissipation)' |
|---|
| 529 | END SELECT |
|---|
| 530 | WRITE(numout,*) ' advection of energy by Rossby waves ln_adv_wav = ', ln_adv_wav |
|---|
| 531 | ENDIF |
|---|
| 532 | ! ! allocate eke arrays |
|---|
| 533 | ALLOCATE( eke_b(jpi,jpj) , eke_n(jpi,jpj) , eke_a(jpi,jpj) , eke_keS(jpi,jpj) , & |
|---|
| 534 | & r1_ekedis(jpi,jpj) , zc1(jpi,jpj) , zc_ros(jpi,jpj) , zadv_wav(jpi,jpj) ) |
|---|
| 535 | ! |
|---|
| 536 | SELECT CASE( nn_eke_dis ) ! Specification of linear dissipation |
|---|
| 537 | ! |
|---|
| 538 | CASE( 0 ) !== constant ==! |
|---|
| 539 | IF(lwp) WRITE(numout,*) ' linear EKE dissipation coef. = constant = ', rn_ekedis, ' days' |
|---|
| 540 | r1_ekedis(:,:) = 1._wp / (rn_ekedis * rday) |
|---|
| 541 | ! |
|---|
| 542 | CASE( -20 ) !== fixed horizontal shape read in file ==! |
|---|
| 543 | IF(lwp) WRITE(numout,*) ' linear EKE dissipation coef. = F(i,j) read in geom_diss_2D.nc file' |
|---|
| 544 | CALL iom_open ( 'geom_diss_2D.nc', inum ) |
|---|
| 545 | CALL iom_get ( inum, jpdom_data, 'r1_ekedis', r1_ekedis ) ! read in as time-scale in days... |
|---|
| 546 | CALL iom_close( inum ) |
|---|
| 547 | DO jj = 1, jpj |
|---|
| 548 | DO ji = 1, jpi |
|---|
| 549 | r1_ekedis(ji,jj) = 1._wp / (r1_ekedis(ji,jj) * rday) ! ...convert rate in per second |
|---|
| 550 | END DO |
|---|
| 551 | END DO |
|---|
| 552 | ! |
|---|
| 553 | CASE DEFAULT |
|---|
| 554 | CALL ctl_stop('ldf_eke_init: wrong choice for nn_eke_dis, the option for linear dissipation in GEOMETRIC') |
|---|
| 555 | END SELECT |
|---|
| 556 | ! |
|---|
| 557 | CALL eke_rst( nit000, 'READ' ) !* read or initialize all required files |
|---|
| 558 | ! |
|---|
| 559 | ! tidy up the initilisation of zc1 |
|---|
| 560 | zc_ros(:,:) = 0._wp |
|---|
| 561 | zadv_wav(:,:) = 0._wp |
|---|
| 562 | IF ( ln_adv_wav ) THEN |
|---|
| 563 | DO jj = 2, jpjm1 |
|---|
| 564 | DO ji = 1, fs_jpim1 |
|---|
| 565 | ! set the ones that are not deep enough to be zero |
|---|
| 566 | IF ( ht_0(ji,jj) < 300._wp ) THEN ! jm: should use something like ht_b really... |
|---|
| 567 | zc1(ji,jj) = 0._wp |
|---|
| 568 | END IF |
|---|
| 569 | END DO |
|---|
| 570 | END DO |
|---|
| 571 | END IF |
|---|
| 572 | ! if using beta_plane, compute beta and f0 locally |
|---|
| 573 | IF ( ln_beta_plane ) THEN |
|---|
| 574 | zf0 = ff_f(1,1) |
|---|
| 575 | zbeta = (ff_f(1,2) - zf0) / e2t(1,1) |
|---|
| 576 | IF(lwp) WRITE(numout,*) ' beta plane option is ln_beta_plane =', ln_beta_plane |
|---|
| 577 | IF(lwp) WRITE(numout,*) ' f0 = ', zf0, ' s-1' |
|---|
| 578 | IF(lwp) WRITE(numout,*) ' beta = ', zbeta, 'm-1 s-1' |
|---|
| 579 | END IF |
|---|
| 580 | ! |
|---|
| 581 | ! |
|---|
| 582 | END SUBROUTINE ldf_eke_init |
|---|
| 583 | |
|---|
| 584 | |
|---|
| 585 | SUBROUTINE eke_rst( kt, cdrw ) |
|---|
| 586 | !!--------------------------------------------------------------------- |
|---|
| 587 | !! *** ROUTINE eke_rst *** |
|---|
| 588 | !! |
|---|
| 589 | !! ** Purpose : Read or write EKE file (eke) in restart file |
|---|
| 590 | !! |
|---|
| 591 | !! ** Method : use of IOM library |
|---|
| 592 | !! if the restart does not contain EKE, eke is set |
|---|
| 593 | !! according to rn_eke_init, and aeiu = aeiv = 10 m s^-2 |
|---|
| 594 | !!---------------------------------------------------------------------- |
|---|
| 595 | INTEGER , INTENT(in) :: kt ! ocean time-step |
|---|
| 596 | CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag |
|---|
| 597 | ! |
|---|
| 598 | INTEGER :: jit, jk ! dummy loop indices |
|---|
| 599 | INTEGER :: id1, id2, id3, id4, id5 ! local integer |
|---|
| 600 | !!---------------------------------------------------------------------- |
|---|
| 601 | ! |
|---|
| 602 | IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise |
|---|
| 603 | ! ! --------------- |
|---|
| 604 | IF( ln_rstart ) THEN !* Read the restart file |
|---|
| 605 | id1 = iom_varid( numror, 'eke_b' , ldstop = .FALSE. ) |
|---|
| 606 | id2 = iom_varid( numror, 'eke_n' , ldstop = .FALSE. ) |
|---|
| 607 | id3 = iom_varid( numror, 'aeiu' , ldstop = .FALSE. ) |
|---|
| 608 | id4 = iom_varid( numror, 'aeiv' , ldstop = .FALSE. ) |
|---|
| 609 | IF(lwp) write(numout,*) ' ldfeke init restart' |
|---|
| 610 | ! |
|---|
| 611 | IF( MIN( id1, id2, id3, id4 ) > 0 ) THEN ! all required arrays exist |
|---|
| 612 | IF(lwp) write(numout,*) ' all files for ldfeke exist, loading ...' |
|---|
| 613 | CALL iom_get( numror, jpdom_autoglo, 'eke_b', eke_b ) |
|---|
| 614 | CALL iom_get( numror, jpdom_autoglo, 'eke_n', eke_n ) |
|---|
| 615 | CALL iom_get( numror, jpdom_autoglo, 'aeiu' , aeiu ) |
|---|
| 616 | CALL iom_get( numror, jpdom_autoglo, 'aeiv' , aeiv ) |
|---|
| 617 | ELSE ! one at least array is missing |
|---|
| 618 | IF(lwp) write(numout,*) ' not all files for ldfeke exist ' |
|---|
| 619 | IF(lwp) write(numout,*) ' --- initialize from namelist' |
|---|
| 620 | eke_b(:,:) = rn_eke_init * ht_n(:,:) * ssmask(:,:) |
|---|
| 621 | eke_n(:,:) = eke_b(:,:) |
|---|
| 622 | aeiu(:,:,:) = 10._wp * umask(:,:,:) ! bottom eddy coeff set to zero at last level |
|---|
| 623 | aeiv(:,:,:) = 10._wp * vmask(:,:,:) |
|---|
| 624 | ENDIF |
|---|
| 625 | ! ln_wav_adv restarts + initialisations |
|---|
| 626 | IF( ln_adv_wav ) THEN |
|---|
| 627 | ! id5 = iom_varid( numror, 'zc1' , ldstop = .FALSE. ) |
|---|
| 628 | ! IF( id5 > 0 ) THEN |
|---|
| 629 | ! IF(lwp) write(numout,*) ' ln_adv_wav on and found zc1 variable, loading ...' |
|---|
| 630 | ! CALL iom_get( numror, jpdom_autoglo, 'zc1' , zc1 ) |
|---|
| 631 | ! ELSE |
|---|
| 632 | IF(lwp) write(numout,*) ' ln_adv_wav on but not found zc1 variable, initialising ...' |
|---|
| 633 | zc1(:,:) = 4._wp ! set it to something everywhere but some are set to zero in ldf_eke_init |
|---|
| 634 | ! END IF |
|---|
| 635 | END IF |
|---|
| 636 | ELSE !* Start from rest with a non zero value (required) |
|---|
| 637 | IF(lwp) write(numout,*) ' ldfeke init restart from namelist' |
|---|
| 638 | eke_b(:,:) = rn_eke_init * ht_n(:,:) * ssmask(:,:) |
|---|
| 639 | eke_n(:,:) = eke_b(:,:) |
|---|
| 640 | aeiu(:,:,:) = 10._wp * umask(:,:,:) ! bottom eddy coeff set to zero at last level |
|---|
| 641 | aeiv(:,:,:) = 10._wp * vmask(:,:,:) |
|---|
| 642 | zc1(:,:) = 4._wp |
|---|
| 643 | ! |
|---|
| 644 | ENDIF |
|---|
| 645 | ! |
|---|
| 646 | ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN ! Create restart file |
|---|
| 647 | ! ! ------------------- |
|---|
| 648 | IF(lwp) WRITE(numout,*) '---- eke-rst ----' |
|---|
| 649 | CALL iom_rstput( kt, nitrst, numrow, 'eke_b', eke_b ) |
|---|
| 650 | CALL iom_rstput( kt, nitrst, numrow, 'eke_n', eke_n ) |
|---|
| 651 | CALL iom_rstput( kt, nitrst, numrow, 'aeiu' , aeiu ) |
|---|
| 652 | CALL iom_rstput( kt, nitrst, numrow, 'aeiv' , aeiv ) |
|---|
| 653 | CALL iom_rstput( kt, nitrst, numrow, 'zc1' , zc1 ) |
|---|
| 654 | ! |
|---|
| 655 | ENDIF |
|---|
| 656 | ! |
|---|
| 657 | END SUBROUTINE eke_rst |
|---|
| 658 | |
|---|
| 659 | !!====================================================================== |
|---|
| 660 | END MODULE ldfeke |
|---|
