- Timestamp:
- 2017-09-01T15:49:35+02:00 (7 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icestp.F90
r8426 r8486 17 17 #if defined key_lim3 18 18 !!---------------------------------------------------------------------- 19 !! 'key_lim3' : LIM 3.0 sea-ice model 20 !!---------------------------------------------------------------------- 21 !! ice_stp : sea-ice model time-stepping and update ocean sbc over ice-covered area 22 !!---------------------------------------------------------------------- 23 USE oce ! ocean dynamics and tracers 24 USE dom_oce ! ocean space and time domain 25 USE ice ! LIM-3: ice variables 26 USE ice1D ! LIM-3: thermodynamical variables 19 !! 'key_lim3' LIM 3.0 sea-ice model 20 !!---------------------------------------------------------------------- 21 !! ice_stp : sea-ice model time-stepping and update ocean surf. boundary cond. over ice-covered area 22 !! ice_init : 23 !! ice_run_init : 24 !!---------------------------------------------------------------------- 25 USE oce ! ocean dynamics and tracers 26 USE dom_oce ! ocean space and time domain 27 USE c1d ! 1D vertical configuration 28 USE ice ! sea-ice: variables 29 USE ice1D ! sea-ice: thermodynamical 1D variables 27 30 ! 28 USE sbc_oce 29 USE sbc_ice 30 USE iceforcing ! Surface boundary condition for sea ice31 USE sbc_oce ! Surface boundary condition: ocean fields 32 USE sbc_ice ! Surface boundary condition: ice fields 33 USE iceforcing ! sea-ice: Surface boundary condition !!gm why not icesbc module name 31 34 ! 32 USE phycst 33 USE eosbn2 34 USE icerhg ! Icerheology35 USE iceadv ! Iceadvection36 USE icethd ! Icethermodynamics37 USE icerdgrft ! Iceridging/rafting38 USE iceupdate ! sea surface boundary condition39 USE icedia ! Icebudget diagnostics40 USE icewri ! Iceoutputs41 USE icerst ! Icerestarts42 USE icecor ! Icecorrections43 USE icevar ! Ice variables switch44 USE icectl !35 USE phycst ! Define parameters for the routines 36 USE eosbn2 ! equation of state 37 USE icerhg ! sea-ice: rheology 38 USE iceadv ! sea-ice: advection 39 USE icethd ! sea-ice: thermodynamics 40 USE icerdgrft ! sea-ice: ridging/rafting 41 USE iceupdate ! sea-ice: sea surface boundary condition update 42 USE icedia ! sea-ice: budget diagnostics 43 USE icewri ! sea-ice: outputs 44 USE icerst ! sea-ice: restarts 45 USE icecor ! sea-ice: corrections 46 USE icevar ! sea-ice: operations 47 USE icectl ! sea-ice: control 45 48 ! MV MP 2016 46 USE limmp 49 USE limmp ! sea-ice: melt ponds 47 50 ! END MV MP 2016 48 USE iceist ! LIMinitial state49 USE icethd_sal ! LIM ice thermodynamics:salinity51 USE iceist ! sea-ice: initial state 52 USE icethd_sal ! sea-ice: thermodynamics and salinity 50 53 ! 51 USE c1d ! 1D vertical configuration 52 USE in_out_manager ! I/O manager 53 USE iom ! I/O manager library 54 USE prtctl ! Print control 55 USE lib_fortran ! 56 USE lbclnk ! lateral boundary condition - MPP link 57 USE lib_mpp ! MPP library 58 USE timing ! Timing 59 60 USE bdy_oce , ONLY: ln_bdy 61 USE bdyice ! unstructured open boundary data 54 USE bdy_oce , ONLY : ln_bdy ! flag for bdy 55 USE bdyice ! unstructured open boundary data for sea-ice 62 56 # if defined key_agrif 63 57 USE agrif_ice … … 65 59 USE agrif_lim3_interp 66 60 # endif 61 ! 62 USE in_out_manager ! I/O manager 63 USE iom ! I/O manager library 64 USE prtctl ! Print control 65 USE lib_fortran ! 66 USE lbclnk ! lateral boundary condition - MPP link 67 USE lib_mpp ! MPP library 68 USE timing ! Timing 67 69 68 70 IMPLICIT NONE 69 71 PRIVATE 70 72 71 PUBLIC ice_stp ! routine called by sbcmod.F90 72 PUBLIC ice_init ! routine called by sbcmod.F90 73 PUBLIC ice_stp ! called by sbcmod.F90 74 PUBLIC ice_init ! called by sbcmod.F90 75 76 INTEGER :: nice_dyn ! choice of the type of advection scheme 77 ! ! associated indices: 78 INTEGER, PARAMETER :: np_dynNO = 0 ! no ice dynamics and ice advection 79 INTEGER, PARAMETER :: np_dynFULL = 1 ! full ice dynamics (rheology + advection + ridging/rafting + correction) 80 INTEGER, PARAMETER :: np_dyn = 2 ! no ridging/rafting (rheology + advection + correction) 81 INTEGER, PARAMETER :: np_dynPURE = 3 ! pure dynamics (rheology + advection) 73 82 74 83 !! * Substitutions 75 84 # include "vectopt_loop_substitute.h90" 76 85 !!---------------------------------------------------------------------- 77 !! NEMO/ OPA 4.0 , UCL NEMO Consortium (2011)86 !! NEMO/ICE 4.0 , NEMO Consortium (2017) 78 87 !! $Id: icestp.F90 8319 2017-07-11 15:00:44Z clem $ 79 88 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 85 94 !! *** ROUTINE ice_stp *** 86 95 !! 87 !! ** Purpose : update the ocean surface boundary condition via the88 !! Louvain la Neuve Sea Ice Model time stepping96 !! ** Purpose : sea-ice model time-stepping and update ocean surface 97 !! boundary condition over ice-covered area 89 98 !! 90 99 !! ** Method : ice model time stepping … … 102 111 !!--------------------------------------------------------------------- 103 112 INTEGER, INTENT(in) :: kt ! ocean time step 104 INTEGER, INTENT(in) :: ksbc ! type of sbc flux ( 1 = user defined formulation, 105 ! 3 = bulk formulation, 106 ! 4 = Pure Coupled formulation) 107 INTEGER :: jl ! dummy loop index 108 !!---------------------------------------------------------------------- 109 113 INTEGER, INTENT(in) :: ksbc ! flux formulation (user defined, bulk, or Pure Coupled) 114 ! 115 INTEGER :: jl ! dummy loop index 116 !!---------------------------------------------------------------------- 117 ! 110 118 IF( nn_timing == 1 ) CALL timing_start('ice_stp') 111 112 !-----------------------! 113 ! --- Ice time step --- ! 114 !-----------------------! 115 IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN 116 117 ! Ice model time step 118 kt_ice = kt 119 119 ! 120 ! !-----------------------! 121 IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN ! --- Ice time step --- ! 122 ! !-----------------------! 123 ! 124 kt_ice = kt ! Ice model time step 125 ! 120 126 # if defined key_agrif 121 127 IF( .NOT. Agrif_Root() ) lim_nbstep = MOD( lim_nbstep, Agrif_irhot() * Agrif_Parent(nn_fsbc) / nn_fsbc ) + 1 122 128 # endif 123 129 124 ! mean surface ocean current at ice velocity point (C-grid dynamics : U- & V-points as the ocean)130 ! ! mean surface ocean current at ice velocity point 125 131 u_oce(:,:) = ssu_m(:,:) * umask(:,:,1) 126 132 v_oce(:,:) = ssv_m(:,:) * vmask(:,:,1) 127 128 ! masked sea surface freezing temperature [Kelvin] (set to rt0 over land) 133 !!gm the umask, vmask above is useless as ssu_m, ssv_m are build from masked un,vn (used t be here for B-grid sea-ice) 134 135 ! ! masked sea surface freezing temperature [Kelvin] (set to rt0 over land) 129 136 CALL eos_fzp( sss_m(:,:) , t_bo(:,:) ) 130 137 t_bo(:,:) = ( t_bo(:,:) + rt0 ) * tmask(:,:,1) + rt0 * ( 1._wp - tmask(:,:,1) ) 131 138 ! 132 CALL ice_bef ! Store previous ice values 139 CALL ice_bef ! Store previous ice values 140 ! 133 141 !------------------------------------------------! 134 142 ! --- Dynamical coupling with the atmosphere --- ! 135 143 !------------------------------------------------! 136 ! it provides:144 ! It provides the following fields used in sea ice model: 137 145 ! utau_ice, vtau_ice = surface ice stress [N/m2] 138 !------------------------------------------------ --139 146 !------------------------------------------------! 147 CALL ice_forcing_tau( kt, ksbc, utau_ice, vtau_ice ) 140 148 141 !-------------------------------------------------------! 142 ! --- ice dynamics and transport (except in 1D case) ---! 143 !-------------------------------------------------------! 144 CALL ice_diag0 ! set diag of mass, heat and salt fluxes to 0 145 CALL ice_rst_opn( kt ) ! Open Ice restart file 146 ! 147 ! --- zap this if no ice dynamics --- ! 148 IF( .NOT. lk_c1d .AND. ln_limdyn ) THEN 149 CALL ice_rhg( kt ) ! -- rheology 150 CALL ice_adv( kt ) ! -- advection 151 IF( nn_limdyn == 2 .AND. nn_monocat /= 2 ) & ! -- ridging/rafting 152 & CALL ice_rdgrft( kt ) 153 IF( nn_limdyn == 2 ) CALL ice_cor( kt , 1 ) ! -- Corrections 154 ! 155 ENDIF 156 ! --- 157 149 !-------------------------------------! 150 ! --- ice dynamics and advection --- ! 151 !-------------------------------------! 152 CALL ice_diag0 ! set diag of mass, heat and salt fluxes to 0 153 CALL ice_rst_opn( kt ) ! Open Ice restart file (if necessary) 154 ! 155 SELECT CASE( nice_dyn ) 156 CASE ( np_dynFULL ) !== all dynamical processes ==! 157 CALL ice_rhg ( kt ) ! -- rheology 158 CALL ice_adv ( kt ) ! -- advection of ice 159 CALL ice_rdgrft( kt ) ! -- ridging/rafting 160 CALL ice_cor ( kt , 1 ) ! -- Corrections 161 CASE ( np_dyn ) !== pure dynamics only ==! (no ridging/rafting) (nono cat. case 2) 162 CALL ice_rhg ( kt ) ! -- rheology 163 CALL ice_adv ( kt ) ! -- advection of ice 164 CALL ice_cor ( kt , 1 ) ! -- Corrections 165 CASE ( np_dynPURE ) !== pure dynamics only ==! (nn_limdyn= 0 or 1 ) 166 CALL ice_rhg ( kt ) ! -- rheology 167 CALL ice_adv ( kt ) ! -- advection of ice 168 END SELECT 169 ! 170 ! !== lateral boundary conditions ==! 158 171 #if defined key_agrif 159 IF( .NOT. Agrif_Root() ) CALL agrif_interp_lim3('T')172 IF( .NOT. Agrif_Root() ) CALL agrif_interp_lim3('T') ! -- AGRIF 160 173 #endif 161 IF( ln_limthd .AND. ln_bdy ) CALL bdy_ice( kt ) ! -- bdy ice thermo 162 ! previous lead fraction and ice volume for flux calculations 163 CALL ice_var_glo2eqv ! ht_i and ht_s for ice albedo calculation 164 CALL ice_var_agg(1) ! at_i for coupling 165 CALL ice_bef 174 IF( ln_limthd .AND. ln_bdy ) CALL bdy_ice( kt ) ! -- bdy ice thermo 175 ! 176 ! 177 ! !== previous lead fraction and ice volume for flux calculations 178 ! 179 CALL ice_var_glo2eqv ! ht_i and ht_s for ice albedo calculation 180 CALL ice_var_agg(1) ! at_i for coupling 181 CALL ice_bef ! Store previous ice values 166 182 167 183 !------------------------------------------------------! … … 169 185 !------------------------------------------------------! 170 186 ! It provides the following fields used in sea ice model: 171 ! fr1_i0 , fr2_i0 = 1sr & 2nd fraction of qsr penetration in ice [%] 172 ! emp_oce , emp_ice = E-P over ocean and sea ice [Kg/m2/s] 173 ! sprecip = solid precipitation [Kg/m2/s] 174 ! evap_ice = sublimation [Kg/m2/s] 175 ! qsr_tot , qns_tot = solar & non solar heat flux (total) [W/m2] 176 ! qsr_ice , qns_ice = solar & non solar heat flux over ice [W/m2] 177 ! dqns_ice = non solar heat sensistivity [W/m2] 178 ! qemp_oce, qemp_ice, qprec_ice, qevap_ice = sensible heat (associated with evap & precip) [W/m2] 179 !------------------------------------------------------------------------------------------------------ 180 CALL ice_forcing_flx( kt, ksbc ) 187 ! fr1_i0 , fr2_i0 = 1sr & 2nd fraction of qsr penetration in ice [%] 188 ! emp_oce , emp_ice = E-P over ocean and sea ice [Kg/m2/s] 189 ! sprecip = solid precipitation [Kg/m2/s] 190 ! evap_ice = sublimation [Kg/m2/s] 191 ! qsr_tot , qns_tot = solar & non solar heat flux (total) [W/m2] 192 ! qsr_ice , qns_ice = solar & non solar heat flux over ice [W/m2] 193 ! dqns_ice = non solar heat sensistivity [W/m2] 194 ! qemp_oce, qemp_ice, = sensible heat (associated with evap & precip) [W/m2] 195 ! qprec_ice, qevap_ice 196 !------------------------------------------------------! 197 CALL ice_forcing_flx( kt, ksbc ) 181 198 182 199 !----------------------------! 183 200 ! --- ice thermodynamics --- ! 184 201 !----------------------------! 185 ! --- zap this if no ice thermo --- ! 186 IF( ln_limthd ) CALL ice_thd( kt ) ! -- Ice thermodynamics 202 IF( ln_limthd ) CALL ice_thd( kt ) ! -- Ice thermodynamics 187 203 188 204 ! MV MP 2016 189 IF ( ln_pnd ) CALL lim_mp( kt )! -- Melt ponds205 IF ( ln_pnd ) CALL lim_mp( kt ) ! -- Melt ponds 190 206 ! END MV MP 2016 191 207 192 IF( ln_limthd ) CALL ice_cor( kt , 2 )! -- Corrections208 IF( ln_limthd ) CALL ice_cor( kt , 2 ) ! -- Corrections 193 209 ! --- 194 210 # if defined key_agrif 195 IF( .NOT. Agrif_Root() ) 211 IF( .NOT. Agrif_Root() ) CALL agrif_update_lim3( kt ) 196 212 # endif 197 CALL ice_var_glo2eqv! necessary calls (at least for coupling)198 CALL ice_var_agg( 2 )! necessary calls (at least for coupling)199 213 CALL ice_var_glo2eqv ! necessary calls (at least for coupling) 214 CALL ice_var_agg( 2 ) ! necessary calls (at least for coupling) 215 ! 200 216 !! clem: one should switch the calculation of the fluxes onto the parent grid but the following calls do not work 201 217 !! moreover it should only be called at the update frequency (as in agrif_lim3_update.F90) … … 203 219 !! IF( .NOT. Agrif_Root() ) CALL Agrif_ChildGrid_To_ParentGrid() 204 220 !!# endif 205 CALL ice_update_flx( kt ) ! -- Update surface oceanmass, heat and salt fluxes221 CALL ice_update_flx( kt ) ! -- Update ocean surface mass, heat and salt fluxes 206 222 !!# if defined key_agrif 207 223 !! IF( .NOT. Agrif_Root() ) CALL Agrif_ParentGrid_To_ChildGrid() 208 224 !!# endif 209 IF( ln_limdiahsb ) CALL ice_dia( kt ) ! -- Diagnostics and outputs225 IF( ln_limdiahsb ) CALL ice_dia( kt ) ! -- Diagnostics and outputs 210 226 ! 211 227 CALL ice_wri( 1 ) ! -- Ice outputs 212 228 ! 213 IF( kt == nit000 .AND. ln_rstart ) & 229 IF( kt == nit000 .AND. ln_rstart ) & !!gm I don't understand the ln_rstart, if needed, add a comment, please 214 230 & CALL iom_close( numrir ) ! close input ice restart file 215 231 ! … … 287 303 IF( ln_limdiahsb) CALL ice_dia_init ! initialization for diags 288 304 ! 289 fr_i (:,:)= at_i(:,:) ! initialisation of sea-ice fraction305 fr_i (:,:) = at_i(:,:) ! initialisation of sea-ice fraction 290 306 tn_ice(:,:,:) = t_su(:,:,:) ! initialisation of surface temp for coupled simu 291 307 ! 292 308 DO jj = 1, jpj 293 309 DO ji = 1, jpi 294 IF( gphit(ji,jj) > 0._wp ) THEN ;rn_amax_2d(ji,jj) = rn_amax_n ! NH295 ELSE ;rn_amax_2d(ji,jj) = rn_amax_s ! SH310 IF( gphit(ji,jj) > 0._wp ) THEN ; rn_amax_2d(ji,jj) = rn_amax_n ! NH 311 ELSE ; rn_amax_2d(ji,jj) = rn_amax_s ! SH 296 312 ENDIF 297 313 END DO … … 318 334 !!------------------------------------------------------------------- 319 335 ! 320 REWIND( numnam_ice_ref ) 336 REWIND( numnam_ice_ref ) ! Namelist namicerun in reference namelist : Parameters for ice 321 337 READ ( numnam_ice_ref, namicerun, IOSTAT = ios, ERR = 901) 322 338 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicerun in reference namelist', lwp ) 323 339 324 REWIND( numnam_ice_cfg ) 340 REWIND( numnam_ice_cfg ) ! Namelist namicerun in configuration namelist : Parameters for ice 325 341 READ ( numnam_ice_cfg, namicerun, IOSTAT = ios, ERR = 902 ) 326 342 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicerun in configuration namelist', lwp ) 327 343 IF(lwm) WRITE ( numoni, namicerun ) 328 344 ! 329 REWIND( numnam_ice_ref ) 345 REWIND( numnam_ice_ref ) ! Namelist namicediag in reference namelist : Parameters for ice 330 346 READ ( numnam_ice_ref, namicediag, IOSTAT = ios, ERR = 903) 331 347 903 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicediag in reference namelist', lwp ) 332 348 333 REWIND( numnam_ice_cfg ) 349 REWIND( numnam_ice_cfg ) ! Namelist namicediag in configuration namelist : Parameters for ice 334 350 READ ( numnam_ice_cfg, namicediag, IOSTAT = ios, ERR = 904 ) 335 351 904 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicediag in configuration namelist', lwp ) 336 352 IF(lwm) WRITE ( numoni, namicediag ) 337 353 ! 338 IF(lwp) THEN 354 IF(lwp) THEN ! control print 339 355 WRITE(numout,*) 340 356 WRITE(numout,*) 'ice_run_init : ice share parameters for dynamics/advection/thermo of sea-ice' 341 357 WRITE(numout,*) ' ~~~~~~' 342 WRITE(numout,*) ' number of ice categories jpl = ', jpl343 WRITE(numout,*) ' number of ice layers nlay_i = ', nlay_i344 WRITE(numout,*) ' number of snow layers nlay_s = ', nlay_s345 WRITE(numout,*) ' virtual ITD mono-category param (1-4) or not (0) nn_monocat = ', nn_monocat346 WRITE(numout,*) ' maximum ice concentration for NH = ', rn_amax_n347 WRITE(numout,*) ' maximum ice concentration for SH = ', rn_amax_s348 WRITE(numout,*) ' Ice thermodynamics (T) or not (F) ln_limthd = ', ln_limthd349 WRITE(numout,*) ' Ice dynamics (T) or not (F) ln_limdyn = ', ln_limdyn350 WRITE(numout,*) ' (ln_limdyn=T) Ice dynamics switch nn_limdyn = ', nn_limdyn351 WRITE(numout,*) ' 2: total'352 WRITE(numout,*) ' 1: advection only (no diffusion, no ridging/rafting)'353 WRITE(numout,*) ' 0: advection only (as 1 + prescribed velocity, bypass rheology)'354 WRITE(numout,*) ' (ln_limdyn=T) prescribed u-vel (case nn_limdyn=0) = ', rn_uice355 WRITE(numout,*) ' (ln_limdyn=T) prescribed v-vel (case nn_limdyn=0) = ', rn_vice358 WRITE(numout,*) ' Namelist namicerun : ' 359 WRITE(numout,*) ' number of ice categories jpl = ', jpl 360 WRITE(numout,*) ' number of ice layers nlay_i = ', nlay_i 361 WRITE(numout,*) ' number of snow layers nlay_s = ', nlay_s 362 WRITE(numout,*) ' virtual ITD mono-category param (1-4) or not (0) nn_monocat = ', nn_monocat 363 WRITE(numout,*) ' maximum ice concentration for NH = ', rn_amax_n 364 WRITE(numout,*) ' maximum ice concentration for SH = ', rn_amax_s 365 WRITE(numout,*) ' Ice thermodynamics (T) or not (F) ln_limthd = ', ln_limthd 366 WRITE(numout,*) ' Ice dynamics (T) or not (F) ln_limdyn = ', ln_limdyn 367 WRITE(numout,*) ' associated switch nn_limdyn = ', nn_limdyn 368 WRITE(numout,*) ' =2 all processes (default option)' 369 WRITE(numout,*) ' =1 advection only (no ridging/rafting)' 370 WRITE(numout,*) ' =0 advection only with prescribed velocity given by ' 371 WRITE(numout,*) ' a uniform field (u,v)_ice = (rn_uice,rn_vice) = (', rn_uice,',', rn_vice,')' 356 372 WRITE(numout,*) 357 WRITE(numout,*) '...and ice diagnostics' 358 WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~' 359 WRITE(numout,*) ' Diagnose online heat/mass/salt budget ln_limdiachk = ', ln_limdiachk 360 WRITE(numout,*) ' Output heat/mass/salt budget ln_limdiahsb = ', ln_limdiahsb 361 WRITE(numout,*) ' control prints in ocean.out for (i,j)=(iiceprt,jiceprt) = ', ln_limctl 362 WRITE(numout,*) ' i-index for control prints (ln_limctl=true) = ', iiceprt 363 WRITE(numout,*) ' j-index for control prints (ln_limctl=true) = ', jiceprt 373 WRITE(numout,*) ' Namelist namicediag : ' 374 WRITE(numout,*) ' Diagnose online heat/mass/salt budget ln_limdiachk = ', ln_limdiachk 375 WRITE(numout,*) ' Output heat/mass/salt budget ln_limdiahsb = ', ln_limdiahsb 376 WRITE(numout,*) ' control prints for a given grid point ln_limctl = ', ln_limctl 377 WRITE(numout,*) ' chosen grid point position (iiceprt,jiceprt) = (', iiceprt,',', jiceprt,')' 364 378 ENDIF 365 379 ! 366 IF ( ( jpl > 1 ) .AND. ( nn_monocat == 1 ) ) THEN 380 ! !--- check consistency 381 IF ( jpl > 1 .AND. nn_monocat == 1 ) THEN 367 382 nn_monocat = 0 368 383 IF(lwp) WRITE(numout,*) 369 384 IF(lwp) WRITE(numout,*) ' nn_monocat forced to 0 as jpl>1, i.e. multi-category case is chosen' 370 385 ENDIF 371 IF ( ( jpl == 1 ) .AND. ( nn_monocat == 0 )) THEN386 IF ( jpl == 1 .AND. nn_monocat == 0 ) THEN 372 387 CALL ctl_stop( 'STOP', 'ice_run_init : if jpl=1 then nn_monocat should be between 1 and 4' ) 373 388 ENDIF 374 389 ! 375 ! sea-ice timestep and inverse 376 rdt_ice = REAL(nn_fsbc) * rdt 377 r1_rdtice = 1._wp / rdt_ice 378 379 ! inverse of nlay_i and nlay_s 380 r1_nlay_i = 1._wp / REAL( nlay_i, wp ) 390 IF( ln_bdy .AND. ln_limdiachk ) CALL ctl_warn('ice_run_init: online conservation check does not work with BDY') 391 ! 392 ! ! set the choice of ice dynamics 393 IF( lk_c1d .OR. .NOT.ln_limdyn ) THEN 394 nice_dyn = np_dynNO !--- no dynamics 395 ELSE 396 SELECT CASE( nn_limdyn ) 397 CASE( 2 ) 398 IF( nn_monocat /= 2 ) THEN !--- full dynamics (rheology + advection + ridging/rafting + correction) 399 nice_dyn = np_dynFULL 400 ELSE 401 nice_dyn = np_dyn !--- dynamics without ridging/rafting 402 ENDIF 403 CASE( 0 , 1 ) !--- dynamics without ridging/rafting and correction 404 nice_dyn = np_dynPURE 405 END SELECT 406 ENDIF 407 ! !--- simple conservative piling, comparable with LIM2 408 l_piling = nn_limdyn == 1 .OR. ( nn_monocat == 2 .AND. jpl == 1 ) 409 ! 410 rdt_ice = REAL(nn_fsbc) * rdt !--- sea-ice timestep and inverse 411 r1_rdtice = 1._wp / rdt_ice 412 IF( lwp ) WRITE(numout,*) ' ice timestep rdt_ice = ', rdt_ice 413 ! 414 r1_nlay_i = 1._wp / REAL( nlay_i, wp ) !--- inverse of nlay_i and nlay_s 381 415 r1_nlay_s = 1._wp / REAL( nlay_s, wp ) 382 !383 IF( lwp .AND. ln_bdy .AND. ln_limdiachk ) &384 & CALL ctl_warn('online conservation check activated but it does not work with BDY')385 !386 IF( lwp ) WRITE(numout,*) ' ice timestep rdt_ice = ', rdt_ice387 416 ! 388 417 END SUBROUTINE ice_run_init … … 397 426 !! ** input : Namelist namiceitd 398 427 !!------------------------------------------------------------------- 399 INTEGER :: ios ! Local integer output status for namelist read 428 INTEGER :: jl ! dummy loop index 429 INTEGER :: ios ! Local integer output status for namelist read 430 REAL(wp) :: zc1, zc2, zc3, zx1 ! local scalars 431 REAL(wp) :: zhmax, znum, zden, zalpha ! - - 432 !! 400 433 NAMELIST/namiceitd/ rn_himean 401 !402 INTEGER :: jl ! dummy loop index403 REAL(wp) :: zc1, zc2, zc3, zx1 ! local scalars404 REAL(wp) :: zhmax, znum, zden, zalpha !405 434 !!------------------------------------------------------------------ 406 435 ! 407 REWIND( numnam_ice_ref ) 436 REWIND( numnam_ice_ref ) ! Namelist namiceitd in reference namelist : Parameters for ice 408 437 READ ( numnam_ice_ref, namiceitd, IOSTAT = ios, ERR = 905) 409 438 905 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceitd in reference namelist', lwp ) 410 439 411 REWIND( numnam_ice_cfg ) 440 REWIND( numnam_ice_cfg ) ! Namelist namiceitd in configuration namelist : Parameters for ice 412 441 READ ( numnam_ice_cfg, namiceitd, IOSTAT = ios, ERR = 906 ) 413 442 906 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceitd in configuration namelist', lwp ) 414 443 IF(lwm) WRITE ( numoni, namiceitd ) 415 444 ! 416 IF(lwp) THEN 445 IF(lwp) THEN ! control print 417 446 WRITE(numout,*) 418 447 WRITE(numout,*) 'ice_itd_init : Initialization of ice cat distribution ' 419 448 WRITE(numout,*) '~~~~~~~~~~~~' 420 WRITE(numout,*) ' mean ice thickness in the domain rn_himean = ', rn_himean 449 WRITE(numout,*) ' Namelist namicerun : ' 450 WRITE(numout,*) ' mean ice thickness in the domain rn_himean = ', rn_himean 421 451 ENDIF 422 452 ! 423 !---------------------------------- 424 !- Thickness categories boundaries 425 !---------------------------------- 426 ! 427 !== h^(-alpha) function ==! 428 zalpha = 0.05_wp 453 !-----------------------------------! 454 ! Thickness categories boundaries ! 455 !-----------------------------------! 456 ! 457 zalpha = 0.05_wp ! max of each category (from h^(-alpha) function) 429 458 zhmax = 3._wp * rn_himean 430 459 DO jl = 1, jpl … … 441 470 ! 442 471 IF(lwp) WRITE(numout,*) 443 IF(lwp) WRITE(numout,*) ' Thickness category boundaries'444 IF(lwp) WRITE(numout,*) ' hi_max', hi_max(0:jpl)472 IF(lwp) WRITE(numout,*) ' ===>>> resulting thickness category boundaries :' 473 IF(lwp) WRITE(numout,*) ' hi_max(:)= ', hi_max(0:jpl) 445 474 ! 446 475 END SUBROUTINE ice_itd_init 476 447 477 448 478 SUBROUTINE ice_bef … … 472 502 END DO 473 503 END DO 474 CALL lbc_lnk_multi( a_i_b, 'T', 1., v_i_b , 'T', 1., v_s_b , 'T', 1., smv_i_b, 'T', 1.,&475 & oa_i_b, 'T', 1., ht_i_b, 'T', 1., ht_s_b, 'T', 1. )504 CALL lbc_lnk_multi( a_i_b, 'T', 1., v_i_b , 'T', 1., ht_i_b, 'T', 1., smv_i_b, 'T', 1., & 505 & oa_i_b, 'T', 1., v_s_b , 'T', 1., ht_s_b, 'T', 1. ) 476 506 CALL lbc_lnk( e_i_b, 'T', 1. ) 477 507 CALL lbc_lnk( e_s_b, 'T', 1. ) 478 508 509 !!gm Question: here , a_i_b, u_ice and v_ice are defined over the whole domain, 510 !!gm so why not just a copy over the whole domain and no lbc_lnk ? 511 !!gm that is some thing like: 512 ! at_i_b(:,:) = SUM( a_i_b(:,:,:), dim=3 ) 513 ! u_ice_b(:,:) = u_ice(:,:) 514 ! v_ice_b(:,:) = v_ice(:,:) 515 ! idem for the loop above 516 !!gm 479 517 ! ice velocities & total concentration 480 518 DO jj = 2, jpjm1 … … 486 524 END DO 487 525 CALL lbc_lnk_multi( at_i_b, 'T', 1., u_ice_b , 'U', -1., v_ice_b , 'V', -1. ) 488 526 ! 489 527 END SUBROUTINE ice_bef 490 528 … … 552 590 !!---------------------------------------------------------------------- 553 591 CONTAINS 554 !555 592 SUBROUTINE ice_stp ( kt, ksbc ) ! Dummy routine 556 593 INTEGER, INTENT(in) :: kt, ksbc 557 594 WRITE(*,*) 'ice_stp: You should not have seen this print! error?', kt 558 595 END SUBROUTINE ice_stp 559 !560 596 SUBROUTINE ice_init ! Dummy routine 561 597 END SUBROUTINE ice_init
Note: See TracChangeset
for help on using the changeset viewer.
