Changeset 2715 for trunk/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90
- Timestamp:
- 2011-03-30T17:58:35+02:00 (13 years ago)
- File:
-
- 1 edited
-
trunk/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90 (modified) (15 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90
r2528 r2715 4 4 !! LIM transport ice model : sea-ice advection/diffusion 5 5 !!====================================================================== 6 !! History : LIM-2 ! 2000-01 (M.A. Morales Maqueda, H. Goosse, and T. Fichefet) Original code 7 !! 3.0 ! 2005-11 (M. Vancoppenolle) Multi-layer sea ice, salinity variations 8 !! 4.0 ! 2011-02 (G. Madec) dynamical allocation 9 !!---------------------------------------------------------------------- 6 10 #if defined key_lim3 7 11 !!---------------------------------------------------------------------- … … 9 13 !!---------------------------------------------------------------------- 10 14 !! lim_trp : advection/diffusion process of sea ice 11 !! lim_trp_init : initialization and namelist read 12 !!---------------------------------------------------------------------- 13 USE phycst 14 USE dom_oce 15 !!---------------------------------------------------------------------- 16 USE phycst ! physical constant 17 USE dom_oce ! ocean domain 18 USE sbc_oce ! ocean surface boundary condition 19 USE par_ice ! LIM-3 parameter 20 USE dom_ice ! LIM-3 domain 21 USE ice ! LIM-3 variables 22 USE limadv ! LIM-3 advection 23 USE limhdf ! LIM-3 horizontal diffusion 15 24 USE in_out_manager ! I/O manager 16 USE sbc_oce ! Surface boundary condition: ocean fields 17 USE dom_ice 18 USE ice 19 USE limadv 20 USE limhdf 21 USE lbclnk 22 USE lib_mpp 23 USE par_ice 25 USE lbclnk ! lateral boundary conditions -- MPP exchanges 26 USE lib_mpp ! MPP library 24 27 USE prtctl ! Print control 25 28 … … 27 30 PRIVATE 28 31 29 !! * Routine accessibility 30 PUBLIC lim_trp ! called by ice_step 31 32 !! * Shared module variables 33 REAL(wp), PUBLIC :: & !: 34 bound = 0.e0 !: boundary condit. (0.0 no-slip, 1.0 free-slip) 35 36 !! * Module variables 37 REAL(wp) :: & ! constant values 38 epsi06 = 1.e-06 , & 39 epsi03 = 1.e-03 , & 40 epsi16 = 1.e-16 , & 41 rzero = 0.e0 , & 42 rone = 1.e0 , & 43 zeps10 = 1.e-10 32 PUBLIC lim_trp ! called by ice_step 33 34 REAL(wp), PUBLIC :: bound = 0._wp !: boundary condit. (0.0 no-slip, 1.0 free-slip) 35 36 REAL(wp) :: epsi06 = 1.e-06_wp ! constant values 37 REAL(wp) :: epsi03 = 1.e-03_wp 38 REAL(wp) :: zeps10 = 1.e-10_wp 39 REAL(wp) :: epsi16 = 1.e-16_wp 40 REAL(wp) :: rzero = 0._wp 41 REAL(wp) :: rone = 1._wp 44 42 45 43 !! * Substitution 46 44 # include "vectopt_loop_substitute.h90" 47 45 !!---------------------------------------------------------------------- 48 !! NEMO/LIM3 3.3 , UCL - NEMO Consortium (2010)46 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 49 47 !! $Id$ 50 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 51 !!---------------------------------------------------------------------- 52 48 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 49 !!---------------------------------------------------------------------- 53 50 CONTAINS 54 51 … … 64 61 !! 65 62 !! ** action : 66 !!67 !! History :68 !! 1.0 ! 00-01 (M.A. Morales Maqueda, H. Goosse, and T. Fichefet) Original code69 !! ! 01-05 (G. Madec, R. Hordoir) opa norm70 !! 2.0 ! 04-01 (G. Madec, C. Ethe) F90, mpp71 !! 3.0 ! 05-11 (M. Vancoppenolle) Multi-layer sea ice, salinity variations72 63 !!--------------------------------------------------------------------- 73 INTEGER, INTENT(in) :: kt ! number of iteration 74 !! * Local Variables 75 INTEGER :: ji, jj, jk, jl, layer, & ! dummy loop indices 76 initad ! number of sub-timestep for the advection 77 INTEGER :: ji_maxu, ji_maxv, jj_maxu, jj_maxv 78 79 REAL(wp) :: & 80 zindb , & 81 zindsn , & 82 zindic , & 83 zusvosn, & 84 zusvoic, & 85 zvbord , & 86 zcfl , & 87 zusnit , & 88 zrtt, zsal, zage, & 89 zbigval, ze, & 90 zmaxu, zmaxv 91 92 REAL(wp), DIMENSION(jpi,jpj) :: & ! temporary workspace 93 zui_u , zvi_v , zsm , & 94 zs0at, zs0ow 95 96 REAL(wp), DIMENSION(jpi,jpj,jpl):: & ! temporary workspace 97 zs0ice, zs0sn, zs0a , & 98 zs0c0 , & 99 zs0sm , zs0oi 100 101 ! MHE Multilayer heat content 102 REAL(wp), DIMENSION(jpi,jpj,jkmax,jpl) :: & ! temporary workspace 103 zs0e 104 105 !--------------------------------------------------------------------- 106 107 IF( numit == nstart ) CALL lim_trp_init ! Initialization (first time-step only) 108 64 USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released 65 USE wrk_nemo, ONLY: zs0at => wrk_2d_1 , zsm => wrk_2d_2 , zs0ow => wrk_2d_3 ! 2D workspace 66 USE wrk_nemo, ONLY: wrk_3d_1, wrk_3d_2, wrk_3d_3, wrk_3d_4, wrk_3d_5, wrk_3d_6 ! 3D workspace 67 USE wrk_nemo, ONLY: wrk_4d_1 ! 4D workspace 68 ! 69 INTEGER, INTENT(in) :: kt ! number of iteration 70 ! 71 INTEGER :: ji, jj, jk, jl, layer ! dummy loop indices 72 INTEGER :: initad ! number of sub-timestep for the advection 73 REAL(wp) :: zindb , zindsn , zindic ! local scalar 74 REAL(wp) :: zusvosn, zusvoic, zbigval ! - - 75 REAL(wp) :: zcfl , zusnit , zrtt ! - - 76 REAL(wp) :: ze , zsal , zage ! - - 77 ! 78 REAL(wp), POINTER, DIMENSION(:,:,:) :: zs0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi ! 3D pointer 79 REAL(wp), POINTER, DIMENSION(:,:,:,:) :: zs0e ! 4D pointer 80 !!--------------------------------------------------------------------- 81 82 IF( wrk_in_use(2, 1,2,3,4,5) ) THEN 83 CALL ctl_stop( 'lim_trp : requested workspace arrays unavailable' ) ; RETURN 84 END IF 85 86 zs0ice => wrk_3d_1(:,:,1:jpl) ; zs0a => wrk_3d_3 ; zs0sm => wrk_3d_3 87 zs0sn => wrk_3d_2(:,:,1:jpl) ; zs0c0 => wrk_3d_3 ; zs0oi => wrk_3d_3 88 zs0e => wrk_4d_1(:,:,1:jkmax,1:jpl) 89 90 91 IF( numit == nstart .AND. lwp ) THEN 92 WRITE(numout,*) 93 IF( ln_limdyn ) THEN ; WRITE(numout,*) 'lim_trp : Ice transport ' 94 ELSE ; WRITE(numout,*) 'lim_trp : No ice advection as ln_limdyn = ', ln_limdyn 95 ENDIF 96 WRITE(numout,*) '~~~~~~~~~~~~' 97 ENDIF 98 109 99 zsm(:,:) = area(:,:) 110 100 111 IF( ln_limdyn ) THEN 112 IF( kt == nit000 .AND. lwp ) THEN 113 WRITE(numout,*) ' lim_trp : Ice Advection' 114 WRITE(numout,*) ' ~~~~~~~' 115 ENDIF 116 117 !-----------------------------------------------------------------------------! 118 ! 1) CFL Test 119 !-----------------------------------------------------------------------------! 120 121 !------------------------------------------ 122 ! ice velocities at ocean U- and V-points 123 !------------------------------------------ 124 125 ! zvbord factor between 1 and 2 to take into account slip or no-slip boundary conditions. 126 zvbord = 1.0 + ( 1.0 - bound ) 127 DO jj = 1, jpjm1 128 DO ji = 1, fs_jpim1 129 zui_u(ji,jj) = u_ice(ji,jj) 130 zvi_v(ji,jj) = v_ice(ji,jj) 131 END DO 132 END DO 133 134 ! Lateral boundary conditions 135 CALL lbc_lnk( zui_u, 'U', -1. ) 136 CALL lbc_lnk( zvi_v, 'V', -1. ) 101 ! !-------------------------------------! 102 IF( ln_limdyn ) THEN ! Advection of sea ice properties ! 103 ! !-------------------------------------! 104 ! 137 105 138 106 !------------------------- 107 ! transported fields 108 !------------------------- 109 ! Snow vol, ice vol, salt and age contents, area 110 zs0ow(:,:) = ato_i(:,:) * area(:,:) ! Open water area 111 DO jl = 1, jpl 112 zs0sn (:,:,jl) = v_s (:,:,jl) * area(:,:) ! Snow volume 113 zs0ice(:,:,jl) = v_i (:,:,jl) * area(:,:) ! Ice volume 114 zs0a (:,:,jl) = a_i (:,:,jl) * area(:,:) ! Ice area 115 zs0sm (:,:,jl) = smv_i(:,:,jl) * area(:,:) ! Salt content 116 zs0oi (:,:,jl) = oa_i (:,:,jl) * area(:,:) ! Age content 117 zs0c0 (:,:,jl) = e_s (:,:,1,jl) ! Snow heat content 118 zs0e (:,:,:,jl) = e_i (:,:,:,jl) ! Ice heat content 119 END DO 120 121 !-------------------------- 122 ! Advection of Ice fields (Prather scheme) 123 !-------------------------- 124 ! If ice drift field is too fast, use an appropriate time step for advection. 139 125 ! CFL test for stability 140 !------------------------- 141 142 zcfl = 0.e0 143 zcfl = MAX( zcfl, MAXVAL( ABS( zui_u(1:jpim1, : ) ) * rdt_ice / e1u(1:jpim1, : ) ) ) 144 zcfl = MAX( zcfl, MAXVAL( ABS( zvi_v( : ,1:jpjm1) ) * rdt_ice / e2v( : ,1:jpjm1) ) ) 145 146 zmaxu = 0.0 147 zmaxv = 0.0 148 DO ji = 1, jpim1 149 DO jj = 1, jpjm1 150 IF ( (ABS(zui_u(ji,jj)) .GT. zmaxu) ) THEN 151 zmaxu = MAX(zui_u(ji,jj), zmaxu ) 152 ji_maxu = ji 153 jj_maxu = jj 154 ENDIF 155 IF ( (ABS(zvi_v(ji,jj)) .GT. zmaxv) ) THEN 156 zmaxv = MAX(zvi_v(ji,jj), zmaxv ) 157 ji_maxv = ji 158 jj_maxv = jj 159 ENDIF 160 END DO 161 END DO 162 163 IF (lk_mpp ) CALL mpp_max(zcfl) 164 165 IF ( zcfl > 0.5 .AND. lwp ) & 166 WRITE(numout,*) 'lim_trp : violation of cfl criterion the ',nday,'th day, cfl = ',zcfl 167 168 !-----------------------------------------------------------------------------! 169 ! 2) Computation of transported fields 170 !-----------------------------------------------------------------------------! 171 172 !------------------------------------------------------ 173 ! 1.1) Snow vol, ice vol, salt and age contents, area 174 !------------------------------------------------------ 175 176 zs0ow (:,:) = ato_i(:,:) * area(:,:) ! Open water area 177 DO jl = 1, jpl !sum over thickness categories 178 ! area -> is the unmasked and masked area of T-grid cell 179 zs0sn (:,:,jl) = v_s(:,:,jl) * area(:,:) ! Snow volume. 180 zs0ice(:,:,jl) = v_i(:,:,jl) * area(:,:) ! Ice volume. 181 zs0a (:,:,jl) = a_i(:,:,jl) * area(:,:) ! Ice area 182 zs0sm (:,:,jl) = smv_i(:,:,jl) * area(:,:) ! Salt content 183 zs0oi (:,:,jl) = oa_i (:,:,jl) * area(:,:) ! Age content 184 185 !---------------------------------- 186 ! 1.2) Ice and snow heat contents 187 !---------------------------------- 188 189 zs0c0 (:,:,jl) = e_s(:,:,1,jl) ! Snow heat cont. 190 DO jk = 1, nlay_i 191 zs0e(:,:,jk,jl) = e_i(:,:,jk,jl) ! Ice heat content 192 END DO 193 END DO 194 195 !-----------------------------------------------------------------------------! 196 ! 3) Advection of Ice fields 197 !-----------------------------------------------------------------------------! 198 199 ! If ice drift field is too fast, use an appropriate time step for advection. 126 zcfl = MAXVAL( ABS( u_ice(:,:) ) * rdt_ice / e1u(:,:) ) 127 zcfl = MAX( zcfl, MAXVAL( ABS( v_ice(:,:) ) * rdt_ice / e2v(:,:) ) ) 128 IF(lk_mpp ) CALL mpp_max( zcfl ) 129 !!gm more readability: 130 ! IF( zcfl > 0.5 ) THEN ; initad = 2 ; zusnit = 0.5_wp 131 ! ELSE ; initad = 1 ; zusnit = 1.0_wp 132 ! ENDIF 133 !!gm end 200 134 initad = 1 + INT( MAX( rzero, SIGN( rone, zcfl-0.5 ) ) ) 201 135 zusnit = 1.0 / REAL( initad ) 202 203 IF( MOD( ( kt - 1) / nn_fsbc , 2 ) == 0) THEN !== odd ice time step: adv_x then adv_y ==! 136 IF( zcfl > 0.5 .AND. lwp ) & 137 WRITE(numout,*) 'lim_trp_2 : CFL violation at day ', nday, ', cfl = ', zcfl, & 138 & ': the ice time stepping is split in two' 139 140 IF( MOD( ( kt - 1) / nn_fsbc , 2 ) == 0 ) THEN !== odd ice time step: adv_x then adv_y ==! 204 141 DO jk = 1,initad 205 !--- ice open water area 206 CALL lim_adv_x( zusnit, zui_u, rone , zsm, zs0ow(:,:) , sxopw(:,:) , & 207 sxxopw(:,:), syopw(:,:) , & 208 syyopw(:,:), sxyopw(:,:) ) 209 CALL lim_adv_y( zusnit, zvi_v, rzero, zsm, zs0ow(:,:) , sxopw (:,:) , & 210 sxxopw(:,:), syopw (:,:) , & 211 syyopw(:,:), sxyopw(:,:) ) 142 CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0ow (:,:), sxopw(:,:), & !--- ice open water area 143 & sxxopw(:,:), syopw(:,:), syyopw(:,:), sxyopw(:,:) ) 144 CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0ow (:,:), sxopw(:,:), & 145 & sxxopw(:,:), syopw(:,:), syyopw(:,:), sxyopw(:,:) ) 212 146 DO jl = 1, jpl 213 !--- ice volume --- 214 CALL lim_adv_x( zusnit, zui_u, rone , zsm, zs0ice(:,:,jl) , sxice (:,:,jl) , & 215 sxxice(:,:,jl) , syice (:,:,jl) , & 216 syyice(:,:,jl) , sxyice(:,:,jl) ) 217 CALL lim_adv_y( zusnit, zvi_v, rzero, zsm, zs0ice(:,:,jl) , sxice (:,:,jl) , & 218 sxxice(:,:,jl) , syice (:,:,jl) , & 219 syyice(:,:,jl) , sxyice(:,:,jl) ) 220 !--- snow volume --- 221 CALL lim_adv_x( zusnit, zui_u, rone , zsm, zs0sn (:,:,jl) , sxsn (:,:,jl) , & 222 sxxsn (:,:,jl) , sysn (:,:,jl) , & 223 syysn (:,:,jl) , sxysn (:,:,jl) ) 224 CALL lim_adv_y( zusnit, zvi_v, rzero, zsm, zs0sn (:,:,jl) , sxsn (:,:,jl) , & 225 sxxsn (:,:,jl) , sysn (:,:,jl) , & 226 syysn (:,:,jl) , sxysn (:,:,jl) ) 227 !--- ice salinity --- 228 CALL lim_adv_x( zusnit, zui_u, rone , zsm, zs0sm (:,:,jl) , sxsal (:,:,jl) , & 229 sxxsal(:,:,jl) , sysal (:,:,jl) , & 230 syysal(:,:,jl) , sxysal(:,:,jl) ) 231 CALL lim_adv_y( zusnit, zvi_v, rzero, zsm, zs0sm (:,:,jl) , sxsal (:,:,jl) , & 232 sxxsal(:,:,jl) , sysal (:,:,jl) , & 233 syysal(:,:,jl) , sxysal(:,:,jl) ) 234 !--- ice age --- 235 CALL lim_adv_x( zusnit, zui_u, rone , zsm, zs0oi (:,:,jl) , sxage (:,:,jl) , & 236 sxxage(:,:,jl) , syage (:,:,jl) , & 237 syyage(:,:,jl) , sxyage(:,:,jl) ) 238 CALL lim_adv_y( zusnit, zvi_v, rzero, zsm, zs0oi (:,:,jl) , sxage (:,:,jl) , & 239 sxxage(:,:,jl) , syage (:,:,jl) , & 240 syyage(:,:,jl) , sxyage(:,:,jl) ) 241 !--- ice concentrations --- 242 CALL lim_adv_x( zusnit, zui_u, rone , zsm, zs0a (:,:,jl) , sxa (:,:,jl) , & 243 sxxa (:,:,jl) , sya (:,:,jl) , & 244 syya (:,:,jl) , sxya (:,:,jl) ) 245 CALL lim_adv_y( zusnit, zvi_v, rzero, zsm, zs0a (:,:,jl) , sxa (:,:,jl) , & 246 sxxa (:,:,jl) , sya (:,:,jl) , & 247 syya (:,:,jl) , sxya (:,:,jl) ) 248 !--- ice / snow thermal energetic contents --- 249 CALL lim_adv_x( zusnit, zui_u, rone , zsm, zs0c0 (:,:,jl) , sxc0 (:,:,jl) , & 250 sxxc0 (:,:,jl) , syc0 (:,:,jl) , & 251 syyc0 (:,:,jl) , sxyc0 (:,:,jl) ) 252 CALL lim_adv_y( zusnit, zvi_v, rzero, zsm, zs0c0 (:,:,jl) , sxc0 (:,:,jl) , & 253 sxxc0 (:,:,jl) , syc0 (:,:,jl) , & 254 syyc0 (:,:,jl) , sxyc0 (:,:,jl) ) 255 DO layer = 1, nlay_i 256 CALL lim_adv_x( zusnit, zui_u, rone , zsm, & 257 zs0e(:,:,layer,jl) , sxe (:,:,layer,jl) , & 258 sxxe(:,:,layer,jl) , sye (:,:,layer,jl) , & 259 syye(:,:,layer,jl) , sxye(:,:,layer,jl) ) 260 CALL lim_adv_y( zusnit, zvi_v, rzero, zsm, & 261 zs0e(:,:,layer,jl) , sxe (:,:,layer,jl) , & 262 sxxe(:,:,layer,jl) , sye (:,:,layer,jl) , & 263 syye(:,:,layer,jl) , sxye(:,:,layer,jl) ) 147 CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0ice(:,:,jl), sxice(:,:,jl), & !--- ice volume --- 148 & sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl) ) 149 CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0ice(:,:,jl), sxice(:,:,jl), & 150 & sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl) ) 151 CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0sn (:,:,jl), sxsn (:,:,jl), & !--- snow volume --- 152 & sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl) ) 153 CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0sn (:,:,jl), sxsn (:,:,jl), & 154 & sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl) ) 155 CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0sm (:,:,jl), sxsal(:,:,jl), & !--- ice salinity --- 156 & sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl) ) 157 CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0sm (:,:,jl), sxsal(:,:,jl), & 158 & sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl) ) 159 CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0oi (:,:,jl), sxage(:,:,jl), & !--- ice age --- 160 & sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl) ) 161 CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0oi (:,:,jl), sxage(:,:,jl), & 162 & sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl) ) 163 CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0a (:,:,jl), sxa (:,:,jl), & !--- ice concentrations --- 164 & sxxa (:,:,jl), sya (:,:,jl), syya (:,:,jl), sxya (:,:,jl) ) 165 CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0a (:,:,jl), sxa (:,:,jl), & 166 & sxxa (:,:,jl), sya (:,:,jl), syya (:,:,jl), sxya (:,:,jl) ) 167 CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl), & !--- snow heat contents --- 168 & sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl) ) 169 CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl), & 170 & sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl) ) 171 DO layer = 1, nlay_i !--- ice heat contents --- 172 CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0e(:,:,layer,jl), sxe (:,:,layer,jl), & 173 & sxxe(:,:,layer,jl), sye (:,:,layer,jl), & 174 & syye(:,:,layer,jl), sxye(:,:,layer,jl) ) 175 CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0e(:,:,layer,jl), sxe (:,:,layer,jl), & 176 & sxxe(:,:,layer,jl), sye (:,:,layer,jl), & 177 & syye(:,:,layer,jl), sxye(:,:,layer,jl) ) 264 178 END DO 265 179 END DO … … 267 181 ELSE 268 182 DO jk = 1, initad 269 !--- ice volume --- 270 CALL lim_adv_y( zusnit, zvi_v, rone , zsm, zs0ow (:,:) , sxopw (:,:) , & 271 sxxopw(:,:) , syopw (:,:) , & 272 syyopw(:,:) , sxyopw(:,:) ) 273 CALL lim_adv_x( zusnit, zui_u, rzero, zsm, zs0ow (:,:) , sxopw (:,:) , & 274 sxxopw(:,:) , syopw (:,:) , & 275 syyopw(:,:) , sxyopw(:,:) ) 183 CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0ow (:,:), sxopw(:,:), & !--- ice open water area 184 & sxxopw(:,:), syopw(:,:), syyopw(:,:), sxyopw(:,:) ) 185 CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0ow (:,:), sxopw(:,:), & 186 & sxxopw(:,:), syopw(:,:), syyopw(:,:), sxyopw(:,:) ) 276 187 DO jl = 1, jpl 277 !--- ice volume --- 278 CALL lim_adv_y( zusnit, zvi_v, rone , zsm, zs0ice(:,:,jl) , sxice (:,:,jl) , & 279 sxxice(:,:,jl) , syice (:,:,jl) , & 280 syyice(:,:,jl) , sxyice(:,:,jl) ) 281 CALL lim_adv_x( zusnit, zui_u, rzero, zsm, zs0ice(:,:,jl) , sxice (:,:,jl) , & 282 sxxice(:,:,jl) , syice (:,:,jl) , & 283 syyice(:,:,jl) , sxyice(:,:,jl) ) 284 !--- snow volume --- 285 CALL lim_adv_y( zusnit, zvi_v, rone , zsm, zs0sn (:,:,jl) , sxsn (:,:,jl) , & 286 sxxsn (:,:,jl) , sysn (:,:,jl) , & 287 syysn (:,:,jl) , sxysn (:,:,jl) ) 288 CALL lim_adv_x( zusnit, zui_u, rzero, zsm, zs0sn (:,:,jl) , sxsn (:,:,jl) , & 289 sxxsn (:,:,jl) , sysn (:,:,jl) , & 290 syysn (:,:,jl) , sxysn (:,:,jl) ) 291 !--- ice salinity --- 292 CALL lim_adv_y( zusnit, zvi_v, rone , zsm, zs0sm (:,:,jl) , sxsal (:,:,jl) , & 293 sxxsal(:,:,jl) , sysal (:,:,jl) , & 294 syysal(:,:,jl) , sxysal(:,:,jl) ) 295 CALL lim_adv_x( zusnit, zui_u, rzero, zsm, zs0sm (:,:,jl) , sxsal (:,:,jl) , & 296 sxxsal(:,:,jl) , sysal (:,:,jl) , & 297 syysal(:,:,jl) , sxysal(:,:,jl) ) 298 !--- ice age --- 299 CALL lim_adv_y( zusnit, zvi_v, rone , zsm, zs0oi (:,:,jl) , sxage (:,:,jl) , & 300 sxxage(:,:,jl) , syage (:,:,jl) , & 301 syyage(:,:,jl) , sxyage(:,:,jl) ) 302 CALL lim_adv_x( zusnit, zui_u, rzero, zsm, zs0oi (:,:,jl) , sxage (:,:,jl) , & 303 sxxage(:,:,jl) , syage (:,:,jl) , & 304 syyage(:,:,jl) , sxyage(:,:,jl) ) 305 !--- ice concentration --- 306 CALL lim_adv_y( zusnit, zvi_v, rone , zsm, zs0a (:,:,jl) , sxa (:,:,jl) , & 307 sxxa (:,:,jl) , sya (:,:,jl) , & 308 syya (:,:,jl) , sxya (:,:,jl) ) 309 CALL lim_adv_x( zusnit, zui_u, rzero, zsm, zs0a (:,:,jl) , sxa (:,:,jl) , & 310 sxxa (:,:,jl) , sya (:,:,jl) , & 311 syya (:,:,jl) , sxya (:,:,jl) ) 312 !--- ice / snow thermal energetic contents --- 313 CALL lim_adv_y( zusnit, zvi_v, rone , zsm, zs0c0 (:,:,jl) , sxc0 (:,:,jl) , & 314 sxxc0 (:,:,jl) , syc0 (:,:,jl) , & 315 syyc0 (:,:,jl) , sxyc0 (:,:,jl) ) 316 CALL lim_adv_x( zusnit, zui_u, rzero, zsm, zs0c0 (:,:,jl) , sxc0 (:,:,jl) , & 317 sxxc0 (:,:,jl) , syc0 (:,:,jl) , & 318 syyc0 (:,:,jl) , sxyc0 (:,:,jl) ) 319 DO layer = 1, nlay_i 320 CALL lim_adv_y( zusnit, zvi_v, rone , zsm, zs0e(:,:,layer,jl) , & 321 sxe (:,:,layer,jl) , sxxe (:,:,layer,jl) , sye (:,:,layer,jl) , & 322 syye (:,:,layer,jl), sxye (:,:,layer,jl) ) 323 CALL lim_adv_x( zusnit, zui_u, rzero, zsm, zs0e(:,:,layer,jl) , & 324 sxe (:,:,layer,jl) , sxxe (:,:,layer,jl) , sye (:,:,layer,jl) , & 325 syye (:,:,layer,jl), sxye (:,:,layer,jl) ) 188 CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0ice(:,:,jl), sxice(:,:,jl), & !--- ice volume --- 189 & sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl) ) 190 CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0ice(:,:,jl), sxice(:,:,jl), & 191 & sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl) ) 192 CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0sn (:,:,jl), sxsn (:,:,jl), & !--- snow volume --- 193 & sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl) ) 194 CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0sn (:,:,jl), sxsn (:,:,jl), & 195 & sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl) ) 196 CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0sm (:,:,jl), sxsal(:,:,jl), & !--- ice salinity --- 197 & sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl) ) 198 CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0sm (:,:,jl), sxsal(:,:,jl), & 199 & sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl) ) 200 201 CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0oi (:,:,jl), sxage(:,:,jl), & !--- ice age --- 202 & sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl) ) 203 CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0oi (:,:,jl), sxage(:,:,jl), & 204 & sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl) ) 205 CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0a (:,:,jl), sxa (:,:,jl), & !--- ice concentrations --- 206 & sxxa (:,:,jl), sya (:,:,jl), syya (:,:,jl), sxya (:,:,jl) ) 207 CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0a (:,:,jl), sxa (:,:,jl), & 208 & sxxa (:,:,jl), sya (:,:,jl), syya (:,:,jl), sxya (:,:,jl) ) 209 CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl), & !--- snow heat contents --- 210 & sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl) ) 211 CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl), & 212 & sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl) ) 213 DO layer = 1, nlay_i !--- ice heat contents --- 214 CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0e(:,:,layer,jl), sxe (:,:,layer,jl), & 215 & sxxe(:,:,layer,jl), sye (:,:,layer,jl), & 216 & syye(:,:,layer,jl), sxye(:,:,layer,jl) ) 217 CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0e(:,:,layer,jl), sxe (:,:,layer,jl), & 218 & sxxe(:,:,layer,jl), sye (:,:,layer,jl), & 219 & syye(:,:,layer,jl), sxye(:,:,layer,jl) ) 326 220 END DO 327 328 221 END DO 329 222 END DO … … 333 226 ! Recover the properties from their contents 334 227 !------------------------------------------- 335 336 zs0ow (:,:) = zs0ow(:,:) / area(:,:) 228 zs0ow(:,:) = zs0ow(:,:) / area(:,:) 337 229 DO jl = 1, jpl 338 230 zs0ice(:,:,jl) = zs0ice(:,:,jl) / area(:,:) … … 351 243 !------------------------------------------------------------------------------! 352 244 245 !-------------------------------- 246 ! diffusion of open water area 247 !-------------------------------- 248 zs0at(:,:) = zs0a(:,:,1) ! total ice fraction 249 DO jl = 2, jpl 250 zs0at(:,:) = zs0at(:,:) + zs0a(:,:,jl) 251 END DO 252 ! 253 ! ! Masked eddy diffusivity coefficient at ocean U- and V-points 254 DO jj = 1, jpjm1 ! NB: has not to be defined on jpj line and jpi row 255 DO ji = 1 , fs_jpim1 ! vector opt. 256 pahu(ji,jj) = ( 1._wp - MAX( rzero, SIGN( rone, -zs0at(ji ,jj) ) ) ) & 257 & * ( 1._wp - MAX( rzero, SIGN( rone, -zs0at(ji+1,jj) ) ) ) * ahiu(ji,jj) 258 pahv(ji,jj) = ( 1._wp - MAX( rzero, SIGN( rone, -zs0at(ji,jj ) ) ) ) & 259 & * ( 1._wp - MAX( rzero, SIGN( rone,- zs0at(ji,jj+1) ) ) ) * ahiv(ji,jj) 260 END DO 261 END DO 262 ! 263 CALL lim_hdf( zs0ow (:,:) ) ! Diffusion 264 353 265 !------------------------------------ 354 ! 4.1) diffusion of open water area266 ! Diffusion of other ice variables 355 267 !------------------------------------ 356 357 ! Compute total ice fraction 358 zs0at(:,:) = 0.0 359 DO jl = 1, jpl 360 DO jj = 1, jpj 361 DO ji = 1, jpi 362 zs0at (ji,jj) = zs0at(ji,jj) + zs0a(ji,jj,jl) ! 363 END DO 364 END DO 365 END DO 366 367 ! Masked eddy diffusivity coefficient at ocean U- and V-points 368 DO jj = 1, jpjm1 ! NB: has not to be defined on jpj line and jpi row 369 DO ji = 1 , fs_jpim1 ! vector opt. 370 pahu(ji,jj) = ( 1.0 - MAX( rzero, SIGN( rone, -zs0at(ji ,jj) ) ) ) & 371 & * ( 1.0 - MAX( rzero, SIGN( rone, -zs0at(ji+1,jj) ) ) ) * ahiu(ji,jj) 372 pahv(ji,jj) = ( 1.0 - MAX( rzero, SIGN( rone, -zs0at(ji,jj ) ) ) ) & 373 & * ( 1.0 - MAX( rzero, SIGN( rone,- zs0at(ji,jj+1) ) ) ) * ahiv(ji,jj) 374 END DO !jj 375 END DO ! ji 376 377 ! Diffusion 378 CALL lim_hdf( zs0ow (:,:) ) 379 380 !---------------------------------------- 381 ! 4.2) Diffusion of other ice variables 382 !---------------------------------------- 383 DO jl = 1, jpl 384 385 ! Masked eddy diffusivity coefficient at ocean U- and V-points 386 DO jj = 1, jpjm1 ! NB: has not to be defined on jpj line and jpi row 268 DO jl = 1, jpl 269 ! ! Masked eddy diffusivity coefficient at ocean U- and V-points 270 DO jj = 1, jpjm1 ! NB: has not to be defined on jpj line and jpi row 387 271 DO ji = 1 , fs_jpim1 ! vector opt. 388 pahu(ji,jj) = ( 1. 0- MAX( rzero, SIGN( rone, -zs0a(ji ,jj,jl) ) ) ) &389 & * ( 1. 0- MAX( rzero, SIGN( rone, -zs0a(ji+1,jj,jl) ) ) ) * ahiu(ji,jj)390 pahv(ji,jj) = ( 1. 0 - MAX( rzero, SIGN( rone, -zs0a(ji,jj,jl) ) ) ) &391 & * ( 1. 0- MAX( rzero, SIGN( rone,- zs0a(ji,jj+1,jl) ) ) ) * ahiv(ji,jj)392 END DO !jj393 END DO ! ji272 pahu(ji,jj) = ( 1._wp - MAX( rzero, SIGN( rone, -zs0a(ji ,jj,jl) ) ) ) & 273 & * ( 1._wp - MAX( rzero, SIGN( rone, -zs0a(ji+1,jj,jl) ) ) ) * ahiu(ji,jj) 274 pahv(ji,jj) = ( 1._wp - MAX( rzero, SIGN( rone, -zs0a(ji,jj ,jl) ) ) ) & 275 & * ( 1._wp - MAX( rzero, SIGN( rone,- zs0a(ji,jj+1,jl) ) ) ) * ahiv(ji,jj) 276 END DO 277 END DO 394 278 395 279 CALL lim_hdf( zs0ice (:,:,jl) ) … … 401 285 DO jk = 1, nlay_i 402 286 CALL lim_hdf( zs0e (:,:,jk,jl) ) 403 END DO ! jk404 END DO !jl287 END DO 288 END DO 405 289 406 290 !----------------------------------------- 407 ! 4.3)Remultiply everything by ice area291 ! Remultiply everything by ice area 408 292 !----------------------------------------- 409 zs0ow(:,:) = MAX( rzero, zs0ow(:,:) * area(:,:) )293 zs0ow(:,:) = MAX( rzero, zs0ow(:,:) * area(:,:) ) 410 294 DO jl = 1, jpl 411 295 zs0ice(:,:,jl) = MAX( rzero, zs0ice(:,:,jl) * area(:,:) ) !!bug: est-ce utile … … 432 316 DO jj = 1, jpj 433 317 DO ji = 1, jpi 434 zs0e (ji,jj,jk,jl) = & 435 MAX( rzero, zs0e (ji,jj,jk,jl) / area(ji,jj) ) 318 zs0e(ji,jj,jk,jl) = MAX( rzero, zs0e(ji,jj,jk,jl) / area(ji,jj) ) 436 319 END DO 437 320 END DO … … 441 324 DO jj = 1, jpj 442 325 DO ji = 1, jpi 443 zs0ow 444 END DO 445 END DO 446 447 zs0at(:,:) = 0. 0326 zs0ow(ji,jj) = MAX( rzero, zs0ow (ji,jj) / area(ji,jj) ) 327 END DO 328 END DO 329 330 zs0at(:,:) = 0._wp 448 331 DO jl = 1, jpl 449 332 DO jj = 1, jpj … … 463 346 ! 5.2) Snow thickness, Ice thickness, Ice concentrations 464 347 !--------------------------------------------------------- 465 466 348 DO jj = 1, jpj 467 349 DO ji = 1, jpi 468 zindb = MAX( 0.0 , SIGN( 1.0, zs0at(ji,jj) - zeps10) ) 469 zs0ow(ji,jj) = (1.0 - zindb) + zindb*MAX( zs0ow(ji,jj), 0.00 ) 470 ato_i(ji,jj) = zs0ow(ji,jj) 471 END DO 472 END DO 473 474 ! Remove very small areas 475 DO jl = 1, jpl 350 zindb = MAX( 0._wp , SIGN( 1.0, zs0at(ji,jj) - zeps10) ) 351 zs0ow(ji,jj) = ( 1._wp - zindb ) + zindb * MAX( zs0ow(ji,jj), 0._wp ) 352 ato_i(ji,jj) = zs0ow(ji,jj) 353 END DO 354 END DO 355 356 DO jl = 1, jpl ! Remove very small areas 476 357 DO jj = 1, jpj 477 358 DO ji = 1, jpi 478 359 zindb = MAX( 0.0 , SIGN( 1.0, zs0a(ji,jj,jl) - zeps10) ) 479 480 zs0a(ji,jj,jl) 481 v_s(ji,jj,jl) 482 v_i(ji,jj,jl) 483 484 zindsn 485 zindic 486 zindb 487 zs0a 488 a_i 489 v_s 490 v_i 360 ! 361 zs0a(ji,jj,jl) = zindb * MIN( zs0a(ji,jj,jl), 0.99 ) 362 v_s(ji,jj,jl) = zindb * zs0sn (ji,jj,jl) 363 v_i(ji,jj,jl) = zindb * zs0ice(ji,jj,jl) 364 ! 365 zindsn = MAX( rzero, SIGN( rone, v_s(ji,jj,jl) - zeps10 ) ) 366 zindic = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - zeps10 ) ) 367 zindb = MAX( zindsn, zindic ) 368 zs0a(ji,jj,jl) = zindb * zs0a(ji,jj,jl) !ice concentration 369 a_i (ji,jj,jl) = zs0a(ji,jj,jl) 370 v_s (ji,jj,jl) = zindsn * v_s(ji,jj,jl) 371 v_i (ji,jj,jl) = zindic * v_i(ji,jj,jl) 491 372 END DO 492 373 END DO … … 495 376 DO jj = 1, jpj 496 377 DO ji = 1, jpi 497 zs0at(ji,jj) 378 zs0at(ji,jj) = SUM( zs0a(ji,jj,1:jpl) ) 498 379 END DO 499 380 END DO … … 503 384 !---------------------- 504 385 505 zbigval = 1.0d+13386 zbigval = 1.d+13 506 387 507 388 DO jl = 1, jpl … … 518 399 519 400 ! Ice salinity and age 520 zsal = MAX( MIN( (rhoic-rhosn)/rhoic*sss_m(ji,jj) , & 521 zusvoic * zs0sm(ji,jj,jl) ), s_i_min ) * & 522 v_i(ji,jj,jl) 401 zsal = MAX( MIN( (rhoic-rhosn)/rhoic*sss_m(ji,jj) , & 402 zusvoic * zs0sm(ji,jj,jl) ), s_i_min ) * v_i(ji,jj,jl) 523 403 IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) & 524 404 smv_i(ji,jj,jl) = zindic*zsal + (1.0-zindic)*0.0 525 405 526 zage = MAX( MIN( zbigval, zs0oi(ji,jj,jl) / & 527 MAX( a_i(ji,jj,jl), epsi16 ) ), 0.0 ) * & 528 a_i(ji,jj,jl) 406 zage = MAX( MIN( zbigval, zs0oi(ji,jj,jl) / & 407 MAX( a_i(ji,jj,jl), epsi16 ) ), 0.0 ) * a_i(ji,jj,jl) 529 408 oa_i (ji,jj,jl) = zindic*zage 530 409 … … 583 462 END DO 584 463 ENDIF 585 464 ! 465 IF( wrk_not_released(2, 1,2,3,4,5) ) CALL ctl_stop('lim_trp_2 : failed to release workspace arrays') 466 ! 586 467 END SUBROUTINE lim_trp 587 588 589 SUBROUTINE lim_trp_init590 !!-------------------------------------------------------------------591 !! *** ROUTINE lim_trp_init ***592 !!593 !! ** Purpose : initialization of ice advection parameters594 !!595 !! ** Method : Read the namicetrp namelist and check the parameter596 !! values called at the first timestep (nit000)597 !!598 !! ** input : Namelist namicetrp599 !!600 !! history :601 !! 2.0 ! 03-08 (C. Ethe) Original code602 !!-------------------------------------------------------------------603 NAMELIST/namicetrp/ bound604 !!-------------------------------------------------------------------605 606 ! Read Namelist namicetrp607 REWIND ( numnam_ice )608 READ ( numnam_ice , namicetrp )609 IF(lwp) THEN610 WRITE(numout,*)611 WRITE(numout,*) 'lim_trp_init : Ice parameters for advection '612 WRITE(numout,*) '~~~~~~~~~~~~'613 WRITE(numout,*) ' boundary conditions (0.0 no-slip, 1.0 free-slip) bound = ', bound614 WRITE(numout,*)615 ENDIF616 617 END SUBROUTINE lim_trp_init618 468 619 469 #else
Note: See TracChangeset
for help on using the changeset viewer.
