Changeset 7646 for trunk/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90
- Timestamp:
- 2017-02-06T10:25:03+01:00 (7 years ago)
- File:
-
- 1 edited
-
trunk/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90 (modified) (8 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90
r5836 r7646 17 17 USE phycst ! physical constants 18 18 USE dom_oce ! ocean space and time domain 19 USE sbc_oce ! Surface boundary condition: ocean fields20 19 USE sbc_ice ! Surface boundary condition: ice fields 21 20 USE ice ! LIM-3 variables 22 USE dom_ice ! LIM-3 domain23 21 USE limrhg ! LIM-3 rheology 24 22 USE lbclnk ! lateral boundary conditions - MPP exchanges … … 26 24 USE wrk_nemo ! work arrays 27 25 USE in_out_manager ! I/O manager 28 USE prtctl ! Print control29 26 USE lib_fortran ! glob_sum 30 USE timing ! Timing 31 USE limcons ! conservation tests 27 USE timing ! Timing 28 USE limcons ! conservation tests 29 USE limctl ! control prints 32 30 USE limvar 33 31 … … 35 33 PRIVATE 36 34 37 PUBLIC lim_dyn ! routine called by ice_step 35 PUBLIC lim_dyn ! routine called by sbcice_lim.F90 36 PUBLIC lim_dyn_init ! routine called by sbcice_lim.F90 38 37 39 38 !! * Substitutions … … 50 49 !! *** ROUTINE lim_dyn *** 51 50 !! 52 !! ** Purpose : compute ice velocity and ocean-ice stress51 !! ** Purpose : compute ice velocity 53 52 !! 54 53 !! ** Method : … … 56 55 !! ** Action : - Initialisation 57 56 !! - Call of the dynamic routine for each hemisphere 58 !! - computation of the stress at the ocean surface59 !! - treatment of the case if no ice dynamic60 57 !!------------------------------------------------------------------------------------ 61 58 INTEGER, INTENT(in) :: kt ! number of iteration 62 59 !! 63 INTEGER :: ji, jj, jl, ja ! dummy loop indices 64 INTEGER :: i_j1, i_jpj ! Starting/ending j-indices for rheology 65 REAL(wp) :: zcoef ! local scalar 66 REAL(wp), POINTER, DIMENSION(:) :: zswitch ! i-averaged indicator of sea-ice 67 REAL(wp), POINTER, DIMENSION(:) :: zmsk ! i-averaged of tmask 68 REAL(wp), POINTER, DIMENSION(:,:) :: zu_io, zv_io ! ice-ocean velocity 69 ! 60 INTEGER :: jl, jk ! dummy loop indices 70 61 REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 71 62 !!--------------------------------------------------------------------- … … 73 64 IF( nn_timing == 1 ) CALL timing_start('limdyn') 74 65 75 CALL wrk_alloc( jpi, jpj, zu_io, zv_io ) 76 CALL wrk_alloc( jpj, zswitch, zmsk ) 77 78 CALL lim_var_agg(1) ! aggregate ice categories 79 80 IF( kt == nit000 ) CALL lim_dyn_init ! Initialization (first time-step only) 81 82 IF( ln_limdyn ) THEN 83 ! 84 ! conservation test 85 IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limdyn', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 86 87 u_ice_b(:,:) = u_ice(:,:) * umask(:,:,1) 88 v_ice_b(:,:) = v_ice(:,:) * vmask(:,:,1) 89 90 ! Rheology (ice dynamics) 91 ! ======== 92 93 ! Define the j-limits where ice rheology is computed 94 ! --------------------------------------------------- 95 96 IF( lk_mpp .OR. lk_mpp_rep ) THEN ! mpp: compute over the whole domain 97 i_j1 = 1 98 i_jpj = jpj 99 IF(ln_ctl) CALL prt_ctl_info( 'lim_dyn : i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj ) 100 CALL lim_rhg( i_j1, i_jpj ) 101 ELSE ! optimization of the computational area 102 ! 103 DO jj = 1, jpj 104 zswitch(jj) = SUM( 1.0 - at_i(:,jj) ) ! = REAL(jpj) if ocean everywhere on a j-line 105 zmsk (jj) = SUM( tmask(:,jj,1) ) ! = 0 if land everywhere on a j-line 106 END DO 107 108 IF( l_jeq ) THEN ! local domain include both hemisphere 109 ! ! Rheology is computed in each hemisphere 110 ! ! only over the ice cover latitude strip 111 ! Northern hemisphere 112 i_j1 = njeq 113 i_jpj = jpj 114 DO WHILE ( i_j1 <= jpj .AND. zswitch(i_j1) == FLOAT(jpi) .AND. zmsk(i_j1) /=0 ) 115 i_j1 = i_j1 + 1 116 END DO 117 i_j1 = MAX( 1, i_j1-2 ) 118 IF(ln_ctl) CALL prt_ctl_info( 'lim_dyn : NH i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj ) 119 CALL lim_rhg( i_j1, i_jpj ) 120 ! 121 ! Southern hemisphere 122 i_j1 = 1 123 i_jpj = njeq 124 DO WHILE ( i_jpj >= 1 .AND. zswitch(i_jpj) == FLOAT(jpi) .AND. zmsk(i_jpj) /=0 ) 125 i_jpj = i_jpj - 1 126 END DO 127 i_jpj = MIN( jpj, i_jpj+1 ) 128 IF(ln_ctl) CALL prt_ctl_info( 'lim_dyn : SH i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj ) 129 ! 130 CALL lim_rhg( i_j1, i_jpj ) 131 ! 132 ELSE ! local domain extends over one hemisphere only 133 ! ! Rheology is computed only over the ice cover 134 ! ! latitude strip 135 i_j1 = 1 136 DO WHILE ( i_j1 <= jpj .AND. zswitch(i_j1) == FLOAT(jpi) .AND. zmsk(i_j1) /=0 ) 137 i_j1 = i_j1 + 1 138 END DO 139 i_j1 = MAX( 1, i_j1-2 ) 140 141 i_jpj = jpj 142 DO WHILE ( i_jpj >= 1 .AND. zswitch(i_jpj) == FLOAT(jpi) .AND. zmsk(i_jpj) /=0 ) 143 i_jpj = i_jpj - 1 144 END DO 145 i_jpj = MIN( jpj, i_jpj+1) 146 ! 147 IF(ln_ctl) CALL prt_ctl_info( 'lim_dyn : one hemisphere: i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj ) 148 ! 149 CALL lim_rhg( i_j1, i_jpj ) 150 ! 151 ENDIF 152 ! 153 ENDIF 154 155 ! computation of friction velocity 156 ! -------------------------------- 157 ! ice-ocean velocity at U & V-points (u_ice v_ice at U- & V-points ; ssu_m, ssv_m at U- & V-points) 158 zu_io(:,:) = u_ice(:,:) - ssu_m(:,:) 159 zv_io(:,:) = v_ice(:,:) - ssv_m(:,:) 160 ! frictional velocity at T-point 161 zcoef = 0.5_wp * rn_cio 162 DO jj = 2, jpjm1 163 DO ji = fs_2, fs_jpim1 ! vector opt. 164 ust2s(ji,jj) = zcoef * ( zu_io(ji,jj) * zu_io(ji,jj) + zu_io(ji-1,jj) * zu_io(ji-1,jj) & 165 & + zv_io(ji,jj) * zv_io(ji,jj) + zv_io(ji,jj-1) * zv_io(ji,jj-1) ) * tmask(ji,jj,1) 166 END DO 167 END DO 168 ! 169 ! conservation test 170 IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limdyn', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 171 ! 172 ELSE ! no ice dynamics : transmit directly the atmospheric stress to the ocean 173 ! 174 zcoef = SQRT( 0.5_wp ) * r1_rau0 175 DO jj = 2, jpjm1 176 DO ji = fs_2, fs_jpim1 ! vector opt. 177 ust2s(ji,jj) = zcoef * SQRT( utau(ji,jj) * utau(ji,jj) + utau(ji-1,jj) * utau(ji-1,jj) & 178 & + vtau(ji,jj) * vtau(ji,jj) + vtau(ji,jj-1) * vtau(ji,jj-1) ) * tmask(ji,jj,1) 179 END DO 180 END DO 181 ! 182 ENDIF 183 CALL lbc_lnk( ust2s, 'T', 1. ) ! T-point 184 185 IF(ln_ctl) THEN ! Control print 186 CALL prt_ctl_info(' ') 187 CALL prt_ctl_info(' - Cell values : ') 188 CALL prt_ctl_info(' ~~~~~~~~~~~~~ ') 189 CALL prt_ctl(tab2d_1=ust2s , clinfo1=' lim_dyn : ust2s :') 190 CALL prt_ctl(tab2d_1=divu_i , clinfo1=' lim_dyn : divu_i :') 191 CALL prt_ctl(tab2d_1=delta_i , clinfo1=' lim_dyn : delta_i :') 192 CALL prt_ctl(tab2d_1=strength , clinfo1=' lim_dyn : strength :') 193 CALL prt_ctl(tab2d_1=e1e2t , clinfo1=' lim_dyn : cell area :') 194 CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_dyn : at_i :') 195 CALL prt_ctl(tab2d_1=vt_i , clinfo1=' lim_dyn : vt_i :') 196 CALL prt_ctl(tab2d_1=vt_s , clinfo1=' lim_dyn : vt_s :') 197 CALL prt_ctl(tab2d_1=stress1_i , clinfo1=' lim_dyn : stress1_i :') 198 CALL prt_ctl(tab2d_1=stress2_i , clinfo1=' lim_dyn : stress2_i :') 199 CALL prt_ctl(tab2d_1=stress12_i, clinfo1=' lim_dyn : stress12_i:') 66 CALL lim_var_agg(1) ! aggregate ice categories 67 ! 68 ! conservation test 69 IF( ln_limdiachk ) CALL lim_cons_hsm(0, 'limdyn', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 70 71 ! ice velocities before rheology 72 u_ice_b(:,:) = u_ice(:,:) * umask(:,:,1) 73 v_ice_b(:,:) = v_ice(:,:) * vmask(:,:,1) 74 75 ! Landfast ice parameterization: define max bottom friction 76 tau_icebfr(:,:) = 0._wp 77 IF( ln_landfast ) THEN 200 78 DO jl = 1, jpl 201 CALL prt_ctl_info(' ') 202 CALL prt_ctl_info(' - Category : ', ivar1=jl) 203 CALL prt_ctl_info(' ~~~~~~~~~~') 204 CALL prt_ctl(tab2d_1=a_i (:,:,jl) , clinfo1= ' lim_dyn : a_i : ') 205 CALL prt_ctl(tab2d_1=ht_i (:,:,jl) , clinfo1= ' lim_dyn : ht_i : ') 206 CALL prt_ctl(tab2d_1=ht_s (:,:,jl) , clinfo1= ' lim_dyn : ht_s : ') 207 CALL prt_ctl(tab2d_1=v_i (:,:,jl) , clinfo1= ' lim_dyn : v_i : ') 208 CALL prt_ctl(tab2d_1=v_s (:,:,jl) , clinfo1= ' lim_dyn : v_s : ') 209 CALL prt_ctl(tab2d_1=e_s (:,:,1,jl) , clinfo1= ' lim_dyn : e_s : ') 210 CALL prt_ctl(tab2d_1=t_su (:,:,jl) , clinfo1= ' lim_dyn : t_su : ') 211 CALL prt_ctl(tab2d_1=t_s (:,:,1,jl) , clinfo1= ' lim_dyn : t_snow : ') 212 CALL prt_ctl(tab2d_1=sm_i (:,:,jl) , clinfo1= ' lim_dyn : sm_i : ') 213 CALL prt_ctl(tab2d_1=smv_i (:,:,jl) , clinfo1= ' lim_dyn : smv_i : ') 214 DO ja = 1, nlay_i 215 CALL prt_ctl_info(' ') 216 CALL prt_ctl_info(' - Layer : ', ivar1=ja) 217 CALL prt_ctl_info(' ~~~~~~~') 218 CALL prt_ctl(tab2d_1=t_i(:,:,ja,jl) , clinfo1= ' lim_dyn : t_i : ') 219 CALL prt_ctl(tab2d_1=e_i(:,:,ja,jl) , clinfo1= ' lim_dyn : e_i : ') 220 END DO 79 WHERE( ht_i(:,:,jl) > ht_n(:,:) * rn_gamma ) tau_icebfr(:,:) = tau_icebfr(:,:) + a_i(:,:,jl) * rn_icebfr 221 80 END DO 222 81 ENDIF 82 83 ! Rheology (ice dynamics) 84 ! ======== 85 CALL lim_rhg 223 86 ! 224 CALL wrk_dealloc( jpi, jpj, zu_io, zv_io ) 225 CALL wrk_dealloc( jpj, zswitch, zmsk ) 87 ! conservation test 88 IF( ln_limdiachk ) CALL lim_cons_hsm(1, 'limdyn', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 89 90 ! Control prints 91 IF( ln_ctl ) CALL lim_prt3D( 'limdyn' ) 226 92 ! 227 93 IF( nn_timing == 1 ) CALL timing_stop('limdyn') … … 243 109 !!------------------------------------------------------------------- 244 110 INTEGER :: ios ! Local integer output status for namelist read 245 NAMELIST/namicedyn/ nn_icestr, ln_icestr_bvf, rn_pe_rdg, rn_pstar, rn_crhg, rn_cio, rn_creepl, rn_ecc, & 246 & nn_nevp, rn_relast, nn_ahi0, rn_ahi0_ref 247 INTEGER :: ji, jj 248 REAL(wp) :: za00, zd_max 111 NAMELIST/namicedyn/ nn_limadv, nn_limadv_ord, & 112 & nn_icestr, ln_icestr_bvf, rn_pe_rdg, rn_pstar, rn_crhg, rn_cio, rn_creepl, rn_ecc, & 113 & nn_nevp, rn_relast, ln_landfast, rn_gamma, rn_icebfr, rn_lfrelax 249 114 !!------------------------------------------------------------------- 250 115 … … 262 127 WRITE(numout,*) 'lim_dyn_init : ice parameters for ice dynamics ' 263 128 WRITE(numout,*) '~~~~~~~~~~~~' 264 WRITE(numout,*)' ice strength parameterization (0=Hibler 1=Rothrock) nn_icestr = ', nn_icestr 265 WRITE(numout,*)' Including brine volume in ice strength comp. ln_icestr_bvf = ', ln_icestr_bvf 266 WRITE(numout,*)' Ratio of ridging work to PotEner change in ridging rn_pe_rdg = ', rn_pe_rdg 267 WRITE(numout,*) ' drag coefficient for oceanic stress rn_cio = ', rn_cio 268 WRITE(numout,*) ' first bulk-rheology parameter rn_pstar = ', rn_pstar 269 WRITE(numout,*) ' second bulk-rhelogy parameter rn_crhg = ', rn_crhg 270 WRITE(numout,*) ' creep limit rn_creepl = ', rn_creepl 271 WRITE(numout,*) ' eccentricity of the elliptical yield curve rn_ecc = ', rn_ecc 272 WRITE(numout,*) ' number of iterations for subcycling nn_nevp = ', nn_nevp 273 WRITE(numout,*) ' ratio of elastic timescale over ice time step rn_relast = ', rn_relast 274 WRITE(numout,*) ' horizontal diffusivity calculation nn_ahi0 = ', nn_ahi0 275 WRITE(numout,*) ' horizontal diffusivity coeff. (orca2 grid) rn_ahi0_ref = ', rn_ahi0_ref 129 ! limtrp 130 WRITE(numout,*)' choose the advection scheme (-1=Prather, 0=Ulimate-Macho) nn_limadv = ', nn_limadv 131 WRITE(numout,*)' choose the order of the scheme (if ultimate) nn_limadv_ord = ', nn_limadv_ord 132 ! limrhg 133 WRITE(numout,*)' ice strength parameterization (0=Hibler 1=Rothrock) nn_icestr = ', nn_icestr 134 WRITE(numout,*)' Including brine volume in ice strength comp. ln_icestr_bvf = ', ln_icestr_bvf 135 WRITE(numout,*)' Ratio of ridging work to PotEner change in ridging rn_pe_rdg = ', rn_pe_rdg 136 WRITE(numout,*) ' drag coefficient for oceanic stress rn_cio = ', rn_cio 137 WRITE(numout,*) ' first bulk-rheology parameter rn_pstar = ', rn_pstar 138 WRITE(numout,*) ' second bulk-rhelogy parameter rn_crhg = ', rn_crhg 139 WRITE(numout,*) ' creep limit rn_creepl = ', rn_creepl 140 WRITE(numout,*) ' eccentricity of the elliptical yield curve rn_ecc = ', rn_ecc 141 WRITE(numout,*) ' number of iterations for subcycling nn_nevp = ', nn_nevp 142 WRITE(numout,*) ' ratio of elastic timescale over ice time step rn_relast = ', rn_relast 143 WRITE(numout,*) ' Landfast: param (T or F) ln_landfast = ', ln_landfast 144 WRITE(numout,*) ' Landfast: fraction of ocean depth that ice must reach rn_gamma = ', rn_gamma 145 WRITE(numout,*) ' Landfast: maximum bottom stress per unit area of contact rn_icebfr = ', rn_icebfr 146 WRITE(numout,*) ' Landfast: relax time scale (s-1) to reach static friction rn_lfrelax = ', rn_lfrelax 276 147 ENDIF 277 148 ! 278 usecc2 = 1._wp / ( rn_ecc * rn_ecc )279 rhoco = rau0 * rn_cio280 !281 ! Diffusion coefficients282 SELECT CASE( nn_ahi0 )283 284 CASE( 0 )285 ahiu(:,:) = rn_ahi0_ref286 ahiv(:,:) = rn_ahi0_ref287 288 IF(lwp) WRITE(numout,*) ''289 IF(lwp) WRITE(numout,*) ' laplacian operator: ahim constant = rn_ahi0_ref'290 291 CASE( 1 )292 293 zd_max = MAX( MAXVAL( e1t(:,:) ), MAXVAL( e2t(:,:) ) )294 IF( lk_mpp ) CALL mpp_max( zd_max ) ! max over the global domain295 296 ahiu(:,:) = rn_ahi0_ref * zd_max * 1.e-05_wp ! 1.e05 = 100km = max grid space at 60° latitude in orca2297 ! (60° = min latitude for ice cover)298 ahiv(:,:) = rn_ahi0_ref * zd_max * 1.e-05_wp299 300 IF(lwp) WRITE(numout,*) ''301 IF(lwp) WRITE(numout,*) ' laplacian operator: ahim proportional to max of e1 e2 over the domain (', zd_max, ')'302 IF(lwp) WRITE(numout,*) ' value for ahim = ', rn_ahi0_ref * zd_max * 1.e-05_wp303 304 CASE( 2 )305 306 zd_max = MAX( MAXVAL( e1t(:,:) ), MAXVAL( e2t(:,:) ) )307 IF( lk_mpp ) CALL mpp_max( zd_max ) ! max over the global domain308 309 za00 = rn_ahi0_ref * 1.e-05_wp ! 1.e05 = 100km = max grid space at 60° latitude in orca2310 ! (60° = min latitude for ice cover)311 DO jj = 1, jpj312 DO ji = 1, jpi313 ahiu(ji,jj) = za00 * MAX( e1t(ji,jj), e2t(ji,jj) ) * umask(ji,jj,1)314 ahiv(ji,jj) = za00 * MAX( e1f(ji,jj), e2f(ji,jj) ) * vmask(ji,jj,1)315 END DO316 END DO317 !318 IF(lwp) WRITE(numout,*) ''319 IF(lwp) WRITE(numout,*) ' laplacian operator: ahim proportional to e1'320 IF(lwp) WRITE(numout,*) ' maximum grid-spacing = ', zd_max, ' maximum value for ahim = ', za00*zd_max321 322 END SELECT323 324 149 END SUBROUTINE lim_dyn_init 325 150
Note: See TracChangeset
for help on using the changeset viewer.
