- Timestamp:
- 2015-08-12T17:46:45+02:00 (5 years ago)
- Location:
- branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/LIM_SRC_3
- Files:
-
- 3 deleted
- 27 edited
- 1 copied
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/LIM_SRC_3/dom_ice.F90
r4161 r5682 5 5 !!====================================================================== 6 6 !! History : 3.0 ! 2003-08 (M. Vancoppenolle) LIM-3 original code 7 !! 4.0! 2011-02 (G. Madec) dynamical allocation7 !! 3.5 ! 2011-02 (G. Madec) dynamical allocation 8 8 !!---------------------------------------------------------------------- 9 USE par_ice ! LIM-3 parameter10 9 USE in_out_manager ! I/O manager 11 10 USE lib_mpp ! MPP library … … 21 20 INTEGER, PUBLIC :: njeq , njeqm1 !: j-index of the equator if it is inside the domain 22 21 23 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: fcor !: coriolis coefficient 24 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: covrai !: sine of geographic latitude 25 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: area !: surface of grid cell 26 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tms, tmi !: temperature mask, mask for stress 27 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tmu, tmv !: mask at u and v velocity points 28 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tmf !: mask at f-point 29 22 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: fcor !: coriolis coefficient 30 23 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: wght !: weight of the 4 neighbours to compute averages 31 24 … … 44 37 !!------------------------------------------------------------------- 45 38 ! 46 ALLOCATE( fcor(jpi,jpj) , & 47 & covrai(jpi,jpj) , area(jpi,jpj) , & 48 & tms (jpi,jpj) , tmi (jpi,jpj) , & 49 & tmu (jpi,jpj) , tmv (jpi,jpj) , & 50 & tmf (jpi,jpj) , & 51 & wght(jpi,jpj,2,2) , STAT = dom_ice_alloc ) 39 ALLOCATE( fcor(jpi,jpj), wght(jpi,jpj,2,2), STAT = dom_ice_alloc ) 52 40 ! 53 41 IF( dom_ice_alloc /= 0 ) CALL ctl_warn( 'dom_ice_alloc: failed to allocate arrays.' ) -
branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/LIM_SRC_3/ice.F90
r4990 r5682 11 11 !! 'key_lim3' LIM-3 sea-ice model 12 12 !!---------------------------------------------------------------------- 13 USE par_ice ! LIM sea-ice parameters14 13 USE in_out_manager ! I/O manager 15 14 USE lib_mpp ! MPP library … … 18 17 PRIVATE 19 18 20 PUBLIC ice_alloc ! Called in iceini.F9019 PUBLIC ice_alloc ! Called in sbc_lim_init 21 20 22 21 !!====================================================================== … … 110 109 !! smv_i | - | Sea ice salt content | ppt.m | 111 110 !! oa_i ! - ! Sea ice areal age content | day | 112 !! e_i ! - ! Ice enthalpy | 10^9 J|111 !! e_i ! - ! Ice enthalpy | J/m2 | 113 112 !! - ! q_i_1d ! Ice enthalpy per unit vol. | J/m3 | 114 !! e_s ! - ! Snow enthalpy | 10^9 J|113 !! e_s ! - ! Snow enthalpy | J/m2 | 115 114 !! - ! q_s_1d ! Snow enthalpy per unit vol. | J/m3 | 116 115 !! | … … 148 147 !! tm_i | - | Mean sea ice temperature | K | 149 148 !! ot_i ! - ! Sea ice areal age content | day | 150 !! et_i ! - ! Total ice enthalpy | 10^9 J|151 !! et_s ! - ! Total snow enthalpy | 10^9 J|149 !! et_i ! - ! Total ice enthalpy | J/m2 | 150 !! et_s ! - ! Total snow enthalpy | J/m2 | 152 151 !! bv_i ! - ! Mean relative brine volume | ??? | 153 152 !!===================================================================== … … 165 164 REAL(wp), PUBLIC :: r1_rdtice !: = 1. / rdt_ice 166 165 167 ! !!** ice-dynamic namelist (namicedyn) ** 168 INTEGER , PUBLIC :: nevp !: number of iterations for subcycling 169 REAL(wp), PUBLIC :: epsd !: tolerance parameter for dynamic 170 REAL(wp), PUBLIC :: om !: relaxation constant 171 REAL(wp), PUBLIC :: cw !: drag coefficient for oceanic stress 172 REAL(wp), PUBLIC :: pstar !: determines ice strength (N/M), Hibler JPO79 173 REAL(wp), PUBLIC :: c_rhg !: determines changes in ice strength 174 REAL(wp), PUBLIC :: creepl !: creep limit : has to be under 1.0e-9 175 REAL(wp), PUBLIC :: ecc !: eccentricity of the elliptical yield curve 176 REAL(wp), PUBLIC :: ahi0 !: sea-ice hor. eddy diffusivity coeff. (m2/s) 177 REAL(wp), PUBLIC :: telast !: timescale for elastic waves (s) 178 REAL(wp), PUBLIC :: relast !: ratio => telast/rdt_ice (1/3 or 1/9 depending on nb of subcycling nevp) 179 REAL(wp), PUBLIC :: alphaevp !: coeficient of the internal stresses 180 REAL(wp), PUBLIC :: hminrhg !: ice volume (a*h, in m) below which ice velocity is set to ocean velocity 166 ! !!** ice-thickness distribution namelist (namiceitd) ** 167 INTEGER , PUBLIC :: nn_catbnd !: categories distribution following: tanh function (1), or h^(-alpha) function (2) 168 REAL(wp), PUBLIC :: rn_himean !: mean thickness of the domain (used to compute the distribution, nn_itdshp = 2 only) 169 170 ! !!** ice-dynamics namelist (namicedyn) ** 171 LOGICAL , PUBLIC :: ln_icestr_bvf !: use brine volume to diminish ice strength 172 INTEGER , PUBLIC :: nn_icestr !: ice strength parameterization (0=Hibler79 1=Rothrock75) 173 INTEGER , PUBLIC :: nn_nevp !: number of iterations for subcycling 174 INTEGER , PUBLIC :: nn_ahi0 !: sea-ice hor. eddy diffusivity coeff. (3 ways of calculation) 175 REAL(wp), PUBLIC :: rn_pe_rdg !: ridging work divided by pot. energy change in ridging, nn_icestr = 1 176 REAL(wp), PUBLIC :: rn_cio !: drag coefficient for oceanic stress 177 REAL(wp), PUBLIC :: rn_pstar !: determines ice strength (N/M), Hibler JPO79 178 REAL(wp), PUBLIC :: rn_crhg !: determines changes in ice strength 179 REAL(wp), PUBLIC :: rn_creepl !: creep limit : has to be under 1.0e-9 180 REAL(wp), PUBLIC :: rn_ecc !: eccentricity of the elliptical yield curve 181 REAL(wp), PUBLIC :: rn_ahi0_ref !: sea-ice hor. eddy diffusivity coeff. (m2/s) 182 REAL(wp), PUBLIC :: rn_relast !: ratio => telast/rdt_ice (1/3 or 1/9 depending on nb of subcycling nevp) 181 183 182 184 ! !!** ice-salinity namelist (namicesal) ** 183 REAL(wp), PUBLIC :: s_i_max !: maximum ice salinity [PSU] 184 REAL(wp), PUBLIC :: s_i_min !: minimum ice salinity [PSU] 185 REAL(wp), PUBLIC :: s_i_0 !: 1st sal. value for the computation of sal .prof. [PSU] 186 REAL(wp), PUBLIC :: s_i_1 !: 2nd sal. value for the computation of sal .prof. [PSU] 187 REAL(wp), PUBLIC :: sal_G !: restoring salinity for gravity drainage [PSU] 188 REAL(wp), PUBLIC :: sal_F !: restoring salinity for flushing [PSU] 189 REAL(wp), PUBLIC :: time_G !: restoring time constant for gravity drainage (= 20 days) [s] 190 REAL(wp), PUBLIC :: time_F !: restoring time constant for gravity drainage (= 10 days) [s] 191 REAL(wp), PUBLIC :: bulk_sal !: bulk salinity (ppt) in case of constant salinity 185 REAL(wp), PUBLIC :: rn_simax !: maximum ice salinity [PSU] 186 REAL(wp), PUBLIC :: rn_simin !: minimum ice salinity [PSU] 187 REAL(wp), PUBLIC :: rn_sal_gd !: restoring salinity for gravity drainage [PSU] 188 REAL(wp), PUBLIC :: rn_sal_fl !: restoring salinity for flushing [PSU] 189 REAL(wp), PUBLIC :: rn_time_gd !: restoring time constant for gravity drainage (= 20 days) [s] 190 REAL(wp), PUBLIC :: rn_time_fl !: restoring time constant for gravity drainage (= 10 days) [s] 191 REAL(wp), PUBLIC :: rn_icesal !: bulk salinity (ppt) in case of constant salinity 192 192 193 193 ! !!** ice-salinity namelist (namicesal) ** 194 INTEGER , PUBLIC :: n um_sal!: salinity configuration used in the model194 INTEGER , PUBLIC :: nn_icesal !: salinity configuration used in the model 195 195 ! ! 1 - constant salinity in both space and time 196 196 ! ! 2 - prognostic salinity (s(z,t)) 197 197 ! ! 3 - salinity profile, constant in time 198 INTEGER , PUBLIC :: thcon_i_swi !: thermal conductivity: =0 Untersteiner (1964) ; =1 Pringle et al (2007) 198 INTEGER , PUBLIC :: nn_ice_thcon !: thermal conductivity: =0 Untersteiner (1964) ; =1 Pringle et al (2007) 199 INTEGER , PUBLIC :: nn_monocat !: virtual ITD mono-category parameterizations (1) or not (0) 200 LOGICAL , PUBLIC :: ln_it_qnsice !: iterate surface flux with changing surface temperature or not (F) 199 201 200 202 ! !!** ice-mechanical redistribution namelist (namiceitdme) 201 REAL(wp), PUBLIC :: Cs !: fraction of shearing energy contributing to ridging 202 REAL(wp), PUBLIC :: Cf !: ratio of ridging work to PE loss 203 REAL(wp), PUBLIC :: fsnowrdg !: fractional snow loss to the ocean during ridging 204 REAL(wp), PUBLIC :: fsnowrft !: fractional snow loss to the ocean during ridging 205 REAL(wp), PUBLIC :: Gstar !: fractional area of young ice contributing to ridging 206 REAL(wp), PUBLIC :: astar !: equivalent of G* for an exponential participation function 207 REAL(wp), PUBLIC :: Hstar !: thickness that determines the maximal thickness of ridged ice 208 REAL(wp), PUBLIC :: hparmeter !: threshold thickness (m) for rafting / ridging 209 REAL(wp), PUBLIC :: Craft !: coefficient for smoothness of the hyperbolic tangent in rafting 210 REAL(wp), PUBLIC :: ridge_por !: initial porosity of ridges (0.3 regular value) 211 REAL(wp), PUBLIC :: betas !: coef. for partitioning of snowfall between leads and sea ice 212 REAL(wp), PUBLIC :: kappa_i !: coef. for the extinction of radiation Grenfell et al. (2006) [1/m] 213 REAL(wp), PUBLIC :: nconv_i_thd !: maximal number of iterations for heat diffusion 214 REAL(wp), PUBLIC :: maxer_i_thd !: maximal tolerated error (C) for heat diffusion 203 REAL(wp), PUBLIC :: rn_cs !: fraction of shearing energy contributing to ridging 204 REAL(wp), PUBLIC :: rn_fsnowrdg !: fractional snow loss to the ocean during ridging 205 REAL(wp), PUBLIC :: rn_fsnowrft !: fractional snow loss to the ocean during ridging 206 REAL(wp), PUBLIC :: rn_gstar !: fractional area of young ice contributing to ridging 207 REAL(wp), PUBLIC :: rn_astar !: equivalent of G* for an exponential participation function 208 REAL(wp), PUBLIC :: rn_hstar !: thickness that determines the maximal thickness of ridged ice 209 REAL(wp), PUBLIC :: rn_hraft !: threshold thickness (m) for rafting / ridging 210 REAL(wp), PUBLIC :: rn_craft !: coefficient for smoothness of the hyperbolic tangent in rafting 211 REAL(wp), PUBLIC :: rn_por_rdg !: initial porosity of ridges (0.3 regular value) 212 REAL(wp), PUBLIC :: rn_betas !: coef. for partitioning of snowfall between leads and sea ice 213 REAL(wp), PUBLIC :: rn_kappa_i !: coef. for the extinction of radiation Grenfell et al. (2006) [1/m] 214 REAL(wp), PUBLIC :: nn_conv_dif !: maximal number of iterations for heat diffusion 215 REAL(wp), PUBLIC :: rn_terr_dif !: maximal tolerated error (C) for heat diffusion 215 216 216 217 ! !!** ice-mechanical redistribution namelist (namiceitdme) 217 INTEGER , PUBLIC :: ridge_scheme_swi !: scheme used for ice ridging218 INTEGER , PUBLIC :: raft_swi !: rafting of ice or not219 INTEGER , PUBLIC :: partfun_swi !: participation function: =0 Thorndike et al. (1975), =1 Lipscomb et al. (2007) 220 INTEGER , PUBLIC :: brinstren_swi !: use brine volume to diminish ice strength221 222 REAL(wp), PUBLIC :: usecc2 !: = 1.0 / ( ecc * ecc )223 REAL(wp), PUBLIC :: r hoco !: = rau0 * cw224 218 LOGICAL , PUBLIC :: ln_rafting !: rafting of ice or not 219 INTEGER , PUBLIC :: nn_partfun !: participation function: =0 Thorndike et al. (1975), =1 Lipscomb et al. (2007) 220 221 REAL(wp), PUBLIC :: usecc2 !: = 1.0 / ( rn_ecc * rn_ecc ) 222 REAL(wp), PUBLIC :: rhoco !: = rau0 * cio 223 REAL(wp), PUBLIC :: r1_nlay_i !: 1 / nlay_i 224 REAL(wp), PUBLIC :: r1_nlay_s !: 1 / nlay_s 225 ! 225 226 ! !!** switch for presence of ice or not 226 227 REAL(wp), PUBLIC :: rswitch 227 228 ! 228 229 ! !!** define some parameters 229 REAL(wp), PUBLIC, PARAMETER :: unit_fac = 1.e+09_wp !: conversion factor for ice / snow enthalpy230 230 REAL(wp), PUBLIC, PARAMETER :: epsi06 = 1.e-06_wp !: small number 231 231 REAL(wp), PUBLIC, PARAMETER :: epsi10 = 1.e-10_wp !: small number … … 266 266 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: wfx_res !: residual component of wfx_ice [kg/m2] 267 267 268 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: afx_tot !: ice concentration tendency (total) [s-1] 269 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: afx_thd !: ice concentration tendency (thermodynamics) [s-1] 270 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: afx_dyn !: ice concentration tendency (dynamics) [s-1] 271 268 272 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: sfx_bog !: salt flux due to ice growth/melt [PSU/m2/s] 269 273 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: sfx_bom !: salt flux due to ice growth/melt [PSU/m2/s] … … 282 286 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hfx_snw !: heat flux for snow melt 283 287 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hfx_err !: heat flux error after heat diffusion 288 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hfx_err_dif !: heat flux remaining due to change in non-solar flux 284 289 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hfx_err_rem !: heat flux error after heat remapping 285 290 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hfx_in !: heat flux available for thermo transformations … … 296 301 297 302 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ftr_ice !: transmitted solar radiation under ice 298 299 ! temporary arrays for dummy version of the code300 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: dh_i_surf2D, dh_i_bott2D, q_s301 303 302 304 !!-------------------------------------------------------------------------- … … 333 335 334 336 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: t_i !: ice temperatures [K] 335 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: e_i !: ice thermal contents [Giga J]337 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: e_i !: ice thermal contents [J/m2] 336 338 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: s_i !: ice salinities [PSU] 337 339 … … 356 358 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: e_i_b !: ice temperatures 357 359 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: u_ice_b, v_ice_b !: ice velocity 358 359 360 !!-------------------------------------------------------------------------- 361 !! * Increment of global variables 362 !!-------------------------------------------------------------------------- 360 361 !!-------------------------------------------------------------------------- 362 !! * Ice thickness distribution variables 363 !!-------------------------------------------------------------------------- 364 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: hi_max !: Boundary of ice thickness categories in thickness space 365 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: hi_mean !: Mean ice thickness in catgories 366 367 !!-------------------------------------------------------------------------- 368 !! * Ice Run 369 !!-------------------------------------------------------------------------- 370 ! !!: ** Namelist namicerun read in sbc_lim_init ** 371 INTEGER , PUBLIC :: jpl !: number of ice categories 372 INTEGER , PUBLIC :: nlay_i !: number of ice layers 373 INTEGER , PUBLIC :: nlay_s !: number of snow layers 374 CHARACTER(len=32), PUBLIC :: cn_icerst_in !: suffix of ice restart name (input) 375 CHARACTER(len=256), PUBLIC :: cn_icerst_indir !: ice restart input directory 376 CHARACTER(len=32), PUBLIC :: cn_icerst_out !: suffix of ice restart name (output) 377 CHARACTER(len=256), PUBLIC :: cn_icerst_outdir!: ice restart output directory 378 LOGICAL , PUBLIC :: ln_limdyn !: flag for ice dynamics (T) or not (F) 379 LOGICAL , PUBLIC :: ln_icectl !: flag for sea-ice points output (T) or not (F) 380 REAL(wp) , PUBLIC :: rn_amax !: maximum ice concentration 381 INTEGER , PUBLIC :: iiceprt !: debug i-point 382 INTEGER , PUBLIC :: jiceprt !: debug j-point 383 ! 384 !!-------------------------------------------------------------------------- 385 !! * Ice diagnostics 386 !!-------------------------------------------------------------------------- 387 ! Increment of global variables 363 388 ! thd refers to changes induced by thermodynamics 364 389 ! trp '' '' '' advection (transport of ice) 365 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: d_a_i_thd , d_a_i_trp !: icefractions 366 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: d_v_s_thd , d_v_s_trp !: snow volume 367 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: d_v_i_thd , d_v_i_trp !: ice volume 368 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: d_smv_i_thd, d_smv_i_trp !: 369 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: d_sm_i_fl , d_sm_i_gd !: 370 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: d_sm_i_se , d_sm_i_si , d_sm_i_la !: 371 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: d_oa_i_thd , d_oa_i_trp !: 372 373 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: d_e_s_thd , d_e_s_trp !: 374 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: d_e_i_thd , d_e_i_trp !: 375 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: d_u_ice_dyn, d_v_ice_dyn !: ice velocity 376 377 !!-------------------------------------------------------------------------- 378 !! * Ice thickness distribution variables 379 !!-------------------------------------------------------------------------- 380 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: hi_max !: Boundary of ice thickness categories in thickness space 381 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: hi_mean !: Mean ice thickness in catgories 382 383 !!-------------------------------------------------------------------------- 384 !! * Ice Run 385 !!-------------------------------------------------------------------------- 386 ! !!: ** Namelist namicerun read in iceini ** 387 CHARACTER(len=32) , PUBLIC :: cn_icerst_in !: suffix of ice restart name (input) 388 CHARACTER(len=32) , PUBLIC :: cn_icerst_out !: suffix of ice restart name (output) 389 LOGICAL , PUBLIC :: ln_limdyn !: flag for ice dynamics (T) or not (F) 390 LOGICAL , PUBLIC :: ln_nicep !: flag for sea-ice points output (T) or not (F) 391 REAL(wp) , PUBLIC :: cai !: atmospheric drag over sea ice 392 REAL(wp) , PUBLIC :: cao !: atmospheric drag over ocean 393 REAL(wp) , PUBLIC :: amax !: maximum ice concentration 390 LOGICAL , PUBLIC :: ln_limdiahsb !: flag for ice diag (T) or not (F) 391 LOGICAL , PUBLIC :: ln_limdiaout !: flag for ice diag (T) or not (F) 392 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_trp_vi !: transport of ice volume 393 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_trp_vs !: transport of snw volume 394 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_trp_ei !: transport of ice enthalpy (W/m2) 395 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_trp_es !: transport of snw enthalpy (W/m2) 396 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_trp_smv !: transport of salt content 394 397 ! 395 !!-------------------------------------------------------------------------- 396 !! * Ice diagnostics 397 !!-------------------------------------------------------------------------- 398 !! Check if everything down here is necessary 399 LOGICAL , PUBLIC :: ln_limdiahsb !: flag for ice diag (T) or not (F) 400 LOGICAL , PUBLIC :: ln_limdiaout !: flag for ice diag (T) or not (F) 401 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: dv_dt_thd !: thermodynamic growth rates 402 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_trp_vi !: transport of ice volume 403 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_trp_vs !: transport of snw volume 404 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_trp_ei !: transport of ice enthalpy (W/m2) 405 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_trp_es !: transport of snw enthalpy (W/m2) 398 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_heat !: snw/ice heat content variation [W/m2] 399 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_smvi !: ice salt content variation [] 400 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_vice !: ice volume variation [m/s] 401 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_vsnw !: snw volume variation [m/s] 406 402 ! 407 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_heat_dhc !: snw/ice heat content variation [W/m2]408 !409 INTEGER , PUBLIC :: jiindx, jjindx !: indexes of the debugging point410 411 403 !!---------------------------------------------------------------------- 412 404 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2010) … … 422 414 INTEGER :: ice_alloc 423 415 ! 424 INTEGER :: ierr(1 9), ii416 INTEGER :: ierr(17), ii 425 417 !!----------------------------------------------------------------- 426 418 … … 439 431 440 432 ii = ii + 1 441 ALLOCATE( sist (jpi,jpj) , icethi (jpi,jpj) , t_bo (jpi,jpj) , &442 & frld (jpi,jpj) , pfrld (jpi,jpj) , phicif (jpi,jpj) , &443 & wfx_snw(jpi,jpj) , wfx_ice(jpi,jpj) , wfx_sub(jpi,jpj) , &433 ALLOCATE( sist (jpi,jpj) , icethi (jpi,jpj) , t_bo (jpi,jpj) , & 434 & frld (jpi,jpj) , pfrld (jpi,jpj) , phicif (jpi,jpj) , & 435 & wfx_snw(jpi,jpj) , wfx_ice(jpi,jpj) , wfx_sub(jpi,jpj) , & 444 436 & wfx_bog(jpi,jpj) , wfx_dyn(jpi,jpj) , wfx_bom(jpi,jpj) , wfx_sum(jpi,jpj) , & 445 & wfx_res(jpi,jpj) , wfx_sni(jpi,jpj) , wfx_opw(jpi,jpj) , wfx_spr(jpi,jpj) , qlead (jpi,jpj) , & 446 & fhtur (jpi,jpj) , ftr_ice(jpi,jpj,jpl) , & 447 & sfx_res(jpi,jpj) , sfx_bri(jpi,jpj) , sfx_dyn(jpi,jpj) , & 448 & sfx_bog(jpi,jpj) , sfx_bom(jpi,jpj) , sfx_sum(jpi,jpj) , sfx_sni(jpi,jpj) , sfx_opw(jpi,jpj) , & 449 & hfx_res(jpi,jpj) , hfx_snw(jpi,jpj) , hfx_sub(jpi,jpj) , hfx_err(jpi,jpj) , hfx_err_rem(jpi,jpj), & 450 & hfx_in (jpi,jpj) , hfx_out(jpi,jpj) , fhld(jpi,jpj) , & 451 & hfx_sum(jpi,jpj) , hfx_bom(jpi,jpj) , hfx_bog(jpi,jpj) , hfx_dif(jpi,jpj) , hfx_opw(jpi,jpj) , & 437 & wfx_res(jpi,jpj) , wfx_sni(jpi,jpj) , wfx_opw(jpi,jpj) , wfx_spr(jpi,jpj) , & 438 & afx_tot(jpi,jpj) , afx_thd(jpi,jpj), afx_dyn(jpi,jpj) , & 439 & fhtur (jpi,jpj) , ftr_ice(jpi,jpj,jpl), qlead (jpi,jpj) , & 440 & sfx_res(jpi,jpj) , sfx_bri(jpi,jpj) , sfx_dyn(jpi,jpj) , & 441 & sfx_bog(jpi,jpj) , sfx_bom(jpi,jpj) , sfx_sum(jpi,jpj) , sfx_sni(jpi,jpj) , sfx_opw(jpi,jpj) , & 442 & hfx_res(jpi,jpj) , hfx_snw(jpi,jpj) , hfx_sub(jpi,jpj) , hfx_err(jpi,jpj) , & 443 & hfx_err_dif(jpi,jpj) , hfx_err_rem(jpi,jpj) , & 444 & hfx_in (jpi,jpj) , hfx_out(jpi,jpj) , fhld(jpi,jpj) , & 445 & hfx_sum(jpi,jpj) , hfx_bom(jpi,jpj) , hfx_bog(jpi,jpj) , hfx_dif(jpi,jpj) , hfx_opw(jpi,jpj) , & 452 446 & hfx_thd(jpi,jpj) , hfx_dyn(jpi,jpj) , hfx_spr(jpi,jpj) , STAT=ierr(ii) ) 453 447 … … 464 458 & bv_i (jpi,jpj) , smt_i(jpi,jpj) , STAT=ierr(ii) ) 465 459 ii = ii + 1 466 ALLOCATE( t_s(jpi,jpj,nlay_s,jpl) , & 467 & e_s(jpi,jpj,nlay_s,jpl) , STAT=ierr(ii) ) 468 ii = ii + 1 469 ALLOCATE( t_i(jpi,jpj,nlay_i+1,jpl) , e_i(jpi,jpj,nlay_i+1,jpl) , s_i(jpi,jpj,nlay_i+1,jpl) , STAT=ierr(ii) ) 460 ALLOCATE( t_s(jpi,jpj,nlay_s,jpl) , e_s(jpi,jpj,nlay_s,jpl) , STAT=ierr(ii) ) 461 ii = ii + 1 462 ALLOCATE( t_i(jpi,jpj,nlay_i,jpl) , e_i(jpi,jpj,nlay_i,jpl) , s_i(jpi,jpj,nlay_i,jpl) , STAT=ierr(ii) ) 470 463 471 464 ! * Moments for advection … … 483 476 & STAT=ierr(ii) ) 484 477 ii = ii + 1 485 ALLOCATE( sxe (jpi,jpj,nlay_i +1,jpl) , sye (jpi,jpj,nlay_i+1,jpl) , sxxe(jpi,jpj,nlay_i+1,jpl) , &486 & syye(jpi,jpj,nlay_i +1,jpl) , sxye(jpi,jpj,nlay_i+1,jpl), STAT=ierr(ii) )478 ALLOCATE( sxe (jpi,jpj,nlay_i,jpl) , sye (jpi,jpj,nlay_i,jpl) , sxxe(jpi,jpj,nlay_i,jpl) , & 479 & syye(jpi,jpj,nlay_i,jpl) , sxye(jpi,jpj,nlay_i,jpl) , STAT=ierr(ii) ) 487 480 488 481 ! * Old values of global variables 489 482 ii = ii + 1 490 483 ALLOCATE( v_s_b (jpi,jpj,jpl) , v_i_b (jpi,jpj,jpl) , e_s_b(jpi,jpj,nlay_s,jpl) , & 491 & a_i_b (jpi,jpj,jpl) , smv_i_b(jpi,jpj,jpl) , e_i_b(jpi,jpj,nlay_i+1 ,jpl) , & 492 & oa_i_b (jpi,jpj,jpl) , & 493 & u_ice_b(jpi,jpj) , v_ice_b(jpi,jpj) , STAT=ierr(ii) ) 494 495 ! * Increment of global variables 496 ii = ii + 1 497 ALLOCATE( d_a_i_thd(jpi,jpj,jpl) , d_a_i_trp (jpi,jpj,jpl) , d_v_s_thd (jpi,jpj,jpl) , d_v_s_trp (jpi,jpj,jpl) , & 498 & d_v_i_thd(jpi,jpj,jpl) , d_v_i_trp (jpi,jpj,jpl) , d_smv_i_thd(jpi,jpj,jpl) , d_smv_i_trp(jpi,jpj,jpl) , & 499 & d_sm_i_fl(jpi,jpj,jpl) , d_sm_i_gd (jpi,jpj,jpl) , d_sm_i_se (jpi,jpj,jpl) , d_sm_i_si (jpi,jpj,jpl) , & 500 & d_sm_i_la(jpi,jpj,jpl) , d_oa_i_thd(jpi,jpj,jpl) , d_oa_i_trp (jpi,jpj,jpl) , & 501 & STAT=ierr(ii) ) 502 ii = ii + 1 503 ALLOCATE( d_e_s_thd(jpi,jpj,nlay_s,jpl) , d_e_i_thd(jpi,jpj,nlay_i,jpl) , d_u_ice_dyn(jpi,jpj) , & 504 & d_e_s_trp(jpi,jpj,nlay_s,jpl) , d_e_i_trp(jpi,jpj,nlay_i,jpl) , d_v_ice_dyn(jpi,jpj) , STAT=ierr(ii) ) 484 & a_i_b (jpi,jpj,jpl) , smv_i_b(jpi,jpj,jpl) , e_i_b(jpi,jpj,nlay_i,jpl) , & 485 & oa_i_b (jpi,jpj,jpl) , u_ice_b(jpi,jpj) , v_ice_b(jpi,jpj) , STAT=ierr(ii) ) 505 486 506 487 ! * Ice thickness distribution variables … … 510 491 ! * Ice diagnostics 511 492 ii = ii + 1 512 ALLOCATE( d v_dt_thd(jpi,jpj,jpl), &513 & diag_trp_ vi(jpi,jpj), diag_trp_vs (jpi,jpj), diag_trp_ei(jpi,jpj), &514 & diag_ trp_es(jpi,jpj), diag_heat_dhc(jpi,jpj),STAT=ierr(ii) )493 ALLOCATE( diag_trp_vi(jpi,jpj), diag_trp_vs (jpi,jpj), diag_trp_ei(jpi,jpj), & 494 & diag_trp_es(jpi,jpj), diag_trp_smv(jpi,jpj), diag_heat (jpi,jpj), & 495 & diag_smvi (jpi,jpj), diag_vice (jpi,jpj), diag_vsnw (jpi,jpj), STAT=ierr(ii) ) 515 496 516 497 ice_alloc = MAXVAL( ierr(:) ) -
branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/LIM_SRC_3/limadv.F90
r4990 r5682 63 63 !! 64 64 INTEGER :: ji, jj ! dummy loop indices 65 REAL(wp) :: zs1max, zrdt, zslpmax, ztemp , zin0! local scalars65 REAL(wp) :: zs1max, zrdt, zslpmax, ztemp ! local scalars 66 66 REAL(wp) :: zs1new, zalf , zalfq , zbt ! - - 67 67 REAL(wp) :: zs2new, zalf1, zalf1q, zbt1 ! - - … … 85 85 zs2new = MIN( 2.0 * zslpmax - 0.3334 * ABS( zs1new ), & 86 86 & MAX( ABS( zs1new ) - zslpmax, psxx(ji,jj) ) ) 87 zin0 = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tms(ji,jj) ! Case of empty boxes & Apply mask87 rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1) ! Case of empty boxes & Apply mask 88 88 89 89 ps0 (ji,jj) = zslpmax 90 psx (ji,jj) = zs1new * zin091 psxx(ji,jj) = zs2new * zin092 psy (ji,jj) = psy (ji,jj) * zin093 psyy(ji,jj) = psyy(ji,jj) * zin094 psxy(ji,jj) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj) ) ) * zin090 psx (ji,jj) = zs1new * rswitch 91 psxx(ji,jj) = zs2new * rswitch 92 psy (ji,jj) = psy (ji,jj) * rswitch 93 psyy(ji,jj) = psyy(ji,jj) * rswitch 94 psxy(ji,jj) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj) ) ) * rswitch 95 95 END DO 96 96 END DO 97 97 98 98 ! Initialize volumes of boxes (=area if adv_x first called, =psm otherwise) 99 psm (:,:) = MAX( pcrh * area(:,:) + ( 1.0 - pcrh ) * psm(:,:) , epsi20 )99 psm (:,:) = MAX( pcrh * e12t(:,:) + ( 1.0 - pcrh ) * psm(:,:) , epsi20 ) 100 100 101 101 ! Calculate fluxes and moments between boxes i<-->i+1 … … 207 207 208 208 !-- Lateral boundary conditions 209 CALL lbc_lnk ( psm , 'T', 1. ) ; CALL lbc_lnk( ps0 , 'T', 1. )210 CALL lbc_lnk( psx , 'T', -1. ) ; CALL lbc_lnk( psy , 'T', -1. )! caution gradient ==> the sign changes211 CALL lbc_lnk( psxx, 'T', 1. ) ; CALL lbc_lnk( psyy, 'T', 1. )212 CALL lbc_lnk(psxy, 'T', 1. )209 CALL lbc_lnk_multi( psm , 'T', 1., ps0 , 'T', 1. & 210 & , psx , 'T', -1., psy , 'T', -1. & ! caution gradient ==> the sign changes 211 & , psxx, 'T', 1., psyy, 'T', 1. & 212 & , psxy, 'T', 1. ) 213 213 214 214 IF(ln_ctl) THEN … … 248 248 !! 249 249 INTEGER :: ji, jj ! dummy loop indices 250 REAL(wp) :: zs1max, zrdt, zslpmax, ztemp , zin0! temporary scalars250 REAL(wp) :: zs1max, zrdt, zslpmax, ztemp ! temporary scalars 251 251 REAL(wp) :: zs1new, zalf , zalfq , zbt ! - - 252 252 REAL(wp) :: zs2new, zalf1, zalf1q, zbt1 ! - - … … 270 270 zs2new = MIN( ( 2.0 * zslpmax - 0.3334 * ABS( zs1new ) ), & 271 271 & MAX( ABS( zs1new )-zslpmax, psyy(ji,jj) ) ) 272 zin0 = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tms(ji,jj) ! Case of empty boxes & Apply mask272 rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1) ! Case of empty boxes & Apply mask 273 273 ! 274 274 ps0 (ji,jj) = zslpmax 275 psx (ji,jj) = psx (ji,jj) * zin0276 psxx(ji,jj) = psxx(ji,jj) * zin0277 psy (ji,jj) = zs1new * zin0278 psyy(ji,jj) = zs2new * zin0279 psxy(ji,jj) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj) ) ) * zin0275 psx (ji,jj) = psx (ji,jj) * rswitch 276 psxx(ji,jj) = psxx(ji,jj) * rswitch 277 psy (ji,jj) = zs1new * rswitch 278 psyy(ji,jj) = zs2new * rswitch 279 psxy(ji,jj) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj) ) ) * rswitch 280 280 END DO 281 281 END DO 282 282 283 283 ! Initialize volumes of boxes (=area if adv_x first called, =psm otherwise) 284 psm(:,:) = MAX( pcrh * area(:,:) + ( 1.0 - pcrh ) * psm(:,:) , epsi20 )284 psm(:,:) = MAX( pcrh * e12t(:,:) + ( 1.0 - pcrh ) * psm(:,:) , epsi20 ) 285 285 286 286 ! Calculate fluxes and moments between boxes j<-->j+1 … … 393 393 394 394 !-- Lateral boundary conditions 395 CALL lbc_lnk ( psm , 'T', 1. ) ; CALL lbc_lnk( ps0 , 'T', 1. )396 CALL lbc_lnk( psx , 'T', -1. ) ; CALL lbc_lnk( psy , 'T', -1. )! caution gradient ==> the sign changes397 CALL lbc_lnk( psxx, 'T', 1. ) ; CALL lbc_lnk( psyy, 'T', 1. )398 CALL lbc_lnk(psxy, 'T', 1. )395 CALL lbc_lnk_multi( psm , 'T', 1., ps0 , 'T', 1. & 396 & , psx , 'T', -1., psy , 'T', -1. & ! caution gradient ==> the sign changes 397 & , psxx, 'T', 1., psyy, 'T', 1. & 398 & , psxy, 'T', 1. ) 399 399 400 400 IF(ln_ctl) THEN -
branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/LIM_SRC_3/limcons.F90
r4873 r5682 6 6 !! History : - ! Original code from William H. Lipscomb, LANL 7 7 !! 3.0 ! 2004-06 (M. Vancoppenolle) Energy Conservation 8 !! 4.0! 2011-02 (G. Madec) add mpp considerations8 !! 3.5 ! 2011-02 (G. Madec) add mpp considerations 9 9 !! - ! 2014-05 (C. Rousset) add lim_cons_hsm 10 !! - ! 2015-03 (C. Rousset) add lim_cons_final 10 11 !!---------------------------------------------------------------------- 11 12 #if defined key_lim3 … … 16 17 !!---------------------------------------------------------------------- 17 18 USE phycst ! physical constants 18 USE par_ice ! LIM-3 parameter19 19 USE ice ! LIM-3 variables 20 20 USE dom_ice ! LIM-3 domain … … 23 23 USE lib_mpp ! MPP library 24 24 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 25 USE sbc_oce , ONLY : sfx ! Surface boundary condition: ocean fields 25 26 26 27 IMPLICIT NONE … … 31 32 PUBLIC lim_cons_check 32 33 PUBLIC lim_cons_hsm 34 PUBLIC lim_cons_final 33 35 34 36 !!---------------------------------------------------------------------- … … 73 75 !! ** Method : Arithmetics 74 76 !!--------------------------------------------------------------------- 75 INTEGER 76 INTEGER 77 REAL(wp), DIMENSION(jpi,jpj,nlay_i +1,jpl), INTENT(in ) :: pin!: input field78 REAL(wp), DIMENSION(jpi,jpj) 77 INTEGER , INTENT(in ) :: ksum !: number of categories 78 INTEGER , INTENT(in ) :: klay !: number of vertical layers 79 REAL(wp), DIMENSION(jpi,jpj,nlay_i,jpl), INTENT(in ) :: pin !: input field 80 REAL(wp), DIMENSION(jpi,jpj) , INTENT( out) :: pout !: output field 79 81 ! 80 82 INTEGER :: jk, jl ! dummy loop indices … … 156 158 157 159 SUBROUTINE lim_cons_hsm( icount, cd_routine, zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b ) 158 !!------------------------------------------------------------------- 159 !! *** ROUTINE lim_cons_hsm *** 160 !! 161 !! ** Purpose : Test the conservation of heat, salt and mass for each routine 162 !! 163 !! ** Method : 164 !!--------------------------------------------------------------------- 165 INTEGER , INTENT(in) :: icount ! determine wether this is the beggining of the routine (0) or the end (1) 166 CHARACTER(len=*), INTENT(in) :: cd_routine ! name of the routine 160 !!-------------------------------------------------------------------------------------------------------- 161 !! *** ROUTINE lim_cons_hsm *** 162 !! 163 !! ** Purpose : Test the conservation of heat, salt and mass for each ice routine 164 !! + test if ice concentration and volume are > 0 165 !! 166 !! ** Method : This is an online diagnostics which can be activated with ln_limdiahsb=true 167 !! It prints in ocean.output if there is a violation of conservation at each time-step 168 !! The thresholds (zv_sill, zs_sill, zh_sill) which determine violations are set to 169 !! a minimum of 1 mm of ice (over the ice area) that is lost/gained spuriously during 100 years. 170 !! For salt and heat thresholds, ice is considered to have a salinity of 10 171 !! and a heat content of 3e5 J/kg (=latent heat of fusion) 172 !!-------------------------------------------------------------------------------------------------------- 173 INTEGER , INTENT(in) :: icount ! determine wether this is the beggining of the routine (0) or the end (1) 174 CHARACTER(len=*), INTENT(in) :: cd_routine ! name of the routine 167 175 REAL(wp) , INTENT(inout) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 168 176 REAL(wp) :: zvi, zsmv, zei, zfs, zfw, zft 169 177 REAL(wp) :: zvmin, zamin, zamax 178 REAL(wp) :: zvtrp, zetrp 179 REAL(wp) :: zarea, zv_sill, zs_sill, zh_sill 180 REAL(wp), PARAMETER :: zconv = 1.e-9 ! convert W to GW and kg to Mt 170 181 171 182 IF( icount == 0 ) THEN 172 183 173 zvi_b = glob_sum( SUM( v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 ) * area(:,:) * tms(:,:) ) 174 zsmv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 175 zei_b = glob_sum( SUM( e_i(:,:,1:nlay_i,:), dim=3 ) + SUM( e_s(:,:,1:nlay_s,:), dim=3 ) ) 176 zfw_b = glob_sum( - ( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) + & 177 & wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) + wfx_sub(:,:) + wfx_spr(:,:) & 178 & ) * area(:,:) * tms(:,:) ) 179 zfs_b = glob_sum( ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + & 180 & sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) & 181 & ) * area(:,:) * tms(:,:) ) 182 zft_b = glob_sum( ( hfx_sum(:,:) + hfx_bom(:,:) + hfx_bog(:,:) + hfx_dif(:,:) + hfx_opw(:,:) + hfx_snw(:,:) & 183 & - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) - hfx_sub(:,:) - hfx_spr(:,:) & 184 & ) * area(:,:) / unit_fac * tms(:,:) ) 184 ! salt flux 185 zfs_b = glob_sum( ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + & 186 & sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) & 187 & ) * e12t(:,:) * tmask(:,:,1) * zconv ) 188 189 ! water flux 190 zfw_b = glob_sum( -( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) + & 191 & wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) + wfx_sub(:,:) + wfx_spr(:,:) & 192 & ) * e12t(:,:) * tmask(:,:,1) * zconv ) 193 194 ! heat flux 195 zft_b = glob_sum( ( hfx_sum(:,:) + hfx_bom(:,:) + hfx_bog(:,:) + hfx_dif(:,:) + hfx_opw(:,:) + hfx_snw(:,:) & 196 & - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) - hfx_sub(:,:) - hfx_spr(:,:) & 197 & ) * e12t(:,:) * tmask(:,:,1) * zconv ) 198 199 zvi_b = glob_sum( SUM( v_i * rhoic + v_s * rhosn, dim=3 ) * e12t * tmask(:,:,1) * zconv ) 200 201 zsmv_b = glob_sum( SUM( smv_i * rhoic , dim=3 ) * e12t * tmask(:,:,1) * zconv ) 202 203 zei_b = glob_sum( ( SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 ) + & 204 & SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 ) & 205 ) * e12t * tmask(:,:,1) * zconv ) 185 206 186 207 ELSEIF( icount == 1 ) THEN 187 208 188 zfs = glob_sum( ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + & 189 & sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) & 190 & ) * area(:,:) * tms(:,:) ) - zfs_b 191 zfw = glob_sum( - ( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) + & 192 & wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) + wfx_sub(:,:) + wfx_spr(:,:) & 193 & ) * area(:,:) * tms(:,:) ) - zfw_b 194 zft = glob_sum( ( hfx_sum(:,:) + hfx_bom(:,:) + hfx_bog(:,:) + hfx_dif(:,:) + hfx_opw(:,:) + hfx_snw(:,:) & 195 & - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) - hfx_sub(:,:) - hfx_spr(:,:) & 196 & ) * area(:,:) / unit_fac * tms(:,:) ) - zft_b 209 ! salt flux 210 zfs = glob_sum( ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + & 211 & sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) & 212 & ) * e12t(:,:) * tmask(:,:,1) * zconv ) - zfs_b 213 214 ! water flux 215 zfw = glob_sum( -( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) + & 216 & wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) + wfx_sub(:,:) + wfx_spr(:,:) & 217 & ) * e12t(:,:) * tmask(:,:,1) * zconv ) - zfw_b 218 219 ! heat flux 220 zft = glob_sum( ( hfx_sum(:,:) + hfx_bom(:,:) + hfx_bog(:,:) + hfx_dif(:,:) + hfx_opw(:,:) + hfx_snw(:,:) & 221 & - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) - hfx_sub(:,:) - hfx_spr(:,:) & 222 & ) * e12t(:,:) * tmask(:,:,1) * zconv ) - zft_b 197 223 198 zvi = ( glob_sum( SUM( v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 ) * area(:,:) * tms(:,:) ) - zvi_b ) * r1_rdtice - zfw 199 zsmv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zsmv_b ) * r1_rdtice + ( zfs / rhoic ) 200 zei = glob_sum( SUM( e_i(:,:,1:nlay_i,:), dim=3 ) + SUM( e_s(:,:,1:nlay_s,:), dim=3 ) ) * r1_rdtice - zei_b * r1_rdtice + zft 201 202 zvmin = glob_min(v_i) 203 zamax = glob_max(SUM(a_i,dim=3)) 204 zamin = glob_min(a_i) 205 224 ! outputs 225 zvi = ( ( glob_sum( SUM( v_i * rhoic + v_s * rhosn, dim=3 ) & 226 & * e12t * tmask(:,:,1) * zconv ) - zvi_b ) * r1_rdtice - zfw ) * rday 227 228 zsmv = ( ( glob_sum( SUM( smv_i * rhoic , dim=3 ) & 229 & * e12t * tmask(:,:,1) * zconv ) - zsmv_b ) * r1_rdtice + zfs ) * rday 230 231 zei = glob_sum( ( SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 ) + & 232 & SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 ) & 233 & ) * e12t * tmask(:,:,1) * zconv ) * r1_rdtice - zei_b * r1_rdtice + zft 234 235 ! zvtrp and zetrp must be close to 0 if the advection scheme is conservative 236 zvtrp = glob_sum( ( diag_trp_vi * rhoic + diag_trp_vs * rhosn ) * e12t * tmask(:,:,1) * zconv ) * rday 237 zetrp = glob_sum( ( diag_trp_ei + diag_trp_es ) * e12t * tmask(:,:,1) * zconv ) 238 239 zvmin = glob_min( v_i ) 240 zamax = glob_max( SUM( a_i, dim=3 ) ) 241 zamin = glob_min( a_i ) 242 243 ! set threshold values and calculate the ice area (+epsi10 to set a threshold > 0 when there is no ice) 244 zarea = glob_sum( SUM( a_i + epsi10, dim=3 ) * e12t * zconv ) ! in 1.e9 m2 245 zv_sill = zarea * 2.5e-5 246 zs_sill = zarea * 25.e-5 247 zh_sill = zarea * 10.e-5 248 206 249 IF(lwp) THEN 207 IF ( ABS( zvi ) > 1.e-4 ) WRITE(numout,*) 'violation volume [kg/day] (',cd_routine,') = ',(zvi * rday)208 IF ( ABS( zsmv ) > 1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (',cd_routine,') = ',(zsmv * rday)209 IF ( ABS( zei ) > 1. ) WRITE(numout,*) 'violation enthalpy [1e9 J] (',cd_routine,') = ',(zei)210 IF ( zvmin < 0. ) WRITE(numout,*) 'violation v_i<0 [m] (',cd_routine,') = ',(zvmin)211 IF( cd_routine /= 'limtrp' .AND. cd_routine /= 'limitd_me' .AND. zamax > amax+1.e-10 ) THEN212 WRITE(numout,*) 'violation a_i>amax (',cd_routine,') = ',zamax250 IF ( ABS( zvi ) > zv_sill ) WRITE(numout,*) 'violation volume [Mt/day] (',cd_routine,') = ',zvi 251 IF ( ABS( zsmv ) > zs_sill ) WRITE(numout,*) 'violation saline [psu*Mt/day] (',cd_routine,') = ',zsmv 252 IF ( ABS( zei ) > zh_sill ) WRITE(numout,*) 'violation enthalpy [GW] (',cd_routine,') = ',zei 253 IF ( ABS(zvtrp ) > zv_sill .AND. cd_routine == 'limtrp' ) THEN 254 WRITE(numout,*) 'violation vtrp [Mt/day] (',cd_routine,') = ',zvtrp 255 WRITE(numout,*) 'violation etrp [GW] (',cd_routine,') = ',zetrp 213 256 ENDIF 214 IF ( zamin < 0. ) WRITE(numout,*) 'violation a_i<0 (',cd_routine,') = ',zamin 257 IF ( zvmin < -epsi10 ) WRITE(numout,*) 'violation v_i<0 [m] (',cd_routine,') = ',zvmin 258 IF ( zamax > rn_amax+epsi10 .AND. cd_routine /= 'limtrp' .AND. cd_routine /= 'limitd_me' ) THEN 259 WRITE(numout,*) 'violation a_i>amax (',cd_routine,') = ',zamax 260 ENDIF 261 IF ( zamin < -epsi10 ) WRITE(numout,*) 'violation a_i<0 (',cd_routine,') = ',zamin 215 262 ENDIF 216 263 … … 218 265 219 266 END SUBROUTINE lim_cons_hsm 267 268 SUBROUTINE lim_cons_final( cd_routine ) 269 !!--------------------------------------------------------------------------------------------------------- 270 !! *** ROUTINE lim_cons_final *** 271 !! 272 !! ** Purpose : Test the conservation of heat, salt and mass at the end of each ice time-step 273 !! 274 !! ** Method : This is an online diagnostics which can be activated with ln_limdiahsb=true 275 !! It prints in ocean.output if there is a violation of conservation at each time-step 276 !! The thresholds (zv_sill, zs_sill, zh_sill) which determine the violation are set to 277 !! a minimum of 1 mm of ice (over the ice area) that is lost/gained spuriously during 100 years. 278 !! For salt and heat thresholds, ice is considered to have a salinity of 10 279 !! and a heat content of 3e5 J/kg (=latent heat of fusion) 280 !!-------------------------------------------------------------------------------------------------------- 281 CHARACTER(len=*), INTENT(in) :: cd_routine ! name of the routine 282 REAL(wp) :: zhfx, zsfx, zvfx 283 REAL(wp) :: zarea, zv_sill, zs_sill, zh_sill 284 REAL(wp), PARAMETER :: zconv = 1.e-9 ! convert W to GW and kg to Mt 285 286 #if ! defined key_bdy 287 ! heat flux 288 zhfx = glob_sum( ( hfx_in - hfx_out - diag_heat - diag_trp_ei - diag_trp_es - hfx_sub ) * e12t * tmask(:,:,1) * zconv ) 289 ! salt flux 290 zsfx = glob_sum( ( sfx + diag_smvi ) * e12t * tmask(:,:,1) * zconv ) * rday 291 ! water flux 292 zvfx = glob_sum( ( wfx_ice + wfx_snw + wfx_spr + wfx_sub + diag_vice + diag_vsnw ) * e12t * tmask(:,:,1) * zconv ) * rday 293 294 ! set threshold values and calculate the ice area (+epsi10 to set a threshold > 0 when there is no ice) 295 zarea = glob_sum( SUM( a_i + epsi10, dim=3 ) * e12t * zconv ) ! in 1.e9 m2 296 zv_sill = zarea * 2.5e-5 297 zs_sill = zarea * 25.e-5 298 zh_sill = zarea * 10.e-5 299 300 IF( ABS( zvfx ) > zv_sill ) WRITE(numout,*) 'violation vfx [Mt/day] (',cd_routine,') = ',(zvfx) 301 IF( ABS( zsfx ) > zs_sill ) WRITE(numout,*) 'violation sfx [psu*Mt/day] (',cd_routine,') = ',(zsfx) 302 IF( ABS( zhfx ) > zh_sill ) WRITE(numout,*) 'violation hfx [GW] (',cd_routine,') = ',(zhfx) 303 #endif 304 305 END SUBROUTINE lim_cons_final 220 306 221 307 #else -
branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/LIM_SRC_3/limdiahsb.F90
- Property svn:keywords set to Id
r4990 r5682 14 14 !!---------------------------------------------------------------------- 15 15 USE ice ! LIM-3: sea-ice variable 16 USE par_ice ! LIM-3: ice parameters17 16 USE dom_ice ! LIM-3: sea-ice domain 18 17 USE dom_oce ! ocean domain … … 32 31 33 32 PUBLIC lim_diahsb ! routine called by ice_step.F90 34 !!PUBLIC lim_diahsb_init ! routine called by ice_init.F9035 !!PUBLIC lim_diahsb_rst ! routine called by ice_init.F9036 33 37 34 real(wp) :: frc_sal, frc_vol ! global forcing trends … … 43 40 !!---------------------------------------------------------------------- 44 41 !! NEMO/OPA 3.4 , NEMO Consortium (2012) 45 !! $Id : limdiahsb.F90 3294 2012-10-18 16:44:18Z rblod$42 !! $Id$ 46 43 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 47 44 !!---------------------------------------------------------------------- … … 74 71 75 72 ! 1/area 76 z1_area = 1._wp / MAX( glob_sum( area(:,:) * tms(:,:) ), epsi06 )77 78 rswitch = MAX( 0._wp , SIGN( 1._wp , glob_sum( area(:,:) * tms(:,:) ) - epsi06 ) )73 z1_area = 1._wp / MAX( glob_sum( e12t(:,:) * tmask(:,:,1) ), epsi06 ) 74 75 rswitch = MAX( 0._wp , SIGN( 1._wp , glob_sum( e12t(:,:) * tmask(:,:,1) ) - epsi06 ) ) 79 76 ! ----------------------- ! 80 77 ! 1 - Content variations ! 81 78 ! ----------------------- ! 82 zbg_ivo = glob_sum( vt_i(:,:) * area(:,:) * tms(:,:) ) ! volume ice83 zbg_svo = glob_sum( vt_s(:,:) * area(:,:) * tms(:,:) ) ! volume snow84 zbg_are = glob_sum( at_i(:,:) * area(:,:) * tms(:,:) ) ! area85 zbg_sal = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) ! mean salt content86 zbg_tem = glob_sum( ( tm_i(:,:) - rt t ) * vt_i(:,:) * area(:,:) * tms(:,:) ) ! mean temp content87 88 !zbg_ihc = glob_sum( et_i(:,:) * area(:,:) * tms(:,:) ) / MAX( zbg_ivo,epsi06 ) ! ice heat content89 !zbg_shc = glob_sum( et_s(:,:) * area(:,:) * tms(:,:) ) / MAX( zbg_svo,epsi06 ) ! snow heat content79 zbg_ivo = glob_sum( vt_i(:,:) * e12t(:,:) * tmask(:,:,1) ) ! volume ice 80 zbg_svo = glob_sum( vt_s(:,:) * e12t(:,:) * tmask(:,:,1) ) ! volume snow 81 zbg_are = glob_sum( at_i(:,:) * e12t(:,:) * tmask(:,:,1) ) ! area 82 zbg_sal = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * e12t(:,:) * tmask(:,:,1) ) ! mean salt content 83 zbg_tem = glob_sum( ( tm_i(:,:) - rt0 ) * vt_i(:,:) * e12t(:,:) * tmask(:,:,1) ) ! mean temp content 84 85 !zbg_ihc = glob_sum( et_i(:,:) * e12t(:,:) * tmask(:,:,1) ) / MAX( zbg_ivo,epsi06 ) ! ice heat content 86 !zbg_shc = glob_sum( et_s(:,:) * e12t(:,:) * tmask(:,:,1) ) / MAX( zbg_svo,epsi06 ) ! snow heat content 90 87 91 88 ! Volume 92 89 ztmp = rswitch * z1_area * r1_rau0 * rday 93 zbg_vfx = ztmp * glob_sum( emp(:,:) * area(:,:) * tms(:,:) )94 zbg_vfx_bog = ztmp * glob_sum( wfx_bog(:,:) * area(:,:) * tms(:,:) )95 zbg_vfx_opw = ztmp * glob_sum( wfx_opw(:,:) * area(:,:) * tms(:,:) )96 zbg_vfx_sni = ztmp * glob_sum( wfx_sni(:,:) * area(:,:) * tms(:,:) )97 zbg_vfx_dyn = ztmp * glob_sum( wfx_dyn(:,:) * area(:,:) * tms(:,:) )98 zbg_vfx_bom = ztmp * glob_sum( wfx_bom(:,:) * area(:,:) * tms(:,:) )99 zbg_vfx_sum = ztmp * glob_sum( wfx_sum(:,:) * area(:,:) * tms(:,:) )100 zbg_vfx_res = ztmp * glob_sum( wfx_res(:,:) * area(:,:) * tms(:,:) )101 zbg_vfx_spr = ztmp * glob_sum( wfx_spr(:,:) * area(:,:) * tms(:,:) )102 zbg_vfx_snw = ztmp * glob_sum( wfx_snw(:,:) * area(:,:) * tms(:,:) )103 zbg_vfx_sub = ztmp * glob_sum( wfx_sub(:,:) * area(:,:) * tms(:,:) )90 zbg_vfx = ztmp * glob_sum( emp(:,:) * e12t(:,:) * tmask(:,:,1) ) 91 zbg_vfx_bog = ztmp * glob_sum( wfx_bog(:,:) * e12t(:,:) * tmask(:,:,1) ) 92 zbg_vfx_opw = ztmp * glob_sum( wfx_opw(:,:) * e12t(:,:) * tmask(:,:,1) ) 93 zbg_vfx_sni = ztmp * glob_sum( wfx_sni(:,:) * e12t(:,:) * tmask(:,:,1) ) 94 zbg_vfx_dyn = ztmp * glob_sum( wfx_dyn(:,:) * e12t(:,:) * tmask(:,:,1) ) 95 zbg_vfx_bom = ztmp * glob_sum( wfx_bom(:,:) * e12t(:,:) * tmask(:,:,1) ) 96 zbg_vfx_sum = ztmp * glob_sum( wfx_sum(:,:) * e12t(:,:) * tmask(:,:,1) ) 97 zbg_vfx_res = ztmp * glob_sum( wfx_res(:,:) * e12t(:,:) * tmask(:,:,1) ) 98 zbg_vfx_spr = ztmp * glob_sum( wfx_spr(:,:) * e12t(:,:) * tmask(:,:,1) ) 99 zbg_vfx_snw = ztmp * glob_sum( wfx_snw(:,:) * e12t(:,:) * tmask(:,:,1) ) 100 zbg_vfx_sub = ztmp * glob_sum( wfx_sub(:,:) * e12t(:,:) * tmask(:,:,1) ) 104 101 105 102 ! Salt 106 zbg_sfx = ztmp * glob_sum( sfx(:,:) * area(:,:) * tms(:,:) )107 zbg_sfx_bri = ztmp * glob_sum( sfx_bri(:,:) * area(:,:) * tms(:,:) )108 zbg_sfx_res = ztmp * glob_sum( sfx_res(:,:) * area(:,:) * tms(:,:) )109 zbg_sfx_dyn = ztmp * glob_sum( sfx_dyn(:,:) * area(:,:) * tms(:,:) )110 111 zbg_sfx_bog = ztmp * glob_sum( sfx_bog(:,:) * area(:,:) * tms(:,:) )112 zbg_sfx_opw = ztmp * glob_sum( sfx_opw(:,:) * area(:,:) * tms(:,:) )113 zbg_sfx_sni = ztmp * glob_sum( sfx_sni(:,:) * area(:,:) * tms(:,:) )114 zbg_sfx_bom = ztmp * glob_sum( sfx_bom(:,:) * area(:,:) * tms(:,:) )115 zbg_sfx_sum = ztmp * glob_sum( sfx_sum(:,:) * area(:,:) * tms(:,:) )103 zbg_sfx = ztmp * glob_sum( sfx(:,:) * e12t(:,:) * tmask(:,:,1) ) 104 zbg_sfx_bri = ztmp * glob_sum( sfx_bri(:,:) * e12t(:,:) * tmask(:,:,1) ) 105 zbg_sfx_res = ztmp * glob_sum( sfx_res(:,:) * e12t(:,:) * tmask(:,:,1) ) 106 zbg_sfx_dyn = ztmp * glob_sum( sfx_dyn(:,:) * e12t(:,:) * tmask(:,:,1) ) 107 108 zbg_sfx_bog = ztmp * glob_sum( sfx_bog(:,:) * e12t(:,:) * tmask(:,:,1) ) 109 zbg_sfx_opw = ztmp * glob_sum( sfx_opw(:,:) * e12t(:,:) * tmask(:,:,1) ) 110 zbg_sfx_sni = ztmp * glob_sum( sfx_sni(:,:) * e12t(:,:) * tmask(:,:,1) ) 111 zbg_sfx_bom = ztmp * glob_sum( sfx_bom(:,:) * e12t(:,:) * tmask(:,:,1) ) 112 zbg_sfx_sum = ztmp * glob_sum( sfx_sum(:,:) * e12t(:,:) * tmask(:,:,1) ) 116 113 117 114 ! Heat budget 118 zbg_ihc = glob_sum( et_i(:,:) * 1.e-20 ) ! ice heat content [1.e-20 J]119 zbg_shc = glob_sum( et_s(:,:) * 1.e-20 ) ! snow heat content [1.e-20 J]120 zbg_hfx_dhc = glob_sum( diag_heat _dhc(:,:) * area(:,:) * tms(:,:) ) ! [in W]121 zbg_hfx_spr = glob_sum( hfx_spr(:,:) * area(:,:) * tms(:,:) ) ! [in W]122 123 zbg_hfx_thd = glob_sum( hfx_thd(:,:) * area(:,:) * tms(:,:) ) ! [in W]124 zbg_hfx_dyn = glob_sum( hfx_dyn(:,:) * area(:,:) * tms(:,:) ) ! [in W]125 zbg_hfx_res = glob_sum( hfx_res(:,:) * area(:,:) * tms(:,:) ) ! [in W]126 zbg_hfx_sub = glob_sum( hfx_sub(:,:) * area(:,:) * tms(:,:) ) ! [in W]127 zbg_hfx_snw = glob_sum( hfx_snw(:,:) * area(:,:) * tms(:,:) ) ! [in W]128 zbg_hfx_sum = glob_sum( hfx_sum(:,:) * area(:,:) * tms(:,:) ) ! [in W]129 zbg_hfx_bom = glob_sum( hfx_bom(:,:) * area(:,:) * tms(:,:) ) ! [in W]130 zbg_hfx_bog = glob_sum( hfx_bog(:,:) * area(:,:) * tms(:,:) ) ! [in W]131 zbg_hfx_dif = glob_sum( hfx_dif(:,:) * area(:,:) * tms(:,:) ) ! [in W]132 zbg_hfx_opw = glob_sum( hfx_opw(:,:) * area(:,:) * tms(:,:) ) ! [in W]133 zbg_hfx_out = glob_sum( hfx_out(:,:) * area(:,:) * tms(:,:) ) ! [in W]134 zbg_hfx_in = glob_sum( hfx_in(:,:) * area(:,:) * tms(:,:) ) ! [in W]115 zbg_ihc = glob_sum( et_i(:,:) * e12t(:,:) * 1.e-20 ) ! ice heat content [1.e20 J] 116 zbg_shc = glob_sum( et_s(:,:) * e12t(:,:) * 1.e-20 ) ! snow heat content [1.e20 J] 117 zbg_hfx_dhc = glob_sum( diag_heat(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W] 118 zbg_hfx_spr = glob_sum( hfx_spr(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W] 119 120 zbg_hfx_thd = glob_sum( hfx_thd(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W] 121 zbg_hfx_dyn = glob_sum( hfx_dyn(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W] 122 zbg_hfx_res = glob_sum( hfx_res(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W] 123 zbg_hfx_sub = glob_sum( hfx_sub(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W] 124 zbg_hfx_snw = glob_sum( hfx_snw(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W] 125 zbg_hfx_sum = glob_sum( hfx_sum(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W] 126 zbg_hfx_bom = glob_sum( hfx_bom(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W] 127 zbg_hfx_bog = glob_sum( hfx_bog(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W] 128 zbg_hfx_dif = glob_sum( hfx_dif(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W] 129 zbg_hfx_opw = glob_sum( hfx_opw(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W] 130 zbg_hfx_out = glob_sum( hfx_out(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W] 131 zbg_hfx_in = glob_sum( hfx_in(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W] 135 132 136 133 ! --------------------------------------------- ! 137 134 ! 2 - Trends due to forcing and ice growth/melt ! 138 135 ! --------------------------------------------- ! 139 z_frc_vol = r1_rau0 * glob_sum( - emp(:,:) * area(:,:) * tms(:,:) ) ! volume fluxes140 z_frc_sal = r1_rau0 * glob_sum( sfx(:,:) * area(:,:) * tms(:,:) ) ! salt fluxes136 z_frc_vol = r1_rau0 * glob_sum( - emp(:,:) * e12t(:,:) * tmask(:,:,1) ) ! volume fluxes 137 z_frc_sal = r1_rau0 * glob_sum( sfx(:,:) * e12t(:,:) * tmask(:,:,1) ) ! salt fluxes 141 138 z_bg_grme = glob_sum( - ( wfx_bog(:,:) + wfx_opw(:,:) + wfx_sni(:,:) + wfx_dyn(:,:) + & 142 & wfx_bom(:,:) + wfx_sum(:,:) + wfx_res(:,:) + wfx_snw(:,:) + wfx_sub(:,:) ) * area(:,:) * tms(:,:) ) ! volume fluxes 139 & wfx_bom(:,:) + wfx_sum(:,:) + wfx_res(:,:) + wfx_snw(:,:) + & 140 & wfx_sub(:,:) ) * e12t(:,:) * tmask(:,:,1) ) ! volume fluxes 143 141 ! 144 142 frc_vol = frc_vol + z_frc_vol * rdt_ice … … 247 245 WRITE(numout,*) '~~~~~~~~~~~~' 248 246 ENDIF 249 250 ! ---------------------------------- !251 ! 2 - initial conservation variables !252 ! ---------------------------------- !253 !frc_vol = 0._wp ! volume trend due to forcing254 !frc_sal = 0._wp ! salt content - - - -255 !bg_grme = 0._wp ! ice growth + melt volume trend256 247 ! 257 248 CALL lim_diahsb_rst( nstart, 'READ' ) !* read or initialize all required files -
branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90
r4990 r5682 6 6 !! history : 1.0 ! 2002-08 (C. Ethe, G. Madec) original VP code 7 7 !! 3.0 ! 2007-03 (MA Morales Maqueda, S. Bouillon, M. Vancoppenolle) LIM3: EVP-Cgrid 8 !! 4.0! 2011-02 (G. Madec) dynamical allocation8 !! 3.5 ! 2011-02 (G. Madec) dynamical allocation 9 9 !!---------------------------------------------------------------------- 10 10 #if defined key_lim3 … … 20 20 USE sbc_ice ! Surface boundary condition: ice fields 21 21 USE ice ! LIM-3 variables 22 USE par_ice ! LIM-3 parameters23 22 USE dom_ice ! LIM-3 domain 24 23 USE limrhg ! LIM-3 rheology … … 31 30 USE timing ! Timing 32 31 USE limcons ! conservation tests 32 USE limvar 33 33 34 34 IMPLICIT NONE … … 76 76 CALL wrk_alloc( jpj, zswitch, zmsk ) 77 77 78 CALL lim_var_agg(1) ! aggregate ice categories 79 78 80 IF( kt == nit000 ) CALL lim_dyn_init ! Initialization (first time-step only) 79 81 … … 83 85 IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limdyn', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 84 86 85 u_ice_b(:,:) = u_ice(:,:) * tmu(:,:)86 v_ice_b(:,:) = v_ice(:,:) * tmv(:,:)87 u_ice_b(:,:) = u_ice(:,:) * umask(:,:,1) 88 v_ice_b(:,:) = v_ice(:,:) * vmask(:,:,1) 87 89 88 90 ! Rheology (ice dynamics) … … 101 103 DO jj = 1, jpj 102 104 zswitch(jj) = SUM( 1.0 - at_i(:,jj) ) ! = REAL(jpj) if ocean everywhere on a j-line 103 zmsk (jj) = SUM( tmask(:,jj,1) ) ! = 0 if land everywhere on a j-line105 zmsk (jj) = SUM( tmask(:,jj,1) ) ! = 0 if land everywhere on a j-line 104 106 END DO 105 107 … … 157 159 zv_io(:,:) = v_ice(:,:) - ssv_m(:,:) 158 160 ! frictional velocity at T-point 159 zcoef = 0.5_wp * cw161 zcoef = 0.5_wp * rn_cio 160 162 DO jj = 2, jpjm1 161 163 DO ji = fs_2, fs_jpim1 ! vector opt. 162 164 ust2s(ji,jj) = zcoef * ( zu_io(ji,jj) * zu_io(ji,jj) + zu_io(ji-1,jj) * zu_io(ji-1,jj) & 163 & + zv_io(ji,jj) * zv_io(ji,jj) + zv_io(ji,jj-1) * zv_io(ji,jj-1) ) * tms(ji,jj)165 & + zv_io(ji,jj) * zv_io(ji,jj) + zv_io(ji,jj-1) * zv_io(ji,jj-1) ) * tmask(ji,jj,1) 164 166 END DO 165 167 END DO … … 170 172 ELSE ! no ice dynamics : transmit directly the atmospheric stress to the ocean 171 173 ! 172 zcoef = SQRT( 0.5_wp ) /rau0174 zcoef = SQRT( 0.5_wp ) * r1_rau0 173 175 DO jj = 2, jpjm1 174 176 DO ji = fs_2, fs_jpim1 ! vector opt. 175 177 ust2s(ji,jj) = zcoef * SQRT( utau(ji,jj) * utau(ji,jj) + utau(ji-1,jj) * utau(ji-1,jj) & 176 & + vtau(ji,jj) * vtau(ji,jj) + vtau(ji,jj-1) * vtau(ji,jj-1) ) * tms(ji,jj)178 & + vtau(ji,jj) * vtau(ji,jj) + vtau(ji,jj-1) * vtau(ji,jj-1) ) * tmask(ji,jj,1) 177 179 END DO 178 180 END DO … … 189 191 CALL prt_ctl(tab2d_1=delta_i , clinfo1=' lim_dyn : delta_i :') 190 192 CALL prt_ctl(tab2d_1=strength , clinfo1=' lim_dyn : strength :') 191 CALL prt_ctl(tab2d_1= area, clinfo1=' lim_dyn : cell area :')193 CALL prt_ctl(tab2d_1=e12t , clinfo1=' lim_dyn : cell area :') 192 194 CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_dyn : at_i :') 193 195 CALL prt_ctl(tab2d_1=vt_i , clinfo1=' lim_dyn : vt_i :') … … 241 243 !!------------------------------------------------------------------- 242 244 INTEGER :: ios ! Local integer output status for namelist read 243 NAMELIST/namicedyn/ epsd, om, cw, pstar, & 244 & c_rhg, creepl, ecc, ahi0, & 245 & nevp, relast, alphaevp, hminrhg 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 246 249 !!------------------------------------------------------------------- 247 250 … … 259 262 WRITE(numout,*) 'lim_dyn_init : ice parameters for ice dynamics ' 260 263 WRITE(numout,*) '~~~~~~~~~~~~' 261 WRITE(numout,*) ' tolerance parameter epsd = ', epsd262 WRITE(numout,*) ' relaxation constant om = ', om263 WRITE(numout,*) ' drag coefficient for oceanic stress cw = ', cw264 WRITE(numout,*) ' first bulk-rheology parameter pstar = ', pstar265 WRITE(numout,*) ' second bulk-rhelogy parameter c_rhg = ', c_rhg266 WRITE(numout,*) ' creep limit creepl = ', creepl267 WRITE(numout,*) ' eccentricity of the elliptical yield curve ecc = ', ecc268 WRITE(numout,*) ' horizontal diffusivity coeff. for sea-ice ahi0 = ', ahi0269 WRITE(numout,*) ' number of iterations for subcycling nevp = ',nevp270 WRITE(numout,*) ' ratio of elastic timescale over ice time step relast = ',relast271 WRITE(numout,*) ' coefficient for the solution of int. stresses alphaevp = ', alphaevp272 WRITE(numout,*) ' min ice thickness for rheology calculations hminrhg = ', hminrhg264 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 273 276 ENDIF 274 277 ! 275 usecc2 = 1._wp / ( ecc * ecc ) 276 rhoco = rau0 * cw 277 278 ! elastic damping 279 telast = relast * rdt_ice 280 281 ! Diffusion coefficients. 282 ahiu(:,:) = ahi0 * umask(:,:,1) 283 ahiv(:,:) = ahi0 * vmask(:,:,1) 284 ! 278 usecc2 = 1._wp / ( rn_ecc * rn_ecc ) 279 rhoco = rau0 * rn_cio 280 ! 281 ! Diffusion coefficients 282 SELECT CASE( nn_ahi0 ) 283 284 CASE( 0 ) 285 ahiu(:,:) = rn_ahi0_ref 286 ahiv(:,:) = rn_ahi0_ref 287 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 domain 295 296 ahiu(:,:) = rn_ahi0_ref * zd_max * 1.e-05_wp ! 1.e05 = 100km = max grid space at 60° latitude in orca2 297 ! (60° = min latitude for ice cover) 298 ahiv(:,:) = rn_ahi0_ref * zd_max * 1.e-05_wp 299 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_wp 303 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 domain 308 309 za00 = rn_ahi0_ref * 1.e-05_wp ! 1.e05 = 100km = max grid space at 60° latitude in orca2 310 ! (60° = min latitude for ice cover) 311 DO jj = 1, jpj 312 DO ji = 1, jpi 313 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 DO 316 END DO 317 ! 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_max 321 322 END SELECT 323 285 324 END SUBROUTINE lim_dyn_init 286 325 -
branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/LIM_SRC_3/limhdf.F90
r4990 r5682 13 13 !!---------------------------------------------------------------------- 14 14 !! lim_hdf : diffusion trend on sea-ice variable 15 !! lim_hdf_init : initialisation of diffusion trend on sea-ice variable 15 16 !!---------------------------------------------------------------------- 16 17 USE dom_oce ! ocean domain … … 26 27 PRIVATE 27 28 28 PUBLIC lim_hdf ! called by lim_tra 29 30 LOGICAL :: linit = .TRUE. ! initialization flag (set to flase after the 1st call) 31 REAL(wp) :: epsi04 = 1.e-04 ! constant 29 PUBLIC lim_hdf ! called by lim_trp 30 PUBLIC lim_hdf_init ! called by sbc_lim_init 31 32 LOGICAL :: linit = .TRUE. ! initialization flag (set to flase after the 1st call) 33 INTEGER :: nn_convfrq !: convergence check frequency of the Crant-Nicholson scheme 32 34 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: efact ! metric coefficient 33 35 … … 54 56 REAL(wp), DIMENSION(jpi,jpj), INTENT( inout ) :: ptab ! Field on which the diffusion is applied 55 57 ! 56 INTEGER :: ji, jj ! dummy loop indices 57 INTEGER :: its, iter, ierr ! local integers 58 REAL(wp) :: zalfa, zrlxint, zconv ! local scalars 59 REAL(wp), POINTER, DIMENSION(:,:) :: zrlx, zflu, zflv, zdiv0, zdiv, ztab0 60 CHARACTER(lc) :: charout ! local character 58 INTEGER :: ji, jj ! dummy loop indices 59 INTEGER :: iter, ierr ! local integers 60 REAL(wp) :: zrlxint, zconv ! local scalars 61 REAL(wp), POINTER, DIMENSION(:,:) :: zrlx, zflu, zflv, zdiv0, zdiv, ztab0 62 CHARACTER(lc) :: charout ! local character 63 REAL(wp), PARAMETER :: zrelax = 0.5_wp ! relaxation constant for iterative procedure 64 REAL(wp), PARAMETER :: zalfa = 0.5_wp ! =1.0/0.5/0.0 = implicit/Cranck-Nicholson/explicit 65 INTEGER , PARAMETER :: its = 100 ! Maximum number of iteration 61 66 !!------------------------------------------------------------------- 62 67 … … 71 76 DO jj = 2, jpjm1 72 77 DO ji = fs_2 , fs_jpim1 ! vector opt. 73 efact(ji,jj) = ( e2u(ji,jj) + e2u(ji-1,jj) + e1v(ji,jj) + e1v(ji,jj-1) ) / ( e1t(ji,jj) * e2t(ji,jj))78 efact(ji,jj) = ( e2u(ji,jj) + e2u(ji-1,jj) + e1v(ji,jj) + e1v(ji,jj-1) ) * r1_e12t(ji,jj) 74 79 END DO 75 80 END DO … … 77 82 ENDIF 78 83 ! ! Time integration parameters 79 zalfa = 0.5_wp ! =1.0/0.5/0.0 = implicit/Cranck-Nicholson/explicit80 its = 100 ! Maximum number of iteration81 84 ! 82 85 ztab0(:, : ) = ptab(:,:) ! Arrays initialization … … 91 94 iter = 0 92 95 ! 93 DO WHILE( zconv > ( 2._wp * epsi04 ) .AND. iter <= its ) ! Sub-time step loop96 DO WHILE( zconv > ( 2._wp * 1.e-04 ) .AND. iter <= its ) ! Sub-time step loop 94 97 ! 95 98 iter = iter + 1 ! incrementation of the sub-time step number … … 97 100 DO jj = 1, jpjm1 ! diffusive fluxes in U- and V- direction 98 101 DO ji = 1 , fs_jpim1 ! vector opt. 99 zflu(ji,jj) = pahu(ji,jj) * e2u(ji,jj) /e1u(ji,jj) * ( ptab(ji+1,jj) - ptab(ji,jj) )100 zflv(ji,jj) = pahv(ji,jj) * e1v(ji,jj) /e2v(ji,jj) * ( ptab(ji,jj+1) - ptab(ji,jj) )102 zflu(ji,jj) = pahu(ji,jj) * e2u(ji,jj) * r1_e1u(ji,jj) * ( ptab(ji+1,jj) - ptab(ji,jj) ) 103 zflv(ji,jj) = pahv(ji,jj) * e1v(ji,jj) * r1_e2v(ji,jj) * ( ptab(ji,jj+1) - ptab(ji,jj) ) 101 104 END DO 102 105 END DO … … 104 107 DO jj= 2, jpjm1 ! diffusive trend : divergence of the fluxes 105 108 DO ji = fs_2 , fs_jpim1 ! vector opt. 106 zdiv (ji,jj) = ( zflu(ji,jj) - zflu(ji-1,jj ) & 107 & + zflv(ji,jj) - zflv(ji ,jj-1) ) / ( e1t (ji,jj) * e2t (ji,jj) ) 109 zdiv(ji,jj) = ( zflu(ji,jj) - zflu(ji-1,jj) + zflv(ji,jj) - zflv(ji,jj-1) ) * r1_e12t(ji,jj) 108 110 END DO 109 111 END DO … … 115 117 zrlxint = ( ztab0(ji,jj) & 116 118 & + rdt_ice * ( zalfa * ( zdiv(ji,jj) + efact(ji,jj) * ptab(ji,jj) ) & 117 & + ( 1.0 - zalfa ) * zdiv0(ji,jj) ) )&118 & / ( 1.0 + zalfa * rdt_ice * efact(ji,jj) )119 zrlx(ji,jj) = ptab(ji,jj) + om* ( zrlxint - ptab(ji,jj) )119 & + ( 1.0 - zalfa ) * zdiv0(ji,jj) ) & 120 & ) / ( 1.0 + zalfa * rdt_ice * efact(ji,jj) ) 121 zrlx(ji,jj) = ptab(ji,jj) + zrelax * ( zrlxint - ptab(ji,jj) ) 120 122 END DO 121 123 END DO 122 124 CALL lbc_lnk( zrlx, 'T', 1. ) ! lateral boundary condition 123 125 ! 124 zconv = 0._wp ! convergence test 125 DO jj = 2, jpjm1 126 DO ji = fs_2, fs_jpim1 127 zconv = MAX( zconv, ABS( zrlx(ji,jj) - ptab(ji,jj) ) ) 128 END DO 129 END DO 130 IF( lk_mpp ) CALL mpp_max( zconv ) ! max over the global domain 126 IF ( MOD( iter, nn_convfrq ) == 0 ) THEN ! convergence test every nn_convfrq iterations (perf. optimization) 127 zconv = 0._wp 128 DO jj = 2, jpjm1 129 DO ji = fs_2, fs_jpim1 130 zconv = MAX( zconv, ABS( zrlx(ji,jj) - ptab(ji,jj) ) ) 131 END DO 132 END DO 133 IF( lk_mpp ) CALL mpp_max( zconv ) ! max over the global domain 134 ENDIF 131 135 ! 132 136 ptab(:,:) = zrlx(:,:) … … 138 142 DO jj = 1, jpjm1 ! diffusive fluxes in U- and V- direction 139 143 DO ji = 1 , fs_jpim1 ! vector opt. 140 zflu(ji,jj) = pahu(ji,jj) * e2u(ji,jj) /e1u(ji,jj) * ( ptab(ji+1,jj) - ptab(ji,jj) )141 zflv(ji,jj) = pahv(ji,jj) * e1v(ji,jj) /e2v(ji,jj) * ( ptab(ji,jj+1) - ptab(ji,jj) )144 zflu(ji,jj) = pahu(ji,jj) * e2u(ji,jj) * r1_e1u(ji,jj) * ( ptab(ji+1,jj) - ptab(ji,jj) ) 145 zflv(ji,jj) = pahv(ji,jj) * e1v(ji,jj) * r1_e2v(ji,jj) * ( ptab(ji,jj+1) - ptab(ji,jj) ) 142 146 END DO 143 147 END DO … … 145 149 DO jj= 2, jpjm1 ! diffusive trend : divergence of the fluxes 146 150 DO ji = fs_2 , fs_jpim1 ! vector opt. 147 zdiv (ji,jj) = ( zflu(ji,jj) - zflu(ji-1,jj ) & 148 & + zflv(ji,jj) - zflv(ji ,jj-1) ) / ( e1t (ji,jj) * e2t (ji,jj) ) 151 zdiv(ji,jj) = ( zflu(ji,jj) - zflu(ji-1,jj) + zflv(ji,jj) - zflv(ji,jj-1) ) * r1_e12t(ji,jj) 149 152 ptab(ji,jj) = ztab0(ji,jj) + 0.5 * ( zdiv(ji,jj) + zdiv0(ji,jj) ) 150 153 END DO … … 164 167 END SUBROUTINE lim_hdf 165 168 169 170 SUBROUTINE lim_hdf_init 171 !!------------------------------------------------------------------- 172 !! *** ROUTINE lim_hdf_init *** 173 !! 174 !! ** Purpose : Initialisation of horizontal diffusion of sea-ice 175 !! 176 !! ** Method : Read the namicehdf namelist 177 !! 178 !! ** input : Namelist namicehdf 179 !!------------------------------------------------------------------- 180 INTEGER :: ios ! Local integer output status for namelist read 181 NAMELIST/namicehdf/ nn_convfrq 182 !!------------------------------------------------------------------- 183 ! 184 IF(lwp) THEN 185 WRITE(numout,*) 186 WRITE(numout,*) 'lim_hdf : Ice horizontal diffusion' 187 WRITE(numout,*) '~~~~~~~' 188 ENDIF 189 ! 190 REWIND( numnam_ice_ref ) ! Namelist namicehdf in reference namelist : Ice horizontal diffusion 191 READ ( numnam_ice_ref, namicehdf, IOSTAT = ios, ERR = 901) 192 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicehdf in reference namelist', lwp ) 193 194 REWIND( numnam_ice_cfg ) ! Namelist namicehdf in configuration namelist : Ice horizontal diffusion 195 READ ( numnam_ice_cfg, namicehdf, IOSTAT = ios, ERR = 902 ) 196 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicehdf in configuration namelist', lwp ) 197 IF(lwm) WRITE ( numoni, namicehdf ) 198 ! 199 IF(lwp) THEN ! control print 200 WRITE(numout,*) 201 WRITE(numout,*)' Namelist of ice parameters for ice horizontal diffusion computation ' 202 WRITE(numout,*)' convergence check frequency of the Crant-Nicholson scheme nn_convfrq = ', nn_convfrq 203 ENDIF 204 ! 205 END SUBROUTINE lim_hdf_init 166 206 #else 167 207 !!---------------------------------------------------------------------- 168 208 !! Default option Dummy module NO LIM sea-ice model 169 209 !!---------------------------------------------------------------------- 170 CONTAINS171 SUBROUTINE lim_hdf ! Empty routine172 END SUBROUTINE lim_hdf173 210 #endif 174 211 -
branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90
r4990 r5682 22 22 USE eosbn2 ! equation of state 23 23 USE ice ! sea-ice variables 24 USE par_ice ! ice parameters25 24 USE par_oce ! ocean parameters 26 25 USE dom_ice ! sea-ice domain … … 36 35 37 36 ! !!** init namelist (namiceini) ** 38 REAL(wp) :: thres_sst ! threshold water temperature for initial sea ice39 REAL(wp) :: hts_ini_n ! initial snow thickness in the north40 REAL(wp) :: hts_ini_s ! initial snow thickness in the south41 REAL(wp) :: hti_ini_n ! initial ice thickness in the north42 REAL(wp) :: hti_ini_s ! initial ice thickness in the south43 REAL(wp) :: ati_ini_n ! initial leads area in the north44 REAL(wp) :: ati_ini_s ! initial leads area in the south45 REAL(wp) :: smi_ini_n ! initial salinity46 REAL(wp) :: smi_ini_s ! initial salinity47 REAL(wp) :: tmi_ini_n ! initial temperature48 REAL(wp) :: tmi_ini_s ! initial temperature49 50 LOGICAL :: ln_ limini ! initialization or not37 REAL(wp) :: rn_thres_sst ! threshold water temperature for initial sea ice 38 REAL(wp) :: rn_hts_ini_n ! initial snow thickness in the north 39 REAL(wp) :: rn_hts_ini_s ! initial snow thickness in the south 40 REAL(wp) :: rn_hti_ini_n ! initial ice thickness in the north 41 REAL(wp) :: rn_hti_ini_s ! initial ice thickness in the south 42 REAL(wp) :: rn_ati_ini_n ! initial leads area in the north 43 REAL(wp) :: rn_ati_ini_s ! initial leads area in the south 44 REAL(wp) :: rn_smi_ini_n ! initial salinity 45 REAL(wp) :: rn_smi_ini_s ! initial salinity 46 REAL(wp) :: rn_tmi_ini_n ! initial temperature 47 REAL(wp) :: rn_tmi_ini_s ! initial temperature 48 49 LOGICAL :: ln_iceini ! initialization or not 51 50 !!---------------------------------------------------------------------- 52 51 !! LIM 3.0, UCL-LOCEAN-IPSL (2008) … … 87 86 !! * Local variables 88 87 INTEGER :: ji, jj, jk, jl ! dummy loop indices 89 REAL(wp) :: epsi20,ztmelts, zdh88 REAL(wp) :: ztmelts, zdh 90 89 INTEGER :: i_hemis, i_fill, jl0 91 90 REAL(wp) :: ztest_1, ztest_2, ztest_3, ztest_4, ztests, zsigma, zarg, zA, zV, zA_cons, zV_cons, zconv … … 101 100 CALL wrk_alloc( 2, zht_i_ini, zat_i_ini, zvt_i_ini, zht_s_ini, zsm_i_ini, ztm_i_ini ) 102 101 103 epsi20 = 1.e-20_wp104 105 102 IF(lwp) WRITE(numout,*) 106 103 IF(lwp) WRITE(numout,*) 'lim_istate : Ice initialization ' … … 115 112 ! surface temperature 116 113 DO jl = 1, jpl ! loop over categories 117 t_su (:,:,jl) = rt t * tms(:,:)118 tn_ice(:,:,jl) = rt t * tms(:,:)114 t_su (:,:,jl) = rt0 * tmask(:,:,1) 115 tn_ice(:,:,jl) = rt0 * tmask(:,:,1) 119 116 END DO 120 117 121 118 ! basal temperature (considered at freezing point) 122 t_bo(:,:) = ( eos_fzp( tsn(:,:,1,jp_sal) ) + rt0 ) * tms(:,:) 123 124 IF( ln_limini ) THEN 119 CALL eos_fzp( sss_m(:,:), t_bo(:,:) ) 120 t_bo(:,:) = ( t_bo(:,:) + rt0 ) * tmask(:,:,1) 121 122 IF( ln_iceini ) THEN 125 123 126 124 !-------------------------------------------------------------------- … … 130 128 DO jj = 1, jpj ! ice if sst <= t-freez + ttest 131 129 DO ji = 1, jpi 132 IF( ( tsn(ji,jj,1,jp_tem) - ( t_bo(ji,jj) - rt0 ) ) * tms(ji,jj) >=thres_sst ) THEN133 zswitch(ji,jj) = 0._wp * tm s(ji,jj) ! no ice130 IF( ( sst_m(ji,jj) - ( t_bo(ji,jj) - rt0 ) ) * tmask(ji,jj,1) >= rn_thres_sst ) THEN 131 zswitch(ji,jj) = 0._wp * tmask(ji,jj,1) ! no ice 134 132 ELSE 135 zswitch(ji,jj) = 1._wp * tm s(ji,jj) ! ice133 zswitch(ji,jj) = 1._wp * tmask(ji,jj,1) ! ice 136 134 ENDIF 137 135 END DO … … 158 156 !----------------------------- 159 157 ! assign initial thickness, concentration, snow depth and salinity to an hemisphere-dependent array 160 zht_i_ini(1) = hti_ini_n ; zht_i_ini(2) =hti_ini_s ! ice thickness161 zht_s_ini(1) = hts_ini_n ; zht_s_ini(2) =hts_ini_s ! snow depth162 zat_i_ini(1) = ati_ini_n ; zat_i_ini(2) =ati_ini_s ! ice concentration163 zsm_i_ini(1) = smi_ini_n ; zsm_i_ini(2) =smi_ini_s ! bulk ice salinity164 ztm_i_ini(1) = tmi_ini_n ; ztm_i_ini(2) =tmi_ini_s ! temperature (ice and snow)158 zht_i_ini(1) = rn_hti_ini_n ; zht_i_ini(2) = rn_hti_ini_s ! ice thickness 159 zht_s_ini(1) = rn_hts_ini_n ; zht_s_ini(2) = rn_hts_ini_s ! snow depth 160 zat_i_ini(1) = rn_ati_ini_n ; zat_i_ini(2) = rn_ati_ini_s ! ice concentration 161 zsm_i_ini(1) = rn_smi_ini_n ; zsm_i_ini(2) = rn_smi_ini_s ! bulk ice salinity 162 ztm_i_ini(1) = rn_tmi_ini_n ; ztm_i_ini(2) = rn_tmi_ini_s ! temperature (ice and snow) 165 163 166 164 zvt_i_ini(:) = zht_i_ini(:) * zat_i_ini(:) ! ice volume … … 197 195 !--- Ice thicknesses in the i_fill - 1 first categories 198 196 DO jl = 1, i_fill - 1 199 zh_i_ini(jl,i_hemis) = 0.5 * ( hi_max(jl) + hi_max(jl-1))197 zh_i_ini(jl,i_hemis) = hi_mean(jl) 200 198 END DO 201 199 202 200 !--- jl0: most likely index where cc will be maximum 203 201 DO jl = 1, jpl 204 IF ( ( zht_i_ini(i_hemis) .GT.hi_max(jl-1) ) .AND. &205 ( zht_i_ini(i_hemis) .LE.hi_max(jl) ) ) THEN202 IF ( ( zht_i_ini(i_hemis) > hi_max(jl-1) ) .AND. & 203 & ( zht_i_ini(i_hemis) <= hi_max(jl) ) ) THEN 206 204 jl0 = jl 207 205 ENDIF … … 267 265 268 266 ! Test 3: thickness of the last category is in-bounds ? 269 IF ( zh_i_ini(i_fill, i_hemis) .GT.hi_max(i_fill-1) ) THEN267 IF ( zh_i_ini(i_fill, i_hemis) > hi_max(i_fill-1) ) THEN 270 268 ztest_3 = 1 271 269 ELSE … … 317 315 DO ji = 1, jpi 318 316 a_i(ji,jj,jl) = zswitch(ji,jj) * za_i_ini (jl,zhemis(ji,jj)) ! concentration 319 ht_i(ji,jj,jl) = zswitch(ji,jj) * zh_i_ini(jl,zhemis(ji,jj)) ! ice thickness317 ht_i(ji,jj,jl) = zswitch(ji,jj) * zh_i_ini(jl,zhemis(ji,jj)) ! ice thickness 320 318 ht_s(ji,jj,jl) = ht_i(ji,jj,jl) * ( zht_s_ini( zhemis(ji,jj) ) / zht_i_ini( zhemis(ji,jj) ) ) ! snow depth 321 sm_i(ji,jj,jl) = zswitch(ji,jj) * zsm_i_ini(zhemis(ji,jj)) !+ ( 1._wp - zswitch(ji,jj) ) * s_i_min! salinity322 o_i(ji,jj,jl) = zswitch(ji,jj) * 1._wp + ( 1._wp - zswitch(ji,jj) ) ! age323 t_su(ji,jj,jl) = zswitch(ji,jj) * ztm_i_ini(zhemis(ji,jj)) + ( 1._wp - zswitch(ji,jj) ) * rt t! surf temp319 sm_i(ji,jj,jl) = zswitch(ji,jj) * zsm_i_ini(zhemis(ji,jj)) ! salinity 320 o_i(ji,jj,jl) = zswitch(ji,jj) * 1._wp ! age (1 day) 321 t_su(ji,jj,jl) = zswitch(ji,jj) * ztm_i_ini(zhemis(ji,jj)) + ( 1._wp - zswitch(ji,jj) ) * rt0 ! surf temp 324 322 325 323 ! This case below should not be used if (ht_s/ht_i) is ok in namelist … … 329 327 ! recompute ht_i, ht_s avoiding out of bounds values 330 328 ht_i(ji,jj,jl) = MIN( hi_max(jl), ht_i(ji,jj,jl) + zdh ) 331 ht_s(ji,jj,jl) = MAX( 0._wp, ht_s(ji,jj,jl) - zdh * rhoic /rhosn )329 ht_s(ji,jj,jl) = MAX( 0._wp, ht_s(ji,jj,jl) - zdh * rhoic * r1_rhosn ) 332 330 333 331 ! ice volume, salt content, age content … … 336 334 smv_i(ji,jj,jl) = MIN( sm_i(ji,jj,jl) , sss_m(ji,jj) ) * v_i(ji,jj,jl) ! salt content 337 335 oa_i(ji,jj,jl) = o_i(ji,jj,jl) * a_i(ji,jj,jl) ! age content 338 END DO ! ji339 END DO ! jj340 END DO ! jl336 END DO 337 END DO 338 END DO 341 339 342 340 ! Snow temperature and heat content … … 345 343 DO jj = 1, jpj 346 344 DO ji = 1, jpi 347 t_s(ji,jj,jk,jl) = zswitch(ji,jj) * ztm_i_ini(zhemis(ji,jj)) + ( 1._wp - zswitch(ji,jj) ) * rt t345 t_s(ji,jj,jk,jl) = zswitch(ji,jj) * ztm_i_ini(zhemis(ji,jj)) + ( 1._wp - zswitch(ji,jj) ) * rt0 348 346 ! Snow energy of melting 349 e_s(ji,jj,jk,jl) = zswitch(ji,jj) * rhosn * ( cpic * ( rtt - t_s(ji,jj,jk,jl) ) + lfus ) 350 ! Change dimensions 351 e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) / unit_fac 352 ! Multiply by volume, so that heat content in Joules 353 e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * area(ji,jj) * v_s(ji,jj,jl) / nlay_s 354 END DO ! ji 355 END DO ! jj 356 END DO ! jl 357 END DO ! jk 347 e_s(ji,jj,jk,jl) = zswitch(ji,jj) * rhosn * ( cpic * ( rt0 - t_s(ji,jj,jk,jl) ) + lfus ) 348 349 ! Mutliply by volume, and divide by number of layers to get heat content in J/m2 350 e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * v_s(ji,jj,jl) * r1_nlay_s 351 END DO 352 END DO 353 END DO 354 END DO 358 355 359 356 ! Ice salinity, temperature and heat content … … 362 359 DO jj = 1, jpj 363 360 DO ji = 1, jpi 364 t_i(ji,jj,jk,jl) = zswitch(ji,jj) * ztm_i_ini(zhemis(ji,jj)) + ( 1._wp - zswitch(ji,jj) ) * rt t365 s_i(ji,jj,jk,jl) = zswitch(ji,jj) * zsm_i_ini(zhemis(ji,jj)) !+ ( 1._wp - zswitch(ji,jj) ) * s_i_min366 ztmelts = - tmut * s_i(ji,jj,jk,jl) + rt t!Melting temperature in K361 t_i(ji,jj,jk,jl) = zswitch(ji,jj) * ztm_i_ini(zhemis(ji,jj)) + ( 1._wp - zswitch(ji,jj) ) * rt0 362 s_i(ji,jj,jk,jl) = zswitch(ji,jj) * zsm_i_ini(zhemis(ji,jj)) !+ ( 1._wp - zswitch(ji,jj) ) * rn_simin 363 ztmelts = - tmut * s_i(ji,jj,jk,jl) + rt0 !Melting temperature in K 367 364 368 365 ! heat content per unit volume 369 366 e_i(ji,jj,jk,jl) = zswitch(ji,jj) * rhoic * ( cpic * ( ztmelts - t_i(ji,jj,jk,jl) ) & 370 + lfus * ( 1._wp - (ztmelts-rtt) / MIN((t_i(ji,jj,jk,jl)-rtt),-epsi20) ) & 371 - rcp * ( ztmelts - rtt ) ) 372 373 ! Correct dimensions to avoid big values 374 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / unit_fac 375 376 ! Mutliply by ice volume, and divide by number of layers to get heat content in J 377 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * area(ji,jj) * v_i(ji,jj,jl) / nlay_i 378 END DO ! ji 379 END DO ! jj 380 END DO ! jl 381 END DO ! jk 367 + lfus * ( 1._wp - (ztmelts-rt0) / MIN((t_i(ji,jj,jk,jl)-rt0),-epsi20) ) & 368 - rcp * ( ztmelts - rt0 ) ) 369 370 ! Mutliply by ice volume, and divide by number of layers to get heat content in J/m2 371 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * v_i(ji,jj,jl) * r1_nlay_i 372 END DO 373 END DO 374 END DO 375 END DO 382 376 383 377 tn_ice (:,:,:) = t_su (:,:,:) 384 378 385 379 ELSE 386 ! if ln_ limini=false380 ! if ln_iceini=false 387 381 a_i (:,:,:) = 0._wp 388 382 v_i (:,:,:) = 0._wp … … 400 394 DO jl = 1, jpl 401 395 DO jk = 1, nlay_i 402 t_i(:,:,jk,jl) = rt t * tms(:,:)396 t_i(:,:,jk,jl) = rt0 * tmask(:,:,1) 403 397 END DO 404 398 DO jk = 1, nlay_s 405 t_s(:,:,jk,jl) = rt t * tms(:,:)399 t_s(:,:,jk,jl) = rt0 * tmask(:,:,1) 406 400 END DO 407 401 END DO 408 402 409 ENDIF ! ln_ limini403 ENDIF ! ln_iceini 410 404 411 405 at_i (:,:) = 0.0_wp … … 481 475 !! 8.5 ! 07-11 (M. Vancoppenolle) rewritten initialization 482 476 !!----------------------------------------------------------------------------- 483 NAMELIST/namiceini/ ln_ limini, thres_sst, hts_ini_n, hts_ini_s, hti_ini_n,hti_ini_s, &484 & ati_ini_n, ati_ini_s, smi_ini_n, smi_ini_s, tmi_ini_n,tmi_ini_s477 NAMELIST/namiceini/ ln_iceini, rn_thres_sst, rn_hts_ini_n, rn_hts_ini_s, rn_hti_ini_n, rn_hti_ini_s, & 478 & rn_ati_ini_n, rn_ati_ini_s, rn_smi_ini_n, rn_smi_ini_s, rn_tmi_ini_n, rn_tmi_ini_s 485 479 INTEGER :: ios ! Local integer output status for namelist read 486 480 !!----------------------------------------------------------------------------- … … 502 496 WRITE(numout,*) 'lim_istate_init : ice parameters inititialisation ' 503 497 WRITE(numout,*) '~~~~~~~~~~~~~~~' 504 WRITE(numout,*) ' initialization with ice (T) or not (F) ln_ limini = ', ln_limini505 WRITE(numout,*) ' threshold water temp. for initial sea-ice thres_sst = ',thres_sst506 WRITE(numout,*) ' initial snow thickness in the north hts_ini_n = ',hts_ini_n507 WRITE(numout,*) ' initial snow thickness in the south hts_ini_s = ',hts_ini_s508 WRITE(numout,*) ' initial ice thickness in the north hti_ini_n = ',hti_ini_n509 WRITE(numout,*) ' initial ice thickness in the south hti_ini_s = ',hti_ini_s510 WRITE(numout,*) ' initial ice concentr. in the north ati_ini_n = ',ati_ini_n511 WRITE(numout,*) ' initial ice concentr. in the north ati_ini_s = ',ati_ini_s512 WRITE(numout,*) ' initial ice salinity in the north smi_ini_n = ',smi_ini_n513 WRITE(numout,*) ' initial ice salinity in the south smi_ini_s = ',smi_ini_s514 WRITE(numout,*) ' initial ice/snw temp in the north tmi_ini_n = ',tmi_ini_n515 WRITE(numout,*) ' initial ice/snw temp in the south tmi_ini_s = ',tmi_ini_s498 WRITE(numout,*) ' initialization with ice (T) or not (F) ln_iceini = ', ln_iceini 499 WRITE(numout,*) ' threshold water temp. for initial sea-ice rn_thres_sst = ', rn_thres_sst 500 WRITE(numout,*) ' initial snow thickness in the north rn_hts_ini_n = ', rn_hts_ini_n 501 WRITE(numout,*) ' initial snow thickness in the south rn_hts_ini_s = ', rn_hts_ini_s 502 WRITE(numout,*) ' initial ice thickness in the north rn_hti_ini_n = ', rn_hti_ini_n 503 WRITE(numout,*) ' initial ice thickness in the south rn_hti_ini_s = ', rn_hti_ini_s 504 WRITE(numout,*) ' initial ice concentr. in the north rn_ati_ini_n = ', rn_ati_ini_n 505 WRITE(numout,*) ' initial ice concentr. in the north rn_ati_ini_s = ', rn_ati_ini_s 506 WRITE(numout,*) ' initial ice salinity in the north rn_smi_ini_n = ', rn_smi_ini_n 507 WRITE(numout,*) ' initial ice salinity in the south rn_smi_ini_s = ', rn_smi_ini_s 508 WRITE(numout,*) ' initial ice/snw temp in the north rn_tmi_ini_n = ', rn_tmi_ini_n 509 WRITE(numout,*) ' initial ice/snw temp in the south rn_tmi_ini_s = ', rn_tmi_ini_s 516 510 ENDIF 517 511 -
branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90
r4990 r5682 18 18 USE thd_ice ! LIM thermodynamics 19 19 USE ice ! LIM variables 20 USE par_ice ! LIM parameters21 20 USE dom_ice ! LIM domain 22 USE limthd_lac ! LIM23 21 USE limvar ! LIM 24 USE in_out_manager ! I/O manager25 22 USE lbclnk ! lateral boundary condition - MPP exchanges 26 23 USE lib_mpp ! MPP library 27 24 USE wrk_nemo ! work arrays 28 25 USE prtctl ! Print control 29 ! Check budget (Rousset) 26 27 USE in_out_manager ! I/O manager 30 28 USE iom ! I/O manager 31 29 USE lib_fortran ! glob_sum … … 40 38 PUBLIC lim_itd_me_icestrength 41 39 PUBLIC lim_itd_me_init 42 PUBLIC lim_itd_me_zapsmall 43 PUBLIC lim_itd_me_alloc ! called by iceini.F90 40 PUBLIC lim_itd_me_alloc ! called by sbc_lim_init 44 41 45 42 !----------------------------------------------------------------------- … … 125 122 !! and Elizabeth C. Hunke, LANL are gratefully acknowledged 126 123 !!--------------------------------------------------------------------! 127 INTEGER :: ji, jj, jk, jl ! dummy loop index 128 INTEGER :: niter, nitermax = 20 ! local integer 129 LOGICAL :: asum_error ! flag for asum .ne. 1 124 INTEGER :: ji, jj, jk, jl ! dummy loop index 125 INTEGER :: niter ! local integer 130 126 INTEGER :: iterate_ridging ! if true, repeat the ridging 131 REAL(wp) :: w1, tmpfac! local scalar127 REAL(wp) :: za, zfac ! local scalar 132 128 CHARACTER (len = 15) :: fieldid 133 REAL(wp), POINTER, DIMENSION(:,:) :: closing_net ! net rate at which area is removed (1/s) 134 ! (ridging ice area - area of new ridges) / dt 135 REAL(wp), POINTER, DIMENSION(:,:) :: divu_adv ! divu as implied by transport scheme (1/s) 136 REAL(wp), POINTER, DIMENSION(:,:) :: opning ! rate of opening due to divergence/shear 137 REAL(wp), POINTER, DIMENSION(:,:) :: closing_gross ! rate at which area removed, not counting area of new ridges 138 REAL(wp), POINTER, DIMENSION(:,:) :: msnow_mlt ! mass of snow added to ocean (kg m-2) 139 REAL(wp), POINTER, DIMENSION(:,:) :: esnow_mlt ! energy needed to melt snow in ocean (J m-2) 140 REAL(wp), POINTER, DIMENSION(:,:) :: vt_i_init, vt_i_final ! ice volume summed over categories 129 REAL(wp), POINTER, DIMENSION(:,:) :: closing_net ! net rate at which area is removed (1/s) 130 ! (ridging ice area - area of new ridges) / dt 131 REAL(wp), POINTER, DIMENSION(:,:) :: divu_adv ! divu as implied by transport scheme (1/s) 132 REAL(wp), POINTER, DIMENSION(:,:) :: opning ! rate of opening due to divergence/shear 133 REAL(wp), POINTER, DIMENSION(:,:) :: closing_gross ! rate at which area removed, not counting area of new ridges 134 REAL(wp), POINTER, DIMENSION(:,:) :: msnow_mlt ! mass of snow added to ocean (kg m-2) 135 REAL(wp), POINTER, DIMENSION(:,:) :: esnow_mlt ! energy needed to melt snow in ocean (J m-2) 136 REAL(wp), POINTER, DIMENSION(:,:) :: vt_i_init, vt_i_final ! ice volume summed over categories 137 ! 138 INTEGER, PARAMETER :: nitermax = 20 141 139 ! 142 140 REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b … … 144 142 IF( nn_timing == 1 ) CALL timing_start('limitd_me') 145 143 146 CALL wrk_alloc( jpi, 144 CALL wrk_alloc( jpi,jpj, closing_net, divu_adv, opning, closing_gross, msnow_mlt, esnow_mlt, vt_i_init, vt_i_final ) 147 145 148 146 IF(ln_ctl) THEN … … 156 154 IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limitd_me', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 157 155 156 CALL lim_var_zapsmall 157 CALL lim_var_glo2eqv ! equivalent variables, requested for rafting 158 158 159 !-----------------------------------------------------------------------------! 159 160 ! 1) Thickness categories boundaries, ice / o.w. concentrations, init_ons 160 161 !-----------------------------------------------------------------------------! 161 Cp = 0.5 * grav * (rau0-rhoic) * rhoic / rau0! proport const for PE162 Cp = 0.5 * grav * (rau0-rhoic) * rhoic * r1_rau0 ! proport const for PE 162 163 ! 163 164 CALL lim_itd_me_ridgeprep ! prepare ridging … … 193 194 ! (thick, newly ridged ice). 194 195 195 closing_net(ji,jj) = Cs * 0.5 * ( Delta_i(ji,jj) - ABS( divu_i(ji,jj) ) ) - MIN( divu_i(ji,jj), 0._wp )196 closing_net(ji,jj) = rn_cs * 0.5 * ( delta_i(ji,jj) - ABS( divu_i(ji,jj) ) ) - MIN( divu_i(ji,jj), 0._wp ) 196 197 197 198 ! 2.2 divu_adv … … 237 238 ! Reduce the closing rate if more than 100% of the open water 238 239 ! would be removed. Reduce the opening rate proportionately. 239 IF ( ato_i(ji,jj) .GT. epsi10 .AND. athorn(ji,jj,0) .GT. 0.0 ) THEN 240 w1 = athorn(ji,jj,0) * closing_gross(ji,jj) * rdt_ice 241 IF ( w1 .GT. ato_i(ji,jj)) THEN 242 tmpfac = ato_i(ji,jj) / w1 243 closing_gross(ji,jj) = closing_gross(ji,jj) * tmpfac 244 opning(ji,jj) = opning(ji,jj) * tmpfac 245 ENDIF !w1 246 ENDIF !at0i and athorn 247 248 END DO ! ji 249 END DO ! jj 240 za = athorn(ji,jj,0) * closing_gross(ji,jj) * rdt_ice 241 IF( za > epsi20 ) THEN 242 zfac = MIN( 1._wp, ato_i(ji,jj) / za ) 243 closing_gross(ji,jj) = closing_gross(ji,jj) * zfac 244 opning (ji,jj) = opning (ji,jj) * zfac 245 ENDIF 246 247 END DO 248 END DO 250 249 251 250 ! correction to closing rate / opening if excessive ice removal … … 253 252 ! Reduce the closing rate if more than 100% of any ice category 254 253 ! would be removed. Reduce the opening rate proportionately. 255 256 254 DO jl = 1, jpl 257 255 DO jj = 1, jpj 258 256 DO ji = 1, jpi 259 IF ( a_i(ji,jj,jl) > epsi10 .AND. athorn(ji,jj,jl) > 0._wp )THEN 260 w1 = athorn(ji,jj,jl) * closing_gross(ji,jj) * rdt_ice 261 IF ( w1 > a_i(ji,jj,jl) ) THEN 262 tmpfac = a_i(ji,jj,jl) / w1 263 closing_gross(ji,jj) = closing_gross(ji,jj) * tmpfac 264 opning (ji,jj) = opning (ji,jj) * tmpfac 265 ENDIF 257 za = athorn(ji,jj,jl) * closing_gross(ji,jj) * rdt_ice 258 IF( za > epsi20 ) THEN 259 zfac = MIN( 1._wp, a_i(ji,jj,jl) / za ) 260 closing_gross(ji,jj) = closing_gross(ji,jj) * zfac 261 opning (ji,jj) = opning (ji,jj) * zfac 266 262 ENDIF 267 END DO !ji268 END DO ! jj269 END DO !jl263 END DO 264 END DO 265 END DO 270 266 271 267 ! 3.3 Redistribute area, volume, and energy. … … 276 272 ! 3.4 Compute total area of ice plus open water after ridging. 277 273 !-----------------------------------------------------------------------------! 278 279 CALL lim_itd_me_asumr 274 ! This is in general not equal to one because of divergence during transport 275 asum(:,:) = ato_i(:,:) 276 DO jl = 1, jpl 277 asum(:,:) = asum(:,:) + a_i(:,:,jl) 278 END DO 280 279 281 280 ! 3.5 Do we keep on iterating ??? … … 288 287 DO jj = 1, jpj 289 288 DO ji = 1, jpi 290 IF (ABS(asum(ji,jj) - kamax ) .LT.epsi10) THEN289 IF (ABS(asum(ji,jj) - kamax ) < epsi10) THEN 291 290 closing_net(ji,jj) = 0._wp 292 291 opning (ji,jj) = 0._wp … … 324 323 ! Convert ridging rate diagnostics to correct units. 325 324 ! Update fresh water and heat fluxes due to snow melt. 326 327 asum_error = .false.328 329 325 DO jj = 1, jpj 330 326 DO ji = 1, jpi 331 332 IF(ABS(asum(ji,jj) - kamax) > epsi10 ) asum_error = .true.333 327 334 328 dardg1dt(ji,jj) = dardg1dt(ji,jj) * r1_rdtice … … 341 335 !-----------------------------------------------------------------------------! 342 336 wfx_snw(ji,jj) = wfx_snw(ji,jj) + msnow_mlt(ji,jj) * r1_rdtice ! fresh water source for ocean 343 hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + esnow_mlt(ji,jj) * unit_fac / area(ji,jj) * r1_rdtice! heat sink for ocean (<0, W.m-2)337 hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + esnow_mlt(ji,jj) * r1_rdtice ! heat sink for ocean (<0, W.m-2) 344 338 345 339 END DO … … 347 341 348 342 ! Check if there is a ridging error 349 DO jj = 1, jpj 350 DO ji = 1, jpi 351 IF( ABS( asum(ji,jj) - kamax) > epsi10 ) THEN ! there is a bug 352 WRITE(numout,*) ' ' 353 WRITE(numout,*) ' ALERTE : Ridging error: total area = ', asum(ji,jj) 354 WRITE(numout,*) ' limitd_me ' 355 WRITE(numout,*) ' POINT : ', ji, jj 356 WRITE(numout,*) ' jpl, a_i, athorn ' 357 WRITE(numout,*) 0, ato_i(ji,jj), athorn(ji,jj,0) 358 DO jl = 1, jpl 359 WRITE(numout,*) jl, a_i(ji,jj,jl), athorn(ji,jj,jl) 360 END DO 361 ENDIF ! asum 362 363 END DO !ji 364 END DO !jj 343 IF( lwp ) THEN 344 DO jj = 1, jpj 345 DO ji = 1, jpi 346 IF( ABS( asum(ji,jj) - kamax) > epsi10 ) THEN ! there is a bug 347 WRITE(numout,*) ' ' 348 WRITE(numout,*) ' ALERTE : Ridging error: total area = ', asum(ji,jj) 349 WRITE(numout,*) ' limitd_me ' 350 WRITE(numout,*) ' POINT : ', ji, jj 351 WRITE(numout,*) ' jpl, a_i, athorn ' 352 WRITE(numout,*) 0, ato_i(ji,jj), athorn(ji,jj,0) 353 DO jl = 1, jpl 354 WRITE(numout,*) jl, a_i(ji,jj,jl), athorn(ji,jj,jl) 355 END DO 356 ENDIF 357 END DO 358 END DO 359 END IF 365 360 366 361 ! Conservation check … … 371 366 ENDIF 372 367 368 CALL lim_var_agg( 1 ) 369 373 370 !-----------------------------------------------------------------------------! 374 ! 6) Updating state variables and trend terms (done in limupdate)371 ! control prints 375 372 !-----------------------------------------------------------------------------! 376 CALL lim_var_glo2eqv 377 CALL lim_itd_me_zapsmall 378 379 380 IF(ln_ctl) THEN ! Control print 373 IF(ln_ctl) THEN 374 CALL lim_var_glo2eqv 375 381 376 CALL prt_ctl_info(' ') 382 377 CALL prt_ctl_info(' - Cell values : ') 383 378 CALL prt_ctl_info(' ~~~~~~~~~~~~~ ') 384 CALL prt_ctl(tab2d_1= area, clinfo1=' lim_itd_me : cell area :')379 CALL prt_ctl(tab2d_1=e12t , clinfo1=' lim_itd_me : cell area :') 385 380 CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_itd_me : at_i :') 386 381 CALL prt_ctl(tab2d_1=vt_i , clinfo1=' lim_itd_me : vt_i :') … … 436 431 !!---------------------------------------------------------------------- 437 432 INTEGER, INTENT(in) :: kstrngth ! = 1 for Rothrock formulation, 0 for Hibler (1979) 438 439 INTEGER :: ji,jj, jl ! dummy loop indices 440 INTEGER :: ksmooth ! smoothing the resistance to deformation 441 INTEGER :: numts_rm ! number of time steps for the P smoothing 442 REAL(wp) :: hi, zw1, zp, zdummy, zzc, z1_3 ! local scalars 433 INTEGER :: ji,jj, jl ! dummy loop indices 434 INTEGER :: ksmooth ! smoothing the resistance to deformation 435 INTEGER :: numts_rm ! number of time steps for the P smoothing 436 REAL(wp) :: zhi, zp, z1_3 ! local scalars 443 437 REAL(wp), POINTER, DIMENSION(:,:) :: zworka ! temporary array used here 444 438 !!---------------------------------------------------------------------- … … 466 460 ! 467 461 IF( a_i(ji,jj,jl) > epsi10 .AND. athorn(ji,jj,jl) > 0._wp ) THEN 468 hi = v_i(ji,jj,jl) / a_i(ji,jj,jl)462 zhi = v_i(ji,jj,jl) / a_i(ji,jj,jl) 469 463 !---------------------------- 470 464 ! PE loss from deforming ice 471 465 !---------------------------- 472 strength(ji,jj) = strength(ji,jj) - athorn(ji,jj,jl) * hi *hi466 strength(ji,jj) = strength(ji,jj) - athorn(ji,jj,jl) * zhi * zhi 473 467 474 468 !-------------------------- 475 469 ! PE gain from rafting ice 476 470 !-------------------------- 477 strength(ji,jj) = strength(ji,jj) + 2._wp * araft(ji,jj,jl) * hi *hi471 strength(ji,jj) = strength(ji,jj) + 2._wp * araft(ji,jj,jl) * zhi * zhi 478 472 479 473 !---------------------------- 480 474 ! PE gain from ridging ice 481 475 !---------------------------- 482 strength(ji,jj) = strength(ji,jj) + aridge(ji,jj,jl) /krdg(ji,jj,jl) &483 * z1_3 * ( hrmax(ji,jj,jl)**3 - hrmin(ji,jj,jl)**3) / ( hrmax(ji,jj,jl)-hrmin(ji,jj,jl) )484 !!gm Optimization: (a**3-b**3)/(a-b) = a*a+ab+b*b ==> less costly operations even if a**3 is replaced by a*a*a...485 ENDIF ! aicen > epsi10476 strength(ji,jj) = strength(ji,jj) + aridge(ji,jj,jl) / krdg(ji,jj,jl) & 477 * z1_3 * ( hrmax(ji,jj,jl)**2 + hrmin(ji,jj,jl)**2 + hrmax(ji,jj,jl) * hrmin(ji,jj,jl) ) 478 !!(a**3-b**3)/(a-b) = a*a+ab+b*b 479 ENDIF 486 480 ! 487 END DO ! ji 488 END DO !jj 489 END DO !jl 490 491 zzc = Cf * Cp ! where Cp = (g/2)*(rhow-rhoi)*(rhoi/rhow) and Cf accounts for frictional dissipation 492 strength(:,:) = zzc * strength(:,:) / aksum(:,:) 493 481 END DO 482 END DO 483 END DO 484 485 strength(:,:) = rn_pe_rdg * Cp * strength(:,:) / aksum(:,:) 486 ! where Cp = (g/2)*(rhow-rhoi)*(rhoi/rhow) and rn_pe_rdg accounts for frictional dissipation 494 487 ksmooth = 1 495 488 … … 499 492 ELSE ! kstrngth ne 1: Hibler (1979) form 500 493 ! 501 strength(:,:) = Pstar * vt_i(:,:) * EXP( - C_rhg * ( 1._wp - at_i(:,:) ) )494 strength(:,:) = rn_pstar * vt_i(:,:) * EXP( - rn_crhg * ( 1._wp - at_i(:,:) ) ) 502 495 ! 503 496 ksmooth = 1 … … 511 504 ! CAN BE REMOVED 512 505 ! 513 IF ( brinstren_swi == 1) THEN506 IF( ln_icestr_bvf ) THEN 514 507 515 508 DO jj = 1, jpj 516 509 DO ji = 1, jpi 517 IF ( bv_i(ji,jj) .GT. 0.0 ) THEN518 zdummy = MIN ( bv_i(ji,jj), 0.10 ) * MIN( bv_i(ji,jj), 0.10 )519 ELSE520 zdummy = 0.0521 ENDIF522 510 strength(ji,jj) = strength(ji,jj) * exp(-5.88*SQRT(MAX(bv_i(ji,jj),0.0))) 523 END DO ! j524 END DO ! i511 END DO 512 END DO 525 513 526 514 ENDIF … … 538 526 CALL lbc_lnk( strength, 'T', 1. ) 539 527 540 DO jj = 2, jpj - 1 541 DO ji = 2, jpi - 1 542 IF ( ( asum(ji,jj) - ato_i(ji,jj) ) .GT. epsi10) THEN ! ice is 543 ! present 544 zworka(ji,jj) = 4.0 * strength(ji,jj) & 545 & + strength(ji-1,jj) * tms(ji-1,jj) & 546 & + strength(ji+1,jj) * tms(ji+1,jj) & 547 & + strength(ji,jj-1) * tms(ji,jj-1) & 548 & + strength(ji,jj+1) * tms(ji,jj+1) 549 550 zw1 = 4.0 + tms(ji-1,jj) + tms(ji+1,jj) + tms(ji,jj-1) + tms(ji,jj+1) 551 zworka(ji,jj) = zworka(ji,jj) / zw1 528 DO jj = 2, jpjm1 529 DO ji = 2, jpim1 530 IF ( ( asum(ji,jj) - ato_i(ji,jj) ) > 0._wp) THEN 531 zworka(ji,jj) = ( 4.0 * strength(ji,jj) & 532 & + strength(ji-1,jj) * tmask(ji-1,jj,1) + strength(ji+1,jj) * tmask(ji+1,jj,1) & 533 & + strength(ji,jj-1) * tmask(ji,jj-1,1) + strength(ji,jj+1) * tmask(ji,jj+1,1) & 534 & ) / ( 4.0 + tmask(ji-1,jj,1) + tmask(ji+1,jj,1) + tmask(ji,jj-1,1) + tmask(ji,jj+1,1) ) 552 535 ELSE 553 536 zworka(ji,jj) = 0._wp … … 556 539 END DO 557 540 558 DO jj = 2, jpj -1559 DO ji = 2, jpi -1541 DO jj = 2, jpjm1 542 DO ji = 2, jpim1 560 543 strength(ji,jj) = zworka(ji,jj) 561 544 END DO … … 563 546 CALL lbc_lnk( strength, 'T', 1. ) 564 547 565 ENDIF ! ksmooth548 ENDIF 566 549 567 550 !-------------------- … … 580 563 DO jj = 1, jpj - 1 581 564 DO ji = 1, jpi - 1 582 IF ( ( asum(ji,jj) - ato_i(ji,jj) ) .GT. epsi10) THEN ! ice is present565 IF ( ( asum(ji,jj) - ato_i(ji,jj) ) > 0._wp) THEN 583 566 numts_rm = 1 ! number of time steps for the running mean 584 IF ( strp1(ji,jj) .GT.0.0 ) numts_rm = numts_rm + 1585 IF ( strp2(ji,jj) .GT.0.0 ) numts_rm = numts_rm + 1567 IF ( strp1(ji,jj) > 0.0 ) numts_rm = numts_rm + 1 568 IF ( strp2(ji,jj) > 0.0 ) numts_rm = numts_rm + 1 586 569 zp = ( strength(ji,jj) + strp1(ji,jj) + strp2(ji,jj) ) / numts_rm 587 570 strp2(ji,jj) = strp1(ji,jj) … … 612 595 !!---------------------------------------------------------------------! 613 596 INTEGER :: ji,jj, jl ! dummy loop indices 614 REAL(wp) :: Gstari, astari, hi, hrmean, zdummy ! local scalar597 REAL(wp) :: Gstari, astari, zhi, hrmean, zdummy ! local scalar 615 598 REAL(wp), POINTER, DIMENSION(:,:) :: zworka ! temporary array used here 616 599 REAL(wp), POINTER, DIMENSION(:,:,:) :: Gsum ! Gsum(n) = sum of areas in categories 0 to n … … 620 603 CALL wrk_alloc( jpi,jpj,jpl+2, Gsum, kkstart = -1 ) 621 604 622 Gstari = 1.0/ Gstar623 astari = 1.0/ astar605 Gstari = 1.0/rn_gstar 606 astari = 1.0/rn_astar 624 607 aksum(:,:) = 0.0 625 608 athorn(:,:,:) = 0.0 … … 632 615 633 616 ! ! Zero out categories with very small areas 634 CALL lim_ itd_me_zapsmall617 CALL lim_var_zapsmall 635 618 636 619 !------------------------------------------------------------------------------! … … 639 622 640 623 ! Compute total area of ice plus open water. 641 CALL lim_itd_me_asumr 642 ! This is in general not equal to one 643 ! because of divergence during transport 624 ! This is in general not equal to one because of divergence during transport 625 asum(:,:) = ato_i(:,:) 626 DO jl = 1, jpl 627 asum(:,:) = asum(:,:) + a_i(:,:,jl) 628 END DO 644 629 645 630 ! Compute cumulative thickness distribution function … … 649 634 650 635 Gsum(:,:,-1) = 0._wp 651 652 DO jj = 1, jpj 653 DO ji = 1, jpi 654 IF( ato_i(ji,jj) > epsi10 ) THEN ; Gsum(ji,jj,0) = ato_i(ji,jj) 655 ELSE ; Gsum(ji,jj,0) = 0._wp 656 ENDIF 657 END DO 658 END DO 636 Gsum(:,:,0 ) = ato_i(:,:) 659 637 660 638 ! for each value of h, you have to add ice concentration then 661 639 DO jl = 1, jpl 662 DO jj = 1, jpj 663 DO ji = 1, jpi 664 IF( a_i(ji,jj,jl) .GT. epsi10 ) THEN ; Gsum(ji,jj,jl) = Gsum(ji,jj,jl-1) + a_i(ji,jj,jl) 665 ELSE ; Gsum(ji,jj,jl) = Gsum(ji,jj,jl-1) 666 ENDIF 667 END DO 668 END DO 640 Gsum(:,:,jl) = Gsum(:,:,jl-1) + a_i(:,:,jl) 669 641 END DO 670 642 … … 687 659 !----------------------------------------------------------------- 688 660 689 IF( partfun_swi== 0 ) THEN !--- Linear formulation (Thorndike et al., 1975)661 IF( nn_partfun == 0 ) THEN !--- Linear formulation (Thorndike et al., 1975) 690 662 DO jl = 0, jpl 691 663 DO jj = 1, jpj 692 664 DO ji = 1, jpi 693 IF( Gsum(ji,jj,jl) < Gstar) THEN694 athorn(ji,jj,jl) = Gstari * ( Gsum(ji,jj,jl)-Gsum(ji,jj,jl-1)) * &695 (2.0 - (Gsum(ji,jj,jl-1)+Gsum(ji,jj,jl))*Gstari)696 ELSEIF (Gsum(ji,jj,jl-1) < Gstar) THEN697 athorn(ji,jj,jl) = Gstari * ( Gstar-Gsum(ji,jj,jl-1)) * &698 (2.0 - (Gsum(ji,jj,jl-1)+Gstar)*Gstari)665 IF( Gsum(ji,jj,jl) < rn_gstar) THEN 666 athorn(ji,jj,jl) = Gstari * ( Gsum(ji,jj,jl) - Gsum(ji,jj,jl-1) ) * & 667 & ( 2.0 - (Gsum(ji,jj,jl-1) + Gsum(ji,jj,jl) ) * Gstari ) 668 ELSEIF (Gsum(ji,jj,jl-1) < rn_gstar) THEN 669 athorn(ji,jj,jl) = Gstari * ( rn_gstar - Gsum(ji,jj,jl-1) ) * & 670 & ( 2.0 - ( Gsum(ji,jj,jl-1) + rn_gstar ) * Gstari ) 699 671 ELSE 700 672 athorn(ji,jj,jl) = 0.0 701 673 ENDIF 702 END DO ! ji703 END DO ! jj704 END DO ! jl674 END DO 675 END DO 676 END DO 705 677 706 678 ELSE !--- Exponential, more stable formulation (Lipscomb et al, 2007) 707 679 ! 708 680 zdummy = 1._wp / ( 1._wp - EXP(-astari) ) ! precompute exponential terms using Gsum as a work array 709 710 681 DO jl = -1, jpl 711 682 Gsum(:,:,jl) = EXP( -Gsum(:,:,jl) * astari ) * zdummy 712 END DO !jl683 END DO 713 684 DO jl = 0, jpl 714 685 athorn(:,:,jl) = Gsum(:,:,jl-1) - Gsum(:,:,jl) 715 686 END DO 716 687 ! 717 ENDIF ! partfun_swi718 719 IF( raft_swi == 1) THEN ! Ridging and rafting ice participation functions688 ENDIF 689 690 IF( ln_rafting ) THEN ! Ridging and rafting ice participation functions 720 691 ! 721 692 DO jl = 1, jpl 722 693 DO jj = 1, jpj 723 694 DO ji = 1, jpi 724 IF ( athorn(ji,jj,jl) .GT.0._wp ) THEN695 IF ( athorn(ji,jj,jl) > 0._wp ) THEN 725 696 !!gm TANH( -X ) = - TANH( X ) so can be computed only 1 time.... 726 aridge(ji,jj,jl) = ( TANH ( Craft * ( ht_i(ji,jj,jl) - hparmeter) ) + 1.0 ) * 0.5 * athorn(ji,jj,jl)727 araft (ji,jj,jl) = ( TANH ( - Craft * ( ht_i(ji,jj,jl) - hparmeter) ) + 1.0 ) * 0.5 * athorn(ji,jj,jl)697 aridge(ji,jj,jl) = ( TANH ( rn_craft * ( ht_i(ji,jj,jl) - rn_hraft ) ) + 1.0 ) * 0.5 * athorn(ji,jj,jl) 698 araft (ji,jj,jl) = ( TANH ( -rn_craft * ( ht_i(ji,jj,jl) - rn_hraft ) ) + 1.0 ) * 0.5 * athorn(ji,jj,jl) 728 699 IF ( araft(ji,jj,jl) < epsi06 ) araft(ji,jj,jl) = 0._wp 729 700 aridge(ji,jj,jl) = MAX( athorn(ji,jj,jl) - araft(ji,jj,jl), 0.0 ) 730 ENDIF ! athorn731 END DO ! ji732 END DO ! jj733 END DO ! jl734 735 ELSE ! raft_swi = 0701 ENDIF 702 END DO 703 END DO 704 END DO 705 706 ELSE 736 707 ! 737 708 DO jl = 1, jpl … … 741 712 ENDIF 742 713 743 IF ( raft_swi == 1) THEN744 745 IF( MAXVAL(aridge + araft - athorn(:,:,1:jpl)) .GT. epsi10) THEN714 IF( ln_rafting ) THEN 715 716 IF( MAXVAL(aridge + araft - athorn(:,:,1:jpl)) > epsi10 .AND. lwp ) THEN 746 717 DO jl = 1, jpl 747 718 DO jj = 1, jpj 748 719 DO ji = 1, jpi 749 IF ( aridge(ji,jj,jl) + araft(ji,jj,jl) - athorn(ji,jj,jl) .GT.epsi10 ) THEN720 IF ( aridge(ji,jj,jl) + araft(ji,jj,jl) - athorn(ji,jj,jl) > epsi10 ) THEN 750 721 WRITE(numout,*) ' ALERTE 96 : wrong participation function ... ' 751 722 WRITE(numout,*) ' ji, jj, jl : ', ji, jj, jl … … 793 764 DO ji = 1, jpi 794 765 795 IF (a_i(ji,jj,jl) .GT. epsi10 .AND. athorn(ji,jj,jl) .GT.0.0 ) THEN796 hi = v_i(ji,jj,jl) / a_i(ji,jj,jl)797 hrmean = MAX(SQRT( Hstar*hi),hi*krdgmin)798 hrmin(ji,jj,jl) = MIN(2.0* hi, 0.5*(hrmean +hi))766 IF (a_i(ji,jj,jl) > epsi10 .AND. athorn(ji,jj,jl) > 0.0 ) THEN 767 zhi = v_i(ji,jj,jl) / a_i(ji,jj,jl) 768 hrmean = MAX(SQRT(rn_hstar*zhi), zhi*krdgmin) 769 hrmin(ji,jj,jl) = MIN(2.0*zhi, 0.5*(hrmean + zhi)) 799 770 hrmax(ji,jj,jl) = 2.0*hrmean - hrmin(ji,jj,jl) 800 hraft(ji,jj,jl) = kraft* hi801 krdg(ji,jj,jl) = hrmean / hi771 hraft(ji,jj,jl) = kraft*zhi 772 krdg(ji,jj,jl) = hrmean / zhi 802 773 ELSE 803 774 hraft(ji,jj,jl) = 0.0 … … 807 778 ENDIF 808 779 809 END DO ! ji810 END DO ! jj811 END DO ! jl780 END DO 781 END DO 782 END DO 812 783 813 784 ! Normalization factor : aksum, ensures mass conservation … … 841 812 LOGICAL, PARAMETER :: l_conservation_check = .true. ! if true, check conservation (useful for debugging) 842 813 ! 843 LOGICAL :: neg_ato_i ! flag for ato_i(i,j) < -puny844 LOGICAL :: large_afrac ! flag for afrac > 1845 LOGICAL :: large_afrft ! flag for afrac > 1846 814 INTEGER :: ji, jj, jl, jl1, jl2, jk ! dummy loop indices 847 815 INTEGER :: ij ! horizontal index, combines i and j loops 848 816 INTEGER :: icells ! number of cells with aicen > puny 849 REAL(wp) :: hL, hR, farea, zdummy, zdummy0, ztmelts ! left and right limits of integration 850 REAL(wp) :: zsstK ! SST in Kelvin 817 REAL(wp) :: hL, hR, farea, ztmelts ! left and right limits of integration 851 818 852 819 INTEGER , POINTER, DIMENSION(:) :: indxi, indxj ! compressed indices … … 864 831 REAL(wp), POINTER, DIMENSION(:,:) :: ardg1 , ardg2 ! area of ice ridged & new ridges 865 832 REAL(wp), POINTER, DIMENSION(:,:) :: vsrdg , esrdg ! snow volume & energy of ridging ice 866 REAL(wp), POINTER, DIMENSION(:,:) :: oirdg1, oirdg2 ! areal age content of ridged & rifging ice867 833 REAL(wp), POINTER, DIMENSION(:,:) :: dhr , dhr2 ! hrmax - hrmin & hrmax^2 - hrmin^2 868 834 … … 873 839 REAL(wp), POINTER, DIMENSION(:,:) :: srdg2 ! sal*volume of new ridges 874 840 REAL(wp), POINTER, DIMENSION(:,:) :: smsw ! sal*volume of water trapped into ridges 841 REAL(wp), POINTER, DIMENSION(:,:) :: oirdg1, oirdg2 ! ice age of ice ridged 875 842 876 843 REAL(wp), POINTER, DIMENSION(:,:) :: afrft ! fraction of category area rafted … … 878 845 REAL(wp), POINTER, DIMENSION(:,:) :: virft , vsrft ! ice & snow volume of rafting ice 879 846 REAL(wp), POINTER, DIMENSION(:,:) :: esrft , smrft ! snow energy & salinity of rafting ice 880 REAL(wp), POINTER, DIMENSION(:,:) :: oirft1, oirft2 ! areal age content of rafted ice & rafting ice847 REAL(wp), POINTER, DIMENSION(:,:) :: oirft1, oirft2 ! ice age of ice rafted 881 848 882 849 REAL(wp), POINTER, DIMENSION(:,:,:) :: eirft ! ice energy of rafting ice … … 886 853 !!---------------------------------------------------------------------- 887 854 888 CALL wrk_alloc( (jpi+1)*(jpj+1), indxi, indxj )889 CALL wrk_alloc( jpi, jpj, vice_init, vice_final, eice_init, eice_final )890 CALL wrk_alloc( jpi, jpj, afrac, fvol , ardg1, ardg2, vsrdg, esrdg, oirdg1, oirdg2, dhr, dhr2 )891 CALL wrk_alloc( jpi, jpj, vrdg1, vrdg2, vsw , srdg1, srdg2, smsw)892 CALL wrk_alloc( jpi, jpj, afrft, arft1, arft2, virft, vsrft, esrft, smrft, oirft1, oirft2 )893 CALL wrk_alloc( jpi, jpj, jpl, aicen_init, vicen_init, vsnwn_init, esnwn_init, smv_i_init, oa_i_init )894 CALL wrk_alloc( jpi, jpj, nlay_i +1, eirft, erdg1, erdg2, ersw )895 CALL wrk_alloc( jpi, jpj, nlay_i +1, jpl, eicen_init )855 CALL wrk_alloc( (jpi+1)*(jpj+1), indxi, indxj ) 856 CALL wrk_alloc( jpi, jpj, vice_init, vice_final, eice_init, eice_final ) 857 CALL wrk_alloc( jpi, jpj, afrac, fvol , ardg1, ardg2, vsrdg, esrdg, dhr, dhr2 ) 858 CALL wrk_alloc( jpi, jpj, vrdg1, vrdg2, vsw , srdg1, srdg2, smsw, oirdg1, oirdg2 ) 859 CALL wrk_alloc( jpi, jpj, afrft, arft1, arft2, virft, vsrft, esrft, smrft, oirft1, oirft2 ) 860 CALL wrk_alloc( jpi, jpj, jpl, aicen_init, vicen_init, vsnwn_init, esnwn_init, smv_i_init, oa_i_init ) 861 CALL wrk_alloc( jpi, jpj, nlay_i, eirft, erdg1, erdg2, ersw ) 862 CALL wrk_alloc( jpi, jpj, nlay_i, jpl, eicen_init ) 896 863 897 864 ! Conservation check … … 901 868 CALL lim_column_sum (jpl, v_i, vice_init ) 902 869 CALL lim_column_sum_energy (jpl, nlay_i, e_i, eice_init ) 903 DO ji = mi0( jiindx), mi1(jiindx)904 DO jj = mj0(j jindx), mj1(jjindx)870 DO ji = mi0(iiceprt), mi1(iiceprt) 871 DO jj = mj0(jiceprt), mj1(jiceprt) 905 872 WRITE(numout,*) ' vice_init : ', vice_init(ji,jj) 906 873 WRITE(numout,*) ' eice_init : ', eice_init(ji,jj) … … 912 879 ! 1) Compute change in open water area due to closing and opening. 913 880 !------------------------------------------------------------------------------- 914 915 neg_ato_i = .false.916 917 881 DO jj = 1, jpj 918 882 DO ji = 1, jpi 919 883 ato_i(ji,jj) = ato_i(ji,jj) - athorn(ji,jj,0) * closing_gross(ji,jj) * rdt_ice & 920 884 & + opning(ji,jj) * rdt_ice 921 IF ( ato_i(ji,jj) < -epsi10 ) THEN922 neg_ato_i = .TRUE.923 ELSEIF( ato_i(ji,jj) < 0._wp ) THEN ! roundoff error885 IF ( ato_i(ji,jj) < -epsi10 ) THEN ! there is a bug 886 IF(lwp) WRITE(numout,*) 'Ridging error: ato_i < 0 -- ato_i : ',ato_i(ji,jj) 887 ELSEIF( ato_i(ji,jj) < 0._wp ) THEN ! roundoff error 924 888 ato_i(ji,jj) = 0._wp 925 889 ENDIF 926 END DO !jj 927 END DO !ji 928 929 ! if negative open water area alert it 930 IF( neg_ato_i ) THEN ! there is a bug 931 DO jj = 1, jpj 932 DO ji = 1, jpi 933 IF( ato_i(ji,jj) < -epsi10 ) THEN 934 WRITE(numout,*) '' 935 WRITE(numout,*) 'Ridging error: ato_i < 0' 936 WRITE(numout,*) 'ato_i : ', ato_i(ji,jj) 937 ENDIF ! ato_i < -epsi10 938 END DO 939 END DO 940 ENDIF 890 END DO 891 END DO 941 892 942 893 !----------------------------------------------------------------- 943 894 ! 2) Save initial state variables 944 895 !----------------------------------------------------------------- 945 946 DO jl = 1, jpl 947 aicen_init(:,:,jl) = a_i(:,:,jl) 948 vicen_init(:,:,jl) = v_i(:,:,jl) 949 vsnwn_init(:,:,jl) = v_s(:,:,jl) 950 ! 951 smv_i_init(:,:,jl) = smv_i(:,:,jl) 952 oa_i_init (:,:,jl) = oa_i (:,:,jl) 953 END DO !jl 954 955 esnwn_init(:,:,:) = e_s(:,:,1,:) 956 957 DO jl = 1, jpl 958 DO jk = 1, nlay_i 959 eicen_init(:,:,jk,jl) = e_i(:,:,jk,jl) 960 END DO 961 END DO 896 aicen_init(:,:,:) = a_i (:,:,:) 897 vicen_init(:,:,:) = v_i (:,:,:) 898 vsnwn_init(:,:,:) = v_s (:,:,:) 899 smv_i_init(:,:,:) = smv_i(:,:,:) 900 esnwn_init(:,:,:) = e_s (:,:,1,:) 901 eicen_init(:,:,:,:) = e_i (:,:,:,:) 902 oa_i_init (:,:,:) = oa_i (:,:,:) 962 903 963 904 ! … … 982 923 indxi(icells) = ji 983 924 indxj(icells) = jj 984 ENDIF ! test on a_icen_init 985 END DO ! ji 986 END DO ! jj 987 988 large_afrac = .false. 989 large_afrft = .false. 990 991 !CDIR NODEP 925 ENDIF 926 END DO 927 END DO 928 992 929 DO ij = 1, icells 993 930 ji = indxi(ij) … … 1003 940 arft2(ji,jj) = arft1(ji,jj) / kraft 1004 941 1005 oirdg1(ji,jj)= aridge(ji,jj,jl1)*closing_gross(ji,jj)*rdt_ice1006 oirft1(ji,jj)= araft (ji,jj,jl1)*closing_gross(ji,jj)*rdt_ice1007 oirdg2(ji,jj)= oirdg1(ji,jj) / krdg(ji,jj,jl1)1008 oirft2(ji,jj)= oirft1(ji,jj) / kraft1009 1010 942 !--------------------------------------------------------------- 1011 943 ! 3.3) Compute ridging /rafting fractions, make sure afrac <=1 … … 1015 947 afrft(ji,jj) = arft1(ji,jj) / aicen_init(ji,jj,jl1) !rafting 1016 948 1017 IF (afrac(ji,jj) > kamax + epsi10) THEN !riging1018 large_afrac = .true.1019 ELSEIF (afrac(ji,jj) > kamax) THEN! roundoff error949 IF( afrac(ji,jj) > kamax + epsi10 ) THEN ! there is a bug 950 IF(lwp) WRITE(numout,*) ' ardg > a_i -- ardg, aicen_init : ', ardg1(ji,jj), aicen_init(ji,jj,jl1) 951 ELSEIF( afrac(ji,jj) > kamax ) THEN ! roundoff error 1020 952 afrac(ji,jj) = kamax 1021 953 ENDIF 1022 IF (afrft(ji,jj) > kamax + epsi10) THEN !rafting 1023 large_afrft = .true. 1024 ELSEIF (afrft(ji,jj) > kamax) THEN ! roundoff error 954 955 IF( afrft(ji,jj) > kamax + epsi10 ) THEN ! there is a bug 956 IF(lwp) WRITE(numout,*) ' arft > a_i -- arft, aicen_init : ', arft1(ji,jj), aicen_init(ji,jj,jl1) 957 ELSEIF( afrft(ji,jj) > kamax) THEN ! roundoff error 1025 958 afrft(ji,jj) = kamax 1026 959 ENDIF … … 1031 964 !-------------------------------------------------------------------------- 1032 965 vrdg1(ji,jj) = vicen_init(ji,jj,jl1) * afrac(ji,jj) 1033 vrdg2(ji,jj) = vrdg1(ji,jj) * ( 1. + ridge_por ) 1034 vsw (ji,jj) = vrdg1(ji,jj) * ridge_por 1035 1036 vsrdg(ji,jj) = vsnwn_init(ji,jj,jl1) * afrac(ji,jj) 1037 esrdg(ji,jj) = esnwn_init(ji,jj,jl1) * afrac(ji,jj) 1038 srdg1(ji,jj) = smv_i_init(ji,jj,jl1) * afrac(ji,jj) 1039 srdg2(ji,jj) = smv_i_init(ji,jj,jl1) * afrac(ji,jj) !! MV HC 2014 this line seems useless 966 vrdg2(ji,jj) = vrdg1(ji,jj) * ( 1. + rn_por_rdg ) 967 vsw (ji,jj) = vrdg1(ji,jj) * rn_por_rdg 968 969 vsrdg (ji,jj) = vsnwn_init(ji,jj,jl1) * afrac(ji,jj) 970 esrdg (ji,jj) = esnwn_init(ji,jj,jl1) * afrac(ji,jj) 971 srdg1 (ji,jj) = smv_i_init(ji,jj,jl1) * afrac(ji,jj) 972 oirdg1(ji,jj) = oa_i_init (ji,jj,jl1) * afrac(ji,jj) 973 oirdg2(ji,jj) = oa_i_init (ji,jj,jl1) * afrac(ji,jj) / krdg(ji,jj,jl1) 1040 974 1041 975 ! rafting volumes, heat contents ... 1042 virft(ji,jj) = vicen_init(ji,jj,jl1) * afrft(ji,jj) 1043 vsrft(ji,jj) = vsnwn_init(ji,jj,jl1) * afrft(ji,jj) 1044 esrft(ji,jj) = esnwn_init(ji,jj,jl1) * afrft(ji,jj) 1045 smrft(ji,jj) = smv_i_init(ji,jj,jl1) * afrft(ji,jj) 976 virft (ji,jj) = vicen_init(ji,jj,jl1) * afrft(ji,jj) 977 vsrft (ji,jj) = vsnwn_init(ji,jj,jl1) * afrft(ji,jj) 978 esrft (ji,jj) = esnwn_init(ji,jj,jl1) * afrft(ji,jj) 979 smrft (ji,jj) = smv_i_init(ji,jj,jl1) * afrft(ji,jj) 980 oirft1(ji,jj) = oa_i_init (ji,jj,jl1) * afrft(ji,jj) 981 oirft2(ji,jj) = oa_i_init (ji,jj,jl1) * afrft(ji,jj) / kraft 1046 982 1047 983 ! substract everything 1048 a_i(ji,jj,jl1) = a_i(ji,jj,jl1) - ardg1(ji,jj) - arft1(ji,jj) 1049 v_i(ji,jj,jl1) = v_i(ji,jj,jl1) - vrdg1(ji,jj) - virft(ji,jj) 1050 v_s(ji,jj,jl1) = v_s(ji,jj,jl1) - vsrdg(ji,jj) - vsrft(ji,jj) 1051 e_s(ji,jj,1,jl1) = e_s(ji,jj,1,jl1) - esrdg(ji,jj) - esrft(ji,jj) 984 a_i(ji,jj,jl1) = a_i(ji,jj,jl1) - ardg1 (ji,jj) - arft1 (ji,jj) 985 v_i(ji,jj,jl1) = v_i(ji,jj,jl1) - vrdg1 (ji,jj) - virft (ji,jj) 986 v_s(ji,jj,jl1) = v_s(ji,jj,jl1) - vsrdg (ji,jj) - vsrft (ji,jj) 987 e_s(ji,jj,1,jl1) = e_s(ji,jj,1,jl1) - esrdg (ji,jj) - esrft (ji,jj) 988 smv_i(ji,jj,jl1) = smv_i(ji,jj,jl1) - srdg1 (ji,jj) - smrft (ji,jj) 1052 989 oa_i(ji,jj,jl1) = oa_i(ji,jj,jl1) - oirdg1(ji,jj) - oirft1(ji,jj) 1053 smv_i(ji,jj,jl1) = smv_i(ji,jj,jl1) - srdg1(ji,jj) - smrft(ji,jj)1054 990 1055 991 !----------------------------------------------------------------- 1056 992 ! 3.5) Compute properties of new ridges 1057 993 !----------------------------------------------------------------- 1058 !--------- ----994 !--------- 1059 995 ! Salinity 1060 !--------- ----996 !--------- 1061 997 smsw(ji,jj) = vsw(ji,jj) * sss_m(ji,jj) ! salt content of seawater frozen in voids !! MV HC2014 1062 998 srdg2(ji,jj) = srdg1(ji,jj) + smsw(ji,jj) ! salt content of new ridge 1063 999 1064 !srdg2(ji,jj) = MIN( s_i_max * vrdg2(ji,jj) , zsrdg2 ) ! impose a maximum salinity1000 !srdg2(ji,jj) = MIN( rn_simax * vrdg2(ji,jj) , zsrdg2 ) ! impose a maximum salinity 1065 1001 1066 1002 sfx_dyn(ji,jj) = sfx_dyn(ji,jj) - smsw(ji,jj) * rhoic * r1_rdtice 1067 wfx_dyn(ji,jj) = wfx_dyn(ji,jj) - vsw (ji,jj) * rhoic * r1_rdtice ! gurvan:increase in ice volume du to seawater frozen in voids1003 wfx_dyn(ji,jj) = wfx_dyn(ji,jj) - vsw (ji,jj) * rhoic * r1_rdtice ! increase in ice volume du to seawater frozen in voids 1068 1004 1069 1005 !------------------------------------ … … 1091 1027 ! ij looping 1-icells 1092 1028 1093 msnow_mlt(ji,jj) = msnow_mlt(ji,jj) + rhosn*vsrdg(ji,jj)*(1.0- fsnowrdg) & ! rafting included1094 & + rhosn*vsrft(ji,jj)*(1.0- fsnowrft)1095 1096 ! in 1e-9 Joules(same as e_s)1097 esnow_mlt(ji,jj) = esnow_mlt(ji,jj) - esrdg(ji,jj)*(1.0- fsnowrdg) & !rafting included1098 & - esrft(ji,jj)*(1.0- fsnowrft)1029 msnow_mlt(ji,jj) = msnow_mlt(ji,jj) + rhosn*vsrdg(ji,jj)*(1.0-rn_fsnowrdg) & ! rafting included 1030 & + rhosn*vsrft(ji,jj)*(1.0-rn_fsnowrft) 1031 1032 ! in J/m2 (same as e_s) 1033 esnow_mlt(ji,jj) = esnow_mlt(ji,jj) - esrdg(ji,jj)*(1.0-rn_fsnowrdg) & !rafting included 1034 & - esrft(ji,jj)*(1.0-rn_fsnowrft) 1099 1035 1100 1036 !----------------------------------------------------------------- … … 1109 1045 dhr2(ji,jj) = hrmax(ji,jj,jl1) * hrmax(ji,jj,jl1) - hrmin(ji,jj,jl1) * hrmin(ji,jj,jl1) 1110 1046 1111 END DO ! ij1047 END DO 1112 1048 1113 1049 !-------------------------------------------------------------------- … … 1116 1052 !-------------------------------------------------------------------- 1117 1053 DO jk = 1, nlay_i 1118 !CDIR NODEP1119 1054 DO ij = 1, icells 1120 1055 ji = indxi(ij) … … 1128 1063 ! enthalpy of the trapped seawater (J/m2, >0) 1129 1064 ! clem: if sst>0, then ersw <0 (is that possible?) 1130 zsstK = sst_m(ji,jj) + rt0 1131 ersw(ji,jj,jk) = - rhoic * vsw(ji,jj) * rcp * ( zsstK - rt0 ) / REAL( nlay_i ) 1065 ersw(ji,jj,jk) = - rhoic * vsw(ji,jj) * rcp * sst_m(ji,jj) * r1_nlay_i 1132 1066 1133 1067 ! heat flux to the ocean 1134 1068 hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + ersw(ji,jj,jk) * r1_rdtice ! > 0 [W.m-2] ocean->ice flux 1135 1069 1136 ! Correct dimensions to avoid big values 1137 ersw(ji,jj,jk) = ersw(ji,jj,jk) / unit_fac 1138 1139 ! Mutliply by ice volume, and divide by number of layers to get heat content in 1.e9 J 1140 ! it is added to sea ice because the sign convention is the opposite of the sign convention for the ocean 1141 !! MV HC 2014 1142 ersw (ji,jj,jk) = ersw(ji,jj,jk) * area(ji,jj) 1143 1070 ! it is added to sea ice because the sign convention is the opposite of the sign convention for the ocean 1144 1071 erdg2(ji,jj,jk) = erdg1(ji,jj,jk) + ersw(ji,jj,jk) 1145 1072 1146 END DO ! ij1147 END DO !jk1073 END DO 1074 END DO 1148 1075 1149 1076 1150 1077 IF( con_i ) THEN 1151 1078 DO jk = 1, nlay_i 1152 !CDIR NODEP1153 1079 DO ij = 1, icells 1154 1080 ji = indxi(ij) 1155 1081 jj = indxj(ij) 1156 1082 eice_init(ji,jj) = eice_init(ji,jj) + erdg2(ji,jj,jk) - erdg1(ji,jj,jk) 1157 END DO ! ij 1158 END DO !jk 1159 ENDIF 1160 1161 IF( large_afrac ) THEN ! there is a bug 1162 !CDIR NODEP 1163 DO ij = 1, icells 1164 ji = indxi(ij) 1165 jj = indxj(ij) 1166 IF( afrac(ji,jj) > kamax + epsi10 ) THEN 1167 WRITE(numout,*) '' 1168 WRITE(numout,*) ' ardg > a_i' 1169 WRITE(numout,*) ' ardg, aicen_init : ', ardg1(ji,jj), aicen_init(ji,jj,jl1) 1170 ENDIF 1171 END DO 1172 ENDIF 1173 IF( large_afrft ) THEN ! there is a bug 1174 !CDIR NODEP 1175 DO ij = 1, icells 1176 ji = indxi(ij) 1177 jj = indxj(ij) 1178 IF( afrft(ji,jj) > kamax + epsi10 ) THEN 1179 WRITE(numout,*) '' 1180 WRITE(numout,*) ' arft > a_i' 1181 WRITE(numout,*) ' arft, aicen_init : ', arft1(ji,jj), aicen_init(ji,jj,jl1) 1182 ENDIF 1083 END DO 1183 1084 END DO 1184 1085 ENDIF … … 1190 1091 DO jl2 = 1, jpl 1191 1092 ! over categories to which ridged ice is transferred 1192 !CDIR NODEP1193 1093 DO ij = 1, icells 1194 1094 ji = indxi(ij) … … 1199 1099 ! Transfer area, volume, and energy accordingly. 1200 1100 1201 IF( hrmin(ji,jj,jl1) >= hi_max(jl2) .OR. & 1202 hrmax(ji,jj,jl1) <= hi_max(jl2-1) ) THEN 1101 IF( hrmin(ji,jj,jl1) >= hi_max(jl2) .OR. hrmax(ji,jj,jl1) <= hi_max(jl2-1) ) THEN 1203 1102 hL = 0._wp 1204 1103 hR = 0._wp … … 1214 1113 a_i (ji,jj ,jl2) = a_i (ji,jj ,jl2) + ardg2 (ji,jj) * farea 1215 1114 v_i (ji,jj ,jl2) = v_i (ji,jj ,jl2) + vrdg2 (ji,jj) * fvol(ji,jj) 1216 v_s (ji,jj ,jl2) = v_s (ji,jj ,jl2) + vsrdg (ji,jj) * fvol(ji,jj) * fsnowrdg1217 e_s (ji,jj,1,jl2) = e_s (ji,jj,1,jl2) + esrdg (ji,jj) * fvol(ji,jj) * fsnowrdg1115 v_s (ji,jj ,jl2) = v_s (ji,jj ,jl2) + vsrdg (ji,jj) * fvol(ji,jj) * rn_fsnowrdg 1116 e_s (ji,jj,1,jl2) = e_s (ji,jj,1,jl2) + esrdg (ji,jj) * fvol(ji,jj) * rn_fsnowrdg 1218 1117 smv_i(ji,jj ,jl2) = smv_i(ji,jj ,jl2) + srdg2 (ji,jj) * fvol(ji,jj) 1219 1118 oa_i (ji,jj ,jl2) = oa_i (ji,jj ,jl2) + oirdg2(ji,jj) * farea 1220 1119 1221 END DO ! ij1120 END DO 1222 1121 1223 1122 ! Transfer ice energy to category jl2 by ridging 1224 1123 DO jk = 1, nlay_i 1225 !CDIR NODEP1226 1124 DO ij = 1, icells 1227 1125 ji = indxi(ij) 1228 1126 jj = indxj(ij) 1229 e_i(ji,jj,jk,jl2) = e_i(ji,jj,jk,jl2) + fvol(ji,jj) *erdg2(ji,jj,jk)1127 e_i(ji,jj,jk,jl2) = e_i(ji,jj,jk,jl2) + fvol(ji,jj) * erdg2(ji,jj,jk) 1230 1128 END DO 1231 1129 END DO … … 1235 1133 DO jl2 = 1, jpl 1236 1134 1237 !CDIR NODEP1238 1135 DO ij = 1, icells 1239 1136 ji = indxi(ij) … … 1242 1139 ! thickness category jl2, transfer area, volume, and energy accordingly. 1243 1140 ! 1244 IF( hraft(ji,jj,jl1) <= hi_max(jl2) .AND. & 1245 hraft(ji,jj,jl1) > hi_max(jl2-1) ) THEN 1141 IF( hraft(ji,jj,jl1) <= hi_max(jl2) .AND. hraft(ji,jj,jl1) > hi_max(jl2-1) ) THEN 1246 1142 a_i (ji,jj ,jl2) = a_i (ji,jj ,jl2) + arft2 (ji,jj) 1247 1143 v_i (ji,jj ,jl2) = v_i (ji,jj ,jl2) + virft (ji,jj) 1248 v_s (ji,jj ,jl2) = v_s (ji,jj ,jl2) + vsrft (ji,jj) * fsnowrft1249 e_s (ji,jj,1,jl2) = e_s (ji,jj,1,jl2) + esrft (ji,jj) * fsnowrft1144 v_s (ji,jj ,jl2) = v_s (ji,jj ,jl2) + vsrft (ji,jj) * rn_fsnowrft 1145 e_s (ji,jj,1,jl2) = e_s (ji,jj,1,jl2) + esrft (ji,jj) * rn_fsnowrft 1250 1146 smv_i(ji,jj ,jl2) = smv_i(ji,jj ,jl2) + smrft (ji,jj) 1251 oa_i (ji,jj ,jl2) = oa_i (ji,jj ,jl2) + oirft2(ji,jj) 1252 ENDIF ! hraft1147 oa_i (ji,jj ,jl2) = oa_i (ji,jj ,jl2) + oirft2(ji,jj) 1148 ENDIF 1253 1149 ! 1254 END DO ! ij1150 END DO 1255 1151 1256 1152 ! Transfer rafted ice energy to category jl2 1257 1153 DO jk = 1, nlay_i 1258 !CDIR NODEP1259 1154 DO ij = 1, icells 1260 1155 ji = indxi(ij) 1261 1156 jj = indxj(ij) 1262 IF( hraft(ji,jj,jl1) <= hi_max(jl2) .AND. & 1263 hraft(ji,jj,jl1) > hi_max(jl2-1) ) THEN 1157 IF( hraft(ji,jj,jl1) <= hi_max(jl2) .AND. hraft(ji,jj,jl1) > hi_max(jl2-1) ) THEN 1264 1158 e_i(ji,jj,jk,jl2) = e_i(ji,jj,jk,jl2) + eirft(ji,jj,jk) 1265 1159 ENDIF 1266 END DO ! ij1267 END DO !jk1268 1269 END DO ! jl21160 END DO 1161 END DO 1162 1163 END DO 1270 1164 1271 1165 END DO ! jl1 (deforming categories) … … 1281 1175 CALL lim_cons_check (eice_init, eice_final, 1.0e-2, fieldid) 1282 1176 1283 DO ji = mi0( jiindx), mi1(jiindx)1284 DO jj = mj0(j jindx), mj1(jjindx)1177 DO ji = mi0(iiceprt), mi1(iiceprt) 1178 DO jj = mj0(jiceprt), mj1(jiceprt) 1285 1179 WRITE(numout,*) ' vice_init : ', vice_init (ji,jj) 1286 1180 WRITE(numout,*) ' vice_final : ', vice_final(ji,jj) … … 1291 1185 ENDIF 1292 1186 ! 1293 CALL wrk_dealloc( (jpi+1)*(jpj+1), indxi, indxj )1294 CALL wrk_dealloc( jpi, jpj, vice_init, vice_final, eice_init, eice_final )1295 CALL wrk_dealloc( jpi, jpj, afrac, fvol , ardg1, ardg2, vsrdg, esrdg, oirdg1, oirdg2, dhr, dhr2 )1296 CALL wrk_dealloc( jpi, jpj, vrdg1, vrdg2, vsw , srdg1, srdg2, smsw)1297 CALL wrk_dealloc( jpi, jpj, afrft, arft1, arft2, virft, vsrft, esrft, smrft, oirft1, oirft2 )1298 CALL wrk_dealloc( jpi, jpj, jpl, aicen_init, vicen_init, vsnwn_init, esnwn_init, smv_i_init, oa_i_init )1299 CALL wrk_dealloc( jpi, jpj, nlay_i +1,eirft, erdg1, erdg2, ersw )1300 CALL wrk_dealloc( jpi, jpj, nlay_i +1, jpl,eicen_init )1187 CALL wrk_dealloc( (jpi+1)*(jpj+1), indxi, indxj ) 1188 CALL wrk_dealloc( jpi, jpj, vice_init, vice_final, eice_init, eice_final ) 1189 CALL wrk_dealloc( jpi, jpj, afrac, fvol , ardg1, ardg2, vsrdg, esrdg, dhr, dhr2 ) 1190 CALL wrk_dealloc( jpi, jpj, vrdg1, vrdg2, vsw , srdg1, srdg2, smsw, oirdg1, oirdg2 ) 1191 CALL wrk_dealloc( jpi, jpj, afrft, arft1, arft2, virft, vsrft, esrft, smrft, oirft1, oirft2 ) 1192 CALL wrk_dealloc( jpi, jpj, jpl, aicen_init, vicen_init, vsnwn_init, esnwn_init, smv_i_init, oa_i_init ) 1193 CALL wrk_dealloc( jpi, jpj, nlay_i, eirft, erdg1, erdg2, ersw ) 1194 CALL wrk_dealloc( jpi, jpj, nlay_i, jpl, eicen_init ) 1301 1195 ! 1302 1196 END SUBROUTINE lim_itd_me_ridgeshift 1303 1304 1305 SUBROUTINE lim_itd_me_asumr1306 !!-----------------------------------------------------------------------------1307 !! *** ROUTINE lim_itd_me_asumr ***1308 !!1309 !! ** Purpose : finds total fractional area1310 !!1311 !! ** Method : Find the total area of ice plus open water in each grid cell.1312 !! This is similar to the aggregate_area subroutine except that the1313 !! total area can be greater than 1, so the open water area is1314 !! included in the sum instead of being computed as a residual.1315 !!-----------------------------------------------------------------------------1316 INTEGER :: jl ! dummy loop index1317 !!-----------------------------------------------------------------------------1318 !1319 asum(:,:) = ato_i(:,:) ! open water1320 DO jl = 1, jpl ! ice categories1321 asum(:,:) = asum(:,:) + a_i(:,:,jl)1322 END DO1323 !1324 END SUBROUTINE lim_itd_me_asumr1325 1326 1197 1327 1198 SUBROUTINE lim_itd_me_init … … 1339 1210 !!------------------------------------------------------------------- 1340 1211 INTEGER :: ios ! Local integer output status for namelist read 1341 NAMELIST/namiceitdme/ r idge_scheme_swi, Cs, Cf, fsnowrdg,fsnowrft, &1342 & Gstar, astar, Hstar, raft_swi, hparmeter, Craft, ridge_por, &1343 & partfun_swi, brinstren_swi1212 NAMELIST/namiceitdme/ rn_cs, rn_fsnowrdg, rn_fsnowrft, & 1213 & rn_gstar, rn_astar, rn_hstar, ln_rafting, rn_hraft, rn_craft, rn_por_rdg, & 1214 & nn_partfun 1344 1215 !!------------------------------------------------------------------- 1345 1216 ! … … 1357 1228 WRITE(numout,*)' lim_itd_me_init : ice parameters for mechanical ice redistribution ' 1358 1229 WRITE(numout,*)' ~~~~~~~~~~~~~~~' 1359 WRITE(numout,*)' Switch choosing the ice redistribution scheme ridge_scheme_swi', ridge_scheme_swi 1360 WRITE(numout,*)' Fraction of shear energy contributing to ridging Cs ', Cs 1361 WRITE(numout,*)' Ratio of ridging work to PotEner change in ridging Cf ', Cf 1362 WRITE(numout,*)' Fraction of snow volume conserved during ridging fsnowrdg ', fsnowrdg 1363 WRITE(numout,*)' Fraction of snow volume conserved during ridging fsnowrft ', fsnowrft 1364 WRITE(numout,*)' Fraction of total ice coverage contributing to ridging Gstar ', Gstar 1365 WRITE(numout,*)' Equivalent to G* for an exponential part function astar ', astar 1366 WRITE(numout,*)' Quantity playing a role in max ridged ice thickness Hstar ', Hstar 1367 WRITE(numout,*)' Rafting of ice sheets or not raft_swi ', raft_swi 1368 WRITE(numout,*)' Parmeter thickness (threshold between ridge-raft) hparmeter ', hparmeter 1369 WRITE(numout,*)' Rafting hyperbolic tangent coefficient Craft ', Craft 1370 WRITE(numout,*)' Initial porosity of ridges ridge_por ', ridge_por 1371 WRITE(numout,*)' Switch for part. function (0) linear (1) exponential partfun_swi ', partfun_swi 1372 WRITE(numout,*)' Switch for including brine volume in ice strength comp. brinstren_swi ', brinstren_swi 1230 WRITE(numout,*)' Fraction of shear energy contributing to ridging rn_cs = ', rn_cs 1231 WRITE(numout,*)' Fraction of snow volume conserved during ridging rn_fsnowrdg = ', rn_fsnowrdg 1232 WRITE(numout,*)' Fraction of snow volume conserved during ridging rn_fsnowrft = ', rn_fsnowrft 1233 WRITE(numout,*)' Fraction of total ice coverage contributing to ridging rn_gstar = ', rn_gstar 1234 WRITE(numout,*)' Equivalent to G* for an exponential part function rn_astar = ', rn_astar 1235 WRITE(numout,*)' Quantity playing a role in max ridged ice thickness rn_hstar = ', rn_hstar 1236 WRITE(numout,*)' Rafting of ice sheets or not ln_rafting = ', ln_rafting 1237 WRITE(numout,*)' Parmeter thickness (threshold between ridge-raft) rn_hraft = ', rn_hraft 1238 WRITE(numout,*)' Rafting hyperbolic tangent coefficient rn_craft = ', rn_craft 1239 WRITE(numout,*)' Initial porosity of ridges rn_por_rdg = ', rn_por_rdg 1240 WRITE(numout,*)' Switch for part. function (0) linear (1) exponential nn_partfun = ', nn_partfun 1373 1241 ENDIF 1374 1242 ! 1375 1243 END SUBROUTINE lim_itd_me_init 1376 1377 1378 SUBROUTINE lim_itd_me_zapsmall1379 !!-------------------------------------------------------------------1380 !! *** ROUTINE lim_itd_me_zapsmall ***1381 !!1382 !! ** Purpose : Remove too small sea ice areas and correct salt fluxes1383 !!1384 !! history :1385 !! author: William H. Lipscomb, LANL1386 !! Nov 2003: Modified by Julie Schramm to conserve volume and energy1387 !! Sept 2004: Modified by William Lipscomb; replaced normalize_state with1388 !! additions to local freshwater, salt, and heat fluxes1389 !! 9.0, LIM3.0 - 02-2006 (M. Vancoppenolle) original code1390 !!-------------------------------------------------------------------1391 INTEGER :: ji, jj, jl, jk ! dummy loop indices1392 INTEGER :: icells ! number of cells with ice to zap1393 1394 REAL(wp), POINTER, DIMENSION(:,:) :: zmask ! 2D workspace1395 REAL(wp) :: zmask_glo, zsal, zvi, zvs, zei, zes1396 !!gm REAL(wp) :: xtmp ! temporary variable1397 !!-------------------------------------------------------------------1398 1399 CALL wrk_alloc( jpi, jpj, zmask )1400 1401 ! to be sure that at_i is the sum of a_i(jl)1402 at_i(:,:) = SUM( a_i(:,:,:), dim=3 )1403 1404 DO jl = 1, jpl1405 !-----------------------------------------------------------------1406 ! Count categories to be zapped.1407 !-----------------------------------------------------------------1408 icells = 01409 zmask(:,:) = 0._wp1410 DO jj = 1, jpj1411 DO ji = 1, jpi1412 IF( a_i(ji,jj,jl) <= epsi10 .OR. v_i(ji,jj,jl) <= epsi10 .OR. at_i(ji,jj) <= epsi10 ) THEN1413 zmask(ji,jj) = 1._wp1414 ENDIF1415 END DO1416 END DO1417 !zmask_glo = glob_sum(zmask)1418 !IF( ln_nicep .AND. lwp ) WRITE(numout,*) zmask_glo, ' cells of ice zapped in the ocean '1419 1420 !-----------------------------------------------------------------1421 ! Zap ice energy and use ocean heat to melt ice1422 !-----------------------------------------------------------------1423 1424 DO jk = 1, nlay_i1425 DO jj = 1 , jpj1426 DO ji = 1 , jpi1427 zei = e_i(ji,jj,jk,jl)1428 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * ( 1._wp - zmask(ji,jj) )1429 t_i(ji,jj,jk,jl) = t_i(ji,jj,jk,jl) * ( 1._wp - zmask(ji,jj) ) + rtt * zmask(ji,jj)1430 ! update exchanges with ocean1431 hfx_res(ji,jj) = hfx_res(ji,jj) + ( e_i(ji,jj,jk,jl) - zei ) * unit_fac / area(ji,jj) * r1_rdtice ! W.m-2 <01432 END DO1433 END DO1434 END DO1435 1436 DO jj = 1 , jpj1437 DO ji = 1 , jpi1438 1439 zsal = smv_i(ji,jj,jl)1440 zvi = v_i(ji,jj,jl)1441 zvs = v_s(ji,jj,jl)1442 zes = e_s(ji,jj,1,jl)1443 !-----------------------------------------------------------------1444 ! Zap snow energy and use ocean heat to melt snow1445 !-----------------------------------------------------------------1446 ! xtmp = esnon(i,j,n) / dt ! < 01447 ! fhnet(i,j) = fhnet(i,j) + xtmp1448 ! fhnet_hist(i,j) = fhnet_hist(i,j) + xtmp1449 ! xtmp is greater than 01450 ! fluxes are positive to the ocean1451 ! here the flux has to be negative for the ocean1452 t_s(ji,jj,1,jl) = rtt * zmask(ji,jj) + t_s(ji,jj,1,jl) * ( 1._wp - zmask(ji,jj) )1453 1454 !-----------------------------------------------------------------1455 ! zap ice and snow volume, add water and salt to ocean1456 !-----------------------------------------------------------------1457 ato_i(ji,jj) = a_i (ji,jj,jl) * zmask(ji,jj) + ato_i(ji,jj)1458 a_i (ji,jj,jl) = a_i (ji,jj,jl) * ( 1._wp - zmask(ji,jj) )1459 v_i (ji,jj,jl) = v_i (ji,jj,jl) * ( 1._wp - zmask(ji,jj) )1460 v_s (ji,jj,jl) = v_s (ji,jj,jl) * ( 1._wp - zmask(ji,jj) )1461 t_su (ji,jj,jl) = t_su (ji,jj,jl) * ( 1._wp - zmask(ji,jj) ) + t_bo(ji,jj) * zmask(ji,jj)1462 oa_i (ji,jj,jl) = oa_i (ji,jj,jl) * ( 1._wp - zmask(ji,jj) )1463 smv_i(ji,jj,jl) = smv_i(ji,jj,jl) * ( 1._wp - zmask(ji,jj) )1464 e_s(ji,jj,1,jl) = e_s(ji,jj,1,jl) * ( 1._wp - zmask(ji,jj) )1465 ! additional condition1466 IF( v_s(ji,jj,jl) <= epsi10 ) THEN1467 v_s(ji,jj,jl) = 0._wp1468 e_s(ji,jj,1,jl) = 0._wp1469 ENDIF1470 ! update exchanges with ocean1471 sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsal ) * rhoic * r1_rdtice1472 wfx_res(ji,jj) = wfx_res(ji,jj) - ( v_i(ji,jj,jl) - zvi ) * rhoic * r1_rdtice1473 wfx_snw(ji,jj) = wfx_snw(ji,jj) - ( v_s(ji,jj,jl) - zvs ) * rhosn * r1_rdtice1474 hfx_res(ji,jj) = hfx_res(ji,jj) + ( e_s(ji,jj,1,jl) - zes ) * unit_fac / area(ji,jj) * r1_rdtice ! W.m-2 <01475 END DO1476 END DO1477 END DO ! jl1478 1479 ! to be sure that at_i is the sum of a_i(jl)1480 at_i(:,:) = SUM( a_i(:,:,:), dim=3 )1481 !1482 CALL wrk_dealloc( jpi, jpj, zmask )1483 !1484 END SUBROUTINE lim_itd_me_zapsmall1485 1244 1486 1245 #else … … 1493 1252 SUBROUTINE lim_itd_me_icestrength 1494 1253 END SUBROUTINE lim_itd_me_icestrength 1495 SUBROUTINE lim_itd_me_sort1496 END SUBROUTINE lim_itd_me_sort1497 1254 SUBROUTINE lim_itd_me_init 1498 1255 END SUBROUTINE lim_itd_me_init 1499 SUBROUTINE lim_itd_me_zapsmall1500 END SUBROUTINE lim_itd_me_zapsmall1501 1256 #endif 1502 1257 !!====================================================================== -
branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90
r4990 r5682 13 13 !! 'key_lim3' : LIM3 sea-ice model 14 14 !!---------------------------------------------------------------------- 15 !! lim_itd_th : thermodynamics of ice thickness distribution16 15 !! lim_itd_th_rem : 17 16 !! lim_itd_th_reb : … … 25 24 USE thd_ice ! LIM-3 thermodynamic variables 26 25 USE ice ! LIM-3 variables 27 USE par_ice ! LIM-3 parameters28 USE limthd_lac ! LIM-3 lateral accretion29 26 USE limvar ! LIM-3 variables 30 USE limcons ! LIM-3 conservation31 27 USE prtctl ! Print control 32 28 USE in_out_manager ! I/O manager … … 34 30 USE wrk_nemo ! work arrays 35 31 USE lib_fortran ! to use key_nosignedzero 36 USE timing ! Timing 37 USE limcons ! conservation tests 32 USE limcons ! conservation tests 38 33 39 34 IMPLICIT NONE 40 35 PRIVATE 41 36 42 PUBLIC lim_itd_th ! called by ice_stp43 37 PUBLIC lim_itd_th_rem 44 38 PUBLIC lim_itd_th_reb 45 PUBLIC lim_itd_fitline46 PUBLIC lim_itd_shiftice47 39 48 40 !!---------------------------------------------------------------------- … … 53 45 CONTAINS 54 46 55 SUBROUTINE lim_itd_th( kt )56 !!------------------------------------------------------------------57 !! *** ROUTINE lim_itd_th ***58 !!59 !! ** Purpose : computes the thermodynamics of ice thickness distribution60 !!61 !! ** Method :62 !!------------------------------------------------------------------63 INTEGER, INTENT(in) :: kt ! time step index64 !65 INTEGER :: ji, jj, jk, jl ! dummy loop index66 !67 REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b68 !!------------------------------------------------------------------69 IF( nn_timing == 1 ) CALL timing_start('limitd_th')70 71 ! conservation test72 IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limitd_th', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)73 74 IF( kt == nit000 .AND. lwp ) THEN75 WRITE(numout,*)76 WRITE(numout,*) 'lim_itd_th : Thermodynamics of the ice thickness distribution'77 WRITE(numout,*) '~~~~~~~~~~~'78 ENDIF79 80 !------------------------------------------------------------------------------|81 ! 1) Transport of ice between thickness categories. |82 !------------------------------------------------------------------------------|83 ! Given thermodynamic growth rates, transport ice between84 ! thickness categories.85 IF( jpl > 1 ) CALL lim_itd_th_rem( 1, jpl, kt )86 !87 CALL lim_var_glo2eqv ! only for info88 CALL lim_var_agg(1)89 90 !------------------------------------------------------------------------------|91 ! 3) Add frazil ice growing in leads.92 !------------------------------------------------------------------------------|93 CALL lim_thd_lac94 CALL lim_var_glo2eqv ! only for info95 96 IF(ln_ctl) THEN ! Control print97 CALL prt_ctl_info(' ')98 CALL prt_ctl_info(' - Cell values : ')99 CALL prt_ctl_info(' ~~~~~~~~~~~~~ ')100 CALL prt_ctl(tab2d_1=area , clinfo1=' lim_itd_th : cell area :')101 CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_itd_th : at_i :')102 CALL prt_ctl(tab2d_1=vt_i , clinfo1=' lim_itd_th : vt_i :')103 CALL prt_ctl(tab2d_1=vt_s , clinfo1=' lim_itd_th : vt_s :')104 DO jl = 1, jpl105 CALL prt_ctl_info(' ')106 CALL prt_ctl_info(' - Category : ', ivar1=jl)107 CALL prt_ctl_info(' ~~~~~~~~~~')108 CALL prt_ctl(tab2d_1=a_i (:,:,jl) , clinfo1= ' lim_itd_th : a_i : ')109 CALL prt_ctl(tab2d_1=ht_i (:,:,jl) , clinfo1= ' lim_itd_th : ht_i : ')110 CALL prt_ctl(tab2d_1=ht_s (:,:,jl) , clinfo1= ' lim_itd_th : ht_s : ')111 CALL prt_ctl(tab2d_1=v_i (:,:,jl) , clinfo1= ' lim_itd_th : v_i : ')112 CALL prt_ctl(tab2d_1=v_s (:,:,jl) , clinfo1= ' lim_itd_th : v_s : ')113 CALL prt_ctl(tab2d_1=e_s (:,:,1,jl) , clinfo1= ' lim_itd_th : e_s : ')114 CALL prt_ctl(tab2d_1=t_su (:,:,jl) , clinfo1= ' lim_itd_th : t_su : ')115 CALL prt_ctl(tab2d_1=t_s (:,:,1,jl) , clinfo1= ' lim_itd_th : t_snow : ')116 CALL prt_ctl(tab2d_1=sm_i (:,:,jl) , clinfo1= ' lim_itd_th : sm_i : ')117 CALL prt_ctl(tab2d_1=smv_i (:,:,jl) , clinfo1= ' lim_itd_th : smv_i : ')118 DO jk = 1, nlay_i119 CALL prt_ctl_info(' ')120 CALL prt_ctl_info(' - Layer : ', ivar1=jk)121 CALL prt_ctl_info(' ~~~~~~~')122 CALL prt_ctl(tab2d_1=t_i(:,:,jk,jl) , clinfo1= ' lim_itd_th : t_i : ')123 CALL prt_ctl(tab2d_1=e_i(:,:,jk,jl) , clinfo1= ' lim_itd_th : e_i : ')124 END DO125 END DO126 ENDIF127 !128 ! conservation test129 IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limitd_th', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)130 !131 IF( nn_timing == 1 ) CALL timing_stop('limitd_th')132 END SUBROUTINE lim_itd_th133 !134 135 47 SUBROUTINE lim_itd_th_rem( klbnd, kubnd, kt ) 136 48 !!------------------------------------------------------------------ … … 153 65 REAL(wp) :: zx1, zwk1, zdh0, zetamin, zdamax ! local scalars 154 66 REAL(wp) :: zx2, zwk2, zda0, zetamax ! - - 155 REAL(wp) :: zx3 , zareamin ! - -67 REAL(wp) :: zx3 156 68 CHARACTER (len = 15) :: fieldid 157 69 … … 179 91 !!------------------------------------------------------------------ 180 92 181 CALL wrk_alloc( jpi,jpj, zremap_flag ) ! integer182 CALL wrk_alloc( jpi,jpj,jpl-1, zdonor ) ! integer93 CALL wrk_alloc( jpi,jpj, zremap_flag ) 94 CALL wrk_alloc( jpi,jpj,jpl-1, zdonor ) 183 95 CALL wrk_alloc( jpi,jpj,jpl, zdhice, g0, g1, hL, hR, zht_i_b, dummy_es ) 184 96 CALL wrk_alloc( jpi,jpj,jpl-1, zdaice, zdvice ) 185 97 CALL wrk_alloc( jpi,jpj,jpl+1, zhbnew, kkstart = 0 ) 186 98 CALL wrk_alloc( (jpi+1)*(jpj+1), zvetamin, zvetamax ) 187 CALL wrk_alloc( (jpi+1)*(jpj+1), nind_i, nind_j ) ! integer99 CALL wrk_alloc( (jpi+1)*(jpj+1), nind_i, nind_j ) 188 100 CALL wrk_alloc( jpi,jpj, zhb0,zhb1,vt_i_init,vt_i_final,vt_s_init,vt_s_final,et_i_init,et_i_final,et_s_init,et_s_final ) 189 190 zareamin = epsi10 !minimum area in thickness categories tolerated by the conceptors of the model191 101 192 102 !!---------------------------------------------------------------------------------------------- … … 216 126 DO jj = 1, jpj 217 127 DO ji = 1, jpi 218 rswitch = 1.0 - MAX( 0.0, SIGN( 1.0, - a_i(ji,jj,jl) +epsi10 ) ) !0 if no ice and 1 if yes128 rswitch = MAX( 0.0, SIGN( 1.0, a_i(ji,jj,jl) - epsi10 ) ) !0 if no ice and 1 if yes 219 129 ht_i(ji,jj,jl) = v_i(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi10 ) * rswitch 220 rswitch = 1.0 - MAX( 0.0, SIGN( 1.0, - a_i_b(ji,jj,jl) + epsi10) ) !0 if no ice and 1 if yes130 rswitch = MAX( 0.0, SIGN( 1.0, a_i_b(ji,jj,jl) - epsi10) ) 221 131 zht_i_b(ji,jj,jl) = v_i_b(ji,jj,jl) / MAX( a_i_b(ji,jj,jl), epsi10 ) * rswitch 222 IF( a_i(ji,jj,jl) > epsi10 ) zdhice(ji,jj,jl) = ht_i(ji,jj,jl) - zht_i_b(ji,jj,jl) 132 IF( a_i(ji,jj,jl) > epsi10 ) zdhice(ji,jj,jl) = ht_i(ji,jj,jl) - zht_i_b(ji,jj,jl) ! clem: useless IF statement? 223 133 END DO 224 134 END DO … … 239 149 DO jj = 1, jpj 240 150 DO ji = 1, jpi 241 IF ( at_i(ji,jj) .gt. zareamin) THEN151 IF ( at_i(ji,jj) > epsi10 ) THEN 242 152 nbrem = nbrem + 1 243 153 nind_i(nbrem) = ji … … 247 157 zremap_flag(ji,jj) = 0 248 158 ENDIF 249 END DO !ji250 END DO !jj159 END DO 160 END DO 251 161 252 162 !----------------------------------------------------------------------------------------------- … … 254 164 !----------------------------------------------------------------------------------------------- 255 165 !- 4.1 Compute category boundaries 256 ! Tricky trick see limitd_me.F90257 ! will be soon removed, CT258 ! hi_max(kubnd) = 99.259 166 zhbnew(:,:,:) = 0._wp 260 167 … … 265 172 ! 266 173 zhbnew(ii,ij,jl) = hi_max(jl) 267 IF ( a_i_b(ii,ij,jl) > epsi10 .AND. a_i_b(ii,ij,jl+1) > epsi10 ) THEN174 IF ( a_i_b(ii,ij,jl) > epsi10 .AND. a_i_b(ii,ij,jl+1) > epsi10 ) THEN 268 175 !interpolate between adjacent category growth rates 269 176 zslope = ( zdhice(ii,ij,jl+1) - zdhice(ii,ij,jl) ) / ( zht_i_b(ii,ij,jl+1) - zht_i_b(ii,ij,jl) ) 270 177 zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl) + zslope * ( hi_max(jl) - zht_i_b(ii,ij,jl) ) 271 ELSEIF 178 ELSEIF( a_i_b(ii,ij,jl) > epsi10) THEN 272 179 zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl) 273 ELSEIF 180 ELSEIF( a_i_b(ii,ij,jl+1) > epsi10) THEN 274 181 zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl+1) 275 182 ENDIF … … 280 187 ii = nind_i(ji) 281 188 ij = nind_j(ji) 282 IF( a_i(ii,ij,jl) > epsi10 .AND. ht_i(ii,ij,jl) >= zhbnew(ii,ij,jl) ) THEN 189 190 ! clem: we do not want ht_i to be too close to either HR or HL otherwise a division by nearly 0 is possible 191 ! in lim_itd_fitline in the case (HR-HL) = 3(Hice - HL) or = 3(HR - Hice) 192 IF ( a_i(ii,ij,jl ) > epsi10 .AND. ht_i(ii,ij,jl ) > ( zhbnew(ii,ij,jl) - epsi10 ) ) THEN 283 193 zremap_flag(ii,ij) = 0 284 ELSEIF( a_i(ii,ij,jl+1) > epsi10 .AND. ht_i(ii,ij,jl+1) < = zhbnew(ii,ij,jl) ) THEN194 ELSEIF( a_i(ii,ij,jl+1) > epsi10 .AND. ht_i(ii,ij,jl+1) < ( zhbnew(ii,ij,jl) + epsi10 ) ) THEN 285 195 zremap_flag(ii,ij) = 0 286 196 ENDIF 287 197 288 198 !- 4.3 Check that each zhbnew does not exceed maximal values hi_max 199 IF( zhbnew(ii,ij,jl) < hi_max(jl-1) ) zremap_flag(ii,ij) = 0 289 200 IF( zhbnew(ii,ij,jl) > hi_max(jl+1) ) zremap_flag(ii,ij) = 0 290 IF( zhbnew(ii,ij,jl) < hi_max(jl-1) ) zremap_flag(ii,ij) = 0 291 END DO 292 293 END DO !jl 201 ! clem bug: why is not the following instead? 202 !!IF( zhbnew(ii,ij,jl) < hi_max(jl-1) ) zremap_flag(ii,ij) = 0 203 !!IF( zhbnew(ii,ij,jl) > hi_max(jl ) ) zremap_flag(ii,ij) = 0 204 205 END DO 206 207 END DO 294 208 295 209 !----------------------------------------------------------------------------------------------- … … 312 226 DO jj = 1, jpj 313 227 DO ji = 1, jpi 314 zhb0(ji,jj) = hi_max(0) ! 0eme 315 zhb1(ji,jj) = hi_max(1) ! 1er 316 317 zhbnew(ji,jj,klbnd-1) = 0._wp 228 zhb0(ji,jj) = hi_max(0) 229 zhb1(ji,jj) = hi_max(1) 318 230 319 231 IF( a_i(ji,jj,kubnd) > epsi10 ) THEN 320 zhbnew(ji,jj,kubnd) = 3._wp * ht_i(ji,jj,kubnd) - 2._wp * zhbnew(ji,jj,kubnd-1)232 zhbnew(ji,jj,kubnd) = MAX( hi_max(kubnd-1), 3._wp * ht_i(ji,jj,kubnd) - 2._wp * zhbnew(ji,jj,kubnd-1) ) 321 233 ELSE 322 zhbnew(ji,jj,kubnd) = hi_max(kubnd) 323 !!? clem bug: since hi_max(jpl)=99, this limit is very high 324 !!? but I think it is erased in fitline subroutine 325 ENDIF 326 327 IF( zhbnew(ji,jj,kubnd) < hi_max(kubnd-1) ) zhbnew(ji,jj,kubnd) = hi_max(kubnd-1) 328 329 END DO !jj 330 END DO !jj 234 !clem bug zhbnew(ji,jj,kubnd) = hi_max(kubnd) 235 zhbnew(ji,jj,kubnd) = hi_max(kubnd-1) ! not used anyway 236 ENDIF 237 238 ! clem: we do not want ht_i_b to be too close to either HR or HL otherwise a division by nearly 0 is possible 239 ! in lim_itd_fitline in the case (HR-HL) = 3(Hice - HL) or = 3(HR - Hice) 240 IF ( zht_i_b(ji,jj,klbnd) < ( zhb0(ji,jj) + epsi10 ) ) THEN 241 zremap_flag(ji,jj) = 0 242 ELSEIF( zht_i_b(ji,jj,klbnd) > ( zhb1(ji,jj) - epsi10 ) ) THEN 243 zremap_flag(ji,jj) = 0 244 ENDIF 245 246 END DO 247 END DO 331 248 332 249 !----------------------------------------------------------------------------------------------- … … 334 251 !----------------------------------------------------------------------------------------------- 335 252 !- 7.1 g(h) for category 1 at start of time step 336 CALL lim_itd_fitline( klbnd, zhb0, zhb1, zht_i_b(:,:,klbnd), & 337 & g0(:,:,klbnd), g1(:,:,klbnd), hL(:,:,klbnd), & 253 CALL lim_itd_fitline( klbnd, zhb0, zhb1, zht_i_b(:,:,klbnd), g0(:,:,klbnd), g1(:,:,klbnd), hL(:,:,klbnd), & 338 254 & hR(:,:,klbnd), zremap_flag ) 339 255 … … 343 259 ij = nind_j(ji) 344 260 345 !ji346 IF (a_i(ii,ij,klbnd) .gt. epsi10) THEN 261 IF( a_i(ii,ij,klbnd) > epsi10 ) THEN 262 347 263 zdh0 = zdhice(ii,ij,klbnd) !decrease of ice thickness in the lower category 348 ! ji, a_i > epsi10 349 IF (zdh0 .lt. 0.0) THEN !remove area from category 1 350 ! ji, a_i > epsi10; zdh0 < 0 351 zdh0 = MIN(-zdh0,hi_max(klbnd)) 352 264 265 IF( zdh0 < 0.0 ) THEN !remove area from category 1 266 zdh0 = MIN( -zdh0, hi_max(klbnd) ) 353 267 !Integrate g(1) from 0 to dh0 to estimate area melted 354 zetamax = MIN(zdh0,hR(ii,ij,klbnd)) - hL(ii,ij,klbnd) 355 IF (zetamax.gt.0.0) THEN 356 zx1 = zetamax 357 zx2 = 0.5 * zetamax*zetamax 358 zda0 = g1(ii,ij,klbnd) * zx2 + g0(ii,ij,klbnd) * zx1 !ice area removed 359 ! Constrain new thickness <= ht_i 360 zdamax = a_i(ii,ij,klbnd) * & 361 (1.0 - ht_i(ii,ij,klbnd)/zht_i_b(ii,ij,klbnd)) ! zdamax > 0 362 !ice area lost due to melting of thin ice 363 zda0 = MIN(zda0, zdamax) 364 268 zetamax = MIN( zdh0, hR(ii,ij,klbnd) ) - hL(ii,ij,klbnd) 269 270 IF( zetamax > 0.0 ) THEN 271 zx1 = zetamax 272 zx2 = 0.5 * zetamax * zetamax 273 zda0 = g1(ii,ij,klbnd) * zx2 + g0(ii,ij,klbnd) * zx1 ! ice area removed 274 zdamax = a_i(ii,ij,klbnd) * (1.0 - ht_i(ii,ij,klbnd) / zht_i_b(ii,ij,klbnd) ) ! Constrain new thickness <= ht_i 275 zda0 = MIN( zda0, zdamax ) ! ice area lost due to melting 276 ! of thin ice (zdamax > 0) 365 277 ! Remove area, conserving volume 366 ht_i(ii,ij,klbnd) = ht_i(ii,ij,klbnd) & 367 * a_i(ii,ij,klbnd) / ( a_i(ii,ij,klbnd) - zda0 ) 278 ht_i(ii,ij,klbnd) = ht_i(ii,ij,klbnd) * a_i(ii,ij,klbnd) / ( a_i(ii,ij,klbnd) - zda0 ) 368 279 a_i(ii,ij,klbnd) = a_i(ii,ij,klbnd) - zda0 369 v_i(ii,ij,klbnd) = a_i(ii,ij,klbnd)*ht_i(ii,ij,klbnd) ! clem-useless ? 370 ENDIF ! zetamax > 0 371 ! ji, a_i > epsi10 372 373 ELSE ! if ice accretion 374 ! ji, a_i > epsi10; zdh0 > 0 375 zhbnew(ii,ij,klbnd-1) = MIN(zdh0,hi_max(klbnd)) 376 ! zhbnew was 0, and is shifted to the right to account for thin ice 377 ! growth in openwater (F0 = f1) 378 ENDIF ! zdh0 379 380 ! a_i > epsi10 381 ENDIF ! a_i > epsi10 382 383 END DO ! ji 280 v_i(ii,ij,klbnd) = a_i(ii,ij,klbnd) * ht_i(ii,ij,klbnd) ! clem-useless ? 281 ENDIF 282 283 ELSE ! if ice accretion zdh0 > 0 284 ! zhbnew was 0, and is shifted to the right to account for thin ice growth in openwater (F0 = f1) 285 zhbnew(ii,ij,klbnd-1) = MIN( zdh0, hi_max(klbnd) ) 286 ENDIF 287 288 ENDIF 289 290 END DO 384 291 385 292 !- 7.3 g(h) for each thickness category 386 293 DO jl = klbnd, kubnd 387 CALL lim_itd_fitline( jl, zhbnew(:,:,jl-1), zhbnew(:,:,jl), ht_i(:,:,jl),&388 g0(:,:,jl), g1(:,:,jl), hL(:,:,jl), hR(:,:,jl), zremap_flag)294 CALL lim_itd_fitline( jl, zhbnew(:,:,jl-1), zhbnew(:,:,jl), ht_i(:,:,jl), & 295 & g0(:,:,jl), g1(:,:,jl), hL(:,:,jl), hR(:,:,jl), zremap_flag ) 389 296 END DO 390 297 … … 406 313 ij = nind_j(ji) 407 314 408 IF (zhbnew(ii,ij,jl) .gt. hi_max(jl)) THEN ! transfer from jl to jl+1 409 315 IF (zhbnew(ii,ij,jl) > hi_max(jl)) THEN ! transfer from jl to jl+1 410 316 ! left and right integration limits in eta space 411 zvetamin(ji) = MAX( hi_max(jl), hL(ii,ij,jl)) - hL(ii,ij,jl)412 zvetamax(ji) = MIN( zhbnew(ii,ij,jl), hR(ii,ij,jl)) - hL(ii,ij,jl)317 zvetamin(ji) = MAX( hi_max(jl), hL(ii,ij,jl) ) - hL(ii,ij,jl) 318 zvetamax(ji) = MIN( zhbnew(ii,ij,jl), hR(ii,ij,jl) ) - hL(ii,ij,jl) 413 319 zdonor(ii,ij,jl) = jl 414 320 415 ELSE ! zhbnew(jl) <= hi_max(jl) ; transfer from jl+1 to jl 416 321 ELSE ! zhbnew(jl) <= hi_max(jl) ; transfer from jl+1 to jl 417 322 ! left and right integration limits in eta space 418 323 zvetamin(ji) = 0.0 419 zvetamax(ji) = MIN( hi_max(jl), hR(ii,ij,jl+1)) - hL(ii,ij,jl+1)324 zvetamax(ji) = MIN( hi_max(jl), hR(ii,ij,jl+1) ) - hL(ii,ij,jl+1) 420 325 zdonor(ii,ij,jl) = jl + 1 421 326 422 ENDIF ! zhbnew(jl) > hi_max(jl)423 424 zetamax = MAX( zvetamax(ji), zvetamin(ji)) ! no transfer if etamax < etamin327 ENDIF 328 329 zetamax = MAX( zvetamax(ji), zvetamin(ji) ) ! no transfer if etamax < etamin 425 330 zetamin = zvetamin(ji) 426 331 427 332 zx1 = zetamax - zetamin 428 zwk1 = zetamin *zetamin429 zwk2 = zetamax *zetamax430 zx2 = 0.5 * ( zwk2 - zwk1)333 zwk1 = zetamin * zetamin 334 zwk2 = zetamax * zetamax 335 zx2 = 0.5 * ( zwk2 - zwk1 ) 431 336 zwk1 = zwk1 * zetamin 432 337 zwk2 = zwk2 * zetamax 433 zx3 = 1.0 /3.0 * (zwk2 - zwk1)338 zx3 = 1.0 / 3.0 * ( zwk2 - zwk1 ) 434 339 nd = zdonor(ii,ij,jl) 435 340 zdaice(ii,ij,jl) = g1(ii,ij,nd)*zx2 + g0(ii,ij,nd)*zx1 436 341 zdvice(ii,ij,jl) = g1(ii,ij,nd)*zx3 + g0(ii,ij,nd)*zx2 + zdaice(ii,ij,jl)*hL(ii,ij,nd) 437 342 438 END DO ! ji439 END DO ! jl klbnd -> kubnd - 1343 END DO 344 END DO 440 345 441 346 !!---------------------------------------------------------------------------------------------- … … 451 356 ii = nind_i(ji) 452 357 ij = nind_j(ji) 453 IF ( a_i(ii,ij,1) > epsi10 .AND. ht_i(ii,ij,1) < hiclim) THEN454 a_i (ii,ij,1) = a_i(ii,ij,1) * ht_i(ii,ij,1) / hiclim455 ht_i(ii,ij,1) = hiclim358 IF ( a_i(ii,ij,1) > epsi10 .AND. ht_i(ii,ij,1) < rn_himin ) THEN 359 a_i (ii,ij,1) = a_i(ii,ij,1) * ht_i(ii,ij,1) / rn_himin 360 ht_i(ii,ij,1) = rn_himin 456 361 ENDIF 457 END DO !ji362 END DO 458 363 459 364 !!---------------------------------------------------------------------------------------------- … … 479 384 ENDIF 480 385 481 CALL wrk_dealloc( jpi,jpj, zremap_flag ) ! integer482 CALL wrk_dealloc( jpi,jpj,jpl-1, zdonor ) ! integer386 CALL wrk_dealloc( jpi,jpj, zremap_flag ) 387 CALL wrk_dealloc( jpi,jpj,jpl-1, zdonor ) 483 388 CALL wrk_dealloc( jpi,jpj,jpl, zdhice, g0, g1, hL, hR, zht_i_b, dummy_es ) 484 389 CALL wrk_dealloc( jpi,jpj,jpl-1, zdaice, zdvice ) 485 390 CALL wrk_dealloc( jpi,jpj,jpl+1, zhbnew, kkstart = 0 ) 486 391 CALL wrk_dealloc( (jpi+1)*(jpj+1), zvetamin, zvetamax ) 487 CALL wrk_dealloc( (jpi+1)*(jpj+1), nind_i, nind_j ) ! integer392 CALL wrk_dealloc( (jpi+1)*(jpj+1), nind_i, nind_j ) 488 393 CALL wrk_dealloc( jpi,jpj, zhb0,zhb1,vt_i_init,vt_i_final,vt_s_init,vt_s_final,et_i_init,et_i_final,et_s_init,et_s_final ) 489 394 … … 491 396 492 397 493 SUBROUTINE lim_itd_fitline( num_cat, HbL, Hbr, hice, & 494 & g0, g1, hL, hR, zremap_flag ) 398 SUBROUTINE lim_itd_fitline( num_cat, HbL, Hbr, hice, g0, g1, hL, hR, zremap_flag ) 495 399 !!------------------------------------------------------------------ 496 400 !! *** ROUTINE lim_itd_fitline *** … … 511 415 INTEGER , DIMENSION(jpi,jpj), INTENT(in ) :: zremap_flag ! 512 416 ! 513 INTEGER :: ji,jj! horizontal indices417 INTEGER :: ji,jj ! horizontal indices 514 418 REAL(wp) :: zh13 ! HbL + 1/3 * (HbR - HbL) 515 419 REAL(wp) :: zh23 ! HbL + 2/3 * (HbR - HbL) … … 518 422 !!------------------------------------------------------------------ 519 423 ! 520 !521 424 DO jj = 1, jpj 522 425 DO ji = 1, jpi 523 426 ! 524 427 IF( zremap_flag(ji,jj) == 1 .AND. a_i(ji,jj,num_cat) > epsi10 & 525 & .AND. hice(ji,jj) > 0._wp )THEN428 & .AND. hice(ji,jj) > 0._wp ) THEN 526 429 527 430 ! Initialize hL and hR 528 529 431 hL(ji,jj) = HbL(ji,jj) 530 432 hR(ji,jj) = HbR(ji,jj) 531 433 532 434 ! Change hL or hR if hice falls outside central third of range 533 534 zh13 = 1.0/3.0 * (2.0*hL(ji,jj) + hR(ji,jj)) 535 zh23 = 1.0/3.0 * (hL(ji,jj) + 2.0*hR(ji,jj)) 435 zh13 = 1.0 / 3.0 * ( 2.0 * hL(ji,jj) + hR(ji,jj) ) 436 zh23 = 1.0 / 3.0 * ( hL(ji,jj) + 2.0 * hR(ji,jj) ) 536 437 537 438 IF ( hice(ji,jj) < zh13 ) THEN ; hR(ji,jj) = 3._wp * hice(ji,jj) - 2._wp * hL(ji,jj) … … 540 441 541 442 ! Compute coefficients of g(eta) = g0 + g1*eta 542 543 443 zdhr = 1._wp / (hR(ji,jj) - hL(ji,jj)) 544 444 zwk1 = 6._wp * a_i(ji,jj,num_cat) * zdhr 545 445 zwk2 = ( hice(ji,jj) - hL(ji,jj) ) * zdhr 546 g0(ji,jj) = zwk1 * ( 2._wp /3._wp - zwk2 )547 g1(ji,jj) = 2._wp * zdhr * zwk1 * ( zwk2 - 0.5)446 g0(ji,jj) = zwk1 * ( 2._wp / 3._wp - zwk2 ) 447 g1(ji,jj) = 2._wp * zdhr * zwk1 * ( zwk2 - 0.5 ) 548 448 ! 549 ELSE 449 ELSE ! remap_flag = .false. or a_i < epsi10 550 450 hL(ji,jj) = 0._wp 551 451 hR(ji,jj) = 0._wp 552 452 g0(ji,jj) = 0._wp 553 453 g1(ji,jj) = 0._wp 554 ENDIF ! a_i > epsi10454 ENDIF 555 455 ! 556 456 END DO … … 576 476 577 477 INTEGER :: ji, jj, jl, jl2, jl1, jk ! dummy loop indices 578 INTEGER :: ii, ij ! indices when changing from 2D-1D is done478 INTEGER :: ii, ij ! indices when changing from 2D-1D is done 579 479 580 480 REAL(wp), POINTER, DIMENSION(:,:,:) :: zaTsfn … … 589 489 INTEGER, POINTER, DIMENSION(:) :: nind_i, nind_j ! compressed indices for i/j directions 590 490 591 INTEGER :: nbrem ! number of cells with ice to transfer 592 593 LOGICAL :: zdaice_negative ! true if daice < -puny 594 LOGICAL :: zdvice_negative ! true if dvice < -puny 595 LOGICAL :: zdaice_greater_aicen ! true if daice > aicen 596 LOGICAL :: zdvice_greater_vicen ! true if dvice > vicen 491 INTEGER :: nbrem ! number of cells with ice to transfer 597 492 !!------------------------------------------------------------------ 598 493 599 494 CALL wrk_alloc( jpi,jpj,jpl, zaTsfn ) 600 495 CALL wrk_alloc( jpi,jpj, zworka ) 601 CALL wrk_alloc( (jpi+1)*(jpj+1), nind_i, nind_j ) ! integer496 CALL wrk_alloc( (jpi+1)*(jpj+1), nind_i, nind_j ) 602 497 603 498 !---------------------------------------------------------------------------------------------- … … 606 501 607 502 DO jl = klbnd, kubnd 608 zaTsfn(:,:,jl) = a_i(:,:,jl)*t_su(:,:,jl) 609 END DO 610 611 !---------------------------------------------------------------------------------------------- 612 ! 2) Check for daice or dvice out of range, allowing for roundoff error 613 !---------------------------------------------------------------------------------------------- 614 ! Note: zdaice < 0 or zdvice < 0 usually happens when category jl 615 ! has a small area, with h(n) very close to a boundary. Then 616 ! the coefficients of g(h) are large, and the computed daice and 617 ! dvice can be in error. If this happens, it is best to transfer 618 ! either the entire category or nothing at all, depending on which 619 ! side of the boundary hice(n) lies. 620 !----------------------------------------------------------------- 621 DO jl = klbnd, kubnd-1 622 623 zdaice_negative = .false. 624 zdvice_negative = .false. 625 zdaice_greater_aicen = .false. 626 zdvice_greater_vicen = .false. 627 628 DO jj = 1, jpj 629 DO ji = 1, jpi 630 631 IF (zdonor(ji,jj,jl) .GT. 0) THEN 632 jl1 = zdonor(ji,jj,jl) 633 634 IF (zdaice(ji,jj,jl) .LT. 0.0) THEN 635 IF (zdaice(ji,jj,jl) .GT. -epsi10) THEN 636 IF ( ( jl1.EQ.jl .AND. ht_i(ji,jj,jl1) .GT. hi_max(jl) ) & 637 .OR. & 638 ( jl1.EQ.jl+1 .AND. ht_i(ji,jj,jl1) .LE. hi_max(jl) ) & 639 ) THEN 640 zdaice(ji,jj,jl) = a_i(ji,jj,jl1) ! shift entire category 641 zdvice(ji,jj,jl) = v_i(ji,jj,jl1) 642 ELSE 643 zdaice(ji,jj,jl) = 0.0 ! shift no ice 644 zdvice(ji,jj,jl) = 0.0 645 ENDIF 646 ELSE 647 zdaice_negative = .true. 648 ENDIF 649 ENDIF 650 651 IF (zdvice(ji,jj,jl) .LT. 0.0) THEN 652 IF (zdvice(ji,jj,jl) .GT. -epsi10 ) THEN 653 IF ( ( jl1.EQ.jl .AND. ht_i(ji,jj,jl1).GT.hi_max(jl) ) & 654 .OR. & 655 ( jl1.EQ.jl+1 .AND. ht_i(ji,jj,jl1) .LE. hi_max(jl) ) & 656 ) THEN 657 zdaice(ji,jj,jl) = a_i(ji,jj,jl1) ! shift entire category 658 zdvice(ji,jj,jl) = v_i(ji,jj,jl1) 659 ELSE 660 zdaice(ji,jj,jl) = 0.0 ! shift no ice 661 zdvice(ji,jj,jl) = 0.0 662 ENDIF 663 ELSE 664 zdvice_negative = .true. 665 ENDIF 666 ENDIF 667 668 ! If daice is close to aicen, set daice = aicen. 669 IF (zdaice(ji,jj,jl) .GT. a_i(ji,jj,jl1) - epsi10 ) THEN 670 IF (zdaice(ji,jj,jl) .LT. a_i(ji,jj,jl1)+epsi10) THEN 671 zdaice(ji,jj,jl) = a_i(ji,jj,jl1) 672 zdvice(ji,jj,jl) = v_i(ji,jj,jl1) 673 ELSE 674 zdaice_greater_aicen = .true. 675 ENDIF 676 ENDIF 677 678 IF (zdvice(ji,jj,jl) .GT. v_i(ji,jj,jl1)-epsi10) THEN 679 IF (zdvice(ji,jj,jl) .LT. v_i(ji,jj,jl1)+epsi10) THEN 680 zdaice(ji,jj,jl) = a_i(ji,jj,jl1) 681 zdvice(ji,jj,jl) = v_i(ji,jj,jl1) 682 ELSE 683 zdvice_greater_vicen = .true. 684 ENDIF 685 ENDIF 686 687 ENDIF ! donor > 0 688 END DO ! i 689 END DO ! j 690 691 END DO !jl 503 zaTsfn(:,:,jl) = a_i(:,:,jl) * t_su(:,:,jl) 504 END DO 692 505 693 506 !------------------------------------------------------------------------------- 694 ! 3) Transfer volume and energy between categories507 ! 2) Transfer volume and energy between categories 695 508 !------------------------------------------------------------------------------- 696 509 … … 699 512 DO jj = 1, jpj 700 513 DO ji = 1, jpi 701 IF (zdaice(ji,jj,jl) .GT.0.0 ) THEN ! daice(n) can be < puny514 IF (zdaice(ji,jj,jl) > 0.0 ) THEN ! daice(n) can be < puny 702 515 nbrem = nbrem + 1 703 516 nind_i(nbrem) = ji 704 517 nind_j(nbrem) = jj 705 ENDIF ! tmask518 ENDIF 706 519 END DO 707 520 END DO … … 712 525 713 526 jl1 = zdonor(ii,ij,jl) 714 rswitch = MAX( 0.0 , SIGN( 1.0, v_i(ii,ij,jl1) - epsi10 ) )715 zworka(ii,ij) = zdvice(ii,ij,jl) / MAX(v_i(ii,ij,jl1),epsi10) * rswitch527 rswitch = MAX( 0._wp , SIGN( 1._wp , v_i(ii,ij,jl1) - epsi10 ) ) 528 zworka(ii,ij) = zdvice(ii,ij,jl) / MAX( v_i(ii,ij,jl1), epsi10 ) * rswitch 716 529 IF( jl1 == jl) THEN ; jl2 = jl1+1 717 ELSE 530 ELSE ; jl2 = jl 718 531 ENDIF 719 532 … … 721 534 ! Ice areas 722 535 !-------------- 723 724 536 a_i(ii,ij,jl1) = a_i(ii,ij,jl1) - zdaice(ii,ij,jl) 725 537 a_i(ii,ij,jl2) = a_i(ii,ij,jl2) + zdaice(ii,ij,jl) … … 728 540 ! Ice volumes 729 541 !-------------- 730 731 542 v_i(ii,ij,jl1) = v_i(ii,ij,jl1) - zdvice(ii,ij,jl) 732 543 v_i(ii,ij,jl2) = v_i(ii,ij,jl2) + zdvice(ii,ij,jl) … … 735 546 ! Snow volumes 736 547 !-------------- 737 738 548 zdvsnow = v_s(ii,ij,jl1) * zworka(ii,ij) 739 549 v_s(ii,ij,jl1) = v_s(ii,ij,jl1) - zdvsnow … … 743 553 ! Snow heat content 744 554 !-------------------- 745 746 555 zdesnow = e_s(ii,ij,1,jl1) * zworka(ii,ij) 747 556 e_s(ii,ij,1,jl1) = e_s(ii,ij,1,jl1) - zdesnow … … 751 560 ! Ice age 752 561 !-------------- 753 754 562 zdo_aice = oa_i(ii,ij,jl1) * zdaice(ii,ij,jl) 755 563 oa_i(ii,ij,jl1) = oa_i(ii,ij,jl1) - zdo_aice … … 759 567 ! Ice salinity 760 568 !-------------- 761 762 569 zdsm_vice = smv_i(ii,ij,jl1) * zworka(ii,ij) 763 570 smv_i(ii,ij,jl1) = smv_i(ii,ij,jl1) - zdsm_vice … … 767 574 ! Surface temperature 768 575 !--------------------- 769 770 576 zdaTsf = t_su(ii,ij,jl1) * zdaice(ii,ij,jl) 771 577 zaTsfn(ii,ij,jl1) = zaTsfn(ii,ij,jl1) - zdaTsf 772 578 zaTsfn(ii,ij,jl2) = zaTsfn(ii,ij,jl2) + zdaTsf 773 579 774 END DO ! ji580 END DO 775 581 776 582 !------------------ … … 779 585 780 586 DO jk = 1, nlay_i 781 !CDIR NODEP782 587 DO ji = 1, nbrem 783 588 ii = nind_i(ji) … … 785 590 786 591 jl1 = zdonor(ii,ij,jl) 787 IF (jl1 .EQ.jl) THEN592 IF (jl1 == jl) THEN 788 593 jl2 = jl+1 789 594 ELSE ! n1 = n+1 … … 794 599 e_i(ii,ij,jk,jl1) = e_i(ii,ij,jk,jl1) - zdeice 795 600 e_i(ii,ij,jk,jl2) = e_i(ii,ij,jk,jl2) + zdeice 796 END DO ! ji797 END DO ! jk601 END DO 602 END DO 798 603 799 604 END DO ! boundaries, 1 to ncat-1 … … 809 614 ht_i(ji,jj,jl) = v_i (ji,jj,jl) / a_i(ji,jj,jl) 810 615 t_su(ji,jj,jl) = zaTsfn(ji,jj,jl) / a_i(ji,jj,jl) 811 rswitch = 1.0 - MAX(0.0,SIGN(1.0,-v_s(ji,jj,jl)+epsi10)) !0 if no ice and 1 if yes812 616 ELSE 813 617 ht_i(ji,jj,jl) = 0._wp 814 t_su(ji,jj,jl) = rt t618 t_su(ji,jj,jl) = rt0 815 619 ENDIF 816 END DO ! ji817 END DO ! jj818 END DO ! jl620 END DO 621 END DO 622 END DO 819 623 ! 820 624 CALL wrk_dealloc( jpi,jpj,jpl, zaTsfn ) 821 625 CALL wrk_dealloc( jpi,jpj, zworka ) 822 CALL wrk_dealloc( (jpi+1)*(jpj+1), nind_i, nind_j ) ! integer626 CALL wrk_dealloc( (jpi+1)*(jpj+1), nind_i, nind_j ) 823 627 ! 824 628 END SUBROUTINE lim_itd_shiftice … … 846 650 REAL(wp), POINTER, DIMENSION(:,:) :: vt_s_init, vt_s_final ! snow volume summed over categories 847 651 !!------------------------------------------------------------------ 848 !! clem 2014/04: be carefull, rebining does not conserve salt(maybe?) => the difference is taken into account in limupdate849 652 850 653 CALL wrk_alloc( jpi,jpj,jpl, zdonor ) ! interger … … 864 667 DO jj = 1, jpj 865 668 DO ji = 1, jpi 866 IF( a_i(ji,jj,jl) > epsi10 ) THEN 867 ht_i(ji,jj,jl) = v_i(ji,jj,jl) / a_i(ji,jj,jl) 868 ELSE 869 ht_i(ji,jj,jl) = 0._wp 870 ENDIF 669 rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi10 ) ) 670 ht_i(ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi10 ) * rswitch 871 671 END DO 872 672 END DO … … 874 674 875 675 !------------------------------------------------------------------------------ 876 ! 2) Make sure thickness of cat klbnd is at least hi_max(klbnd) 877 !------------------------------------------------------------------------------ 878 DO jj = 1, jpj 879 DO ji = 1, jpi 880 IF( a_i(ji,jj,klbnd) > epsi10 ) THEN 881 IF( ht_i(ji,jj,klbnd) <= hi_max(0) .AND. hi_max(0) > 0._wp ) THEN 882 a_i(ji,jj,klbnd) = v_i(ji,jj,klbnd) / hi_max(0) 883 ht_i(ji,jj,klbnd) = hi_max(0) 884 ENDIF 885 ENDIF 886 END DO 887 END DO 888 889 !------------------------------------------------------------------------------ 890 ! 3) If a category thickness is not in bounds, shift the 676 ! 2) If a category thickness is not in bounds, shift the 891 677 ! entire area, volume, and energy to the neighboring category 892 678 !------------------------------------------------------------------------------ … … 917 703 zdonor(ji,jj,jl) = jl 918 704 ! begin TECLIM change 919 !zdaice(ji,jj,jl) = a_i(ji,jj,jl)920 !zdvice(ji,jj,jl) = v_i(ji,jj,jl)921 705 !zdaice(ji,jj,jl) = a_i(ji,jj,jl) * 0.5_wp 922 706 !zdvice(ji,jj,jl) = v_i(ji,jj,jl)-zdaice(ji,jj,jl)*(hi_max(jl)+hi_max(jl-1)) * 0.5_wp 923 707 ! end TECLIM change 924 708 ! clem: how much of a_i you send in cat sup is somewhat arbitrary 925 zdaice(ji,jj,jl) = a_i(ji,jj,jl) * ( ht_i(ji,jj,jl) - hi_max(jl) + epsi 10 ) / ht_i(ji,jj,jl)926 zdvice(ji,jj,jl) = v_i(ji,jj,jl) - ( a_i(ji,jj,jl) - zdaice(ji,jj,jl) ) * ( hi_max(jl) - epsi 10 )709 zdaice(ji,jj,jl) = a_i(ji,jj,jl) * ( ht_i(ji,jj,jl) - hi_max(jl) + epsi20 ) / ht_i(ji,jj,jl) 710 zdvice(ji,jj,jl) = v_i(ji,jj,jl) - ( a_i(ji,jj,jl) - zdaice(ji,jj,jl) ) * ( hi_max(jl) - epsi20 ) 927 711 ENDIF 928 END DO ! ji929 END DO ! jj712 END DO 713 END DO 930 714 IF(lk_mpp) CALL mpp_max( zshiftflag ) 931 715 … … 938 722 ENDIF 939 723 ! 940 END DO ! jl724 END DO 941 725 942 726 !---------------------------- … … 951 735 zshiftflag = 0 952 736 953 !clem-change954 737 DO jj = 1, jpj 955 738 DO ji = 1, jpi … … 961 744 zdvice(ji,jj,jl) = v_i(ji,jj,jl+1) 962 745 ENDIF 963 END DO ! ji964 END DO ! jj746 END DO 747 END DO 965 748 966 749 IF(lk_mpp) CALL mpp_max( zshiftflag ) … … 973 756 zdvice(:,:,jl) = 0._wp 974 757 ENDIF 975 !clem-change 976 977 ! ! clem-change begin: why not doing that? 978 ! DO jj = 1, jpj 979 ! DO ji = 1, jpi 980 ! IF( a_i(ji,jj,jl+1) > epsi10 .AND. ht_i(ji,jj,jl+1) <= hi_max(jl) ) THEN 981 ! ht_i(ji,jj,jl+1) = hi_max(jl) + epsi10 982 ! a_i (ji,jj,jl+1) = v_i(ji,jj,jl+1) / ht_i(ji,jj,jl+1) 983 ! ENDIF 984 ! END DO ! ji 985 ! END DO ! jj 986 ! clem-change end 987 988 END DO ! jl 758 759 END DO 989 760 990 761 !------------------------------------------------------------------------------ 991 ! 4) Conservation check762 ! 3) Conservation check 992 763 !------------------------------------------------------------------------------ 993 764 … … 1002 773 ENDIF 1003 774 ! 1004 CALL wrk_dealloc( jpi,jpj,jpl, zdonor ) ! interger775 CALL wrk_dealloc( jpi,jpj,jpl, zdonor ) 1005 776 CALL wrk_dealloc( jpi,jpj,jpl, zdaice, zdvice ) 1006 777 CALL wrk_dealloc( jpi,jpj, vt_i_init, vt_i_final, vt_s_init, vt_s_final ) … … 1013 784 !!---------------------------------------------------------------------- 1014 785 CONTAINS 1015 SUBROUTINE lim_itd_th ! Empty routines1016 END SUBROUTINE lim_itd_th1017 SUBROUTINE lim_itd_th_ini1018 END SUBROUTINE lim_itd_th_ini1019 786 SUBROUTINE lim_itd_th_rem 1020 787 END SUBROUTINE lim_itd_th_rem -
branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/LIM_SRC_3/limmsh.F90
r4161 r5682 23 23 PRIVATE 24 24 25 PUBLIC lim_msh ! routine called by ice_ini.F9025 PUBLIC lim_msh ! routine called by sbcice_lim.F90 26 26 27 27 !!---------------------------------------------------------------------- … … 41 41 !! - Definition of some constants linked with the grid 42 42 !! - Definition of the metric coef. for the sea/ice 43 !! - Initialization of the ice masks (tmsk, umsk)44 43 !! 45 44 !! Reference : Deleersnijder et al. Ocean Modelling 100, 7-10 … … 103 102 !!gm end 104 103 105 ! !== ice masks ==!106 tms(:,:) = tmask(:,:,1) ! ice T-point : use surface tmask107 tmu(:,:) = umask(:,:,1) ! ice U-point : use surface umask (C-grid EVP)108 tmv(:,:) = vmask(:,:,1) ! ice V-point : use surface vmask (C-grid EVP)109 DO jj = 1, jpjm1 ! ice F-point : recompute fmask (due to nn_shlat)110 DO ji = 1 , jpim1 ! NO vector opt.111 tmf(ji,jj) = tms(ji,jj) * tms(ji+1,jj) * tms(ji,jj+1) * tms(ji+1,jj+1)112 END DO113 END DO114 CALL lbc_lnk( tmf(:,:), 'F', 1. ) ! lateral boundary conditions115 116 ! !== unmasked and masked area of T-grid cell117 area(:,:) = e1t(:,:) * e2t(:,:)118 104 ! 119 105 END SUBROUTINE lim_msh -
branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90
r4990 r5682 102 102 !! and charge ellipse. 103 103 !! The user should make sure that the parameters 104 !! n evp, telast andcreepl maintain stress state104 !! nn_nevp, elastic time scale and rn_creepl maintain stress state 105 105 !! on the charge ellipse for plastic flow 106 106 !! e.g. in the Canadian Archipelago … … 108 108 !! References : Hunke and Dukowicz, JPO97 109 109 !! Bouillon et al., Ocean Modelling 2009 110 !! Vancoppenolle et al., Ocean Modelling 2008111 110 !!------------------------------------------------------------------- 112 111 INTEGER, INTENT(in) :: k_j1 ! southern j-index for ice computation … … 117 116 CHARACTER (len=50) :: charout 118 117 REAL(wp) :: zt11, zt12, zt21, zt22, ztagnx, ztagny, delta ! 119 REAL(wp) :: za, zstms , zmask! local scalars120 REAL(wp) :: zc1, zc2, zc3 121 122 REAL(wp) :: dtevp ! time step for subcycling123 REAL(wp) :: dtotel, ecc2, ecci ! square of yield ellipse eccenticity124 REAL(wp) :: z0, zr, zcca, zccb ! temporary scalars125 REAL(wp) :: zu_ice2, zv_ice1 !126 REAL(wp) :: zddc, zdtc ! delta on corners and on centre127 REAL(wp) :: zdst ! shear at the center of the grid point128 REAL(wp) :: zdsshx, zdsshy ! term for the gradient of ocean surface129 REAL(wp) :: sigma1, sigma2 ! internal ice stress118 REAL(wp) :: za, zstms ! local scalars 119 REAL(wp) :: zc1, zc2, zc3 ! ice mass 120 121 REAL(wp) :: dtevp , z1_dtevp ! time step for subcycling 122 REAL(wp) :: dtotel, z1_dtotel, ecc2, ecci ! square of yield ellipse eccenticity 123 REAL(wp) :: z0, zr, zcca, zccb ! temporary scalars 124 REAL(wp) :: zu_ice2, zv_ice1 ! 125 REAL(wp) :: zddc, zdtc ! delta on corners and on centre 126 REAL(wp) :: zdst ! shear at the center of the grid point 127 REAL(wp) :: zdsshx, zdsshy ! term for the gradient of ocean surface 128 REAL(wp) :: sigma1, sigma2 ! internal ice stress 130 129 131 130 REAL(wp) :: zresm ! Maximal error on ice velocity 132 REAL(wp) :: zdummy ! dummy argument133 131 REAL(wp) :: zintb, zintn ! dummy argument 134 132 … … 139 137 REAL(wp), POINTER, DIMENSION(:,:) :: zcorl1, zcorl2 ! coriolis parameter on U/V points 140 138 REAL(wp), POINTER, DIMENSION(:,:) :: za1ct , za2ct ! temporary arrays 141 REAL(wp), POINTER, DIMENSION(:,:) :: u_oce1, v_oce1! ocean u/v component on U points142 REAL(wp), POINTER, DIMENSION(:,:) :: u_oce2 , v_oce2! ocean u/v component on V points139 REAL(wp), POINTER, DIMENSION(:,:) :: v_oce1 ! ocean u/v component on U points 140 REAL(wp), POINTER, DIMENSION(:,:) :: u_oce2 ! ocean u/v component on V points 143 141 REAL(wp), POINTER, DIMENSION(:,:) :: u_ice2, v_ice1 ! ice u/v component on V/U point 144 142 REAL(wp), POINTER, DIMENSION(:,:) :: zf1 , zf2 ! arrays for internal stresses 143 REAL(wp), POINTER, DIMENSION(:,:) :: zmask ! mask ocean grid points 145 144 146 145 REAL(wp), POINTER, DIMENSION(:,:) :: zdt ! tension at centre of grid cells … … 152 151 ! ocean surface (ssh_m) if ice is not embedded 153 152 ! ice top surface if ice is embedded 153 154 REAL(wp), PARAMETER :: zepsi = 1.0e-20_wp ! tolerance parameter 155 REAL(wp), PARAMETER :: zvmin = 1.0e-03_wp ! ice volume below which ice velocity equals ocean velocity 154 156 !!------------------------------------------------------------------- 155 157 156 158 CALL wrk_alloc( jpi,jpj, zpresh, zfrld1, zmass1, zcorl1, za1ct , zpreshc, zfrld2, zmass2, zcorl2, za2ct ) 157 CALL wrk_alloc( jpi,jpj, u_oce 1, u_oce2, u_ice2, v_oce1 , v_oce2, v_ice1)159 CALL wrk_alloc( jpi,jpj, u_oce2, u_ice2, v_oce1 , v_ice1 , zmask ) 158 160 CALL wrk_alloc( jpi,jpj, zf1 , zu_ice, zf2 , zv_ice , zdt , zds ) 159 161 CALL wrk_alloc( jpi,jpj, zdt , zds , zs1 , zs2 , zs12 , zresr , zpice ) … … 161 163 #if defined key_lim2 && ! defined key_lim2_vp 162 164 # if defined key_agrif 163 USE ice_2, vt_s => hsnm164 USE ice_2, vt_i => hicm165 USE ice_2, vt_s => hsnm 166 USE ice_2, vt_i => hicm 165 167 # else 166 vt_s => hsnm167 vt_i => hicm168 vt_s => hsnm 169 vt_i => hicm 168 170 # endif 169 at_i(:,:) = 1. - frld(:,:)171 at_i(:,:) = 1. - frld(:,:) 170 172 #endif 171 173 #if defined key_agrif && defined key_lim2 172 CALL agrif_rhg_lim2_load ! First interpolation of coarse values174 CALL agrif_rhg_lim2_load ! First interpolation of coarse values 173 175 #endif 174 176 ! … … 186 188 187 189 #if defined key_lim3 188 CALL lim_itd_me_icestrength( ridge_scheme_swi ) ! LIM-3: Ice strength on T-points 189 #endif 190 191 !CDIR NOVERRCHK 190 CALL lim_itd_me_icestrength( nn_icestr ) ! LIM-3: Ice strength on T-points 191 #endif 192 192 193 DO jj = k_j1 , k_jpj ! Ice mass and temp variables 193 !CDIR NOVERRCHK194 194 DO ji = 1 , jpi 195 195 #if defined key_lim3 196 zpresh(ji,jj) = tm s(ji,jj) * strength(ji,jj)196 zpresh(ji,jj) = tmask(ji,jj,1) * strength(ji,jj) 197 197 #endif 198 198 #if defined key_lim2 199 zpresh(ji,jj) = tm s(ji,jj) * pstar * vt_i(ji,jj) * EXP( -c_rhg * (1. - at_i(ji,jj) ) )200 #endif 201 ! tmi= 1 where there is ice or on land202 tmi(ji,jj) = 1._wp - ( 1._wp - MAX( 0._wp , SIGN ( 1._wp , vt_i(ji,jj) - epsd ) ) ) * tms(ji,jj)199 zpresh(ji,jj) = tmask(ji,jj,1) * pstar * vt_i(ji,jj) * EXP( -c_rhg * (1. - at_i(ji,jj) ) ) 200 #endif 201 ! zmask = 1 where there is ice or on land 202 zmask(ji,jj) = 1._wp - ( 1._wp - MAX( 0._wp , SIGN ( 1._wp , vt_i(ji,jj) - zepsi ) ) ) * tmask(ji,jj,1) 203 203 END DO 204 204 END DO … … 206 206 ! Ice strength on grid cell corners (zpreshc) 207 207 ! needed for calculation of shear stress 208 !CDIR NOVERRCHK209 208 DO jj = k_j1+1, k_jpj-1 210 !CDIR NOVERRCHK211 209 DO ji = 2, jpim1 !RB caution no fs_ (ji+1,jj+1) 212 zstms = tms(ji+1,jj+1) * wght(ji+1,jj+1,2,2) + & 213 & tms(ji,jj+1) * wght(ji+1,jj+1,1,2) + & 214 & tms(ji+1,jj) * wght(ji+1,jj+1,2,1) + & 215 & tms(ji,jj) * wght(ji+1,jj+1,1,1) 216 zpreshc(ji,jj) = ( zpresh(ji+1,jj+1) * wght(ji+1,jj+1,2,2) + & 217 & zpresh(ji,jj+1) * wght(ji+1,jj+1,1,2) + & 218 & zpresh(ji+1,jj) * wght(ji+1,jj+1,2,1) + & 219 & zpresh(ji,jj) * wght(ji+1,jj+1,1,1) & 220 & ) / MAX( zstms, epsd ) 210 zstms = tmask(ji+1,jj+1,1) * wght(ji+1,jj+1,2,2) + tmask(ji,jj+1,1) * wght(ji+1,jj+1,1,2) + & 211 & tmask(ji+1,jj,1) * wght(ji+1,jj+1,2,1) + tmask(ji,jj,1) * wght(ji+1,jj+1,1,1) 212 zpreshc(ji,jj) = ( zpresh(ji+1,jj+1) * wght(ji+1,jj+1,2,2) + zpresh(ji,jj+1) * wght(ji+1,jj+1,1,2) + & 213 & zpresh(ji+1,jj) * wght(ji+1,jj+1,2,1) + zpresh(ji,jj) * wght(ji+1,jj+1,1,1) & 214 & ) / MAX( zstms, zepsi ) 221 215 END DO 222 216 END DO 223 224 217 CALL lbc_lnk( zpreshc(:,:), 'F', 1. ) 225 218 ! … … 236 229 ! zcorl2: Coriolis parameter on V-points 237 230 ! (ztagnx,ztagny): wind stress on U/V points 238 ! u_oce1: ocean u component on u points239 231 ! v_oce1: ocean v component on u points 240 232 ! u_oce2: ocean u component on v points 241 ! v_oce2: ocean v component on v points242 233 243 234 IF( nn_ice_embd == 2 ) THEN !== embedded sea ice: compute representative ice top surface ==! 244 245 246 235 ! 236 ! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[n/nn_fsbc], n=0,nn_fsbc-1} 237 ! = (1/nn_fsbc)^2 * {SUM[n], n=0,nn_fsbc-1} 247 238 zintn = REAL( nn_fsbc - 1 ) / REAL( nn_fsbc ) * 0.5_wp 248 249 250 239 ! 240 ! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[1-n/nn_fsbc], n=0,nn_fsbc-1} 241 ! = (1/nn_fsbc)^2 * (nn_fsbc^2 - {SUM[n], n=0,nn_fsbc-1}) 251 242 zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp 252 243 ! 253 244 zpice(:,:) = ssh_m(:,:) + ( zintn * snwice_mass(:,:) + zintb * snwice_mass_b(:,:) ) * r1_rau0 254 245 ! 255 246 ELSE !== non-embedded sea ice: use ocean surface for slope calculation ==! 256 247 zpice(:,:) = ssh_m(:,:) … … 260 251 DO ji = fs_2, fs_jpim1 261 252 262 zc1 = tm s(ji ,jj) * ( rhosn * vt_s(ji ,jj ) + rhoic * vt_i(ji ,jj ) )263 zc2 = tm s(ji+1,jj) * ( rhosn * vt_s(ji+1,jj ) + rhoic * vt_i(ji+1,jj ) )264 zc3 = tm s(ji ,jj+1) * ( rhosn * vt_s(ji ,jj+1) + rhoic * vt_i(ji ,jj+1) )265 266 zt11 = tm s(ji ,jj) * e1t(ji ,jj)267 zt12 = tm s(ji+1,jj) * e1t(ji+1,jj)268 zt21 = tm s(ji,jj) * e2t(ji,jj )269 zt22 = tm s(ji,jj+1) * e2t(ji,jj+1)253 zc1 = tmask(ji ,jj ,1) * ( rhosn * vt_s(ji ,jj ) + rhoic * vt_i(ji ,jj ) ) 254 zc2 = tmask(ji+1,jj ,1) * ( rhosn * vt_s(ji+1,jj ) + rhoic * vt_i(ji+1,jj ) ) 255 zc3 = tmask(ji ,jj+1,1) * ( rhosn * vt_s(ji ,jj+1) + rhoic * vt_i(ji ,jj+1) ) 256 257 zt11 = tmask(ji ,jj,1) * e1t(ji ,jj) 258 zt12 = tmask(ji+1,jj,1) * e1t(ji+1,jj) 259 zt21 = tmask(ji,jj ,1) * e2t(ji,jj ) 260 zt22 = tmask(ji,jj+1,1) * e2t(ji,jj+1) 270 261 271 262 ! Leads area. 272 zfrld1(ji,jj) = ( zt12 * ( 1.0 - at_i(ji,jj) ) + zt11 * ( 1.0 - at_i(ji+1,jj) ) ) / ( zt11 + zt12 + epsd)273 zfrld2(ji,jj) = ( zt22 * ( 1.0 - at_i(ji,jj) ) + zt21 * ( 1.0 - at_i(ji,jj+1) ) ) / ( zt21 + zt22 + epsd)263 zfrld1(ji,jj) = ( zt12 * ( 1.0 - at_i(ji,jj) ) + zt11 * ( 1.0 - at_i(ji+1,jj) ) ) / ( zt11 + zt12 + zepsi ) 264 zfrld2(ji,jj) = ( zt22 * ( 1.0 - at_i(ji,jj) ) + zt21 * ( 1.0 - at_i(ji,jj+1) ) ) / ( zt21 + zt22 + zepsi ) 274 265 275 266 ! Mass, coriolis coeff. and currents 276 zmass1(ji,jj) = ( zt12 *zc1 + zt11*zc2 ) / (zt11+zt12+epsd)277 zmass2(ji,jj) = ( zt22 *zc1 + zt21*zc3 ) / (zt21+zt22+epsd)278 zcorl1(ji,jj) = zmass1(ji,jj) * ( e1t(ji+1,jj) *fcor(ji,jj) + e1t(ji,jj)*fcor(ji+1,jj) ) &279 & / ( e1t(ji,jj) + e1t(ji+1,jj) + epsd)280 zcorl2(ji,jj) = zmass2(ji,jj) * ( e2t(ji,jj+1) *fcor(ji,jj) + e2t(ji,jj)*fcor(ji,jj+1) ) &281 & / ( e2t(ji,jj+1) + e2t(ji,jj) + epsd)267 zmass1(ji,jj) = ( zt12 * zc1 + zt11 * zc2 ) / ( zt11 + zt12 + zepsi ) 268 zmass2(ji,jj) = ( zt22 * zc1 + zt21 * zc3 ) / ( zt21 + zt22 + zepsi ) 269 zcorl1(ji,jj) = zmass1(ji,jj) * ( e1t(ji+1,jj) * fcor(ji,jj) + e1t(ji,jj) * fcor(ji+1,jj) ) & 270 & / ( e1t(ji,jj) + e1t(ji+1,jj) + zepsi ) 271 zcorl2(ji,jj) = zmass2(ji,jj) * ( e2t(ji,jj+1) * fcor(ji,jj) + e2t(ji,jj) * fcor(ji,jj+1) ) & 272 & / ( e2t(ji,jj+1) + e2t(ji,jj) + zepsi ) 282 273 ! 283 u_oce1(ji,jj) = u_oce(ji,jj)284 v_oce2(ji,jj) = v_oce(ji,jj)285 286 274 ! Ocean has no slip boundary condition 287 v_oce1(ji,jj) = 0.5 *( (v_oce(ji,jj)+v_oce(ji,jj-1))*e1t(ji,jj)&288 & +(v_oce(ji+1,jj)+v_oce(ji+1,jj-1))*e1t(ji+1,jj))&289 & /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj)290 291 u_oce2(ji,jj) = 0.5 *((u_oce(ji,jj)+u_oce(ji-1,jj))*e2t(ji,jj)&292 & +(u_oce(ji,jj+1)+u_oce(ji-1,jj+1))*e2t(ji,jj+1))&293 & / (e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj)275 v_oce1(ji,jj) = 0.5 * ( ( v_oce(ji ,jj) + v_oce(ji ,jj-1) ) * e1t(ji,jj) & 276 & + ( v_oce(ji+1,jj) + v_oce(ji+1,jj-1) ) * e1t(ji+1,jj) ) & 277 & / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1) 278 279 u_oce2(ji,jj) = 0.5 * ( ( u_oce(ji,jj ) + u_oce(ji-1,jj ) ) * e2t(ji,jj) & 280 & + ( u_oce(ji,jj+1) + u_oce(ji-1,jj+1) ) * e2t(ji,jj+1) ) & 281 & / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1) 294 282 295 283 ! Wind stress at U,V-point … … 303 291 ! include it later 304 292 305 zdsshx = ( zpice(ji+1,jj) - zpice(ji,jj) ) /e1u(ji,jj)306 zdsshy = ( zpice(ji,jj+1) - zpice(ji,jj) ) /e2v(ji,jj)293 zdsshx = ( zpice(ji+1,jj) - zpice(ji,jj) ) * r1_e1u(ji,jj) 294 zdsshy = ( zpice(ji,jj+1) - zpice(ji,jj) ) * r1_e2v(ji,jj) 307 295 308 296 za1ct(ji,jj) = ztagnx - zmass1(ji,jj) * grav * zdsshx … … 318 306 ! 319 307 ! Time step for subcycling 320 dtevp = rdt_ice / nevp 308 dtevp = rdt_ice / nn_nevp 309 #if defined key_lim3 310 dtotel = dtevp / ( 2._wp * rn_relast * rdt_ice ) 311 #else 321 312 dtotel = dtevp / ( 2._wp * telast ) 322 313 #endif 314 z1_dtotel = 1._wp / ( 1._wp + dtotel ) 315 z1_dtevp = 1._wp / dtevp 323 316 !-ecc2: square of yield ellipse eccenticrity (reminder: must become a namelist parameter) 324 ecc2 = ecc *ecc317 ecc2 = rn_ecc * rn_ecc 325 318 ecci = 1. / ecc2 326 319 … … 331 324 332 325 ! !----------------------! 333 DO jter = 1 , n evp! loop over jter !326 DO jter = 1 , nn_nevp ! loop over jter ! 334 327 ! !----------------------! 335 328 DO jj = k_j1, k_jpj-1 … … 339 332 340 333 DO jj = k_j1+1, k_jpj-1 341 DO ji = fs_2, jpim1 !RB bug no vect opt due to tmi334 DO ji = fs_2, fs_jpim1 !RB bug no vect opt due to zmask 342 335 343 336 ! … … 360 353 ! 361 354 ! 362 divu_i(ji,jj) = ( e2u(ji,jj)*u_ice(ji,jj) & 363 & -e2u(ji-1,jj)*u_ice(ji-1,jj) & 364 & +e1v(ji,jj)*v_ice(ji,jj) & 365 & -e1v(ji,jj-1)*v_ice(ji,jj-1) & 366 & ) & 367 & / area(ji,jj) 368 369 zdt(ji,jj) = ( ( u_ice(ji,jj)/e2u(ji,jj) & 370 & -u_ice(ji-1,jj)/e2u(ji-1,jj) & 371 & )*e2t(ji,jj)*e2t(ji,jj) & 372 & -( v_ice(ji,jj)/e1v(ji,jj) & 373 & -v_ice(ji,jj-1)/e1v(ji,jj-1) & 374 & )*e1t(ji,jj)*e1t(ji,jj) & 375 & ) & 376 & / area(ji,jj) 355 divu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj) & 356 & + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) & 357 & ) * r1_e12t(ji,jj) 358 359 zdt(ji,jj) = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj) & 360 & - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj) & 361 & ) * r1_e12t(ji,jj) 377 362 378 363 ! 379 zds(ji,jj) = ( ( u_ice(ji,jj+1)/e1u(ji,jj+1) & 380 & -u_ice(ji,jj)/e1u(ji,jj) & 381 & )*e1f(ji,jj)*e1f(ji,jj) & 382 & +( v_ice(ji+1,jj)/e2v(ji+1,jj) & 383 & -v_ice(ji,jj)/e2v(ji,jj) & 384 & )*e2f(ji,jj)*e2f(ji,jj) & 385 & ) & 386 & / ( e1f(ji,jj) * e2f(ji,jj) ) * ( 2.0 - tmf(ji,jj) ) & 387 & * tmi(ji,jj) * tmi(ji,jj+1) & 388 & * tmi(ji+1,jj) * tmi(ji+1,jj+1) 389 390 391 v_ice1(ji,jj) = 0.5*( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji+1,jj) & 392 & +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji,jj)) & 393 & /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj) 394 395 u_ice2(ji,jj) = 0.5*( (u_ice(ji,jj)+u_ice(ji-1,jj))*e2t(ji,jj+1) & 396 & +(u_ice(ji,jj+1)+u_ice(ji-1,jj+1))*e2t(ji,jj)) & 397 & /(e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj) 398 399 END DO 400 END DO 401 CALL lbc_lnk( v_ice1, 'U', -1. ) ; CALL lbc_lnk( u_ice2, 'V', -1. ) ! lateral boundary cond. 402 403 !CDIR NOVERRCHK 364 zds(ji,jj) = ( ( u_ice(ji,jj+1) * r1_e1u(ji,jj+1) - u_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj) & 365 & + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) & 366 & ) * r1_e12f(ji,jj) * ( 2._wp - fmask(ji,jj,1) ) & 367 & * zmask(ji,jj) * zmask(ji,jj+1) * zmask(ji+1,jj) * zmask(ji+1,jj+1) 368 369 370 v_ice1(ji,jj) = 0.5_wp * ( ( v_ice(ji ,jj) + v_ice(ji ,jj-1) ) * e1t(ji+1,jj) & 371 & + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji ,jj) ) & 372 & / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1) 373 374 u_ice2(ji,jj) = 0.5_wp * ( ( u_ice(ji,jj ) + u_ice(ji-1,jj ) ) * e2t(ji,jj+1) & 375 & + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj ) ) & 376 & / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1) 377 END DO 378 END DO 379 380 CALL lbc_lnk_multi( v_ice1, 'U', -1., u_ice2, 'V', -1. ) ! lateral boundary cond. 381 404 382 DO jj = k_j1+1, k_jpj-1 405 !CDIR NOVERRCHK406 383 DO ji = fs_2, fs_jpim1 407 384 408 385 !- Calculate Delta at centre of grid cells 409 zdst = ( e2u(ji , jj) * v_ice1(ji ,jj) & 410 & - e2u(ji-1, jj) * v_ice1(ji-1,jj) & 411 & + e1v(ji, jj ) * u_ice2(ji,jj ) & 412 & - e1v(ji, jj-1) * u_ice2(ji,jj-1) & 413 & ) & 414 & / area(ji,jj) 415 416 delta = SQRT( divu_i(ji,jj)*divu_i(ji,jj) + ( zdt(ji,jj)*zdt(ji,jj) + zdst*zdst ) * usecc2 ) 417 delta_i(ji,jj) = delta + creepl 418 !-Calculate stress tensor components zs1 and zs2 419 !-at centre of grid cells (see section 3.5 of CICE user's guide). 420 zs1(ji,jj) = ( zs1(ji,jj) + dtotel * ( ( divu_i(ji,jj) / delta_i(ji,jj) - delta / delta_i(ji,jj) ) & 421 & * zpresh(ji,jj) ) ) / ( 1._wp + dtotel ) 422 zs2(ji,jj) = ( zs2(ji,jj) + dtotel * ( ecci * zdt(ji,jj) / delta_i(ji,jj) * zpresh(ji,jj) ) ) & 423 & / ( 1._wp + dtotel ) 424 425 END DO 426 END DO 427 428 CALL lbc_lnk( zs1(:,:), 'T', 1. ) 429 CALL lbc_lnk( zs2(:,:), 'T', 1. ) 430 431 !CDIR NOVERRCHK 432 DO jj = k_j1+1, k_jpj-1 433 !CDIR NOVERRCHK 434 DO ji = fs_2, fs_jpim1 386 zdst = ( e2u(ji,jj) * v_ice1(ji,jj) - e2u(ji-1,jj ) * v_ice1(ji-1,jj ) & 387 & + e1v(ji,jj) * u_ice2(ji,jj) - e1v(ji ,jj-1) * u_ice2(ji ,jj-1) & 388 & ) * r1_e12t(ji,jj) 389 390 delta = SQRT( divu_i(ji,jj)**2 + ( zdt(ji,jj)**2 + zdst**2 ) * usecc2 ) 391 delta_i(ji,jj) = delta + rn_creepl 392 435 393 !- Calculate Delta on corners 436 zddc = ( ( v_ice1(ji,jj+1)/e1u(ji,jj+1) & 437 & -v_ice1(ji,jj)/e1u(ji,jj) & 438 & )*e1f(ji,jj)*e1f(ji,jj) & 439 & +( u_ice2(ji+1,jj)/e2v(ji+1,jj) & 440 & -u_ice2(ji,jj)/e2v(ji,jj) & 441 & )*e2f(ji,jj)*e2f(ji,jj) & 442 & ) & 443 & / ( e1f(ji,jj) * e2f(ji,jj) ) 444 445 zdtc = (-( v_ice1(ji,jj+1)/e1u(ji,jj+1) & 446 & -v_ice1(ji,jj)/e1u(ji,jj) & 447 & )*e1f(ji,jj)*e1f(ji,jj) & 448 & +( u_ice2(ji+1,jj)/e2v(ji+1,jj) & 449 & -u_ice2(ji,jj)/e2v(ji,jj) & 450 & )*e2f(ji,jj)*e2f(ji,jj) & 451 & ) & 452 & / ( e1f(ji,jj) * e2f(ji,jj) ) 453 454 zddc = SQRT(zddc**2+(zdtc**2+zds(ji,jj)**2)*usecc2) + creepl 455 456 !-Calculate stress tensor component zs12 at corners (see section 3.5 of CICE user's guide). 457 zs12(ji,jj) = ( zs12(ji,jj) + dtotel * & 458 & ( ecci * zds(ji,jj) / ( 2._wp * zddc ) * zpreshc(ji,jj) ) ) & 459 & / ( 1.0 + dtotel ) 460 461 END DO ! ji 462 END DO ! jj 463 464 CALL lbc_lnk( zs12(:,:), 'F', 1. ) 465 394 zddc = ( ( v_ice1(ji,jj+1) * r1_e1u(ji,jj+1) - v_ice1(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj) & 395 & + ( u_ice2(ji+1,jj) * r1_e2v(ji+1,jj) - u_ice2(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) & 396 & ) * r1_e12f(ji,jj) 397 398 zdtc = (- ( v_ice1(ji,jj+1) * r1_e1u(ji,jj+1) - v_ice1(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj) & 399 & + ( u_ice2(ji+1,jj) * r1_e2v(ji+1,jj) - u_ice2(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) & 400 & ) * r1_e12f(ji,jj) 401 402 zddc = SQRT( zddc**2 + ( zdtc**2 + zds(ji,jj)**2 ) * usecc2 ) + rn_creepl 403 404 !-Calculate stress tensor components zs1 and zs2 at centre of grid cells (see section 3.5 of CICE user's guide). 405 zs1(ji,jj) = ( zs1 (ji,jj) + dtotel * ( divu_i(ji,jj) - delta ) / delta_i(ji,jj) * zpresh(ji,jj) & 406 & ) * z1_dtotel 407 zs2(ji,jj) = ( zs2 (ji,jj) + dtotel * ecci * zdt(ji,jj) / delta_i(ji,jj) * zpresh(ji,jj) & 408 & ) * z1_dtotel 409 !-Calculate stress tensor component zs12 at corners 410 zs12(ji,jj) = ( zs12(ji,jj) + dtotel * ecci * zds(ji,jj) / ( 2._wp * zddc ) * zpreshc(ji,jj) & 411 & ) * z1_dtotel 412 413 END DO 414 END DO 415 416 CALL lbc_lnk_multi( zs1 , 'T', 1., zs2, 'T', 1., zs12, 'F', 1. ) 417 466 418 ! Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) 467 419 DO jj = k_j1+1, k_jpj-1 468 420 DO ji = fs_2, fs_jpim1 469 421 !- contribution of zs1, zs2 and zs12 to zf1 470 zf1(ji,jj) = 0.5 *( (zs1(ji+1,jj)-zs1(ji,jj))*e2u(ji,jj)&471 & +(zs2(ji+1,jj)*e2t(ji+1,jj)**2-zs2(ji,jj)*e2t(ji,jj)**2)/e2u(ji,jj)&472 & +2.0*(zs12(ji,jj)*e1f(ji,jj)**2-zs12(ji,jj-1)*e1f(ji,jj-1)**2)/e1u(ji,jj)&473 & ) / ( e1u(ji,jj)*e2u(ji,jj))422 zf1(ji,jj) = 0.5 * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj) & 423 & + ( zs2(ji+1,jj) * e2t(ji+1,jj)**2 - zs2(ji,jj) * e2t(ji,jj)**2 ) * r1_e2u(ji,jj) & 424 & + 2.0 * ( zs12(ji,jj) * e1f(ji,jj)**2 - zs12(ji,jj-1) * e1f(ji,jj-1)**2 ) * r1_e1u(ji,jj) & 425 & ) * r1_e12u(ji,jj) 474 426 ! contribution of zs1, zs2 and zs12 to zf2 475 zf2(ji,jj) = 0.5*( (zs1(ji,jj+1)-zs1(ji,jj))*e1v(ji,jj) & 476 & -(zs2(ji,jj+1)*e1t(ji,jj+1)**2 - zs2(ji,jj)*e1t(ji,jj)**2)/e1v(ji,jj) & 477 & + 2.0*(zs12(ji,jj)*e2f(ji,jj)**2 - & 478 zs12(ji-1,jj)*e2f(ji-1,jj)**2)/e2v(ji,jj) & 479 & ) / ( e1v(ji,jj)*e2v(ji,jj) ) 427 zf2(ji,jj) = 0.5 * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj) & 428 & - ( zs2(ji,jj+1) * e1t(ji,jj+1)**2 - zs2(ji,jj) * e1t(ji,jj)**2 ) * r1_e1v(ji,jj) & 429 & + 2.0 * ( zs12(ji,jj) * e2f(ji,jj)**2 - zs12(ji-1,jj) * e2f(ji-1,jj)**2 ) * r1_e2v(ji,jj) & 430 & ) * r1_e12v(ji,jj) 480 431 END DO 481 432 END DO … … 487 438 IF (MOD(jter,2).eq.0) THEN 488 439 489 !CDIR NOVERRCHK490 440 DO jj = k_j1+1, k_jpj-1 491 !CDIR NOVERRCHK492 441 DO ji = fs_2, fs_jpim1 493 zmask = (1.0-MAX(0._wp,SIGN(1._wp,-zmass1(ji,jj))))*tmu(ji,jj)494 z0 = zmass1(ji,jj) /dtevp442 rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass1(ji,jj) ) ) ) * umask(ji,jj,1) 443 z0 = zmass1(ji,jj) * z1_dtevp 495 444 496 445 ! SB modif because ocean has no slip boundary condition 497 zv_ice1 = 0.5*( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji,jj) & 498 & +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji+1,jj)) & 499 & /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj) 500 za = rhoco*SQRT((u_ice(ji,jj)-u_oce1(ji,jj))**2 + & 501 (zv_ice1-v_oce1(ji,jj))**2) * (1.0-zfrld1(ji,jj)) 502 zr = z0*u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + & 503 za*(u_oce1(ji,jj)) 504 zcca = z0+za 446 zv_ice1 = 0.5 * ( ( v_ice(ji ,jj) + v_ice(ji ,jj-1) ) * e1t(ji ,jj) & 447 & + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji+1,jj) ) & 448 & / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1) 449 za = rhoco * SQRT( ( u_ice(ji,jj) - u_oce(ji,jj) )**2 + & 450 & ( zv_ice1 - v_oce1(ji,jj) )**2 ) * ( 1.0 - zfrld1(ji,jj) ) 451 zr = z0 * u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + za * u_oce(ji,jj) 452 zcca = z0 + za 505 453 zccb = zcorl1(ji,jj) 506 u_ice(ji,jj) = (zr+zccb*zv_ice1)/(zcca+epsd)*zmask 507 454 u_ice(ji,jj) = ( zr + zccb * zv_ice1 ) / ( zcca + zepsi ) * rswitch 508 455 END DO 509 456 END DO … … 511 458 CALL lbc_lnk( u_ice(:,:), 'U', -1. ) 512 459 #if defined key_agrif && defined key_lim2 513 CALL agrif_rhg_lim2( jter, n evp, 'U' )460 CALL agrif_rhg_lim2( jter, nn_nevp, 'U' ) 514 461 #endif 515 462 #if defined key_bdy … … 517 464 #endif 518 465 519 !CDIR NOVERRCHK520 466 DO jj = k_j1+1, k_jpj-1 521 !CDIR NOVERRCHK522 467 DO ji = fs_2, fs_jpim1 523 468 524 zmask = (1.0-MAX(0._wp,SIGN(1._wp,-zmass2(ji,jj))))*tmv(ji,jj)525 z0 = zmass2(ji,jj) /dtevp469 rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass2(ji,jj) ) ) ) * vmask(ji,jj,1) 470 z0 = zmass2(ji,jj) * z1_dtevp 526 471 ! SB modif because ocean has no slip boundary condition 527 zu_ice2 = 0.5*( (u_ice(ji,jj)+u_ice(ji-1,jj))*e2t(ji,jj) & 528 & + (u_ice(ji,jj+1)+u_ice(ji-1,jj+1))*e2t(ji,jj+1)) & 529 & /(e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj) 530 za = rhoco*SQRT((zu_ice2-u_oce2(ji,jj))**2 + & 531 (v_ice(ji,jj)-v_oce2(ji,jj))**2)*(1.0-zfrld2(ji,jj)) 532 zr = z0*v_ice(ji,jj) + zf2(ji,jj) + & 533 za2ct(ji,jj) + za*(v_oce2(ji,jj)) 534 zcca = z0+za 472 zu_ice2 = 0.5 * ( ( u_ice(ji,jj ) + u_ice(ji-1,jj ) ) * e2t(ji,jj) & 473 & + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj+1) ) & 474 & / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1) 475 za = rhoco * SQRT( ( zu_ice2 - u_oce2(ji,jj) )**2 + & 476 & ( v_ice(ji,jj) - v_oce(ji,jj))**2 ) * ( 1.0 - zfrld2(ji,jj) ) 477 zr = z0 * v_ice(ji,jj) + zf2(ji,jj) + za2ct(ji,jj) + za * v_oce(ji,jj) 478 zcca = z0 + za 535 479 zccb = zcorl2(ji,jj) 536 v_ice(ji,jj) = (zr-zccb*zu_ice2)/(zcca+epsd)*zmask 537 480 v_ice(ji,jj) = ( zr - zccb * zu_ice2 ) / ( zcca + zepsi ) * rswitch 538 481 END DO 539 482 END DO … … 541 484 CALL lbc_lnk( v_ice(:,:), 'V', -1. ) 542 485 #if defined key_agrif && defined key_lim2 543 CALL agrif_rhg_lim2( jter, n evp, 'V' )486 CALL agrif_rhg_lim2( jter, nn_nevp, 'V' ) 544 487 #endif 545 488 #if defined key_bdy … … 548 491 549 492 ELSE 550 !CDIR NOVERRCHK551 493 DO jj = k_j1+1, k_jpj-1 552 !CDIR NOVERRCHK553 494 DO ji = fs_2, fs_jpim1 554 zmask = (1.0-MAX(0._wp,SIGN(1._wp,-zmass2(ji,jj))))*tmv(ji,jj)555 z0 = zmass2(ji,jj) /dtevp495 rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass2(ji,jj) ) ) ) * vmask(ji,jj,1) 496 z0 = zmass2(ji,jj) * z1_dtevp 556 497 ! SB modif because ocean has no slip boundary condition 557 zu_ice2 = 0.5*( (u_ice(ji,jj)+u_ice(ji-1,jj))*e2t(ji,jj) & 558 & +(u_ice(ji,jj+1)+u_ice(ji-1,jj+1))*e2t(ji,jj+1)) & 559 & /(e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj) 560 561 za = rhoco*SQRT((zu_ice2-u_oce2(ji,jj))**2 + & 562 (v_ice(ji,jj)-v_oce2(ji,jj))**2)*(1.0-zfrld2(ji,jj)) 563 zr = z0*v_ice(ji,jj) + zf2(ji,jj) + & 564 za2ct(ji,jj) + za*(v_oce2(ji,jj)) 565 zcca = z0+za 498 zu_ice2 = 0.5 * ( ( u_ice(ji,jj ) + u_ice(ji-1,jj ) ) * e2t(ji,jj) & 499 & +( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj+1) ) & 500 & / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1) 501 502 za = rhoco * SQRT( ( zu_ice2 - u_oce2(ji,jj) )**2 + & 503 & ( v_ice(ji,jj) - v_oce(ji,jj) )**2 ) * ( 1.0 - zfrld2(ji,jj) ) 504 zr = z0 * v_ice(ji,jj) + zf2(ji,jj) + za2ct(ji,jj) + za * v_oce(ji,jj) 505 zcca = z0 + za 566 506 zccb = zcorl2(ji,jj) 567 v_ice(ji,jj) = (zr-zccb*zu_ice2)/(zcca+epsd)*zmask 568 507 v_ice(ji,jj) = ( zr - zccb * zu_ice2 ) / ( zcca + zepsi ) * rswitch 569 508 END DO 570 509 END DO … … 572 511 CALL lbc_lnk( v_ice(:,:), 'V', -1. ) 573 512 #if defined key_agrif && defined key_lim2 574 CALL agrif_rhg_lim2( jter, n evp, 'V' )513 CALL agrif_rhg_lim2( jter, nn_nevp, 'V' ) 575 514 #endif 576 515 #if defined key_bdy … … 578 517 #endif 579 518 580 !CDIR NOVERRCHK581 519 DO jj = k_j1+1, k_jpj-1 582 !CDIR NOVERRCHK583 520 DO ji = fs_2, fs_jpim1 584 zmask = (1.0-MAX(0._wp,SIGN(1._wp,-zmass1(ji,jj))))*tmu(ji,jj) 585 z0 = zmass1(ji,jj)/dtevp 586 zv_ice1 = 0.5*( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji,jj) & 587 & +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji+1,jj)) & 588 & /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj) 589 590 za = rhoco*SQRT((u_ice(ji,jj)-u_oce1(ji,jj))**2 + & 591 (zv_ice1-v_oce1(ji,jj))**2)*(1.0-zfrld1(ji,jj)) 592 zr = z0*u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + & 593 za*(u_oce1(ji,jj)) 594 zcca = z0+za 521 rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass1(ji,jj) ) ) ) * umask(ji,jj,1) 522 z0 = zmass1(ji,jj) * z1_dtevp 523 zv_ice1 = 0.5 * ( ( v_ice(ji ,jj) + v_ice(ji ,jj-1) ) * e1t(ji,jj) & 524 & + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji+1,jj) ) & 525 & / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1) 526 527 za = rhoco * SQRT( ( u_ice(ji,jj) - u_oce(ji,jj) )**2 + & 528 & ( zv_ice1 - v_oce1(ji,jj) )**2 ) * ( 1.0 - zfrld1(ji,jj) ) 529 zr = z0 * u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + za * u_oce(ji,jj) 530 zcca = z0 + za 595 531 zccb = zcorl1(ji,jj) 596 u_ice(ji,jj) = ( zr+zccb*zv_ice1)/(zcca+epsd)*zmask597 END DO ! ji598 END DO ! jj532 u_ice(ji,jj) = ( zr + zccb * zv_ice1 ) / ( zcca + zepsi ) * rswitch 533 END DO 534 END DO 599 535 600 536 CALL lbc_lnk( u_ice(:,:), 'U', -1. ) 601 537 #if defined key_agrif && defined key_lim2 602 CALL agrif_rhg_lim2( jter, n evp, 'U' )538 CALL agrif_rhg_lim2( jter, nn_nevp, 'U' ) 603 539 #endif 604 540 #if defined key_bdy … … 611 547 !--- Convergence test. 612 548 DO jj = k_j1+1 , k_jpj-1 613 zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ) , & 614 ABS( v_ice(:,jj) - zv_ice(:,jj) ) ) 615 END DO 616 zresm = MAXVAL( zresr( 1:jpi , k_j1+1:k_jpj-1 ) ) 549 zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ), ABS( v_ice(:,jj) - zv_ice(:,jj) ) ) 550 END DO 551 zresm = MAXVAL( zresr( 1:jpi, k_j1+1:k_jpj-1 ) ) 617 552 IF( lk_mpp ) CALL mpp_max( zresm ) ! max over the global domain 618 553 ENDIF … … 625 560 ! 4) Prevent ice velocities when the ice is thin 626 561 !------------------------------------------------------------------------------! 627 ! If the ice thickness is below hminrhg (5cm) then ice velocity should equal the 628 ! ocean velocity, 629 ! This prevents high velocity when ice is thin 630 !CDIR NOVERRCHK 562 ! If the ice volume is below zvmin then ice velocity should equal the 563 ! ocean velocity. This prevents high velocity when ice is thin 631 564 DO jj = k_j1+1, k_jpj-1 632 !CDIR NOVERRCHK633 565 DO ji = fs_2, fs_jpim1 634 zdummy = vt_i(ji,jj) 635 IF ( zdummy .LE. hminrhg ) THEN 566 IF ( vt_i(ji,jj) <= zvmin ) THEN 636 567 u_ice(ji,jj) = u_oce(ji,jj) 637 568 v_ice(ji,jj) = v_oce(ji,jj) 638 ENDIF ! zdummy569 ENDIF 639 570 END DO 640 571 END DO 641 572 642 CALL lbc_lnk ( u_ice(:,:), 'U', -1. )643 CALL lbc_lnk( v_ice(:,:), 'V', -1. ) 573 CALL lbc_lnk_multi( u_ice(:,:), 'U', -1., v_ice(:,:), 'V', -1. ) 574 644 575 #if defined key_agrif && defined key_lim2 645 CALL agrif_rhg_lim2( n evp ,nevp, 'U' )646 CALL agrif_rhg_lim2( n evp ,nevp, 'V' )576 CALL agrif_rhg_lim2( nn_nevp , nn_nevp, 'U' ) 577 CALL agrif_rhg_lim2( nn_nevp , nn_nevp, 'V' ) 647 578 #endif 648 579 #if defined key_bdy … … 653 584 DO jj = k_j1+1, k_jpj-1 654 585 DO ji = fs_2, fs_jpim1 655 zdummy = vt_i(ji,jj) 656 IF ( zdummy .LE. hminrhg ) THEN 657 v_ice1(ji,jj) = 0.5*( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji+1,jj) & 658 & +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji,jj)) & 659 & /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj) 660 661 u_ice2(ji,jj) = 0.5*( (u_ice(ji,jj)+u_ice(ji-1,jj))*e2t(ji,jj+1) & 662 & +(u_ice(ji,jj+1)+u_ice(ji-1,jj+1))*e2t(ji,jj)) & 663 & /(e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj) 664 ENDIF ! zdummy 586 IF ( vt_i(ji,jj) <= zvmin ) THEN 587 v_ice1(ji,jj) = 0.5_wp * ( ( v_ice(ji ,jj) + v_ice(ji, jj-1) ) * e1t(ji+1,jj) & 588 & + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji ,jj) ) & 589 & / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1) 590 591 u_ice2(ji,jj) = 0.5_wp * ( ( u_ice(ji,jj ) + u_ice(ji-1,jj ) ) * e2t(ji,jj+1) & 592 & + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj ) ) &nbs