Changeset 834 for trunk/NEMO
- Timestamp:
- 2008-03-07T18:11:35+01:00 (16 years ago)
- Location:
- trunk/NEMO/LIM_SRC_3
- Files:
-
- 2 deleted
- 28 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMO/LIM_SRC_3/dom_ice.F90
r825 r834 6 6 !! History : 7 7 !! 2.0 ! 03-08 (C. Ethe) Free form and module 8 !! 3.0 ! 03-08 (M. Vancoppenolle), Ice cats, salinity and EVP-C 8 9 !!---------------------------------------------------------------------- 9 !! LIM 2.0, UCL-LOCEAN-IPSL (2005)10 !! LIM 3.0, UCL-ASTR-LOCEAN-IPSL (2005) 10 11 !! $Header: /home/opalod/NEMOCVSROOT/NEMO/LIM_SRC/dom_ice.F90,v 1.4 2005/03/27 18:34:40 opalod Exp $ 11 12 !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt -
trunk/NEMO/LIM_SRC_3/ice.F90
r825 r834 1 1 MODULE ice 2 !!======================================================================3 !! *** MODULE ice ***4 !! Sea Ice physics: diagnostics variables of ice defined in memory5 !!=====================================================================6 2 #if defined key_lim3 7 3 !!---------------------------------------------------------------------- 8 !! 'key_lim3' : LIM sea-ice model4 !! 'key_lim3' : LIM3 sea-ice model 9 5 !!---------------------------------------------------------------------- 10 6 !! History : 11 7 !! 2.0 ! 03-08 (C. Ethe) F90: Free form and module 8 !! 3.0 ! 08-03 (M. Vancoppenolle) : LIM3 ! 12 9 !!---------------------------------------------------------------------- 13 !! LIM 2.0, UCL-LOCEAN-IPSL (2005)10 !! LIM 3.0, UCL-LOCEAN-IPSL (2005) 14 11 !! $Header: /home/opalod/NEMOCVSROOT/NEMO/LIM_SRC/ice.F90,v 1.4 2005/03/27 18:34:41 opalod Exp $ 15 12 !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt … … 20 17 IMPLICIT NONE 21 18 PRIVATE 19 !! 20 !!====================================================================== 21 !! *** MODULE ice *** 22 !! 23 !! ************** 24 !! * L I M 3.0 * 25 !! ************** 26 !! 27 !! ''in ice we trust'' 28 !! 29 !! This module contains the sea ice 30 !! diagnostics variables of ice defined 31 !! in memory 32 !! 33 !!====================================================================== 34 !! 35 !! LIM3 by the use of sweat, agile fingers and sometimes brain juice, 36 !! was developed in Louvain-la-Neuve by : 37 !! 38 !! * Martin Vancoppenolle (UCL-ASTR, Belgium) 39 !! * Sylvain Bouillon (UCL-ASTR, Belgium) 40 !! * Miguel Angel Morales Maqueda (POL, UK) 41 !! 42 !! Based on extremely valuable earlier work by 43 !! 44 !! * Thierry Fichefet 45 !! * Hugues Goosse 46 !! 47 !! The following persons also contributed to the code in various ways 48 !! 49 !! * Gurvan Madec, Claude Talandier, Christian Ethe 50 !! and Rachid Benshila (LOCEAN-IPSL, France) 51 !! * Xavier Fettweis (UCL-ASTR), Ralph Timmermann (AWI, Germany) 52 !! * Bill Lipscomb (LANL), Cecilia Bitz (UWa) 53 !! and Elisabeth Hunke (LANL), USA. 54 !! 55 !! (c) UCL-ASTR, 2005-2008 56 !! 57 !! For more info, the interested user is kindly invited to consult the 58 !! following references 59 !! For model description and validation : 60 !! * Vancoppenolle et al., Ocean Modelling, 2008a. 61 !! * Vancoppenolle et al., Ocean Modelling, 2008b. 62 !! 63 !! For a specific description of EVP : 64 !! * Bouillon et al., in prep for 2008. 65 !! 66 !! Or the reference manual, that should be available by 2009 67 !! 68 !!====================================================================== 69 !! | 70 !! ***************************************** | 71 !! * * | 72 !! ************ I C E S T A T E V A R I A B L E S **************** | 73 !! * * | 74 !! ***************************************** | 75 !! | 76 !! Introduction : | 77 !! -------------- | 78 !! | 79 !! Every ice-covered grid cell is characterized by a series of state | 80 !! variables. To account for unresolved spatial variability in ice | 81 !! thickness, the ice cover in divided in ice thickness categories. | 82 !! | 83 !! Sea ice state variables depend on the ice thickness category | 84 !! | 85 !! Those variables are divided into two groups | 86 !! * Extensive (or global) variables. | 87 !! These are the variables that are transported by all means | 88 !! * Intensive (or equivalent) variables. | 89 !! These are the variables that are either physically more | 90 !! meaningful and/or used in ice thermodynamics | 91 !! | 92 !! Routines in limvar.F90 perform conversions | 93 !! - lim_var_glo2eqv : from global to equivalent variables | 94 !! - lim_var_eqv2glo : from equivalent to global variables | 95 !! | 96 !! For various purposes, the sea ice state variables have sometimes | 97 !! to be aggregated over all ice thickness categories. This operation | 98 !! is done in : | 99 !! - lim_var_agg | 100 !! | 101 !! in icestp.F90, the routines that compute the changes in the ice | 102 !! state variables are called | 103 !! - lim_dyn : ice dynamics | 104 !! - lim_trp : ice transport | 105 !! - lim_itd_me : mechanical redistribution (ridging and rafting) | 106 !! - lim_thd : ice halo-thermodynamics | 107 !! - lim_itd_th : thermodynamic changes in ice thickness distribution | 108 !! and creation of new ice | 109 !! | 110 !! See the associated routines for more information | 111 !! | 112 !! List of ice state variables : | 113 !! ----------------------------- | 114 !! | 115 !!-------------|-------------|---------------------------------|-------| 116 !! name in | name in | meaning | units | 117 !! 2D routines | 1D routines | | | 118 !!-------------|-------------|---------------------------------|-------| 119 !! | 120 !! ******************************************************************* | 121 !! *** Dynamical variables (prognostic) *** | 122 !! ******************************************************************* | 123 !! | 124 !! u_ice | - | Comp. U of the ice velocity | m/s | 125 !! v_ice | - | Comp. V of the ice velocity | m/s | 126 !! | 127 !! ******************************************************************* | 128 !! *** Category dependent state variables (prognostic) *** | 129 !! ******************************************************************* | 130 !! | 131 !! ** Global variables | 132 !! | 133 !!-------------|-------------|---------------------------------|-------| 134 !! a_i | a_i_b | Ice concentration | | 135 !! v_i | - | Ice volume per unit area | m | 136 !! v_s | - | Snow volume per unit area | m | 137 !! smv_i | - | Sea ice salt content | ppt.m | 138 !! oa_i ! - ! Sea ice areal age content | day | 139 !! e_i ! - ! Ice enthalpy | 10^9 J| 140 !! - ! q_i_b ! Ice enthalpy per unit vol. | J/m3 | 141 !! e_s ! - ! Snow enthalpy | 10^9 J| 142 !! - ! q_s_b ! Snow enthalpy per unit vol. | J/m3 | 143 !! | 144 !!-------------|-------------|---------------------------------|-------| 145 !! | 146 !! ** Equivalent variables | 147 !! | 148 !!-------------|-------------|---------------------------------|-------| 149 !! | 150 !! ht_i | ht_i_b | Ice thickness | m | 151 !! ht_s ! ht_s_b | Snow depth | m | 152 !! sm_i ! sm_i_b | Sea ice bulk salinity ! ppt | 153 !! s_i ! s_i_b | Sea ice salinity profile ! ppt | 154 !! o_i ! - | Sea ice Age ! days | 155 !! t_i ! t_i_b | Sea ice temperature ! K | 156 !! t_s ! t_s_b | Snow temperature ! K | 157 !! t_su ! t_su_b | Sea ice surface temperature ! K | 158 !! | 159 !! notes: the ice model only sees a bulk (i.e., vertically averaged) | 160 !! salinity, except in thermodynamic computations, for which | 161 !! the salinity profile is computed as a function of bulk | 162 !! salinity | 163 !! | 164 !! the sea ice surface temperature is not associated to any | 165 !! heat content. Therefore, it is not a state variable and | 166 !! does not have to be advected. Nevertheless, it has to be | 167 !! computed to determine whether the ice is melting or not | 168 !! | 169 !! ******************************************************************* | 170 !! *** Category-summed state variables (diagnostic) *** | 171 !! ******************************************************************* | 172 !! at_i | at_i_b | Total ice concentration | | 173 !! vt_i | - | Total ice vol. per unit area | m | 174 !! vt_s | - | Total snow vol. per unit ar. | m | 175 !! smt_i | - | Mean sea ice salinity | ppt | 176 !! tm_i | - | Mean sea ice temperature | K | 177 !! ot_i ! - ! Sea ice areal age content | day | 178 !! et_i ! - ! Total ice enthalpy | 10^9 J| 179 !! et_s ! - ! Total snow enthalpy | 10^9 J| 180 !! bv_i ! - ! Mean relative brine volume | ??? | 181 !! | 182 !! | 183 !!===================================================================== 22 184 23 185 LOGICAL, PUBLIC :: & 24 186 con_i = .false. ! switch for conservation test 25 187 188 !!-------------------------------------------------------------------------- 26 189 !! * Share Module variables 27 INTEGER , PUBLIC :: & !!: ** ice-dynamic namelist (namicedyn) ** 28 nbiter = 1 , & !: number of sub-time steps for relaxation 29 nbitdr = 250 , & !: maximum number of iterations for relaxation 30 nevp = 400 , & !: number of iterations for subcycling 31 nlay_i = 5 !: number of layers in the ice 32 33 REAL(wp), PUBLIC :: & !!: ** ice-dynamic namelist (namicedyn) ** 34 epsd = 1.0e-20, & !: tolerance parameter for dynamic 35 alpha = 0.5 , & !: coefficient for semi-implicit coriolis 36 dm = 0.6e+03, & !: diffusion constant for dynamics 37 om = 0.5 , & !: relaxation constant 38 resl = 5.0e-05, & !: maximum value for the residual of relaxation 39 cw = 5.0e-03, & !: drag coefficient for oceanic stress 40 angvg = 0.0 , & !: turning angle for oceanic stress 41 pstar = 1.0e+04, & !: first bulk-rheology parameter 42 c_rhg = 20.0 , & !: second bulk-rhelogy parameter 43 etamn = 0.0e+07, & !: minimun value for viscosity 44 creepl = 2.0e-08, & !: creep limit 190 !!-------------------------------------------------------------------------- 191 INTEGER , PUBLIC :: & !!: ** ice-dynamic namelist (namicedyn) ** 192 nbiter = 1 , & !: number of sub-time steps for relaxation 193 nbitdr = 250 , & !: maximum number of iterations for relaxation 194 nevp = 400 , & !: number of iterations for subcycling 195 nlay_i = 5 !: number of layers in the ice 196 197 REAL(wp), PUBLIC :: & !!: ** ice-dynamic namelist (namicedyn) ** 198 epsd = 1.0e-20, & !: tolerance parameter for dynamic 199 alpha = 0.5 , & !: coefficient for semi-implicit coriolis 200 dm = 0.6e+03, & !: diffusion constant for dynamics 201 om = 0.5 , & !: relaxation constant 202 resl = 5.0e-05, & !: maximum value for the residual of relaxation 203 cw = 5.0e-03, & !: drag coefficient for oceanic stress 204 angvg = 0.0 , & !: turning angle for oceanic stress 205 pstar = 1.0e+04, & !: determines ice strength (N/M), Hibler JPO79 206 c_rhg = 20.0 , & !: determines changes in ice strength 207 etamn = 0.0e+07, & !: minimun value for viscosity : has to be 0 208 creepl = 2.0e-08, & !: creep limit : has to be under 1.0e-9 45 209 ecc = 2.0 , & !: eccentricity of the elliptical yield curve 46 210 ahi0 = 350.e0 , & !: sea-ice hor. eddy diffusivity coeff. (m2/s) 47 ! EVP148 211 telast = 2880.0 , & !: timescale for elastic waves (s) !SB 49 212 alphaevp = 1.0 , & !: coeficient of the internal stresses !SB … … 51 214 52 215 REAL(wp), PUBLIC :: & !!: ** ice-salinity namelist (namicesal) ** 53 s_i_max = 20.0 , & !: maximum ice salinity 54 s_i_min = 0.1 , & !: minimum ice salinity 216 s_i_max = 20.0 , & !: maximum ice salinity (ppt) 217 s_i_min = 0.1 , & !: minimum ice salinity (ppt) 55 218 s_i_0 = 3.5 , & !: 1st sal. value for the computation of sal .prof. 219 !: (ppt) 56 220 s_i_1 = 4.5 , & !: 2nd sal. value for the computation of sal .prof. 221 !: (ppt) 57 222 sal_G = 5.00 , & !: restoring salinity for gravity drainage 223 !: (ppt) 58 224 sal_F = 2.50 , & !: restoring salinity for flushing 59 time_G = 1.728e+06,&!: restoring time constant for gravity drainage 60 time_F = 1.728e+06,&!: restoring time constant for gravity drainage 61 bulk_sal = 4.0 !!: bulk salinity in case of constant salinity 225 !: (ppt) 226 time_G = 1.728e+06,&!: restoring time constant for gravity drainage 227 !: (= 20 days, in s) 228 time_F = 8.640e+05,&!: restoring time constant for gravity drainage 229 !: (= 10 days, in s) 230 bulk_sal = 4.0 !: bulk salinity (ppt) in case of constant salinity 62 231 63 232 INTEGER , PUBLIC :: & !!: ** ice-salinity namelist (namicesal) ** 64 num_sal = 1 , & !: 233 num_sal = 1 , & !: salinity configuration used in the model 234 !: 1 - s constant in space and time 235 !: 2 - prognostic salinity (s(z,t)) 236 !: 3 - salinity profile, constant in time 237 !: 4 - salinity variations affect only ice 238 ! thermodynamics 65 239 sal_prof = 1 , & !: salinity profile or not 66 thcon_i_swi = 1 !: thermal conductivity of Untersteiner (1964) or67 !: Pringle et al (2007) 240 thcon_i_swi = 1 !: thermal conductivity of Untersteiner (1964) (1) or 241 !: Pringle et al (2007) (2) 68 242 69 243 REAL(wp), PUBLIC :: & !!: ** ice-mechanical redistribution namelist (namiceitdme) … … 74 248 Gstar = 0.15 , & !!: fractional area of young ice contributing to ridging 75 249 astar = 0.05 , & !!: equivalent of G* for an exponential participation function 76 Hstar = 100.0 , & !!: fractional area of young ice contributing to ridging 77 hparmeter = 0.75, & !!: threshold thickness for rafting / ridging 78 Craft = 5.0 , & !!: coefficient for the hyperbolic tangent in rafting 79 ridge_por = 0.0 , & !!: initial porosity of ridges 80 sal_max_ridge = 15.0, & !!: maximum ridged ice salinity 250 Hstar = 100.0 , & !!: thickness that determines the maximal thickness of ridged 251 !!: ice 252 hparmeter = 0.75, & !!: threshold thickness (m) for rafting / ridging 253 Craft = 5.0 , & !!: coefficient for smoothness of the hyperbolic tangent in rafting 254 ridge_por = 0.0 , & !!: initial porosity of ridges (0.3 regular value) 255 sal_max_ridge = 15.0, & !!: maximum ridged ice salinity (ppt) 81 256 betas = 1.0 , & !:: coef. for partitioning of snowfall between leads and sea ice 82 257 kappa_i = 1.0 , & !!: coefficient for the extinction of radiation 83 nconv_i_thd = 50 , & !!: maximal number of iterations in heat diffusion 84 maxer_i_thd = 1.0e-4 !!: maximal tolerated error for heat diffusion 258 !!: Grenfell et al. (2006) (m-1) 259 nconv_i_thd = 50 , & !!: maximal number of iterations for heat diffusion 260 maxer_i_thd = 1.0e-4 !!: maximal tolerated error (C) for heat diffusion 85 261 86 262 INTEGER , PUBLIC :: & !!: ** ice-mechanical redistribution namelist (namiceitdme) 87 263 ridge_scheme_swi = 0, & !!: scheme used for ice ridging 88 264 raftswi = 1, & !!: rafting of ice or not 89 partfun_swi = 1, & !!: participation function TH75 (0) or Letal07 (1) 90 transfun_swi = 0, & !!: transfer function of H80 (0) or Letal07 (1) 265 partfun_swi = 1, & !!: participation function Thorndike et al. JGR75 (0) 266 !!: or Lipscomb et al. JGR07 (1) 267 transfun_swi = 0, & !!: transfer function of Hibler, MWR80 (0) 268 !!: or Lipscomb et al., 2007 (1) 91 269 brinstren_swi = 0 !!: use brine volume to diminish ice strength 92 270 … … 162 340 163 341 INTEGER, PUBLIC, DIMENSION(jpi, jpj, jpl) :: & !:: 164 patho_case ! number of the pathological case (if any) 165 166 !--------------------------------------------------------------------------------------------------- 167 ! ) Ice global state variables 168 !--------------------------------------------------------------------------------------------------- 342 patho_case ! number of the pathological case (if any, of course) 343 344 !!-------------------------------------------------------------------------- 345 !! * Ice global state variables 346 !!-------------------------------------------------------------------------- 347 !! Variables defined for each ice category 169 348 REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpl) :: & !: 170 ht_i , & !: Ice thickness 171 a_i , & !: Ice fractional areas 172 v_i , & !: Ice volume per unit area 173 v_s , & !: Snow volume per unit area 174 ht_s , & !: Snow thickness 175 t_su , & !: Sea-Ice Surface Temperature [K] 176 sm_i , & !: Sea-Ice Bulk salinity [ppt] 177 smv_i , & !: Sea-Ice Bulk salinity [ppt] times volume 349 ht_i , & !: Ice thickness (m) 350 a_i , & !: Ice fractional areas (concentration) 351 v_i , & !: Ice volume per unit area (m) 352 v_s , & !: Snow volume per unit area(m) 353 ht_s , & !: Snow thickness (m) 354 t_su , & !: Sea-Ice Surface Temperature (K) 355 sm_i , & !: Sea-Ice Bulk salinity (ppt) 356 smv_i , & !: Sea-Ice Bulk salinity times volume per area (ppt.m) 357 !: this is an extensive variable that has to be transported 178 358 o_i , & !: Sea-Ice Age (days) 179 ov_i , & !: Sea-Ice Age times volume /area (days.m)359 ov_i , & !: Sea-Ice Age times volume per area (days.m) 180 360 oa_i !: Sea-Ice Age times ice area (days) 181 361 362 !! Variables summed over all categories, or associated to 363 !! all the ice in a single grid cell 182 364 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: & !: 183 365 u_ice, v_ice, & !: two components of the ice velocity (m/s) 184 366 tio_u, tio_v, & !: two components of the ice-ocean stress (N/m2) 185 vt_i , & !: ice total volume 186 vt_s , & !: snow total volume 187 at_i , & !: ice total fractional area 188 ato_i , & !: total open water fractional area 367 vt_i , & !: ice total volume per unit area (m) 368 vt_s , & !: snow total volume per unit area (m) 369 at_i , & !: ice total fractional area (ice concentration) 370 ato_i , & !: total open water fractional area (1-at_i) 371 et_i , & !: total ice heat content 372 et_s , & !: total snow heat content 189 373 ot_i , & !: mean age over all categories 190 et_i , & !: total ice heat content 191 bv_i , & !: total ice heat content 192 tm_i , & !: mean ice temperature 193 et_s , & !: total snow heat content 194 smt_i 374 tm_i , & !: mean ice temperature over all categories 375 bv_i , & !: brine volume averaged over all categories 376 smt_i !: mean sea ice salinity averaged over all categories 195 377 196 378 REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpm) :: & !: … … 216 398 reslum !: Relative absorption of solar radiation in each ocean level 217 399 218 !--------------------------------------------------------------------------------------------------- 219 ! ) Moments for Advection 220 !--------------------------------------------------------------------------------------------------- 221 400 !!-------------------------------------------------------------------------- 401 !! * Moments for advection 402 !!-------------------------------------------------------------------------- 222 403 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: & !: 223 404 sxopw, syopw, sxxopw, syyopw, sxyopw !: open water in sea ice … … 234 415 sxe , sye , sxxe , syye , sxye !: ice layers heat content 235 416 236 !---------------------------------------------------------------------------------------------------237 ! )Old values of global variables238 !---------------------------------------------------------------------------------------------------417 !!-------------------------------------------------------------------------- 418 !! * Old values of global variables 419 !!-------------------------------------------------------------------------- 239 420 REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpl) :: & !: 240 421 old_v_s, old_v_i, & !: snow and ice volumes … … 247 428 old_u_ice, old_v_ice 248 429 249 !--------------------------------------------------------------------------------------------------- 250 ! ) Values of increments of global variables 251 !--------------------------------------------------------------------------------------------------- 252 ! thd refers to changes induced by thermodynamics 253 ! trp '' '' '' advection (transport of ice) 254 430 !!-------------------------------------------------------------------------- 431 !! * Increment of global variables 432 !!-------------------------------------------------------------------------- 433 ! thd refers to changes induced by thermodynamics 434 ! trp '' '' '' advection (transport of ice) 255 435 REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpl) :: & !: 256 436 d_a_i_thd , d_a_i_trp , & !: icefractions … … 271 451 d_u_ice_dyn, d_v_ice_dyn 272 452 273 !---------------------------------------------------------------------------------------------------274 ! )Ice thickness distribution variables275 !---------------------------------------------------------------------------------------------------276 453 !!-------------------------------------------------------------------------- 454 !! * Ice thickness distribution variables 455 !!-------------------------------------------------------------------------- 456 ! REMOVE 277 457 INTEGER(wp), PUBLIC, DIMENSION(jpl) :: & !: 278 458 ice_types !: Vector making the connection between types and categories 279 459 280 460 INTEGER(wp), PUBLIC, DIMENSION(jpm,2) :: & !: 281 ice_cat_bounds !: Matrix containing the integer upper and lower boundaries of ice thickness categories 282 461 ice_cat_bounds !: Matrix containing the integer upper and 462 !: lower boundaries of ice thickness categories 463 464 ! REMOVE 283 465 INTEGER(wp), PUBLIC, DIMENSION(jpm) :: & !: 284 466 ice_ncat_types !: Vector containing the number of thickness categories in each ice type … … 290 472 hi_mean !: Mean ice thickness in catgories 291 473 474 ! REMOVE 292 475 REAL(wp), PUBLIC, DIMENSION(0:jpl,jpm) :: & !: 293 hi_max_typ !: Boundary of ice thickness categories in thickness space (same but specific for each ice type) 294 295 !--------------------------------------------------------------------------------------------------- 296 ! ) Ice diagnostics 297 !--------------------------------------------------------------------------------------------------- 476 hi_max_typ !: Boundary of ice thickness categories 477 !:in thickness space (same but specific for each ice type) 478 479 !!-------------------------------------------------------------------------- 480 !! * Ice diagnostics 481 !!-------------------------------------------------------------------------- 482 !! Check if everything down here is necessary 298 483 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: & !: volume of ice formed in the leads 299 484 v_newice 300 485 REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpl) :: & !: thermodynamic growth rates 301 486 dv_dt_thd, & 302 izero, fstroc, fhbricat ! to remove 303 487 izero, fstroc, fhbricat 304 488 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: & 305 489 diag_sni_gr, & ! snow ice growth … … 309 493 diag_bot_me, & ! vertical bottom melt 310 494 diag_sur_me ! vertical surface melt 311 312 495 INTEGER , PUBLIC :: & !: indexes of the debugging 313 496 jiindex, & ! point -
trunk/NEMO/LIM_SRC_3/iceini.F90
r830 r834 22 22 USE par_ice 23 23 USE limvar 24 USE limicepoints25 24 26 25 IMPLICIT NONE … … 51 50 !!---------------------------------------------------------------------- 52 51 !! LIM 2.0, UCL-LOCEAN-IPSL (2005) 52 !! LIM 3.0, UCL-ASTR-LOCEAN-IPSL (2008) 53 53 !! $Header: /home/opalod/NEMOCVSROOT/NEMO/LIM_SRC/iceini.F90,v 1.4 2005/03/27 18:34:41 opalod Exp $ 54 54 !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt … … 64 64 !! 65 65 !! History : 66 !! 8.5 ! 02-08 (G. Madec) F90: Free form and modules 66 !! 2.0 ! 02-08 (G. Madec) F90: Free form and modules 67 !! 3.0 ! 08-03 (M. Vancop) ITD, salinity, EVP-C 67 68 !!---------------------------------------------------------------------- 68 69 INTEGER :: ji,jj,jk,jl, index70 69 71 70 ! Open the namelist file … … 73 72 74 73 CALL ice_run ! read in namelist some run parameters 75 76 CALL lim_icepoints ! define specific checking ice points77 74 78 75 ! Louvain la Neuve Ice model … … 87 84 CALL lim_msh ! ice mesh initialization 88 85 89 CALL lim_itd_ini 86 CALL lim_itd_ini ! initialize the ice thickness 87 ! distribution 90 88 ! Initial sea-ice state 91 89 IF( .NOT.ln_rstart ) THEN 92 90 numit = 0 93 ! martin ajoute94 91 numit = nit000 - 1 95 92 CALL lim_istate ! start from rest: sea-ice deduced from sst 96 CALL lim_var_agg(1) 97 CALL lim_var_glo2eqv 93 CALL lim_var_agg(1) ! aggregate category variables in 94 ! bulk variables 95 CALL lim_var_glo2eqv ! convert global variables in equivalent 96 ! variables 98 97 ELSE 99 98 CALL lim_rst_read( numit ) ! start from a restart file 100 ! martin ajoute101 99 numit = nit000 - 1 102 CALL lim_var_agg(1) 103 CALL lim_var_glo2eqv 100 CALL lim_var_agg(1) ! aggregate ice variables 101 CALL lim_var_glo2eqv ! convert global var in equivalent variables 104 102 ENDIF 105 103 … … 130 128 !! history : 131 129 !! 2.0 ! 03-08 (C. Ethe) Original code 130 !! 3.0 ! 08-03 (M. Vancop) LIM3 132 131 !!------------------------------------------------------------------- 133 132 NAMELIST/namicerun/ ln_limdyn, acrit, hsndif, hicdif, cai, cao, ln_nicep … … 158 157 !! Initializes the ice thickness distribution 159 158 !! ** Method : 160 !! Very simple 159 !! Very simple. Currently there are no ice types in the 160 !! model... 161 161 !! 162 162 !! ** Arguments : … … 172 172 !! ** History : 173 173 !! (12-2005) Martin Vancoppenolle 174 !! Rien n'est jamais acquis175 !! A l'homme ni sa force176 !! Ni sa faiblesse ni son coeur177 !! Et quand il croit178 !! Ouvrir ses bras son ombre179 !! Est celle d'une croix180 174 !! 181 175 !!------------------------------------------------------------------ … … 183 177 184 178 !! * Local variables 185 INTEGER :: ji, & ! spatial dummy loop index 186 jl, & ! ice category dummy loop index 179 INTEGER :: jl, & ! ice category dummy loop index 187 180 jm ! ice types dummy loop index 188 181 … … 198 191 199 192 !!-- End of declarations 200 !!---------------------------------------------------------------------------------------------- 201 202 !--------------------------------------------------------------------------------------------------! 203 ! 1) Ice thickness distribution parameters initialization ! 204 !--------------------------------------------------------------------------------------------------! 205 206 !- Types boundaries in integer thickness space 207 !---------------------------------------------- 208 !- Type 1 (undeformed ice) Boundaries 193 !!------------------------------------------------------------------------------ 194 195 !------------------------------------------------------------------------------! 196 ! 1) Ice thickness distribution parameters initialization 197 !------------------------------------------------------------------------------! 198 199 !- Types boundaries (integer) 200 !---------------------------- 209 201 ice_cat_bounds(1,1) = 1 210 202 ice_cat_bounds(1,2) = jpl 211 ! !- Type 2 (ridged ice ) Boundaries212 ! ice_cat_bounds(2,1) = 4213 ! ice_cat_bounds(2,2) = 5214 ! !- Type 3 (rafted ice ) Boundaries215 ! ice_cat_bounds(3,1) = 6216 ! ice_cat_bounds(3,2) = 6217 203 218 204 !- Number of ice thickness categories in each ice type … … 223 209 !- Make the correspondence between thickness categories and ice types 224 210 !--------------------------------------------------------------------- 225 !- ice_types = 1 -> undeformed ice226 !- ice_types = 2 -> ridged ice227 !- ice_types = 3 -> rafted ice228 211 DO jm = 1, jpm !over types 229 212 DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2) !over thickness categories … … 242 225 WRITE(numout,*) 243 226 244 !- Thickness categories boundaries in thickness space 245 !--------------------------------------------------------------------- 246 247 ! NEW DEFINITION 248 249 hi_max(:) = 0.0 250 hi_max_typ(:,:) = 0.0 251 252 ! new categories 253 ! !- Type 3 - rafted ice 254 ! hi_max(ice_cat_bounds(3,2)) = 5.0 255 256 ! !- Type 2 - ridged ice 257 ! zc1 = 3./REAL(ice_cat_bounds(2,2)-ice_cat_bounds(2,1)+1) 258 ! zc2 = 10.0*zc1 259 ! zc3 = 3.0 260 261 ! DO jl = ice_cat_bounds(2,1),ice_cat_bounds(2,2) 262 ! zx1 = REAL(jl-ice_cat_bounds(2,1)+1) / REAL(ice_cat_bounds(2,2)-ice_cat_bounds(2,1)+1) 263 ! hi_max(jl) = hi_max(jl-1) + zc1 + zc2 * (1.0 + TANH ( zc3 * (zx1 - 1.0 ) ) ) 264 ! END DO 265 266 ! ! force the thickness into the ice categories 267 ! hi_max(ice_cat_bounds(2,1)) = 15.0 268 ! hi_max(ice_cat_bounds(2,2)) = 25.0 227 !- Thickness categories boundaries 228 !---------------------------------- 229 hi_max(:) = 0.0 230 hi_max_typ(:,:) = 0.0 269 231 270 232 !- Type 1 - undeformed ice 271 272 273 274 275 276 277 278 233 zc1 = 3./REAL(ice_cat_bounds(1,2)-ice_cat_bounds(1,1)+1) 234 zc2 = 10.0*zc1 235 zc3 = 3.0 236 237 DO jl = ice_cat_bounds(1,1), ice_cat_bounds(1,2) 238 zx1 = REAL(jl-1) / REAL(ice_cat_bounds(1,2)-ice_cat_bounds(1,1)+1) 239 hi_max(jl) = hi_max(jl-1) + zc1 + zc2 * (1.0 + TANH ( zc3 * (zx1 - 1.0 ) ) ) 240 END DO 279 241 280 242 !- Fill in the hi_max_typ vector, useful in other circumstances 281 DO jl = ice_cat_bounds(1,1), ice_cat_bounds(1,2) 282 hi_max_typ(jl,1) = hi_max(jl) 283 END DO 284 ! DO jl = ice_cat_bounds(2,1), ice_cat_bounds(2,2) 285 ! hi_max_typ(jl-ice_cat_bounds(2,1)+1,2) = hi_max(jl) 286 ! END DO 287 ! DO jl = ice_cat_bounds(3,1), ice_cat_bounds(3,2) 288 ! hi_max_typ(jl-ice_cat_bounds(3,1)+1,3) = hi_max(jl) 289 ! END DO 290 291 WRITE(numout,*) ' Thickness category boundaries independently of ice type ' 292 WRITE(numout,*) ' hi_max ', hi_max(0:jpl) 293 294 WRITE(numout,*) ' Thickness category boundaries inside ice types ' 295 DO jm = 1, jpm 296 WRITE(numout,*) ' Type number ', jm 297 WRITE(numout,*) ' hi_max_typ : ', hi_max_typ(0:ice_ncat_types(jm),jm) 298 END DO 299 300 DO jl = 1, jpl 301 hi_mean(jl) = ( hi_max(jl) + hi_max(jl-1) ) / 2.0 302 END DO 303 304 tn_ice(:,:,:) = t_su(:,:,:) 243 DO jl = ice_cat_bounds(1,1), ice_cat_bounds(1,2) 244 hi_max_typ(jl,1) = hi_max(jl) 245 END DO 246 247 WRITE(numout,*) ' Thickness category boundaries independently of ice type ' 248 WRITE(numout,*) ' hi_max ', hi_max(0:jpl) 249 250 WRITE(numout,*) ' Thickness category boundaries inside ice types ' 251 DO jm = 1, jpm 252 WRITE(numout,*) ' Type number ', jm 253 WRITE(numout,*) ' hi_max_typ : ', hi_max_typ(0:ice_ncat_types(jm),jm) 254 END DO 255 256 DO jl = 1, jpl 257 hi_mean(jl) = ( hi_max(jl) + hi_max(jl-1) ) / 2.0 258 END DO 259 260 tn_ice(:,:,:) = t_su(:,:,:) 305 261 306 262 END SUBROUTINE lim_itd_ini -
trunk/NEMO/LIM_SRC_3/icestp.F90
r825 r834 6 6 #if defined key_lim3 7 7 !!---------------------------------------------------------------------- 8 !! 'key_lim3' : L imsea-ice model8 !! 'key_lim3' : LIM3 sea-ice model 9 9 !!---------------------------------------------------------------------- 10 10 !! ice_stp : sea-ice model time-stepping … … 27 27 USE limdyn 28 28 USE limtrp 29 USE limicepoints30 29 USE limthd 31 30 USE limitd_th … … 102 101 !+++++ 103 102 index = 12 104 jiindex = arc_sp_grid(index,1)105 jjindex = arc_sp_grid(index,2)106 103 jiindex = 44 107 104 jjindex = 140 -
trunk/NEMO/LIM_SRC_3/limadv.F90
r825 r834 6 6 #if defined key_lim3 7 7 !!---------------------------------------------------------------------- 8 !! 'key_lim3' LIM sea-ice model8 !! 'key_lim3' LIM3 sea-ice model 9 9 !!---------------------------------------------------------------------- 10 10 !! lim_adv_x : advection of sea ice on x axis -
trunk/NEMO/LIM_SRC_3/limcons.F90
r825 r834 1 1 MODULE limcons 2 #if defined key_lim3 3 !!---------------------------------------------------------------------- 4 !! 'key_lim3' : LIM3 sea-ice model 5 !!---------------------------------------------------------------------- 6 !! 2 7 !!====================================================================== 3 8 !! *** MODULE limcons *** … … 10 15 !! 11 16 !!====================================================================== 12 !!----------------------------------------------------------------------13 !! 'key_lim3' : LIM sea-ice model14 !!----------------------------------------------------------------------15 17 !! lim_cons : checks whether energy/mass are conserved 16 18 !!---------------------------------------------------------------------- 19 !! 17 20 !! * Modules used 18 21 … … 33 36 34 37 !! * Module variables 35 36 !!---------------------------------------------------------------------- 37 !! LIM 2.1 , UCL-ASTR-LODYC-IPSL (2004) 38 !!---------------------------------------------------------------------- 39 !! LIM 3.0 , UCL-ASTR-LODYC-IPSL (2008) 38 40 !!---------------------------------------------------------------------- 39 41 … … 222 224 END SUBROUTINE lim_cons_check 223 225 224 !============================================================================== 225 226 END MODULE limcons 226 #else 227 !!---------------------------------------------------------------------- 228 !! Default option Empty module NO LIM sea-ice model 229 !!---------------------------------------------------------------------- 230 #endif 231 !!====================================================================== 232 END MODULE limcons -
trunk/NEMO/LIM_SRC_3/limdia.F90
r825 r834 10 10 #if defined key_lim3 11 11 !!---------------------------------------------------------------------- 12 !! 'key_lim3' : LIM sea-ice model12 !! 'key_lim3' : LIM3 sea-ice model 13 13 !!---------------------------------------------------------------------- 14 14 !! lim_dia : computation of the time evolution of keys var. … … 95 95 !! * Local variables 96 96 INTEGER :: jv,ji,jj,jl ! dummy loop indices 97 INTEGER :: nv ! indice of variable98 97 REAL(wp), DIMENSION(jpinfmx) :: & 99 98 vinfor ! temporary working space … … 138 137 vinfor(7) = vinfor(7) + vt_i(ji,jj)*aire(ji,jj) / 1.0e12 !ice volume 139 138 vinfor(9) = vinfor(9) + vt_s(ji,jj)*aire(ji,jj) / 1.0e12 !snow volume 140 ! vinfor(11) = vinfor(11) + vt_i(ji,jj)*aire(ji,jj) / 1.0e12 !undef ice volume141 ! vinfor(13) = vinfor(13) + vt_i(ji,jj)*aire(ji,jj) / 1.0e12 !def ice volume142 139 vinfor(15) = vinfor(15) + ot_i(ji,jj) *vt_i(ji,jj)*aire(ji,jj) / 1.0e12 !mean age 143 140 vinfor(29) = vinfor(29) + smt_i(ji,jj)*vt_i(ji,jj)*aire(ji,jj) / 1.0e12 !mean salinity 141 ! the computation of this diagnostic is not reliable 144 142 vinfor(31) = vinfor(31) + vt_i(ji,jj)*( u_ice(ji,jj)*u_ice(ji,jj) + & 145 v_ice(ji,jj)*v_ice(ji,jj) )*aire(ji,jj)/1.0e12 !ice vel143 v_ice(ji,jj)*v_ice(ji,jj) )*aire(ji,jj)/1.0e12 146 144 vinfor(53) = vinfor(53) + fsalt(ji,jj)*aire(ji,jj) / 1.0e12 !salt flux 147 145 vinfor(55) = vinfor(55) + fsbri(ji,jj)*aire(ji,jj) / 1.0e12 !brine drainage flux 148 146 vinfor(57) = vinfor(57) + fseqv(ji,jj)*aire(ji,jj) / 1.0e12 !equivalent salt flux 149 !slightly modified this150 147 vinfor(59) = vinfor(59) + sst_io(ji,jj)*at_i(ji,jj)*aire(ji,jj) / 1.0e12 !SST 151 148 vinfor(61) = vinfor(61) + sss_io(ji,jj)*at_i(ji,jj)*aire(ji,jj) / 1.0e12 !SSS … … 157 154 vinfor(75) = vinfor(75) + v_i(ji,jj,4)*aire(ji,jj) / 1.0e12 !ice volume 158 155 vinfor(77) = vinfor(77) + v_i(ji,jj,5)*aire(ji,jj) / 1.0e12 !ice volume 159 ! vinfor(79) = vinfor(79) + v_i(ji,jj,6)*aire(ji,jj) / 1.0e12 !ice volume160 156 vinfor(79) = 0.0 161 157 vinfor(81) = vinfor(81) + fmass(ji,jj)*aire(ji,jj) / 1.0e12 ! mass flux 162 ! vinfor(83) = vinfor(83) + ht_i(ji,jj,5)*aire(ji,jj) / 1.0e12 ! mass flux163 158 ENDIF 164 159 END DO … … 173 168 END DO 174 169 175 ! DO jl = ice_cat_bounds(2,1), ice_cat_bounds(3,2)176 ! DO jj = njeq, jpjm1177 ! DO ji = fs_2, fs_jpim1 ! vector opt.178 ! vinfor(13) = vinfor(13) + v_i(ji,jj,jl)*aire(ji,jj) / 1.0e12 !def ice volume (inc. rafted ice)179 ! END DO180 ! END DO181 ! END DO182 170 vinfor(13) = 0.0 183 171 184 172 vinfor(15) = vinfor(15) / vinfor(7) ! these have to be divided by total ice volume to have the 185 173 vinfor(29) = vinfor(29) / vinfor(7) ! right value 186 ! vinfor(31) = vinfor(31) / vinfor(7)187 174 vinfor(31) = SQRT( vinfor(31) / MAX( vinfor(7) , epsi06 ) ) 188 175 vinfor(67) = vinfor(67) / vinfor(7) … … 192 179 vinfor(57) = vinfor(57) / vinfor(5) ! 193 180 vinfor(79) = vinfor(79) / vinfor(5) ! 194 ! vinfor(83) = vinfor(83) / vinfor(5) !195 181 196 182 zindb = 1.0 - MAX(0.0,SIGN(1.0,-vinfor(3))) ! … … 273 259 !! Fram strait in ORCA2 = 5 points 274 260 !! export = -v_ice*e1t*ddtb*at_i or -v_ice*e1t*ddtb*at_i*h_i 275 ! jj = 137276 261 jj = 136 ! C grid 277 262 vinfor(83) = 0.0 278 263 vinfor(84) = 0.0 279 ! WRITE(numout,*) ' Fram Strait Export '280 264 DO ji = 134, 138 281 ! vinfor(83) = vinfor(83) - ( v_ice(ji,jj) + v_ice(ji+1,jj) ) / 2.0 * &282 ! e1t(ji,jj)*at_i(ji,jj)*rdt_ice / 1.0e12283 ! vinfor(84) = vinfor(84) - ( v_ice(ji,jj) + v_ice(ji+1,jj) ) / 2.0 * &284 ! e1t(ji,jj)*vt_i(ji,jj)*rdt_ice / 1.0e12285 ! C grid286 265 vinfor(83) = vinfor(83) - v_ice(ji,jj) * & 287 266 e1t(ji,jj)*at_i(ji,jj)*rdt_ice / 1.0e12 288 267 vinfor(84) = vinfor(84) - v_ice(ji,jj) * & 289 268 e1t(ji,jj)*vt_i(ji,jj)*rdt_ice / 1.0e12 290 ! WRITE(numout,*)291 ! WRITE(numout,*) ' ji , jj : ', ji, jj292 ! WRITE(numout,*) ' v_ice ji -jj : ', v_ice(ji,jj)293 ! WRITE(numout,*) ' v_ice ji+1-jj : ', v_ice(ji+1,jj)294 ! WRITE(numout,*) ' e1t(ji,jj) : ', e1t(ji,jj)295 ! WRITE(numout,*) ' at_i(ji,jj) : ', at_i(ji,jj)296 ! WRITE(numout,*) ' rdt_ice : ', rdt_ice297 ! WRITE(numout,*) ' area export : ', vinfor(83)298 ! WRITE(numout,*) ' vol export : ', vinfor(84)299 269 END DO 300 270 … … 311 281 vinfor(8) = vinfor(8) + vt_i(ji,jj)*aire(ji,jj) / 1.0e12 !ice volume 312 282 vinfor(10) = vinfor(10) + vt_s(ji,jj)*aire(ji,jj) / 1.0e12 !snow volume 313 ! vinfor(12) = vinfor(12) + vt_i(ji,jj)*aire(ji,jj) / 1.0e12 !undef ice volume314 ! vinfor(14) = vinfor(14) + vt_i(ji,jj)*aire(ji,jj) / 1.0e12 !def ice volume315 283 vinfor(16) = vinfor(16) + ot_i(ji,jj)*vt_i(ji,jj)*aire(ji,jj) / 1.0e12 !mean age 316 284 vinfor(30) = vinfor(30) + smt_i(ji,jj)*vt_i(ji,jj)*aire(ji,jj) / 1.0e12 !mean salinity 317 ! vinfor(32) = vinfor(32) + u_ice(ji,jj)*u_ice(ji,jj) + v_ice(ji,jj)*v_ice(ji,jj) !ice vel 285 ! this diagnostic is not well computed (weighted by vol instead 286 ! of area) 318 287 vinfor(32) = vinfor(32) + vt_i(ji,jj)*( u_ice(ji,jj)*u_ice(ji,jj) + & 319 288 v_ice(ji,jj)*v_ice(ji,jj) )*aire(ji,jj)/1.0e12 !ice vel … … 324 293 vinfor(62) = vinfor(62) + sss_io(ji,jj)*at_i(ji,jj)*aire(ji,jj) / 1.0e12 !SSS 325 294 vinfor(66) = vinfor(66) + et_s(ji,jj)/1.0e9*aire(ji,jj) / 1.0e12 ! snow temperature 326 vinfor(68) = vinfor(68) + et_i(ji,jj)/1.0e9*aire(ji,jj) / 1.0e12 ! ice heat content295 vinfor(68) = vinfor(68) + et_i(ji,jj)/1.0e9*aire(ji,jj) / 1.0e12 ! ice enthalpy 327 296 vinfor(70) = vinfor(70) + v_i(ji,jj,1)*aire(ji,jj) / 1.0e12 !ice volume 328 297 vinfor(72) = vinfor(72) + v_i(ji,jj,2)*aire(ji,jj) / 1.0e12 !ice volume … … 330 299 vinfor(76) = vinfor(76) + v_i(ji,jj,4)*aire(ji,jj) / 1.0e12 !ice volume 331 300 vinfor(78) = vinfor(78) + v_i(ji,jj,5)*aire(ji,jj) / 1.0e12 !ice volume 332 ! vinfor(80) = vinfor(80) + v_i(ji,jj,6)*aire(ji,jj) / 1.0e12 ! mass flux333 301 vinfor(80) = 0.0 334 302 vinfor(82) = vinfor(82) + fmass(ji,jj)*aire(ji,jj) / 1.0e12 ! mass flux 335 ! vinfor(84) = vinfor(84) + ht_i(ji,jj,5)*aire(ji,jj) / 1.0e12 ! mass flux336 303 ENDIF 337 304 END DO … … 346 313 END DO 347 314 348 ! DO jl = ice_cat_bounds(2,1), ice_cat_bounds(3,2)349 ! DO jj = 2, njeqm1350 ! DO ji = fs_2, fs_jpim1 ! vector opt.351 ! vinfor(14) = vinfor(14) + v_i(ji,jj,jl)*aire(ji,jj) / 1.0e12 !def ice volume (inc. rafted ice)352 ! END DO353 ! END DO354 ! END DO355 315 vinfor(14) = 0.0 356 316 357 317 zindb = 1.0 - MAX(0.0,SIGN(1.0,-vinfor(8))) 358 vinfor(16) = zindb * vinfor(16) / MAX(vinfor(8),epsi06) ! these have to be divided by total ice volume to have the359 vinfor(30) = zindb * vinfor(30) / MAX(vinfor(8),epsi06) ! right value318 vinfor(16) = zindb * vinfor(16) / MAX(vinfor(8),epsi06) ! these have to be divided by ice vol 319 vinfor(30) = zindb * vinfor(30) / MAX(vinfor(8),epsi06) ! 360 320 vinfor(32) = zindb * SQRT( vinfor(32) / MAX( vinfor(8) , epsi06 ) ) 361 vinfor(68) = zindb * vinfor(68) / MAX(vinfor(8),epsi06) ! right value321 vinfor(68) = zindb * vinfor(68) / MAX(vinfor(8),epsi06) ! 362 322 363 323 zindb = 1.0 - MAX(0.0,SIGN(1.0,-vinfor(6))) 364 vinfor(54) = zindb * vinfor(54) / MAX(vinfor(6),epsi06) ! these have to be divided by total ice extent to have the365 vinfor(56) = zindb * vinfor(56) / MAX(vinfor(6),epsi06) ! right value324 vinfor(54) = zindb * vinfor(54) / MAX(vinfor(6),epsi06) ! these have to be divided by ice extt 325 vinfor(56) = zindb * vinfor(56) / MAX(vinfor(6),epsi06) ! 366 326 vinfor(58) = zindb * vinfor(58) / MAX(vinfor(6),epsi06) ! 367 327 vinfor(80) = zindb * vinfor(80) / MAX(vinfor(6),epsi06) ! … … 408 368 END DO 409 369 zindb = 1.0 - MAX(0.0,SIGN(1.0,-vinfor(4))) ! 410 vinfor(64) = zindb * vinfor(64) / MAX(vinfor(4),epsi06) ! these have to be divided by total ice extent to have the370 vinfor(64) = zindb * vinfor(64) / MAX(vinfor(4),epsi06) ! divide by ice extt 411 371 !! 2.2) Diagnostics dependent on age 412 372 !!------------------------------------ … … 472 432 !! history : 473 433 !! 8.5 ! 03-08 (C. Ethe) original code 434 !! 9.0 ! 08-03 (M. Vancoppenolle) LIM3 474 435 !!------------------------------------------------------------------- 475 436 NAMELIST/namicedia/fmtinf, nfrinf, ninfo, ntmoy … … 479 440 & ndeb , & 480 441 & irecl 481 482 INTEGER :: nv ! indice of variable483 442 484 443 REAL(wp) :: zxx0, zxx1 ! temporary scalars -
trunk/NEMO/LIM_SRC_3/limdyn.F90
r825 r834 6 6 #if defined key_lim3 7 7 !!---------------------------------------------------------------------- 8 !! 'key_lim3' : LIMsea-ice model8 !! 'key_lim3' : LIM3 sea-ice model 9 9 !!---------------------------------------------------------------------- 10 10 !! lim_dyn : computes ice velocities … … 36 36 37 37 !!---------------------------------------------------------------------- 38 !! LIM 2.0, UCL-LOCEAN-IPSL (2005)38 !! LIM 3.0, UCL-ASTR-LOCEAN-IPSL (2008) 39 39 !! $Header: /home/opalod/NEMOCVSROOT/NEMO/LIM_SRC/limdyn.F90,v 1.5 2005/03/27 18:34:41 opalod Exp $ 40 40 !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt … … 58 58 !! 1.0 ! 01-04 (LIM) Original code 59 59 !! 2.0 ! 02-08 (C. Ethe, G. Madec) F90, mpp 60 !! 3.0 ! 2007-03 (M.A. Morales Maqueda, S. Bouillon, M. Vancoppenolle) LIM3, EVP, C-grid 60 !! 3.0 ! 2007-03 (M.A. Morales Maqueda, S. Bouillon, M. Vancoppenolle) 61 !! LIM3, EVP, C-grid 61 62 !!------------------------------------------------------------------------------------ 62 63 !! * Local variables 63 INTEGER :: ji, jj , jl! dummy loop indices64 INTEGER :: ji, jj ! dummy loop indices 64 65 INTEGER :: i_j1, i_jpj ! Starting/ending j-indices for rheology 65 ! nemo modif 66 ! jhemis ! jhemis = 1 (NH) ; jhemis = -1 (SH) 66 67 67 REAL(wp) :: & 68 68 ztairx, ztairy, & ! tempory scalars … … 70 70 ztglx , ztgly , & 71 71 zt11, zt12, zt21, zt22 , & 72 zustm, zsfrld, zsfrldm4,&72 zustm, & 73 73 zsfrldmx2, zsfrldmy2, & 74 74 zu_ice, zv_ice, ztair2 75 ! nemo modif 75 76 76 REAL(wp),DIMENSION(jpj) :: & 77 77 zind, & ! i-averaged indicator of sea-ice … … 86 86 IF ( ln_limdyn ) THEN 87 87 88 ! Mean ice and snow thicknesses. 89 88 ! ocean velocity 90 89 u_oce(:,:) = u_io(:,:) * tmu(:,:) 91 90 v_oce(:,:) = v_io(:,:) * tmv(:,:) … … 94 93 old_v_ice(:,:) = v_ice(:,:) * tmv(:,:) 95 94 96 ! !Rheology (ice dynamics)97 ! !========95 ! Rheology (ice dynamics) 96 ! ======== 98 97 99 98 ! Define the j-limits where ice rheology is computed … … 105 104 IF(ln_ctl) CALL prt_ctl_info( 'lim_dyn : i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj ) 106 105 CALL lim_rhg( i_j1, i_jpj ) 107 108 106 ELSE ! optimization of the computational area 109 107 … … 124 122 i_j1 = MAX( 1, i_j1-1 ) 125 123 IF(ln_ctl) CALL prt_ctl_info( 'lim_dyn : NH i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj ) 126 127 124 CALL lim_rhg( i_j1, i_jpj ) 128 125 … … 169 166 ENDIF 170 167 171 ! !Ice-Ocean stress172 ! !================168 ! Ice-Ocean stress 169 ! ================ 173 170 DO jj = 2, jpjm1 174 ! jhemis = SIGN(1, jj - jeq )175 ! zsang = REAL(jhemis) * sangvg176 !nemo new version modif177 171 zsang = SIGN(1.e0, gphif(1,jj-1) ) * sangvg 178 172 … … 193 187 ztairy = ( 2.0 - at_i(ji,jj) - at_i(ji,jj+1) ) * gtauy(ji,jj) / cai * cao 194 188 195 ! ! test on oscillations196 ! ztairx = 0.0197 ! ztairy = 0.0198 199 189 zsfrldmx2 = at_i(ji,jj) + at_i(ji+1,jj) 200 190 zsfrldmy2 = at_i(ji,jj) + at_i(ji,jj+1) … … 208 198 ztglx = zsfrldmx2 * rhoco * zmod * ( cangvg * zu_ice - zsang * zv_ice ) 209 199 ztgly = zsfrldmy2 * rhoco * zmod * ( cangvg * zv_ice + zsang * zu_ice ) 210 211 ! ! test on oscillations 212 ! ztglx = 0.0 213 ! ztgly = 0.0 200 ! 201 ! ! IMPORTANT 202 ! ! these lignes are bound to prevent numerical oscillations 203 ! ! in the ice-ocean stress 204 ! ! They are physically ill-based. There is a cleaner solution 205 ! ! to try (remember discussion in Paris Gurvan) 206 ! 214 207 ztglx = ztglx * exp( - zmod / 0.5 ) 215 208 ztgly = ztglx * exp( - zmod / 0.5 ) … … 281 274 CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_dyn : at_i :') 282 275 ENDIF 283 284 !OPA9_DYN LIM2.1 2005285 !END OPA9_DYN LIM2.1 2005286 276 287 277 END SUBROUTINE lim_dyn -
trunk/NEMO/LIM_SRC_3/limflx.F90
r825 r834 6 6 #if defined key_lim3 7 7 !!---------------------------------------------------------------------- 8 !! 'key_lim3' LIMsea-ice model8 !! 'key_lim3' LIM3 sea-ice model 9 9 !!---------------------------------------------------------------------- 10 10 !! lim_flx : flux at the ice / ocean interface … … 21 21 USE in_out_manager 22 22 USE albedo 23 USE limicepoints24 23 USE par_ice 25 24 USE prtctl ! Print control … … 69 68 !! * Modules used 70 69 !! * Local variables 71 INTEGER :: ji, jj , jl! dummy loop indices70 INTEGER :: ji, jj ! dummy loop indices 72 71 73 72 INTEGER :: & … … 80 79 !! zfcm2 , & ! non solar heat fluxes 81 80 zfold , & 82 zdfseqv, & ! increment to the ice salt flux83 81 #if defined key_lim_fdd 84 82 zfons, & ! salt exchanges at the ice/ocean interface -
trunk/NEMO/LIM_SRC_3/limhdf.F90
r825 r834 6 6 #if defined key_lim3 7 7 !!---------------------------------------------------------------------- 8 !! 'key_lim3' LIMsea-ice model8 !! 'key_lim3' LIM3 sea-ice model 9 9 !!---------------------------------------------------------------------- 10 10 !! lim_hdf : diffusion trend on sea-ice variable -
trunk/NEMO/LIM_SRC_3/limistate.F90
r825 r834 6 6 #if defined key_lim3 7 7 !!---------------------------------------------------------------------- 8 !! 'key_lim3' : LIMsea-ice model8 !! 'key_lim3' : LIM3 sea-ice model 9 9 !!---------------------------------------------------------------------- 10 10 !! lim_istate : Initialisation of diagnostics ice variables … … 21 21 USE dom_ice 22 22 USE ice 23 USE limicepoints24 23 USE lbclnk 25 24 … … 72 71 73 72 !! * Local variables 74 INTEGER :: ji, jj, jk, jl , jm! dummy loop indices73 INTEGER :: ji, jj, jk, jl ! dummy loop indices 75 74 76 75 REAL(wp) :: zidto, & ! temporary scalar 77 zs0, ztf, zbin, zai, zvi, zvs, zmi, zms,&76 zs0, ztf, zbin, & 78 77 zeps6, zeps, ztmelts, & 79 78 epsi06 -
trunk/NEMO/LIM_SRC_3/limitd_me.F90
r825 r834 2 2 3 3 #if defined key_lim3 4 !!---------------------------------------------------------------------- 5 !! 'key_lim3' : LIM3 sea-ice model 6 !!---------------------------------------------------------------------- 7 !! 4 8 !!====================================================================== 5 9 !! *** MODULE limitd_me *** … … 7 11 !! computation of changes in g(h) 8 12 !!====================================================================== 9 13 !! 10 14 !!---------------------------------------------------------------------- 11 15 !! * Modules used … … 17 21 USE ice_oce ! ice variables 18 22 USE thd_ice 19 USE limicepoints20 23 USE limistate 21 24 USE in_out_manager 22 25 USE ice 23 26 USE par_ice 24 USE limthd_lab25 27 USE limthd_lac 26 28 USE limvar … … 74 76 REAL(wp) :: & 75 77 Cp 76 77 78 ! 78 79 !----------------------------------------------------------------------- 79 ! 80 ! Ridging diagnostic arrays for history files 80 81 !----------------------------------------------------------------------- 81 82 ! … … 88 89 89 90 !!---------------------------------------------------------------------- 90 !! LIM-@ 3.0, UCL-ASTR (200 5-2006)91 !! LIM-@ 3.0, UCL-ASTR (2008) 91 92 !! (c) UCL-ASTR and Martin Vancoppenolle 92 93 !!---------------------------------------------------------------------- … … 114 115 !! ** External : 115 116 !! 117 !! ** Steps : 118 !! 1) Thickness categories boundaries, ice / o.w. concentrations 119 !! Ridge preparation 120 !! 2) Dynamical inputs (closing rate, divu_adv, opning) 121 !! 3) Ridging iteration 122 !! 4) Ridging diagnostics 123 !! 5) Heat, salt and freshwater fluxes 124 !! 6) Compute increments of tate variables and come back to old values 125 !! 116 126 !! ** References : There are a lot of references and can be difficult / 117 127 !! boring to read … … 138 148 !! 139 149 !! ** History : 140 !! This routine was freely inpired from CICE 141 !! authors William H. Lipscomb, LANL 142 !! Elizabeth C. Hunke (LANL) 143 !! the routines were available online 144 !! so, big thanks to the authors and to LANL 145 !! The modified features are the ridging/rafting algorithm 150 !! This routine is based on CICE code 151 !! and authors William H. Lipscomb, 152 !! and Elizabeth C. Hunke, LANL 153 !! are gratefully acknowledged 146 154 !! 147 155 !! (02-2006) Martin Vancoppenolle, UCL-ASTR 148 !! Voila, c'est la derniere etape du codage!!! J'en ai149 !! plein le cul mais je m'en sors mieux qu'avant!!!150 !! je ne sais pas quand ca va marcher151 !! O toi, futur utilisateur, si tu pestes sur ce code152 !! Pardonne moi, sache que j'ai bien peste dessus aussi153 156 !! 154 157 !!--------------------------------------------------------------------! 155 158 !! * Arguments 156 159 157 !! * Local variables 158 INTEGER :: ji, & ! spatial dummy loop index 159 jj, & ! spatial dummy loop index 160 jk, & ! vertical layering dummy loop index 161 jl, & ! ice category dummy loop index 162 jm, & ! ice types dummy loop index 163 index, & ! indexes for ice points debugging 164 jbnd1, & 165 jbnd2, & 166 niter, & ! iteration counter 167 nitermax = 20 ! max number of ridging iterations (it used to 168 ! be 20 169 170 REAL(wp) :: & ! constant values 171 zeps = 1.0e-10, & 172 epsi10 = 1.0e-10, & 173 epsi06 = 1.0e-6 174 175 REAL(wp) :: & ! constant values for ice enthalpy 176 zindb 177 178 REAL(wp), DIMENSION(jpi,jpj) :: & 179 closing_net, & ! net rate at which area is removed (1/s) 180 ! (ridging ice area - area of new ridges) / dt 181 divu_adv , & ! divu as implied by transport scheme (1/s) 182 opning , & ! rate of opening due to divergence/shear 183 closing_gross, & ! rate at which area removed, not counting 184 ! area of new ridges 185 msnow_mlt , & ! mass of snow added to ocean (kg m-2) 186 esnow_mlt ! energy needed to melt snow in ocean (J m-2) 187 188 REAL(wp) :: & 189 w1, & ! temporary variable 190 tmpfac, & ! factor by which opening/closing rates are cut 191 dti ! 1 / dt 192 193 LOGICAL :: & 194 iterate_ridging, & ! if true, repeat the ridging 195 asum_error ! flag for asum .ne. 1 196 197 REAL(wp) :: & 198 big = 1.0e8 199 200 REAL (wp), DIMENSION(jpi,jpj) :: & ! 201 vt_i_init, vt_i_final, & ! ice volume summed over categories 202 vt_s_init, vt_s_final, & ! snow volume summed over categories 203 et_i_init, et_i_final, & ! ice energy summed over categories 204 et_s_init, et_s_final ! snow energy summed over categories 205 206 CHARACTER (len = 15) :: fieldid 160 !! * Local variables 161 INTEGER :: ji, & ! spatial dummy loop index 162 jj, & ! spatial dummy loop index 163 jk, & ! vertical layering dummy loop index 164 jl, & ! ice category dummy loop index 165 niter, & ! iteration counter 166 nitermax = 20 ! max number of ridging iterations 167 168 REAL(wp) :: & ! constant values 169 zeps = 1.0e-10, & 170 epsi10 = 1.0e-10, & 171 epsi06 = 1.0e-6 172 173 REAL(wp), DIMENSION(jpi,jpj) :: & 174 closing_net, & ! net rate at which area is removed (1/s) 175 ! (ridging ice area - area of new ridges) / dt 176 divu_adv , & ! divu as implied by transport scheme (1/s) 177 opning , & ! rate of opening due to divergence/shear 178 closing_gross, & ! rate at which area removed, not counting 179 ! area of new ridges 180 msnow_mlt , & ! mass of snow added to ocean (kg m-2) 181 esnow_mlt ! energy needed to melt snow in ocean (J m-2) 182 183 REAL(wp) :: & 184 w1, & ! temporary variable 185 tmpfac, & ! factor by which opening/closing rates are cut 186 dti ! 1 / dt 187 188 LOGICAL :: & 189 iterate_ridging, & ! if true, repeat the ridging 190 asum_error ! flag for asum .ne. 1 191 192 REAL(wp) :: & 193 big = 1.0e8 194 195 REAL (wp), DIMENSION(jpi,jpj) :: & ! 196 vt_i_init, vt_i_final ! ice volume summed over categories 197 198 CHARACTER (len = 15) :: fieldid 207 199 208 200 !!-- End of declarations … … 223 215 ! Set hi_max(ncat) to a big value to ensure that all ridged ice 224 216 ! is thinner than hi_max(ncat). 217 225 218 Cp = 0.5* grav * (rau0-rhoic)*rhoic/rau0 ! proport const for PE 226 ! BACK TO LIMA_MEC 227 ! hi_max(ice_cat_bounds(2,2)) = big 228 229 CALL lim_itd_me_ridgeprep 230 231 ! ! conservation check 232 ! CALL lim_column_sum (jpl, v_i, vt_i_init) 219 CALL lim_itd_me_ridgeprep ! prepare ridging 220 221 ! conservation check 222 IF ( con_i) CALL lim_column_sum (jpl, v_i, vt_i_init) 233 223 234 224 ! Initialize arrays. … … 288 278 END DO 289 279 290 !+++++ [291 WRITE(numout,*) ' 2. Dynamical inputs '292 WRITE(numout,*) ' closing_net ', closing_net(jiindex,jjindex)293 WRITE(numout,*) ' divu_adv ', divu_adv (jiindex,jjindex)294 WRITE(numout,*) ' opning ', opning (jiindex,jjindex)295 !+++++ ]296 297 298 280 !-----------------------------------------------------------------------------! 299 281 ! 3) Ridging iteration … … 304 286 305 287 DO WHILE ( iterate_ridging .AND. niter < nitermax ) 306 307 !+++++ [308 ! WRITE(numout,*)309 WRITE(numout,*) ' 3. Ridging iteration : ', niter310 WRITE(numout,*) ' ---------------------- '311 ! WRITE(numout,*)312 ! WRITE(numout,*) 'a_i ', a_i(jiindex, jjindex,1:jpl)313 ! WRITE(numout,*) 'v_i ', v_i(jiindex, jjindex,1:jpl)314 !+++++ ]315 316 ! 3.1 ice / o.w. concentrationts, participation and transfer functions317 288 318 289 DO jj = 1, jpj … … 342 313 END DO ! jj 343 314 344 !+++++ [345 WRITE(numout,*) ' correction 1'346 WRITE(numout,*) ' closing_gross ', closing_gross(jiindex,jjindex)347 WRITE(numout,*) ' opning ', opning (jiindex,jjindex)348 !+++++ [349 350 315 ! correction to closing rate / opening if excessive ice removal 351 316 !--------------------------------------------------------------- … … 356 321 DO jj = 1, jpj 357 322 DO ji = 1, jpi 358 ! IF (numit.GE.565) WRITE(numout,*) 'Stop Hier : ',ji,jj,jl359 323 IF ( a_i(ji,jj,jl) .GT. epsi11 .AND. athorn(ji,jj,jl) .GT. 0.0 ) THEN 360 324 w1 = athorn(ji,jj,jl) * closing_gross(ji,jj) * rdt_ice … … 369 333 END DO !jl 370 334 371 !+++++ [372 WRITE(numout,*) ' correction 2'373 WRITE(numout,*) ' closing_gross ', closing_gross(jiindex,jjindex)374 WRITE(numout,*) ' opning ', opning(jiindex,jjindex)375 !+++++ ]376 377 335 ! 3.3 Redistribute area, volume, and energy. 378 336 !-----------------------------------------------------------------------------! … … 385 343 386 344 CALL lim_itd_me_asumr 387 388 !+++++ [389 WRITE(numout,*)390 WRITE(numout,*) ' back from asumr'391 WRITE(numout,*) ' asum, ato_i, at_i : ', asum(jiindex,jjindex), &392 ato_i(jiindex,jjindex), at_i(jiindex, jjindex)393 WRITE(numout,*) ' a_i : ', a_i(jiindex,jjindex,1:jpl)394 WRITE(numout,*) ' v_i : ', v_i(jiindex,jjindex,1:jpl)395 !+++++ ]396 345 397 346 ! 3.5 Do we keep on iterating ??? … … 457 406 !-----------------------------------------------------------------------------! 458 407 ! fresh water source for ocean 459 fmmec(ji,jj) = fmmec(ji,jj) + msnow_mlt(ji,jj)*dti !in kg / m2 /s !!! very important 460 ! fresh_hist(ji,jj) = fresh_hist(ji,jj) + msnow_mlt(ji,jj)*dti 408 fmmec(ji,jj) = fmmec(ji,jj) + msnow_mlt(ji,jj)*dti 461 409 462 410 ! heat sink for ocean 463 411 fhmec(ji,jj) = fhmec(ji,jj) + esnow_mlt(ji,jj)*dti 464 ! fhnet_hist(ji,jj) = fhnet_hist(ji,jj) + esnow_mlt(ji,jj)*dti465 412 466 413 END DO … … 486 433 487 434 ! Conservation check 488 ! CALL lim_column_sum (jpl, v_i, vt_i_final) 489 ! fieldid = ' v_i : limitd_me ' 490 ! CALL lim_cons_check (vt_i_init, vt_i_final, 1.0e-6, fieldid) 435 IF ( con_i ) THEN 436 CALL lim_column_sum (jpl, v_i, vt_i_final) 437 fieldid = ' v_i : limitd_me ' 438 CALL lim_cons_check (vt_i_init, vt_i_final, 1.0e-6, fieldid) 439 ENDIF 491 440 492 441 !-----------------------------------------------------------------------------! … … 494 443 !-----------------------------------------------------------------------------! 495 444 496 !+++++ 497 ! BACK TO LIMA_MEC 498 ! hi_max(ice_cat_bounds(2,2)) = 20.0 499 ! only for debugging 500 CALL lim_var_glo2eqv 501 502 WRITE(numout,*) ' End of lim_itd_me ****** ' 503 WRITE(numout,*) ' u_ice : ', u_ice(jiindex,jjindex) 504 WRITE(numout,*) ' v_ice : ', v_ice(jiindex,jjindex) 505 506 DO jl = 1, jpl 507 WRITE(numout,*) '* - category number ', jl 508 WRITE(numout,*) ' a_i : ', a_i(jiindex,jjindex,jl) 509 WRITE(numout,*) ' ht_i : ', ht_i(jiindex,jjindex,jl) 510 WRITE(numout,*) ' v_i : ', v_i(jiindex,jjindex,jl) 511 WRITE(numout,*) ' v_s : ', v_s(jiindex,jjindex,jl) 512 WRITE(numout,*) ' e_s : ', e_s(jiindex,jjindex,1,jl)/1.0e9 513 WRITE(numout,*) ' e_i1 : ', e_i(jiindex,jjindex,1,jl)/1.0e9 514 WRITE(numout,*) ' e_i2 : ', e_i(jiindex,jjindex,2,jl)/1.0e9 515 WRITE(numout,*) ' smv_i : ', smv_i(jiindex,jjindex,jl) 516 WRITE(numout,*) ' sm_i : ', sm_i(jiindex,jjindex,jl) 517 WRITE(numout,*) ' oa_i : ', oa_i (jiindex,jjindex,jl) 518 WRITE(numout,*) ' o_i : ', o_i (jiindex,jjindex,jl) 519 WRITE(numout,*) ' t_snow : ', t_s(jiindex,jjindex,1,jl) 520 WRITE(numout,*) ' t_su : ', t_su(jiindex,jjindex,jl) 521 WRITE(numout,*) ' t_i1 : ', t_i(jiindex,jjindex,1,jl) 522 WRITE(numout,*) ' t_i2 : ', t_i(jiindex,jjindex,2,jl) 523 WRITE(numout,*) 524 END DO 525 526 ! DO jl = 1, jpl 527 ! DO jj = 1, jpj 528 ! DO ji = 1, jpi 529 ! IF (old_a_i(ji,jj,jl).eq.0.00) THEN 530 ! o_i(ji,jj,jl) = MAX( MIN( rdt_ice*float(numit)/86400.0, o_i(ji,jj,jl) ), 0.0 ) 531 ! ENDIF 532 ! END DO 533 ! END DO 534 ! END DO 535 !+++++ 536 537 ! ! Zero out categories with very small areas 445 CALL lim_var_glo2eqv 538 446 CALL lim_itd_me_zapsmall 539 447 … … 615 523 END DO 616 524 END DO 617 618 619 525 620 526 END SUBROUTINE lim_itd_me … … 624 530 SUBROUTINE lim_itd_me_icestrength (kstrngth) ! (subroutine 2/6) 625 531 626 !!---------------------------------------------------------------------- -------532 !!---------------------------------------------------------------------- 627 533 !! *** ROUTINE lim_itd_me_icestrength *** 628 534 !! ** Purpose : … … 634 540 !! dissipated per unit area removed from the ice pack under compression, 635 541 !! and assumed proportional to the change in potential energy caused 636 !! by ridging. 542 !! by ridging. Note that only Hibler's formulation is stable and that 543 !! ice strength has to be smoothed 637 544 !! 638 545 !! ** Inputs / Ouputs : kstrngth (what kind of ice strength we are using) … … 640 547 !! ** External : 641 548 !! 642 !! ** References : There are a lot of references and they are difficult to read 643 !! so the best thing to do is to accept what is here 644 !! 645 !! ** History : 646 !! This routine was freely inpired from CICE 647 !! authors William H. Lipscomb, LANL 648 !! Elizabeth C. Hunke (LANL) 649 !! the routines were available online 650 !! so, big thanks to the authors and to LANL 651 !! The modified features are the ridging/rafting algorithm 652 !! 653 !! (02-2006) Martin Vancoppenolle, UCL-ASTR 654 !! Voila, c'est la derniere etape du codage!!! J'en ai 655 !! plein le cul mais je m'en sors mieux qu'avant!!! 656 !!------------------------------------------------------------------ 549 !! ** References : 550 !! 551 !!---------------------------------------------------------------------- 657 552 !! * Arguments 658 553 … … 674 569 REAL(wp), DIMENSION(jpi,jpj) :: & 675 570 zworka !: temporary array used here 676 677 WRITE(numout,*) ' lim_itd_me_icestrength : compute ice strength '678 WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~~~ '679 571 680 572 !------------------------------------------------------------------------------! … … 707 599 hi * hi 708 600 709 ! strength(ji,jj) = strength(ji,jj) &710 ! + athorn(ji,jj,jl)/krdg(ji,jj,jl) &711 ! * 1.0/3.0 * (hrmax(ji,jj,jl)**3 - hrmin(ji,jj,jl)**3) &712 ! / (hrmax(ji,jj,jl)-hrmin(ji,jj,jl))713 601 !-------------------------- 714 602 ! PE gain from rafting ice … … 720 608 ! PE gain from ridging ice 721 609 !---------------------------- 722 ! anything from ridge porosity ?723 610 strength(ji,jj) = strength(ji,jj) & 724 611 + aridge(ji,jj,jl)/krdg(ji,jj,jl) & … … 753 640 END DO ! i 754 641 755 WRITE(numout,*) ' Pstar : ', Pstar756 WRITE(numout,*) ' vt_i : ', vt_i(jiindex,jjindex)757 WRITE(numout,*) ' C_rhg : ', C_rhg758 WRITE(numout,*) ' at_i : ', at_i(jiindex,jjindex)759 760 642 ksmooth = 1 761 643 … … 766 648 ! 5) Impact of brine volume 767 649 !------------------------------------------------------------------------------! 650 ! CAN BE REMOVED 768 651 ! 769 652 IF ( brinstren_swi .EQ. 1 ) THEN 770 771 WRITE(numout,*) ' mean brine volume : ', bv_i(jiindex,jjindex)772 653 773 654 DO jj = 1, jpj … … 778 659 zdummy = 0.0 779 660 ENDIF 780 ! zdummy = 0.0781 ! strength(ji,jj) = strength(ji,jj) * ( 1.0 - zdummy )782 ! Timco and O'brien, CRST, 1994, Flexural strength equation for sea ice783 661 strength(ji,jj) = strength(ji,jj) * exp(-5.88*SQRT(MAX(bv_i(ji,jj),0.0))) 784 662 END DO ! j 785 663 END DO ! i 786 664 787 WRITE(numout,*) ' strength after brine volume correction : ', &788 strength(jiindex,jjindex)789 790 665 ENDIF 791 666 … … 861 736 ENDIF ! ksmooth 862 737 863 WRITE(numout,*) ' Ice strength : ', strength(jiindex,jjindex)864 865 738 ! Boundary conditions 866 739 CALL lbc_lnk( strength, 'T', 1. ) … … 885 758 !! ** External : 886 759 !! 887 !! ** History :888 !! This routine was freely inspired from CICE889 !! authors William H. Lipscomb, LANL890 !! Elizabeth C. Hunke (LANL)891 !! the routines were available online892 !! so, big thanks to the authors and to LANL893 !! The modified features are the ridging/rafting algorithm894 !!895 !! (02-2006) Martin Vancoppenolle, UCL-ASTR896 !! J'espere que nos lignes sont courbes897 !! Mais j'ai bien peur qu'elles ne soient droites898 760 !!---------------------------------------------------------------------! 899 761 !! * Arguments … … 901 763 INTEGER :: & 902 764 ji,jj, & ! horizontal indices 903 jl,jl1, & ! thickness category index 904 jm, & ! ice type index 905 index, & 765 jl, & ! thickness category index 906 766 krdg_index ! which participation function using 907 767 … … 913 773 Gsum ! Gsum(n) = sum of areas in categories 0 to n 914 774 915 REAL(wp), DIMENSION(jpi,jpj,-1:jpl) :: &916 Gsumra ! Gsum(n) = sum of areas in categories 0 to n917 918 REAL(wp), DIMENSION(jpi,jpj,jpm) :: &919 ai_typ ! Surface of respective ice types920 921 775 REAL(wp) :: & 922 776 hi, & ! ice thickness for each cat (m) … … 927 781 928 782 REAL(wp) :: & 929 zindb, zdummy, & 930 zweight, & 783 zdummy, & 931 784 epsi06 = 1.0e-6 932 785 933 INTEGER, DIMENSION(jpi,jpj,-1:jpl) :: &934 index_cumul ! temporary array used here935 786 !------------------------------------------------------------------------------! 936 787 … … 958 809 ! because of divergence during transport 959 810 960 !+++++961 ! 0.0 Compute surfaces of relative ice types962 ! BACK TO LIMA_MEC963 ! ai_typ(:,:,:) = 0.0964 ! DO jm = 1, jpm965 ! DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)966 ! DO jj = 1, jpj967 ! DO ji = 1, jpi968 ! ai_typ(ji,jj,jm) = ai_typ(ji,jj,jm) + a_i(ji,jj,jl)969 ! END DO970 ! END DO971 ! END DO972 ! END DO973 ! WRITE(numout,*) ' ai_typ ', ai_typ(jiindex,jjindex,1:3)974 ! BACK TO LIMA_MEC975 976 811 ! Compute cumulative thickness distribution function 977 812 ! Compute the cumulative thickness distribution function Gsum, 978 813 ! where Gsum(n) is the fractional area in categories 0 to n. 979 980 814 ! initial value (in h = 0) equals open water area 981 815 … … 1006 840 1007 841 ! Normalize the cumulative distribution to 1 1008 !former line1009 ! zworka(:,:) = 1.0 / Gsum(:,:,jpl)1010 WRITE(numout,*) ' AVANT : numit '1011 !new lines1012 842 DO jj = 1, jpj 1013 843 DO ji = 1, jpi 1014 ! IF (numit.GE.46365) THEN1015 IF ( (ji .EQ. jiindex) .AND. ( jj.EQ.jjindex ) ) THEN1016 WRITE(numout,*) ' ji ,jj : ',ji,jj, 'Gsum:', &1017 Gsum(ji,jj,jpl)1018 ENDIF1019 844 zworka(ji,jj) = 1.0 / Gsum(ji,jj,jpl) 1020 845 END DO 1021 846 END DO 1022 847 1023 WRITE(numout,*) ' APRES '1024 848 DO jl = 0, jpl 1025 849 DO jj = 1, jpj 1026 850 DO ji = 1, jpi 1027 IF ( (ji .EQ. jiindex) .AND. ( jj.EQ.jjindex ) ) THEN1028 WRITE(numout,*) ' ji ,jj : ',ji,jj,jl, 'Gsum:', &1029 Gsum(ji,jj,jl)1030 WRITE(numout,*) ' zworka ', zworka(ji,jj)1031 ENDIF1032 1033 851 Gsum(ji,jj,jl) = Gsum(ji,jj,jl) * zworka(ji,jj) 1034 852 END DO 1035 853 END DO 1036 854 END DO 1037 1038 ! BACK TO LIMA_MEC1039 ! For rafted ice1040 ! jm = 31041 ! DO jj = 1, jpj1042 ! DO ji = 1, jpi1043 ! Gsumra(ji,jj,ice_cat_bounds(jm,1)-1) = 0.01044 !! zworka(ji,jj) = 1.0 / ai_typ(ji,jj,jm) ! used to normalize Gsum1045 ! END DO1046 ! END DO1047 1048 ! DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)1049 ! DO jj = 1, jpj1050 ! DO ji = 1, jpi1051 ! Gsumra(ji,jj,jl) = Gsumra(ji,jj,jl-1) + a_i(ji,jj,jl)*zworka(ji,jj)1052 ! END DO1053 ! END DO1054 ! END DO1055 1056 ! +++++ [1057 WRITE(numout,*)1058 WRITE(numout,*) ' Gsum : ', Gsum(jiindex,jjindex,-1:jpl)1059 ! +++++ ]1060 1061 855 1062 856 ! 1.3 Compute participation function a(h) = b(h).g(h) (athorn) … … 1098 892 ! precompute exponential terms using Gsum as a work array 1099 893 zdummy = 1.0 / (1.0-EXP(-astari)) 1100 WRITE(numout,*) ' zdummy : ', zdummy1101 894 1102 895 DO jl = -1, jpl … … 1119 912 ENDIF ! krdg_index 1120 913 1121 WRITE(numout,*) ' astari : ', astari1122 WRITE(numout,*) ' athorn : ', athorn(jiindex,jjindex,0:jpl)1123 1124 ! BACK TO LIMA_MEC1125 ! jm = 31126 ! DO jl = ice_cat_bounds(3,1), ice_cat_bounds(3,2) ! only undeformed ice participates1127 ! DO jj = 1, jpj1128 ! DO ji = 1, jpi1129 ! IF (Gsumra(ji,jj,jl) < Gstar) THEN1130 ! athorn(ji,jj,jl) = Gstari * (Gsumra(ji,jj,jl)-Gsumra(ji,jj,jl-1)) * &1131 ! (2.0 - (Gsumra(ji,jj,jl-1)+Gsumra(ji,jj,jl))*Gstari)1132 ! ELSEIF (Gsumra(ji,jj,jl-1) < Gstar) THEN1133 ! athorn(ji,jj,jl) = Gstari * (Gstar-Gsumra(ji,jj,jl-1)) * &1134 ! (2.0 - (Gsumra(ji,jj,jl-1)+Gstar)*Gstari)1135 ! ELSE1136 ! athorn(ji,jj,jl) = 0.01137 ! ENDIF1138 ! END DO ! ji1139 ! END DO ! jj1140 ! END DO ! jl1141 ! BACK TO LIMA_MEC1142 1143 !!!1144 ! BACK TO LIMA_MEC1145 ! Weighting1146 ! Undeformed ice1147 ! DO jl = ice_cat_bounds(1,1), ice_cat_bounds(1,2)1148 ! DO jj = 1, jpj1149 ! DO ji = 1, jpi1150 ! zindb = MAX( 0.0, SIGN( 1.0, a_i(ji,jj,jl) - epsi06 ) )1151 ! zweight = zindb*ai_typ(ji,jj,1) / &1152 ! MAX( ai_typ(ji,jj,1) + ai_typ(ji,jj,3), epsi06 )1153 ! athorn(ji,jj,jl) = zweight*athorn(ji,jj,jl)1154 ! END DO1155 ! END DO1156 ! END DO1157 1158 ! Rafted ice1159 ! jm = 31160 ! DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)1161 ! DO jj = 1, jpj1162 ! DO ji = 1, jpi1163 ! zindb = MAX( 0.0, SIGN( 1.0, a_i(ji,jj,jl) - epsi06 ) )1164 ! zweight = zindb*ai_typ(ji,jj,jm) / &1165 ! MAX( ai_typ(ji,jj,1) + ai_typ(ji,jj,3) , epsi06)1166 ! athorn(ji,jj,jl) = zweight*athorn(ji,jj,jl)1167 ! END DO1168 ! END DO1169 ! END DO1170 1171 !+++++1172 ! BACK TO LIMA_MEC1173 !!!1174 1175 914 ! Ridging and rafting ice participation functions 1176 915 IF ( raftswi .EQ. 1 ) THEN … … 1224 963 1225 964 ENDIF 1226 !+++++1227 WRITE(numout,*)1228 WRITE(numout,*) ' aridge : ', aridge(jiindex, jjindex, 1:jpl)1229 WRITE(numout,*) ' araft : ', araft (jiindex, jjindex, 1:jpl)1230 !+++++1231 965 1232 966 !----------------------------------------------------------------- … … 1290 1024 DO ji = 1, jpi 1291 1025 aksum(ji,jj) = aksum(ji,jj) & 1292 ! + athorn(ji,jj,jl) * (1.0 - 1.0/krdg(ji,jj,jl))1293 1026 + aridge(ji,jj,jl) * (1.0 - 1.0/krdg(ji,jj,jl)) & 1294 1027 + araft (ji,jj,jl) * (1.0 - 1.0/kraft) … … 1297 1030 END DO 1298 1031 1299 WRITE(numout,*) ' krdg : ', krdg(jiindex,jjindex,1:jpl)1300 WRITE(numout,*) ' hraft : ', hraft(jiindex,jjindex,1:jpl)1301 WRITE(numout,*) ' aksum : ', aksum(jiindex,jjindex)1302 1303 1032 END SUBROUTINE lim_itd_me_ridgeprep 1304 1033 … … 1324 1053 !! ** External : 1325 1054 !! 1326 !! ** History :1327 !! This routine was freely inpired from CICE1328 !! authors William H. Lipscomb, LANL1329 !! Elizabeth C. Hunke (LANL)1330 !! the routines were available online1331 !! so, big thanks to the authors and to LANL1332 !! The modified features are the ridging/rafting algorithm1333 !!1334 !! (02-2006) Martin Vancoppenolle, UCL-ASTR1335 !! Dans 5 mois je suis sensé avoir fini ma thèse alors j'ai interet a cavaler1336 !! sinon ca va mal se passer1337 !! J'ai trop chaud, je suis fatigué1338 !! Mon ame se desseche t elle ??1339 !! Suis je en train de changer ???1340 !! Oublie-je l'essentiel de mon histoire ??1341 !!1342 1055 1343 1056 REAL (wp), DIMENSION(jpi,jpj), INTENT(IN) :: & … … 1355 1068 jk, & ! ice layer index 1356 1069 ij, & ! horizontal index, combines i and j loops 1357 icells, & ! number of cells with aicen > puny 1358 index ! index for debug ice points 1070 icells ! number of cells with aicen > puny 1359 1071 1360 1072 INTEGER, DIMENSION(1:(jpi+1)*(jpj+1)) :: & … … 1370 1082 vsnon_init, & ! snow volume before ridging 1371 1083 esnon_init, & ! snow energy before ridging 1372 smv_i_init, 1373 oa_i_init 1084 smv_i_init, & ! ice salinity before ridging 1085 oa_i_init ! ice age before ridging 1374 1086 1375 1087 REAL(wp), DIMENSION(jpi,jpj,jkmax,jpl) :: & … … 1380 1092 ardg1 , & ! area of ice ridged 1381 1093 ardg2 , & ! area of new ridges 1382 virdg , & ! ice volume of ridging ice1383 1094 vsrdg , & ! snow volume of ridging ice 1384 1095 esrdg , & ! snow energy of ridging ice 1385 smrdg , & ! salinity of ridging ice1386 oirdg , & ! volumic age content of ridging ice1387 1096 oirdg1 , & ! areal age content of ridged ice 1388 1097 oirdg2 , & ! areal age content of ridging ice … … 1397 1106 srdg1 , & ! sal*volume of ice ridged 1398 1107 srdg2 , & ! sal*volume of new ridges 1399 smsw , & ! sal*volume of water trapped into ridges 1400 totvrdg1 , & ! total volume of ice ridged 1401 totvrdg2 ! total volume of new ridges 1108 smsw ! sal*volume of water trapped into ridges 1402 1109 1403 1110 REAL(wp), DIMENSION(jpi,jpj) :: & … … 1409 1116 esrft , & ! snow energy of rafting ice 1410 1117 smrft , & ! salinity of rafting ice 1411 oirft , & ! volumic age content of rafting ice1412 1118 oirft1 , & ! areal age content of rafted ice 1413 1119 oirft2 ! areal age content of rafting ice 1414 1120 1415 1121 REAL(wp), DIMENSION(jpi,jpj,jkmax) :: & 1416 eirdg , & ! ice energy of ridging ice1417 1122 eirft , & ! ice energy of rafting ice 1418 1123 erdg1 , & ! enth*volume of ice ridged … … 1477 1182 END DO !jj 1478 1183 END DO !ji 1479 1480 !+++++1481 ! WRITE(numout,*) ' 1. Change in open water area '1482 ! WRITE(numout,*) ' closing_gross', closing_gross(jiindex, jjindex)1483 ! WRITE(numout,*) ' opning ', opning (jiindex, jjindex)1484 ! WRITE(numout,*) ' ato_i ', ato_i (jiindex, jjindex)1485 ! WRITE(numout,*)1486 !+++++1487 1184 1488 1185 ! if negative open water area alert it … … 1527 1224 END DO !jl 1528 1225 1529 !+++++ 1530 WRITE(numout,*) ' State variables before ridging is computed ' 1531 WRITE(numout,*) ' a_i ', a_i(jiindex, jjindex, 1:jpl) 1532 WRITE(numout,*) ' v_i ', v_i(jiindex, jjindex, 1:jpl) 1533 WRITE(numout,*) ' v_s ', v_s(jiindex, jjindex, 1:jpl) 1534 WRITE(numout,*) ' e_s ', e_s(jiindex, jjindex, 1, 1:jpl) 1535 WRITE(numout,*) ' smv_i ', smv_i(jiindex, jjindex, 1:jpl) 1536 WRITE(numout,*) ' oa_i ', oa_i (jiindex, jjindex, 1:jpl) 1537 WRITE(numout,*) ' e_i 1 ', e_i (jiindex, jjindex, 1, 1:jpl) 1538 WRITE(numout,*) ' e_i 2 ', e_i (jiindex, jjindex, 2, 1:jpl) 1539 WRITE(numout,*) 1540 ! !+++++ 1541 1226 ! 1542 1227 !----------------------------------------------------------------- 1543 1228 ! 3) Pump everything from ice which is being ridged / rafted … … 1552 1237 !------------------------------------------------ 1553 1238 1554 !+++++1555 ! WRITE(numout,*)1556 ! WRITE(numout,*) ' *** jl1 *** ', jl11557 !+++++1558 1239 icells = 0 1559 1240 DO jj = 1, jpj … … 1571 1252 large_afrft = .false. 1572 1253 1573 ! !+++++1574 ! WRITE(numout,*) ' 3.1 ridging cells '1575 ! WRITE(numout,*) ' icells ', icells1576 ! WRITE(numout,*)1577 !+++++1578 1579 1254 DO ij = 1, icells 1580 1255 ji = indxi(ij) … … 1594 1269 oirdg2(ji,jj)= oirdg1(ji,jj) / krdg(ji,jj,jl1) 1595 1270 oirft2(ji,jj)= oirft1(ji,jj) / kraft 1596 1597 !+++++1598 ! IF ( ( ji.EQ.jiindex ) .AND. (jj.EQ.jjindex) ) THEN1599 ! WRITE(numout,*) ' 3.2 area of ridging ice and of new ridge '1600 ! WRITE(numout,*) ' aridge ', aridge(ji,jj,jl1)1601 ! WRITE(numout,*) ' araft ', araft (ji,jj,jl1)1602 ! WRITE(numout,*) ' krdg ', krdg (ji,jj,jl1)1603 ! WRITE(numout,*) ' kraft ', kraft1604 ! WRITE(numout,*) ' clogross ', closing_gross(ji,jj)1605 ! WRITE(numout,*) ' ardg1 : ', ardg1(jiindex,jjindex)1606 ! WRITE(numout,*) ' arft1 : ', arft1(jiindex,jjindex)1607 ! WRITE(numout,*) ' ardg2 : ', ardg2(jiindex,jjindex)1608 ! WRITE(numout,*) ' arft2 : ', arft2(jiindex,jjindex)1609 ! WRITE(numout,*)1610 ! ENDIF1611 !+++++1612 1271 1613 1272 !--------------------------------------------------------------- … … 1629 1288 ENDIF 1630 1289 1631 ! !+++++1632 ! IF ( ( ji.EQ.jiindex ) .AND. (jj.EQ.jjindex) ) THEN1633 ! WRITE(numout,*) ' 3.3 Ridging / rafting fractions '1634 ! WRITE(numout,*) ' afrac : ', afrac(ji,jj)1635 ! WRITE(numout,*) ' afrft : ', afrft(ji,jj)1636 ! WRITE(numout,*)1637 ! ENDIF1638 ! !+++++1639 1640 1290 !-------------------------------------------------------------------------- 1641 1291 ! 3.4) Subtract area, volume, and energy from ridging 1642 1292 ! / rafting category n1. 1643 1293 !-------------------------------------------------------------------------- 1644 ! jl1 looping 1-jpl1645 ! ij looping 1-icells1646 1647 ! ridging volumes, heat contents ...1648 ! virdg is the volume of ridged ice, at the end of ridging process1649 ! i.e., it includes 30% porosity1650 ! virdg(ji,jj) = vicen_init(ji,jj,jl1) * afrac(ji,jj)1651 1294 vrdg1(ji,jj) = vicen_init(ji,jj,jl1) * afrac(ji,jj) / & 1652 1295 ( 1.0 + ridge_por ) … … 1659 1302 ( 1. + ridge_por ) 1660 1303 srdg2(ji,jj) = smv_i_init(ji,jj,jl1) * afrac(ji,jj) 1661 ! oirdg(ji,jj) = oa_i_init(ji,jj,jl1) * afrac(ji,jj) * &1662 ! ( 1.0 - ridge_por )1663 1304 1664 1305 ! rafting volumes, heat contents ... … … 1667 1308 esrft(ji,jj) = esnon_init(ji,jj,jl1) * afrft(ji,jj) 1668 1309 smrft(ji,jj) = smv_i_init(ji,jj,jl1) * afrft(ji,jj) 1669 ! oirft(ji,jj) = oa_i_init(ji,jj,jl1) * afrft(ji,jj)1670 1310 1671 1311 ! substract everything … … 1683 1323 ! Salinity 1684 1324 !------------- 1685 ! salinity times volume of the seawater into ridges1686 ! WRITE(numout,*) ' sss_io : ', sss_io(ji,jj)1687 ! WRITE(numout,*) ' vsw : ', vsw(ji,jj)1688 ! WRITE(numout,*) ' ridge_por : ', ridge_por1689 1690 1325 smsw(ji,jj) = sss_io(ji,jj) * vsw(ji,jj) * ridge_por 1691 1326 … … 1704 1339 srdg2(ji,jj) = sm_newridge * vrdg2(ji,jj) 1705 1340 1706 1707 1341 !------------------------------------ 1708 ! 3. 5Increment ridging diagnostics1342 ! 3.6 Increment ridging diagnostics 1709 1343 !------------------------------------ 1710 1344 … … 1720 1354 1721 1355 !------------------------------------------ 1722 ! 3. 6Put the snow somewhere in the ocean1356 ! 3.7 Put the snow somewhere in the ocean 1723 1357 !------------------------------------------ 1724 1358 … … 1742 1376 1743 1377 !----------------------------------------------------------------- 1744 ! 3. 7)Compute quantities used to apportion ice among categories1378 ! 3.8 Compute quantities used to apportion ice among categories 1745 1379 ! in the n2 loop below 1746 1380 !----------------------------------------------------------------- … … 1754 1388 1755 1389 1756 !++++++++1757 IF ( vrdg1(ji,jj) .GT. 0.0 ) THEN1758 IF ( ( vrdg2(ji,jj) / vrdg1(ji,jj) - ( 1. + ridge_por ) .GT. epsi10 ) ) &1759 THEN1760 WRITE(numout,*) ' ALERTE in ridging routine '1761 WRITE(numout,*) ' difference : ', vrdg2(ji,jj) / vrdg1(ji,jj) - ( 1. + ridge_por )1762 WRITE(numout,*) ' vrdg2/vrdg1 ', vrdg2(ji,jj) / &1763 vrdg1(ji,jj)1764 WRITE(numout,*) ' should be equal to 1+p : ', 1 + ridge_por1765 WRITE(numout,*) ' vrdg1 : ', vrdg1(ji,jj)1766 WRITE(numout,*) ' vrdg2 : ', vrdg2(ji,jj)1767 ENDIF1768 ENDIF1769 1770 !++++++++1771 !1772 1773 1390 END DO ! ij 1774 1391 1775 ! !+++++1776 ! !+++++1777 1778 ! !+++++1779 ! WRITE(numout,*)1780 ! WRITE(numout,*) ' Mass of salt of ridging ice 1 : ', srdg1(jiindex,jjindex)1781 ! WRITE(numout,*) ' Mass of salt of ridging ice 2 : ', srdg2(jiindex,jjindex)1782 ! WRITE(numout,*) ' Salinity of ridging ice : ', smrdg(jiindex,jjindex) / &1783 ! virdg(jiindex,jjindex) / (1.0-ridge_por)1784 ! WRITE(numout,*) ' Mass of salt after ridging : ', sm_i(ji,jj,jl1)1785 ! WRITE(numout,*) ' Sal after removal of ridg.ice ', sm_i(ji,jj,jl1) / &1786 ! vicen_init(ji,jj,jl1)1787 !+++++1788 1789 1790 1392 !-------------------------------------------------------------------- 1791 ! 3. 8)Compute ridging ice enthalpy, remove it from ridging ice and1393 ! 3.9 Compute ridging ice enthalpy, remove it from ridging ice and 1792 1394 ! compute ridged ice enthalpy 1793 1395 !-------------------------------------------------------------------- … … 1838 1440 ji = indxi(ij) 1839 1441 jj = indxj(ij) 1840 ! ridging diagnostics1841 ! WRITE(numout,*) ' debut ', eice_init(ji,jj)1842 ! WRITE(numout,*) ' erdg2 ', erdg2(ji,jj,jk)1843 ! WRITE(numout,*) ' erdg1 ', erdg1(ji,jj,jk)1844 ! WRITE(numout,*) ' ersw ', ersw(ji,jj,jk)1845 ! WRITE(numout,*) ' vsw ', vsw(ji,jj)1846 1442 eice_init(ji,jj) = eice_init(ji,jj) + erdg2(ji,jj,jk) - & 1847 1443 erdg1(ji,jj,jk) … … 1849 1445 END DO !jk 1850 1446 ENDIF 1851 1852 !+++++ [1853 ! WRITE(numout,*)1854 ! WRITE(numout,*) ' ridging ice enthalpy ', jl11855 ! WRITE(numout,*) ' eird1 ', erdg1(jiindex,jjindex,1:nlay_i)1856 ! WRITE(numout,*) ' eird2 old ', erdg2(jiindex,jjindex,1:nlay_i)1857 ! WRITE(numout,*) ' eirft ', eirft(jiindex,jjindex,1:2)1858 ! WRITE(numout,*) ' e_i ', e_i (jiindex,jjindex,1:nlay_i,jl1)1859 !+++++ ]1860 1447 1861 1448 IF (large_afrac) THEN ! there is a bug … … 1884 1471 ENDIF ! large_afrft 1885 1472 1886 !+++++ [1887 ! IF ( ( ji.EQ.jiindex ) .AND. (jj.EQ.jjindex) ) THEN1888 ! WRITE(numout,*) ' 3.4 Substract area, volume and energy '1889 ! WRITE(numout,*) ' virdg : ', virdg(ji,jj)1890 ! WRITE(numout,*) ' vsrdg ', vsrdg(ji,jj)1891 ! WRITE(numout,*) ' esrdg ', esrdg(ji,jj)1892 ! WRITE(numout,*) ' smrdg ', smrdg(ji,jj)1893 ! WRITE(numout,*) ' oirdg ', oirdg(ji,jj)1894 ! WRITE(numout,*)1895 ! WRITE(numout,*) ' virft : ', virft(ji,jj)1896 ! WRITE(numout,*) ' vsrft ', vsrft(ji,jj)1897 ! WRITE(numout,*) ' esrft ', esrft(ji,jj)1898 ! WRITE(numout,*) ' smrft ', smrft(ji,jj)1899 ! WRITE(numout,*) ' oirft ', oirft(ji,jj)1900 ! WRITE(numout,*)1901 1902 ! WRITE(numout,*) ' a_i : ', a_i(ji,jj,jl1)1903 ! WRITE(numout,*) ' v_i : ', v_i(ji,jj,jl1)1904 ! WRITE(numout,*) ' v_s ', v_s(ji,jj,jl1)1905 ! WRITE(numout,*) ' e_s ', e_s(ji,jj,1,jl1)1906 ! WRITE(numout,*) ' smv_i ', smv_i(ji,jj,jl1)1907 ! WRITE(numout,*) ' oa_i ', oa_i(ji,jj,jl1)1908 ! ENDIF1909 !+++++ ]1910 1911 1473 !------------------------------------------------------------------------------- 1912 1474 ! 4) Add area, volume, and energy of new ridge to each category jl2 … … 1946 1508 1947 1509 END DO ! ij 1948 ! WRITE(numout,*) ' new ridges being built from ', jl11949 ! WRITE(numout,*) ' added to ', jl21950 ! WRITE(numout,*) ' a_i : ', a_i(jiindex,jjindex,jl2)1951 ! WRITE(numout,*) ' v_i : ', v_i(jiindex,jjindex,jl2)1952 1510 1953 1511 ! Transfer ice energy to category jl2 by ridging … … 1964 1522 END DO ! jl2 (new ridges) 1965 1523 1966 ! WRITE(numout,*) ' new ridges built '1967 ! WRITE(numout,*) ' a_i : ', a_i(jiindex,jjindex,jl1)1968 ! WRITE(numout,*) ' v_i : ', v_i(jiindex,jjindex,jl1)1969 ! WRITE(numout,*) ' e_i : ', e_i(jiindex,jjindex,:,jl1)1970 1971 1524 DO jl2 = ice_cat_bounds(1,1), ice_cat_bounds(1,2) 1972 ! jl2 = ice_cat_bounds(3,1) ! rafted category1973 1525 1974 1526 DO ij = 1, icells … … 1990 1542 oa_i(ji,jj,jl2) = oa_i(ji,jj,jl2) & 1991 1543 + oirft2(ji,jj) 1992 1993 ! IF ( (ji.EQ.jiindex) .AND. (jj.EQ.jjindex) ) THEN1994 ! WRITE(numout,*) ' raft ice from jl1 ', jl1 , &1995 ! ' to jl2 :', jl21996 ! WRITE(numout,*) ' hraft (jl1) ', hraft(ji,jj,jl1)1997 ! WRITE(numout,*) ' hi_maxs ', hi_max(jl2-1), hi_max(jl2)1998 ! ENDIF1999 2000 1544 ENDIF ! hraft 2001 1545 … … 2017 1561 END DO ! jl2 2018 1562 2019 ! WRITE(numout,*) ' new rafted ice built '2020 ! WRITE(numout,*) ' a_i : ', a_i(jiindex,jjindex,jl1)2021 ! WRITE(numout,*) ' v_i : ', v_i(jiindex,jjindex,jl1)2022 ! WRITE(numout,*) ' e_i : ', e_i(jiindex,jjindex,:,jl1)2023 2024 1563 END DO ! jl1 (deforming categories) 2025 1564 … … 2060 1599 2061 1600 INTEGER :: ji, jj, jl 2062 2063 ! WRITE(numout,*) ' lim_itd_me_asumr : sum the ice / open water fractions'2064 ! WRITE(numout,*) ' ~~~~~~~~~~~~~~~~ '2065 1601 2066 1602 !----------------------------------------------------------------- … … 2172 1708 xtmp ! temporary variable 2173 1709 2174 WRITE(numout,*) ' lim_itd_me_zapsmall '2175 WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~ '2176 2177 1710 DO jl = 1, jpl 2178 1711 -
trunk/NEMO/LIM_SRC_3/limitd_th.F90
r825 r834 1 1 MODULE limitd_th 2 2 #if defined key_lim3 3 !!---------------------------------------------------------------------- 4 !! 'key_lim3' : LIM3 sea-ice model 5 !!---------------------------------------------------------------------- 3 6 !!====================================================================== 4 7 !! *** MODULE limitd_th *** … … 15 18 USE ice_oce ! ice variables 16 19 USE thd_ice 17 USE limicepoints18 20 USE limistate 19 21 USE in_out_manager 20 22 USE ice 21 23 USE par_ice 22 USE limthd_lab23 24 USE limthd_lac 24 25 USE limvar … … 44 45 45 46 !!---------------------------------------------------------------------- 46 !! LIM -@ 4.0, UCL-ASTR (2005)47 !! LIM3.0, UCL-ASTR (2008) 47 48 !! (c) UCL-ASTR and Martin Vancoppenolle 48 49 !!---------------------------------------------------------------------- … … 73 74 !! ** History : 74 75 !! (12-2005) Martin Vancoppenolle 75 !! Au moment ou j'ecris ces lignes, je ne me rends pas76 !! compte du boulot que j'entame... un truc de malate comme77 !! on dit ici chez les Belches78 76 !! 79 77 !!------------------------------------------------------------------ … … 81 79 82 80 !! * Local variables 83 INTEGER :: ji, & ! spatial dummy loop index 84 jj, & ! spatial dummy loop index 85 jk, & ! vertical layering dummy loop index 86 jl, & ! ice category dummy loop index 87 jm, & ! ice types dummy loop index 88 index, & 81 INTEGER :: jm, & ! ice types dummy loop index 89 82 jbnd1, & 90 83 jbnd2 … … 94 87 epsi10 = 1.0e-10 95 88 96 REAL(wp) :: & ! constant values for ice enthalpy97 zindb98 99 89 !!-- End of declarations 100 90 !!---------------------------------------------------------------------------------------------- 101 91 102 IF (lwp) THEN92 IF (lwp) THEN 103 93 WRITE(numout,*) 104 WRITE(numout,*) 'lim_itd_th : Thermodynamics of the ice thickness distribution'94 WRITE(numout,*) 'lim_itd_th : Thermodynamics of the ice thickness distribution' 105 95 WRITE(numout,*) '~~~~~~~~~~~' 106 96 ENDIF … … 120 110 CALL lim_var_agg(1) 121 111 122 !+++++123 WRITE(numout,*) ' From limitd_th : '124 WRITE(numout,*) ' at_i : ', at_i(jiindex,jjindex)125 WRITE(numout,*) ' vt_i : ', vt_i(jiindex,jjindex)126 WRITE(numout,*) ' vt_s : ', vt_s(jiindex,jjindex)127 DO jl = 1, jpl128 WRITE(numout,*) '* - category number ', jl129 WRITE(numout,*) ' a_i : ', a_i(jiindex,jjindex,jl)130 WRITE(numout,*) ' v_i : ', v_i(jiindex,jjindex,jl)131 WRITE(numout,*) ' ht_i : ', ht_i(jiindex,jjindex,jl)132 WRITE(numout,*) ' v_s : ', v_s(jiindex,jjindex,jl)133 WRITE(numout,*) ' ht_s : ', ht_s(jiindex,jjindex,jl)134 WRITE(numout,*) ' e_s : ', e_s(jiindex,jjindex,1,jl)/1.0e9135 WRITE(numout,*) ' e_i : ', e_i(jiindex,jjindex,1:nlay_i,jl)/1.0e9136 WRITE(numout,*) ' t_su : ', t_su(jiindex,jjindex,jl)137 WRITE(numout,*) ' t_snow : ', t_s(jiindex,jjindex,1,jl)138 WRITE(numout,*) ' t_i : ', t_i(jiindex,jjindex,1:nlay_i,jl)139 WRITE(numout,*) ' smv_i : ', smv_i(jiindex,jjindex,jl)140 WRITE(numout,*) ' oa_i : ', oa_i(jiindex,jjindex,jl)141 WRITE(numout,*)142 END DO143 !+++++144 145 !------------------------------------------------------------------------------|146 ! 2) Melt ice laterally.147 !------------------------------------------------------------------------------|148 ! DO jm = 1, jpm149 ! CALL lim_thd_lab(ice_cat_bounds(jm,1),ice_cat_bounds(jm,2))150 ! END DO151 ! CALL lim_thd_lab152 153 112 !------------------------------------------------------------------------------| 154 113 ! 3) Add frazil ice growing in leads. … … 158 117 CALL lim_var_glo2eqv ! only for info 159 118 160 !+++++161 WRITE(numout,*) ' limthd_lac, new values ***** '162 DO jl = 1, jpl163 WRITE(numout,*) '* - category number ', jl164 WRITE(numout,*) ' a_i : ', a_i(jiindex,jjindex,jl)165 WRITE(numout,*) ' ht_i : ', ht_i(jiindex,jjindex,jl)166 WRITE(numout,*) ' v_i : ', v_i(jiindex,jjindex,jl)167 WRITE(numout,*) ' v_s : ', v_s(jiindex,jjindex,jl)168 WRITE(numout,*) ' e_i : ', e_i(jiindex,jjindex,1:nlay_i,jl)/1.0e9169 WRITE(numout,*) ' smv_i : ', smv_i(jiindex,jjindex,jl)170 WRITE(numout,*) ' t_su : ', t_su(jiindex,jjindex,jl)171 WRITE(numout,*) ' t_snow : ', t_s(jiindex,jjindex,1,jl)172 WRITE(numout,*) ' t_i : ', t_i(jiindex,jjindex,1:nlay_i,jl)173 WRITE(numout,*)174 END DO175 !+++++176 119 !---------------------------------------------------------------------------------------- 177 120 ! 4) Computation of trend terms and get back to old values … … 230 173 !! (06-2006) Adaptation to include salt, age and types 231 174 !! (04-2007) Mass conservation checked 232 !!233 !! Je suis d'humeur massacrante aujourd'hui, tout234 !! le monde m'embete et m'empeche de coder235 !!236 !! Muere lentamente237 !! quien evita una pasion y su remolino de emociones238 !! justamente estas que regresan el brillo a los ojos239 !! y restauran los corazones destrozados240 !!241 175 !!------------------------------------------------------------------ 242 176 !! * Arguments … … 250 184 INTEGER :: ji, & ! spatial dummy loop index 251 185 jj, & ! spatial dummy loop index 252 jk, & ! vertical layering dummy loop index253 186 jl, & ! ice category dummy loop index 254 index, & ! for ice points255 187 zji, zjj, & ! dummy indices used when changing coordinates 256 188 nd ! used for thickness categories … … 307 239 308 240 REAL(wp) :: & ! constant values for ice enthalpy 309 zdummy, zdummy2, &310 241 zslope ! used to compute local thermodynamic "speeds" 311 242 … … 337 268 !! 1) Compute thickness and changes in each ice category 338 269 !!---------------------------------------------------------------------------------------------- 339 IF (lwp) THEN340 341 WRITE(numout,*) 'lim_itd_th_rem: Remapping the ice thickness distribution'342 343 WRITE(numout,*) 'klbnd : ', klbnd344 WRITE(numout,*) 'kubnd : ', kubnd345 WRITE(numout,*) 'ntyp : ', ntyp270 IF (lwp) THEN 271 WRITE(numout,*) 272 WRITE(numout,*) 'lim_itd_th_rem : Remapping the ice thickness distribution' 273 WRITE(numout,*) '~~~~~~~~~~~~~~~' 274 WRITE(numout,*) ' klbnd : ', klbnd 275 WRITE(numout,*) ' kubnd : ', kubnd 276 WRITE(numout,*) ' ntyp : ', ntyp 346 277 ENDIF 347 348 ! +++++ [349 ! index = 1350 ! jiindex = arc_sp_grid(index,1)351 ! jjindex = arc_sp_grid(index,2)352 ! WRITE(numout,*) '*', arc_sp_acro(index), ' ', arc_sp_name(index)353 ! WRITE(numout,*)354 ! WRITE(numout,*) ' a_i : ', a_i(jiindex,jjindex,klbnd:kubnd)355 ! WRITE(numout,*) ' ht_i : ', ht_i(jiindex,jjindex,klbnd:kubnd)356 ! WRITE(numout,*) ' v_i : ', v_i(jiindex,jjindex,klbnd:kubnd)357 ! +++++ ]358 278 359 279 zdhice(:,:,:) = 0.0 … … 361 281 DO jj = 1, jpj 362 282 DO ji = 1, jpi 363 364 283 zindb = 1.0-MAX(0.0,SIGN(1.0,-a_i(ji,jj,jl))) !0 if no ice and 1 if yes 365 284 ht_i(ji,jj,jl) = v_i(ji,jj,jl) / MAX(a_i(ji,jj,jl),zeps) * zindb … … 369 288 zdhice(ji,jj,jl) = ht_i(ji,jj,jl) - zht_i_o(ji,jj,jl) 370 289 ENDIF 371 372 290 END DO 373 291 END DO 374 292 END DO 375 !+++++376 ! WRITE(numout,*) ' 1 *** '377 ! WRITE(numout,*) ' klbnd -> kubnd ', klbnd, kubnd378 ! WRITE(numout,*) 'ht_i ', ht_i (jiindex,jjindex,klbnd:kubnd)379 ! WRITE(numout,*) 'zht_i_o', zht_i_o(jiindex,jjindex,klbnd:kubnd)380 ! WRITE(numout,*) 'zdhice ', zdhice (jiindex,jjindex,klbnd:kubnd)381 !+++++382 293 383 294 !----------------------------------------------------------------------------------------------- … … 392 303 END DO 393 304 END DO 394 395 !+++++396 ! WRITE(numout,*) ' 2 *** '397 ! WRITE(numout,*) ' klbnd -> kubnd ', klbnd, kubnd398 ! WRITE(numout,*) 'a_i ', a_i (jiindex,jjindex,klbnd:kubnd)399 ! WRITE(numout,*) 'at_i ', at_i (jiindex,jjindex)400 !+++++401 305 402 306 !----------------------------------------------------------------------------------------------- … … 478 382 ! ji 479 383 END DO !jl 480 !+++++481 ! WRITE(numout,*) ' 4 *** '482 ! WRITE(numout,*) ' klbnd -> kubnd - 1 ', klbnd, kubnd - 1483 ! WRITE(numout,*) ' hi_max ', hi_max(klbnd:kubnd-1)484 ! WRITE(numout,*) ' zhbnew ', zhbnew(jiindex,jjindex,klbnd:kubnd-1)485 !+++++486 384 487 385 !----------------------------------------------------------------------------------------------- … … 521 419 END DO !jj 522 420 523 !+++++524 ! WRITE(numout,*) ' 6 *** '525 ! WRITE(numout,*) ' klbnd -1, kubnd ', klbnd - 1, kubnd526 ! WRITE(numout,*) ' zhb0 ', zhb0(jiindex,jjindex)527 ! WRITE(numout,*) ' zhb1 ', zhb1(jiindex,jjindex)528 ! WRITE(numout,*) ' zhbnew klbnd-1 ', zhbnew(jiindex,jjindex,klbnd-1)529 ! WRITE(numout,*) ' zhbnew kubnd ', zhbnew(jiindex,jjindex,klbnd)530 !+++++531 532 421 !----------------------------------------------------------------------------------------------- 533 422 ! 7) Compute g(h) … … 537 426 g0(:,:,klbnd), g1(:,:,klbnd), hL(:,:,klbnd), & 538 427 hR(:,:,klbnd), zremap_flag) 539 540 !+++++541 ! WRITE(numout,*) ' 7a *** klbnd ', klbnd542 ! WRITE(numout,*) ' g0(klbnd) ', g0(jiindex,jjindex,klbnd)543 ! WRITE(numout,*) ' g1(klbnd) ', g1(jiindex,jjindex,klbnd)544 ! WRITE(numout,*) ' hL(klbnd) ', hL(jiindex,jjindex,klbnd)545 ! WRITE(numout,*) ' hR(klbnd) ', hR(jiindex,jjindex,klbnd)546 !+++++547 428 548 429 !- 7.2 Area lost due to melting of thin ice (first category, klbnd) … … 601 482 END DO 602 483 603 !+++++604 ! WRITE(numout,*) ' 7b *** klbnd->kubnd ', klbnd, kubnd605 ! WRITE(numout,*) ' g0 ', g0(jiindex,jjindex,klbnd:kubnd)606 ! WRITE(numout,*) ' g1 ', g1(jiindex,jjindex,klbnd:kubnd)607 ! WRITE(numout,*) ' hL ', hL(jiindex,jjindex,klbnd:kubnd)608 ! WRITE(numout,*) ' hR ', hR(jiindex,jjindex,klbnd:kubnd)609 ! WRITE(numout,*)610 ! WRITE(numout,*) ' ht_i ', ht_i(jiindex,jjindex,klbnd:kubnd)611 ! WRITE(numout,*) ' a_i ', a_i (jiindex,jjindex,klbnd:kubnd)612 ! WRITE(numout,*) ' v_i ', v_i (jiindex,jjindex,klbnd:kubnd)613 !+++++614 615 484 !----------------------------------------------------------------------------------------------- 616 485 ! 8) Compute area and volume to be shifted across each boundary … … 669 538 CALL lim_itd_shiftice ( klbnd, kubnd, zdonor, zdaice, zdvice ) 670 539 671 ! WRITE(numout,*) ' 9 *** '672 ! WRITE(numout,*) ' ht_i ', ht_i(jiindex,jjindex,klbnd:kubnd)673 ! WRITE(numout,*) ' a_i ', a_i (jiindex,jjindex,klbnd:kubnd)674 ! WRITE(numout,*) ' v_i ', v_i (jiindex,jjindex,klbnd:kubnd)675 676 540 !!---------------------------------------------------------------------------------------------- 677 541 !! 10) Make sure ht_i >= minimum ice thickness hi_min … … 741 605 !! authors: William H. Lipscomb, LANL, Elizabeth C. Hunke, LANL 742 606 !! (01-2006) Martin Vancoppenolle 743 !! Au moment ou j'ecris ces lignes, je ne me rends pas744 !! compte du boulot que j'entame... un truc de malate comme745 !! on dit ici chez les Belches746 !! This routine was inspired from CICE (W.H. Lipscomb, E. Hunke, C. M. Bitz)747 !! the sea ice model of LANL, Los Alamos, USA.748 !! Thanks to these guys and their team to put their routines online749 607 !! 750 608 !!------------------------------------------------------------------ … … 783 641 !!-- End of declarations 784 642 !!---------------------------------------------------------------------------------------------- 785 786 ! WRITE(numout,*) ' lim_itd_fitline : linearly fitting the function g(h) '787 ! WRITE(numout,*) ' ~~~~~~~~~~~~~~~ in category number ', num_cat788 643 789 644 DO jj = 1, jpj … … 853 708 !! authors: William H. Lipscomb, LANL, Elizabeth C. Hunke, LANL 854 709 !! (01-2006) Martin Vancoppenolle 855 !! Et quand fatigues de s'etre souvenus856 !! Nos souvenirs fatigues ne seront plus que des loques...857 !! "!!! Il a attrape la ratatinette, l'epouvantable ratatinette"858 !!859 !! This routine was largely inspired from CICE860 !! (W.H. Lipscomb, E. Hunke, C. M. Bitz)861 !! the sea ice model of LANL, Los Alamos, USA.862 !! Merci a eux et a leur equipe de mettre leurs routines en ligne863 !!864 !! J'ajoute aujourd'hui : le secret de ma vie pour l'instant865 !! c'est de trouver la balance entre mon ego et l'absence d'ego866 !! balance entre respect des autres et respect de moi867 !! trouver comment realiser la compréhension de l'autre868 !! sans déclencher ma propre aliénation869 710 !! 870 711 !!------------------------------------------------------------------ … … 887 728 jl1, & ! donor category 888 729 jk, & ! ice layer index 889 zji, zjj, & ! indices when changing from 2D-1D is done 890 index ! for ice points referencing 730 zji, zjj ! indices when changing from 2D-1D is done 891 731 892 732 REAL(wp), DIMENSION(jpi,jpj,jpl) :: & 893 zaTsfn, & 894 zqsnow 733 zaTsfn 895 734 896 735 REAL(wp), DIMENSION(jpi,jpj) :: & … … 902 741 zdeice , & ! ice energy transferred 903 742 zdsm_vice , & ! ice salinity times volume transferred 904 zsm_v1 , & ! ice salinity times volume905 zsm_v2 , & ! ice salinity times volume906 743 zdo_aice , & ! ice age times volume transferred 907 zo_v1 , & ! ice age times volume908 zo_v2 , & ! ice age times volume909 744 zdaTsf , & ! aicen*Tsfcn transferred 910 745 zindsn , & ! snow or not … … 928 763 929 764 !!-- End of declarations 930 ! WRITE(numout,*) ' lim_itd_shiftice : shifting ice between categories ' 931 ! WRITE(numout,*) ' ~~~~~~~~~~~~~~~~ ' 765 932 766 !---------------------------------------------------------------------------------------------- 933 767 ! 1) Define a variable equal to a_i*T_su … … 1024 858 END DO !jl 1025 859 1026 !-----------------------------------------------------------------1027 ! error messages1028 !-----------------------------------------------------------------1029 1030 ! if (daice_negative) then1031 ! do j = jlo,jhi1032 ! do i = ilo,ihi1033 ! if (donor(i,j,n) > 0 .and. daice(i,j,n) <= -puny) then1034 ! write(nu_diag,*) my_task,':',i,j,1035 ! & 'ITD Neg daice =',daice(i,j,n),' boundary',n1036 ! call abort_ice ('ice: ITD Neg daice')1037 ! endif1038 ! enddo1039 ! enddo1040 ! endif1041 1042 ! if (dvice_negative) then1043 ! do j = jlo,jhi1044 ! do i = ilo,ihi1045 ! if (donor(i,j,n) > 0 .and. dvice(i,j,n) <= -puny) then1046 ! write(nu_diag,*) my_task,':',i,j,1047 ! & 'ITD Neg dvice =',dvice(i,j,n),' boundary',n1048 ! call abort_ice ('ice: ITD Neg dvice')1049 ! endif1050 ! enddo1051 ! enddo1052 ! endif1053 1054 ! if (daice_greater_aicen) then1055 ! do j = jlo,jhi1056 ! do i = ilo,ihi1057 ! if (donor(i,j,n) > 0) then1058 ! n1 = donor(i,j,n)1059 ! if (daice(i,j,n) >= aicen(i,j,n1)+puny) then1060 ! write(nu_diag,*) my_task,':',i,j,1061 ! & 'ITD daice > aicen, cat',n11062 ! write(nu_diag,*) my_task,':',i,j,1063 ! & 'daice =', daice(i,j,n),1064 ! & 'aicen =', aicen(i,j,n1)1065 ! call abort_ice ('ice: ITD daice > aicen')1066 ! endif1067 ! endif1068 ! enddo1069 ! enddo1070 ! endif1071 1072 ! if (dvice_greater_vicen) then1073 ! do j = jlo,jhi1074 ! do i = ilo,ihi1075 ! if (donor(i,j,n) > 0) then1076 ! n1 = donor(i,j,n)1077 ! if (dvice(i,j,n) >= vicen(i,j,n1)+puny) then1078 ! write(nu_diag,*) my_task,':',i,j,1079 ! & 'ITD dvice > vicen, cat',n11080 ! write(nu_diag,*) my_task,':',i,j,1081 ! & 'dvice =', dvice(i,j,n),1082 ! & 'vicen =', vicen(i,j,n1)1083 ! call abort_ice ('ice: ITD dvice > vicen')1084 ! endif1085 ! endif1086 ! enddo1087 ! enddo1088 ! endif1089 1090 ! enddo ! boundaries 1 to ncat-11091 1092 860 !------------------------------------------------------------------------------- 1093 861 ! 3) Transfer volume and energy between categories … … 1210 978 t_su(ji,jj,jl) = zaTsfn(ji,jj,jl) / a_i(ji,jj,jl) 1211 979 zindsn = 1.0 - MAX(0.0,SIGN(1.0,-v_s(ji,jj,jl))) !0 if no ice and 1 if yes 1212 ! t_s(ji,jj,1,jl) = - zQsnow(ji,jj,jl) / MAX(v_s(ji,jj,jl),zeps) * zindsn + rtt1213 980 ELSE 1214 981 ht_i(ji,jj,jl) = 0.0 1215 982 t_su(ji,jj,jl) = rtt 1216 ! t_s(ji,jj,1,jl) = rtt1217 983 ENDIF 1218 984 END DO ! ji … … 1246 1012 !! authors: William H. Lipscomb, LANL, Elizabeth C. Hunke, LANL 1247 1013 !! (01-2006) Martin Vancoppenolle (adaptation) 1248 !!1249 !! Quel etre encore que celui-ci! Le Jugement Dernier sera la1250 !! avant qu'il vous fasse jamais une avance sur votre mois,1251 !! Seigneur! Tu peux supplier, te mettre en quatre,1252 !! meme si tu es dans la misere, il ne te donnera rien,1253 !! le vieux demon! Et quant on pense que, chez lui,1254 !! sa cuisiniere lui donne des gifles! Je ne vois pas l'interet1255 !! qu'il y a a travailler dans un ministere. Cela ne rapporte1256 !! absolument rien.1257 !!1258 !! This routine was largely inspired from CICE1259 !! (W.H. Lipscomb, E. Hunke, C. M. Bitz)1260 !! the sea ice model of LANL, Los Alamos, USA.1261 1014 !! 1262 1015 !!------------------------------------------------------------------ … … 1283 1036 REAL(wp) :: & ! constant values 1284 1037 zeps = 1.0e-10, & 1285 epsi10 = 1.0e-10, & 1286 zindb 1038 epsi10 = 1.0e-10 1287 1039 1288 1040 REAL (wp), DIMENSION(jpi,jpj) :: & ! … … 1294 1046 !!-- End of declarations 1295 1047 !------------------------------------------------------------------------------ 1296 ! WRITE(numout,*) ' lim_itd_th_reb '1297 ! WRITE(numout,*) ' ~~~~~~~~~~~~~~ '1298 ! WRITE(numout,*) ' Ice Type no ', ntyp1299 ! WRITE(numout,*) ' bounds of categories ', klbnd, kubnd1300 1048 1301 1049 ! ! conservation check 1302 ! CALL lim_column_sum (jpl, v_i, vt_i_init) 1303 ! CALL lim_column_sum (jpl, v_s, vt_s_init) 1050 IF ( con_i ) THEN 1051 CALL lim_column_sum (jpl, v_i, vt_i_init) 1052 CALL lim_column_sum (jpl, v_s, vt_s_init) 1053 ENDIF 1304 1054 1305 1055 ! … … 1307 1057 ! 1) Compute ice thickness. 1308 1058 !------------------------------------------------------------------------------ 1309 ! nothing to do1310 1059 DO jl = klbnd, kubnd 1311 1060 DO jj = 1, jpj … … 1447 1196 !------------------------------------------------------------------------------ 1448 1197 1449 ! CALL lim_column_sum (jpl, v_i, vt_i_final) 1450 ! fieldid = ' v_i : limitd_reb ' 1451 ! CALL lim_cons_check (vt_i_init, vt_i_final, 1.0e-6, fieldid) 1452 1453 ! CALL lim_column_sum (jpl, v_s, vt_s_final) 1454 ! fieldid = ' v_s : limitd_reb ' 1455 ! CALL lim_cons_check (vt_s_init, vt_s_final, 1.0e-6, fieldid) 1198 IF ( con_i ) THEN 1199 CALL lim_column_sum (jpl, v_i, vt_i_final) 1200 fieldid = ' v_i : limitd_reb ' 1201 CALL lim_cons_check (vt_i_init, vt_i_final, 1.0e-6, fieldid) 1202 1203 CALL lim_column_sum (jpl, v_s, vt_s_final) 1204 fieldid = ' v_s : limitd_reb ' 1205 CALL lim_cons_check (vt_s_init, vt_s_final, 1.0e-6, fieldid) 1206 ENDIF 1456 1207 1457 1208 END SUBROUTINE lim_itd_th_reb -
trunk/NEMO/LIM_SRC_3/limmsh.F90
r825 r834 6 6 #if defined key_lim3 7 7 !!---------------------------------------------------------------------- 8 !! 'key_lim3' LIMsea-ice model8 !! 'key_lim3' LIM3 sea-ice model 9 9 !!---------------------------------------------------------------------- 10 10 !! lim_msh : definition of the ice mesh … … 24 24 25 25 !!---------------------------------------------------------------------- 26 !! LIM 2.0, UCL-LOCEAN-IPSL (2005)26 !! LIM 3.0, UCL-ASTR-LOCEAN-IPSL (2008) 27 27 !! $Header: /home/opalod/NEMOCVSROOT/NEMO/LIM_SRC/limmsh.F90,v 1.5 2005/03/27 18:34:42 opalod Exp $ 28 28 !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt … … 75 75 fcor(:,:) = 2. * omega * SIN( gphit(:,:) * rad ) ! coriolis factor 76 76 77 !i DO jj = 1, jpj78 !i zmsk(jj) = SUM( tmask(:,jj,:) ) ! = 0 if land everywhere on a j-line79 !!ii write(numout,*) jj, zind(jj)80 !i END DO81 82 77 IF( fcor(1,1) * fcor(1,nlcj) < 0.e0 ) THEN ! local domain include both hemisphere 83 78 l_jeq = .TRUE. … … 118 113 tmv(:,:) = 0.0e0 ! CGrid EVP 119 114 !!i 120 121 122 115 ! metric coefficients for sea ice dynamic 123 116 !---------------------------------------- … … 169 162 & + e2t(ji-1,jj-1) * wght(ji,jj,1,1) 170 163 171 ! better written but change the last digit and thus solver in less than 100 timestep172 ! zh1p = e1t(ji-1,jj ) * wght(ji,jj,1,2) + e1t(ji,jj ) * wght(ji,jj,2,2) &173 ! & + e1t(ji-1,jj-1) * wght(ji,jj,1,1) + e1t(ji,jj-1) * wght(ji,jj,2,1)174 175 ! zh2p = e2t(ji-1,jj ) * wght(ji,jj,1,2) + e2t(ji,jj ) * wght(ji,jj,2,2) &176 ! & + e2t(ji-1,jj-1) * wght(ji,jj,1,1) + e2t(ji,jj-1) * wght(ji,jj,2,1)177 178 !!ibug =0 zusden = 1.0 / ( zh1p * zh2p * 4.e0 )179 164 zusden = 1.0 / MAX( zh1p * zh2p * 4.e0 , 1.e-20 ) 180 165 zusden2 = zusden * 2.0 … … 231 216 tms(:,:) = tmask(:,:,1) ! ice T-point : use surface tmask 232 217 233 !i here we can use umask with a i and j shift of -1,-1234 218 tmu(:,1) = 0.e0 235 219 tmu(1,:) = 0.e0 236 220 tmv(:,1) = 0.e0 237 221 tmv(1,:) = 0.e0 238 ! DO jj = 2, jpj ! ice U.V-point: computed from ice T-point mask 239 ! DO ji = 2, jpim1 240 ! bug 222 241 223 DO jj = 1, jpj - 1 242 224 DO ji = 2 , jpi - 1 -
trunk/NEMO/LIM_SRC_3/limrhg.F90
r825 r834 2 2 !!====================================================================== 3 3 !! *** MODULE limrhg *** 4 !! Ice rheology : performssea ice rheology4 !! Ice rheology : sea ice rheology 5 5 !!====================================================================== 6 6 #if defined key_lim3 7 7 !!---------------------------------------------------------------------- 8 !! 'key_lim3' LIMsea-ice model8 !! 'key_lim3' LIM3 sea-ice model 9 9 !!---------------------------------------------------------------------- 10 10 !! lim_rhg : computes ice velocities … … 36 36 rone = 1.e0 37 37 !!---------------------------------------------------------------------- 38 !! LIM 2.0, UCL-LOCEAN-IPSL (2005)38 !! LIM 3.0, UCL-LOCEAN-IPSL (2008) 39 39 !! $Header: /home/opalod/NEMOCVSROOT/NEMO/LIM_SRC/limrhg.F90,v 1.5 2005/03/27 18:34:42 opalod Exp $ 40 40 !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt … … 44 44 45 45 SUBROUTINE lim_rhg( k_j1, k_jpj ) 46 46 47 !!------------------------------------------------------------------- 47 !! *** SUBROUTINR lim_rhg *** 48 !! 49 !! ** purpose : determines sea ice drift from wind stress, ice-ocean 48 !! *** SUBROUTINE lim_rhg *** 49 !! EVP-C-grid 50 !! 51 !! ** purpose : determines sea ice drift from wind stress, ice-ocean 50 52 !! stress and sea-surface slope. Ice-ice interaction is described by 51 !! a non-linear elasto-viscous-plastic law including shear strength 52 !! and a bulk rheology (Hunke and Dukowicz, 2002). 53 !! 54 !! Grid B: 55 !! the routine calculates 4 estimates of the divergence, tension, 56 !! and shear on each corner of a given scalar grid box. Likewise, 57 !! four estimates of the components of the stress tensor are 58 !! calculated on each corner. The internal forces on a corner are 59 !! calculated as averages of the four stress tensor contributions. 60 !! 61 !! ** Action : - compute u_ice, v_ice the sea-ice velocity 62 !! 63 !! ** NOTE : for the ice dynamics, the labeling convention is 64 !! that the velocity point (i,j) is located on the southwest corner 65 !! of a scalar (i,j) grid box. This choice is somewhat unfortunate, 66 !! as most models locate the velocity point (i,j) on the northeast 67 !! corner... 68 !! 69 !! ** MORE IMPORTANT NOTES : related to the note above, it is crucial 70 !! to make sure that all variables and coefficients are located on 71 !! right points of the grid. Verify everywhere! 53 !! a non-linear elasto-viscous-plastic (EVP) law including shear 54 !! strength and a bulk rheology (Hunke and Dukowicz, 2002). 55 !! 56 !! The points in the C-grid look like this, dear reader 57 !! 58 !! (ji,jj) 59 !! | 60 !! | 61 !! (ji-1,jj) | (ji,jj) 62 !! --------- 63 !! | | 64 !! | (ji,jj) |------(ji,jj) 65 !! | | 66 !! --------- 67 !! (ji-1,jj-1) (ji,jj-1) 68 !! 69 !! ** Inputs : - wind forcing (stress), oceanic currents 70 !! ice total volume (vt_i) per unit area 71 !! snow total volume (vt_s) per unit area 72 !! 73 !! ** Action : - compute u_ice, v_ice : the components of the 74 !! sea-ice velocity vector 75 !! - compute delta_i, shear_i, divu_i, which are inputs 76 !! of the ice thickness distribution 77 !! 78 !! ** Steps : 1) Compute ice snow mass, ice strength 79 !! 2) Compute wind, oceanic stresses, mass terms and 80 !! coriolis terms of the momentum equation 81 !! 3) Solve the momentum equation (iterative procedure) 82 !! 4) Prevent high velocities if the ice is thin 83 !! 5) Recompute invariants of the strain rate tensor 84 !! which are inputs of the ITD, store stress 85 !! for the next time step 86 !! 6) Control prints of residual (convergence) 87 !! and charge ellipse. 88 !! The user should make sure that the parameters 89 !! nevp, telast and creepl maintain stress state 90 !! on the charge ellipse for plastic flow 91 !! e.g. in the Canadian Archipelago 92 !! 93 !! ** References : Hunke and Dukowicz, JPO97 94 !! Bouillon et al., 08, in prep (update this when 95 !! published) 96 !! Vancoppenolle et al., OM08 72 97 !! 73 98 !! History : 74 !! 1.0 ! 07-03 (M.A. Morales Maqueda, S. Bouillon , M. Vancoppenolle)75 !! EVP, C grid,LIM399 !! 1.0 ! 07-03 (M.A. Morales Maqueda, S. Bouillon) 100 !! 2.0 ! 08-03 M. Vancoppenolle : LIM3 76 101 !! 77 102 !!------------------------------------------------------------------- … … 79 104 ! 80 105 INTEGER, INTENT(in) :: & 81 k_j1 , &!: southern j-index for ice computation82 k_jpj !: northern j-index for ice computation106 k_j1 , & !: southern j-index for ice computation 107 k_jpj !: northern j-index for ice computation 83 108 84 109 ! * Local variables 85 INTEGER :: ji, jj , jl!: dummy loop indices110 INTEGER :: ji, jj !: dummy loop indices 86 111 87 112 INTEGER :: & 88 iim1, ijm1, iip1 , ijp1 , & ! temporary integers 89 iter, jter ! " " 113 iter, jter !: temporary integers 90 114 91 115 CHARACTER (len=50) :: charout 92 116 93 117 REAL(wp) :: & 94 zt11 , zt12 , zt21 , zt22 , & ! temporary scalars 95 ztagnx, ztagny , delta ! " " 118 zt11, zt12, zt21, zt22, & !: temporary scalars 119 ztagnx, ztagny, & !: wind stress on U/V points 120 delta ! 96 121 97 122 REAL(wp) :: & 98 za, zstms, zsang, zmask 99 100 REAL(wp),DIMENSION(jpi,jpj) :: & 101 zpresh, zpreshc, zfrld1, zfrld2, zmass1, zmass2, zcorl1, zcorl2, & 102 za1ct, za2ct, zc1, zusw , & 103 u_oce1, v_oce1, u_oce2, v_oce2, u_ice2, v_ice1 123 za, & !: 124 zstms, & !: temporary scalar for ice strength 125 zsang, & !: temporary scalar for coriolis term 126 zmask !: mask for the computation of ice mass 127 128 REAL(wp),DIMENSION(jpi,jpj) :: & 129 zpresh , & !: temporary array for ice strength 130 zpreshc , & !: Ice strength on grid cell corners (zpreshc) 131 zfrld1, zfrld2, & !: lead fraction on U/V points 132 zmass1, zmass2, & !: ice/snow mass on U/V points 133 zcorl1, zcorl2, & !: coriolis parameter on U/V points 134 za1ct, za2ct , & !: temporary arrays 135 zc1 , & !: ice mass 136 zusw , & !: temporary weight for the computation 137 !: of ice strength 138 u_oce1, v_oce1, & !: ocean u/v component on U points 139 u_oce2, v_oce2, & !: ocean u/v component on V points 140 u_ice2, & !: ice u component on V point 141 v_ice1 !: ice v component on U point 104 142 105 143 REAL(wp) :: & 106 dtevp,dtotel,ecc2,z0,zr,zcca,zccb,denr, & 107 zu_ice2,zv_ice1, zddc, zdtc, zdst, zdsshx, zdsshy, & 108 sigma1, sigma2 109 110 REAL(wp),DIMENSION(jpi,jpj) :: & 111 zf1, zf2 112 113 REAL(wp),DIMENSION(jpi,jpj) :: & 114 zdd, & !: Divergence [ 115 zdt, & !: 116 zds, & !: 117 deltat, & 118 deltac, & 119 zs1, & 120 zs2, & 121 zs12 122 123 REAL :: & 124 zumax !: Maximal tolerated ice velocity 125 126 REAL(wp) :: & 127 zresm !: Maximal error on ice velocity 128 129 REAL(wp),DIMENSION(jpi,jpj) :: & 130 zu_ice , & !: Ice velocity on previous time step 131 zv_ice , & 132 zresr !: Local error on velocity 144 dtevp, & !: time step for subcycling 145 dtotel, & !: 146 ecc2, & !: square of yield ellipse eccenticity 147 z0, & !: temporary scalar 148 zr, & !: temporary scalar 149 zcca, zccb, & !: temporary scalars 150 zu_ice2, & !: 151 zv_ice1, & !: 152 zddc, zdtc, & !: temporary array for delta on corners 153 zdst, & !: temporary array for delta on centre 154 zdsshx, zdsshy, & !: term for the gradient of ocean surface 155 sigma1, sigma2 !: internal ice stress 156 157 REAL(wp),DIMENSION(jpi,jpj) :: & 158 zf1, zf2 !: arrays for internal stresses 159 160 REAL(wp),DIMENSION(jpi,jpj) :: & 161 zdd, zdt, & !: Divergence and tension at centre of grid cells 162 zds, & !: Shear on northeast corner of grid cells 163 deltat, & !: Delta at centre of grid cells 164 deltac, & !: Delta on corners 165 zs1, zs2, & !: Diagonal stress tensor components zs1 and zs2 166 zs12 !: Non-diagonal stress tensor component zs12 133 167 134 168 REAL(wp) :: & 135 zindb , & !: ice (1) or not (0) 136 zdummy !: ice computation 169 zresm , & !: Maximal error on ice velocity 170 zindb , & !: ice (1) or not (0) 171 zdummy !: dummy argument 172 173 REAL(wp),DIMENSION(jpi,jpj) :: & 174 zu_ice , & !: Ice velocity on previous time step 175 zv_ice , & 176 zresr !: Local error on velocity 137 177 138 178 INTEGER :: & 139 179 stress_mean_swi = 1 140 141 !!---------------------------------------------------------------------- 142 143 ! 144 !------------------------------------------------------------------------------! 145 ! 0) Ice-Snow mass (zc1), ice strength (zpresh) ! 146 !------------------------------------------------------------------------------! 147 ! 148 ! Maximal tolerated velocity 149 ! zumax = 1.0 150 180 ! 181 !------------------------------------------------------------------------------! 182 ! 1) Ice-Snow mass (zc1), ice strength (zpresh) ! 183 !------------------------------------------------------------------------------! 184 ! 151 185 ! Put every vector to 0 152 186 zpresh(:,:) = 0.0 ; zpreshc(:,:) = 0.0 ; zfrld1(:,:) = 0.0 … … 154 188 za1ct(:,:) = 0.0 ; za2ct(:,:) = 0.0 155 189 zc1(:,:) = 0.0 ; zusw(:,:) = 0.0 156 157 190 u_oce1(:,:) = 0.0 ; v_oce1(:,:) = 0.0 ; u_oce2(:,:) = 0.0 158 191 v_oce2(:,:) = 0.0 ; u_ice2(:,:) = 0.0 ; v_ice1(:,:) = 0.0 159 160 zf1(:,:) = 0.0 ; zf2(:,:) = 0.0 161 192 zf1(:,:) = 0.0 ; zf2(:,:) = 0.0 162 193 zdd(:,:) = 0.0 ; zdt(:,:) = 0.0 ; zds(:,:) = 0.0 163 194 deltat(:,:) = 0.0 ; deltac(:,:) = 0.0 164 165 ! 166 !------------------------------------------------------------------------------! 167 ! 1) Ice-Snow mass (zc1), ice strength (zpresh) ! 168 !------------------------------------------------------------------------------! 169 ! 170 ! compute ice strength 195 zpresh(:,:) = 0.0 196 197 ! Ice strength on T-points 171 198 CALL lim_itd_me_icestrength(ridge_scheme_swi) 172 199 173 zpresh(:,:) = 0.0 174 200 ! Ice mass and temp variables 175 201 DO jj = k_j1 , k_jpj-1 176 202 DO ji = 1 , jpi … … 178 204 zpresh(ji,jj) = tms(ji,jj) * strength(ji,jj) / 2. 179 205 ! tmi = 1 where ther is ice or on land 206 ! Claude sees a bug, next line : tmi(ji,jj) 180 207 tmi = 1.0 - ( 1.0 - MAX( 0.0 , SIGN ( 1.0 , vt_i(ji,jj) - & 181 208 epsd ) ) ) * tms(ji,jj) … … 186 213 CALL lbc_lnk( zpresh(:,:), 'T', 1. ) 187 214 188 ! Ice strength (zpreshc) on grid cell corners (needed for 189 ! calculation of shear stress 215 ! Claude sees a bug 216 ! CALL lbc_lnk( tmi(:,:), 'T', 1. ) 217 218 ! Ice strength on grid cell corners (zpreshc) 219 ! needed for calculation of shear stress 190 220 DO jj = k_j1+1, k_jpj-1 191 221 DO ji = 2, jpim1 192 zstms = tms(ji+1,jj+1) * wght(ji+1,jj+1,2,2) + tms(ji,jj+1) * wght(ji+1,jj+1,1,2) & 193 & + tms(ji+1,jj ) * wght(ji+1,jj+1,2,1) + tms(ji,jj ) * wght(ji+1,jj+1,1,1) 222 zstms = tms(ji+1,jj+1) * wght(ji+1,jj+1,2,2) + & 223 & tms(ji,jj+1) * wght(ji+1,jj+1,1,2) + & 224 & tms(ji+1,jj) * wght(ji+1,jj+1,2,1) + & 225 & tms(ji,jj) * wght(ji+1,jj+1,1,1) 194 226 zusw(ji,jj) = 1.0 / MAX( zstms, epsd ) 195 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) & 196 & + zpresh(ji+1,jj ) * wght(ji+1,jj+1,2,1) + zpresh(ji,jj ) * wght(ji+1,jj+1,1,1) ) & 197 & * zusw(ji,jj) 227 zpreshc(ji,jj) = ( zpresh(ji+1,jj+1) * wght(ji+1,jj+1,2,2) + & 228 & zpresh(ji,jj+1) * wght(ji+1,jj+1,1,2) + & 229 & zpresh(ji+1,jj) * wght(ji+1,jj+1,2,1) + & 230 & zpresh(ji,jj) * wght(ji+1,jj+1,1,1) & 231 & ) * zusw(ji,jj) 198 232 END DO 199 233 END DO … … 233 267 234 268 ! Mass, coriolis coeff. and currents 235 zmass1(ji,jj) = (zt12*zc1(ji,jj)+zt11*zc1(ji+1,jj))/(zt11+zt12+epsd) 236 zmass2(ji,jj) = (zt22*zc1(ji,jj)+zt21*zc1(ji,jj+1))/(zt21+zt22+epsd) 237 238 zcorl1(ji,jj) = zmass1(ji,jj)*(e1t(ji+1,jj)*fcor(ji,jj)+e1t(ji,jj)*fcor(ji+1,jj))/(e1t(ji,jj)+e1t(ji+1,jj)+epsd) 239 zcorl2(ji,jj) = zmass2(ji,jj)*(e2t(ji,jj+1)*fcor(ji,jj)+e2t(ji,jj)*fcor(ji,jj+1))/(e2t(ji,jj+1)+e2t(ji,jj)+epsd) 240 ! 241 ! Reminder: since this is a C grid, the location of ocean currents 242 ! calculated by OPA should coincide with ice model vector points: 243 ! no need for interpolation once the definition of u_io and v_io 244 ! has been changed in module "icestp". Arrays u_oce1 and v_oce2 could 245 ! then be replaced by u_oce and v_oce, respectively. 269 zmass1(ji,jj) = ( zt12*zc1(ji,jj) + zt11*zc1(ji+1,jj) ) / (zt11+zt12+epsd) 270 zmass2(ji,jj) = ( zt22*zc1(ji,jj) + zt21*zc1(ji,jj+1) ) / (zt21+zt22+epsd) 271 zcorl1(ji,jj) = zmass1(ji,jj) * ( e1t(ji+1,jj)*fcor(ji,jj) + & 272 e1t(ji,jj)*fcor(ji+1,jj) ) & 273 / (e1t(ji,jj) + e1t(ji+1,jj) + epsd ) 274 zcorl2(ji,jj) = zmass2(ji,jj) * ( e2t(ji,jj+1)*fcor(ji,jj) + & 275 e2t(ji,jj)*fcor(ji,jj+1) ) & 276 / ( e2t(ji,jj+1) + e2t(ji,jj) + epsd ) 246 277 ! 247 278 u_oce1(ji,jj) = u_io(ji,jj) 248 249 ! SB modif because ocean has no slip boundary condition 279 v_oce2(ji,jj) = v_io(ji,jj) 280 281 ! Ocean has no slip boundary condition 250 282 v_oce1(ji,jj) = 0.5*( (v_io(ji,jj)+v_io(ji,jj-1))*e1t(ji+1,jj) & 251 283 & +(v_io(ji+1,jj)+v_io(ji+1,jj-1))*e1t(ji,jj)) & 252 284 & /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj) 253 285 254 u_oce2(ji,jj) = 0.5*( (u_io(ji,jj)+u_io(ji-1,jj))*e2t(ji,jj+1)&286 u_oce2(ji,jj) = 0.5*((u_io(ji,jj)+u_io(ji-1,jj))*e2t(ji,jj+1) & 255 287 & +(u_io(ji,jj+1)+u_io(ji-1,jj+1))*e2t(ji,jj)) & 256 & /(e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj) 257 258 v_oce2(ji,jj) = v_io(ji,jj) 288 & / (e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj) 259 289 260 290 ! Wind stress. … … 262 292 ztagny = ( 1. - zfrld2(ji,jj) ) * gtauy(ji,jj) 263 293 264 ! Computation of the velocity field taking into account the ice -ice interaction.294 ! Computation of the velocity field taking into account the ice internal interaction. 265 295 ! Terms that are independent of the velocity field. 266 ! Reminder: does zcorl*(-v_oce,u_oce) represent the surface pressure gradient? Why is it still267 ! calculated in this way? An ocean model with a free surface should provide the gradient268 ! of surface elevation directly...269 296 270 297 ! SB On utilise maintenant le gradient de la pente de l'ocean … … 272 299 ! zdsshx = (ssh_io(ji+1,jj) - ssh_io(ji,jj))/e1u(ji,jj) 273 300 ! zdsshy = (ssh_io(ji,jj+1) - ssh_io(ji,jj))/e2v(ji,jj) 301 274 302 zdsshx = 0.0 275 303 zdsshy = 0.0 … … 283 311 ! 284 312 !------------------------------------------------------------------------------! 285 ! 3) Solution of the momentum equation 313 ! 3) Solution of the momentum equation, iterative procedure 286 314 !------------------------------------------------------------------------------! 287 315 ! … … 319 347 END DO 320 348 321 322 349 DO jj = k_j1+1, k_jpj-1 350 DO ji = 2, jpim1 323 351 324 352 ! … … 327 355 !- zds(:,:): shear on northeast corner of grid cells 328 356 ! 329 !- IMPORTANT REMINDER: note that, the way these terms are coded, there are many repeated 330 ! calculations. Speed could be improved by regrouping terms. For 357 !- IMPORTANT REMINDER: Dear Gurvan, note that, the way these terms are coded, 358 ! there are many repeated calculations. 359 ! Speed could be improved by regrouping terms. For 331 360 ! the moment, however, the stress is on clarity of coding to avoid 332 ! bugs (mamm). 333 !- ALSO: arrays zdd, zdt, zds and delta could be removed in the future to minimise memory demand. 361 ! bugs (Martin, for Miguel). 362 ! 363 !- ALSO: arrays zdd, zdt, zds and delta could 364 ! be removed in the future to minimise memory demand. 365 ! 334 366 !- MORE NOTES: Note that we are calculating deformation rates and stresses on the corners of 335 367 ! grid cells, exactly as in the B grid case. For simplicity, the indexation on … … 337 369 ! 338 370 ! 339 ! (ji,jj)340 ! |341 ! |342 ! (ji-1,jj) | (ji,jj)343 ! ---------344 ! | |345 ! | (ji,jj) |------(ji,jj)346 ! | |347 ! ---------348 !(ji-1,jj-1) (ji,jj-1)349 !350 351 371 zdd(ji,jj) = ( e2u(ji,jj)*u_ice(ji,jj) & 352 372 & -e2u(ji-1,jj)*u_ice(ji-1,jj) & … … 366 386 367 387 ! 368 ! SB modif because ocean has no slip boundary condition369 388 zds(ji,jj) = ( ( u_ice(ji,jj+1)/e1u(ji,jj+1) & 370 389 & -u_ice(ji,jj)/e1u(ji,jj) & … … 379 398 380 399 381 ! SB modif because need to compute zddc and zdtc correctly382 400 v_ice1(ji,jj) = 0.5*( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji+1,jj) & 383 401 & +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji,jj)) & … … 406 424 & / area(ji,jj) 407 425 408 ! Old lines409 ! deltat(ji,jj) = SQRT( zdd(ji,jj)*zdd(ji,jj) + &410 ! & ( zdt(ji,jj)*zdt(ji,jj) + zdst*zdst ) * usecc2 &411 ! & ) + creepl412 ! New lines413 ! Regularisation of dP/dx414 426 delta = SQRT( zdd(ji,jj)*zdd(ji,jj) + & 415 427 & ( zdt(ji,jj)*zdt(ji,jj) + zdst*zdst ) * usecc2 ) 416 428 deltat(ji,jj) = MAX( SQRT(zdd(ji,jj)**2 + & 417 429 (zdt(ji,jj)**2 + zdst**2)*usecc2), creepl ) 418 ! End of new lines419 430 420 431 !-Calculate stress tensor components zs1 and zs2 421 432 !-at centre of grid cells (see section 3.5 of CICE user's guide). 422 ! Old lines423 ! zs1(ji,jj) = ( zs1(ji,jj) &424 ! & - dtotel*((1.0-alphaevp)*zs1(ji,jj)+(1.0-zdd(ji,jj)/deltat(ji,jj))*zpresh(ji,jj))) &425 ! & /(1.0+alphaevp*dtotel)426 ! New lines427 433 zs1(ji,jj) = ( zs1(ji,jj) & 428 434 & - dtotel*( ( 1.0 - alphaevp) * zs1(ji,jj) + & 429 & ( delta / deltat(ji,jj) - zdd(ji,jj) / deltat(ji,jj) ) * zpresh(ji,jj) ) ) & 435 & ( delta / deltat(ji,jj) - zdd(ji,jj) / deltat(ji,jj) ) & 436 * zpresh(ji,jj) ) ) & 430 437 & / ( 1.0 + alphaevp * dtotel ) 431 ! End of new lines 432 433 zs2(ji,jj) = ( zs2(ji,jj)&434 & - dtotel*((1.0-alphaevp)*ecc2*zs2(ji,jj)-zdt(ji,jj)/deltat(ji,jj)*zpresh(ji,jj)))&435 & / (1.0+alphaevp*ecc2*dtotel)438 439 zs2(ji,jj) = ( zs2(ji,jj) & 440 & - dtotel*((1.0-alphaevp)*ecc2*zs2(ji,jj) - & 441 zdt(ji,jj)/deltat(ji,jj)*zpresh(ji,jj)) ) & 442 & / ( 1.0 + alphaevp*ecc2*dtotel ) 436 443 437 444 END DO … … 443 450 DO jj = k_j1+1, k_jpj-1 444 451 DO ji = 2, jpim1 445 446 452 !- Calculate Delta on corners 447 448 453 zddc = ( ( v_ice1(ji,jj+1)/e1u(ji,jj+1) & 449 454 & -v_ice1(ji,jj)/e1u(ji,jj) & … … 455 460 & / ( e1f(ji,jj) * e2f(ji,jj) ) 456 461 457 458 462 zdtc = (-( v_ice1(ji,jj+1)/e1u(ji,jj+1) & 459 463 & -v_ice1(ji,jj)/e1u(ji,jj) & … … 465 469 & / ( e1f(ji,jj) * e2f(ji,jj) ) 466 470 467 deltac(ji,jj) = sqrt(zddc**2+(zdtc**2+zds(ji,jj)**2)*usecc2) + creepl471 deltac(ji,jj) = SQRT(zddc**2+(zdtc**2+zds(ji,jj)**2)*usecc2) + creepl 468 472 469 473 !-Calculate stress tensor component zs12 at corners (see section 3.5 of CICE user's guide). 470 471 zs12(ji,jj) = ( zs12(ji,jj)&472 & -dtotel*((1.0-alphaevp)*ecc2*zs12(ji,jj)-zds(ji,jj)/(2.0*deltac(ji,jj))*zpreshc(ji,jj)))&473 & / (1.0+alphaevp*ecc2*dtotel)474 475 END DO 476 END DO 474 zs12(ji,jj) = ( zs12(ji,jj) & 475 & - dtotel*( (1.0-alphaevp)*ecc2*zs12(ji,jj) - zds(ji,jj) / & 476 & ( 2.0*deltac(ji,jj) ) * zpreshc(ji,jj))) & 477 & / ( 1.0 + alphaevp*ecc2*dtotel ) 478 479 END DO ! ji 480 END DO ! jj 477 481 478 482 CALL lbc_lnk( zs12(:,:), 'F', 1. ) … … 481 485 DO jj = k_j1+1, k_jpj-1 482 486 DO ji = 2, jpim1 483 484 ! 485 !- contribution of zs1, zs2 and zs12 to zf1 486 ! 487 ! (ji,jj) 488 ! | 489 ! | 490 ! (ji-1,jj) | (ji,jj) 491 ! --------- 492 ! | | 493 ! | (ji,jj) |------(ji,jj) 494 ! | | 495 ! --------- 496 !(ji-1,jj-1) (ji,jj-1) 497 ! 498 499 zf1(ji,jj) = 0.5*( (zs1(ji+1,jj)-zs1(ji,jj))*e2u(ji,jj) & 500 & +(zs2(ji+1,jj)*e2t(ji+1,jj)**2-zs2(ji,jj)*e2t(ji,jj)**2)/e2u(ji,jj) & 501 & +2.0*(zs12(ji,jj)*e1f(ji,jj)**2-zs12(ji,jj-1)*e1f(ji,jj-1)**2)/e1u(ji,jj) & 502 & ) & 503 & /(e1u(ji,jj)*e2u(ji,jj)) 504 505 ! 506 !contribution of zs1, zs2 and zs12 to zf2 507 ! 508 ! (ji,jj) 509 ! | 510 ! | 511 ! (ji-1,jj) | (ji,jj) 512 ! --------- 513 ! | | 514 ! | (ji,jj) |------(ji,jj) 515 ! | | 516 ! --------- 517 !(ji-1,jj-1) (ji,jj-1) 518 ! 519 520 zf2(ji,jj) = 0.5*( (zs1(ji,jj+1)-zs1(ji,jj))*e1v(ji,jj) & 521 & -(zs2(ji,jj+1)*e1t(ji,jj+1)**2-zs2(ji,jj)*e1t(ji,jj)**2)/e1v(ji,jj) & 522 & +2.0*(zs12(ji,jj)*e2f(ji,jj)**2-zs12(ji-1,jj)*e2f(ji-1,jj)**2)/e2v(ji,jj) & 523 & ) & 524 & /(e1v(ji,jj)*e2v(ji,jj)) 525 487 !- contribution of zs1, zs2 and zs12 to zf1 488 zf1(ji,jj) = 0.5*( (zs1(ji+1,jj)-zs1(ji,jj))*e2u(ji,jj) & 489 & +(zs2(ji+1,jj)*e2t(ji+1,jj)**2-zs2(ji,jj)*e2t(ji,jj)**2)/e2u(ji,jj) & 490 & +2.0*(zs12(ji,jj)*e1f(ji,jj)**2-zs12(ji,jj-1)*e1f(ji,jj-1)**2)/e1u(ji,jj) & 491 & ) / ( e1u(ji,jj)*e2u(ji,jj) ) 492 ! contribution of zs1, zs2 and zs12 to zf2 493 zf2(ji,jj) = 0.5*( (zs1(ji,jj+1)-zs1(ji,jj))*e1v(ji,jj) & 494 & -(zs2(ji,jj+1)*e1t(ji,jj+1)**2 - zs2(ji,jj)*e1t(ji,jj)**2)/e1v(ji,jj) & 495 & + 2.0*(zs12(ji,jj)*e2f(ji,jj)**2 - & 496 zs12(ji-1,jj)*e2f(ji-1,jj)**2)/e2v(ji,jj) & 497 & ) / ( e1v(ji,jj)*e2v(ji,jj) ) 526 498 END DO 527 499 END DO … … 531 503 ! Both the Coriolis term and the ice-ocean drag are solved semi-implicitly. 532 504 ! 533 IF ( mod(iter,2).eq.0) THEN505 IF (MOD(iter,2).eq.0) THEN 534 506 535 507 DO jj = k_j1+1, k_jpj-1 536 508 DO ji = 2, jpim1 537 ! 538 ! (ji,jj) 539 ! | 540 ! | 541 ! (ji-1,jj) | (ji,jj) 542 ! --------- 543 ! | | 544 ! | (ji,jj) |------(ji,jj) 545 ! | | 546 ! --------- 547 !(ji-1,jj-1) (ji,jj-1) 548 ! 549 zmask = (1.0-max(rzero,sign(rone,-zmass1(ji,jj))))*tmu(ji,jj) 509 zmask = (1.0-MAX(rzero,SIGN(rone,-zmass1(ji,jj))))*tmu(ji,jj) 550 510 zsang = SIGN ( 1.0 , fcor(ji,jj) ) * sangvg 551 511 z0 = zmass1(ji,jj)/dtevp … … 555 515 & +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji,jj)) & 556 516 & /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj) 557 za = rhoco*sqrt((u_ice(ji,jj)-u_oce1(ji,jj))**2+(zv_ice1-v_oce1(ji,jj))**2)*(1.0-zfrld1(ji,jj)) 558 zr = z0*u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + za*(cangvg*u_oce1(ji,jj)-zsang*v_oce1(ji,jj)) 517 za = rhoco*SQRT((u_ice(ji,jj)-u_oce1(ji,jj))**2 + & 518 (zv_ice1-v_oce1(ji,jj))**2) * (1.0-zfrld1(ji,jj)) 519 zr = z0*u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + & 520 za*(cangvg*u_oce1(ji,jj)-zsang*v_oce1(ji,jj)) 559 521 zcca = z0+za*cangvg 560 522 zccb = zcorl1(ji,jj)+za*zsang 561 523 u_ice(ji,jj) = (zr+zccb*zv_ice1)/(zcca+epsd)*zmask 562 ! u_ice(ji,jj) = MAX( -1.0 , MIN( u_ice(ji,jj), 1.0 ) )563 524 564 525 END DO … … 569 530 DO jj = k_j1+1, k_jpj-1 570 531 DO ji = 2, jpim1 571 ! 572 ! (ji,jj) 573 ! | 574 ! | 575 ! (ji-1,jj) | (ji,jj) 576 ! --------- 577 ! | | 578 ! | (ji,jj) |------(ji,jj) 579 ! | | 580 ! --------- 581 !(ji-1,jj-1) (ji,jj-1) 582 ! 532 583 533 zmask = (1.0-max(rzero,sign(rone,-zmass2(ji,jj))))*tmv(ji,jj) 584 534 zsang = sign(1.0,fcor(ji,jj))*sangvg … … 588 538 & + (u_ice(ji,jj+1)+u_ice(ji-1,jj+1))*e2t(ji,jj)) & 589 539 & /(e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj) 590 za = rhoco*sqrt((zu_ice2-u_oce2(ji,jj))**2+(v_ice(ji,jj)-v_oce2(ji,jj))**2)*(1.0-zfrld2(ji,jj)) 591 zr = z0*v_ice(ji,jj) + zf2(ji,jj) + za2ct(ji,jj) + za*(cangvg*v_oce2(ji,jj)+zsang*u_oce2(ji,jj)) 540 za = rhoco*SQRT((zu_ice2-u_oce2(ji,jj))**2 + & 541 (v_ice(ji,jj)-v_oce2(ji,jj))**2)*(1.0-zfrld2(ji,jj)) 542 zr = z0*v_ice(ji,jj) + zf2(ji,jj) + & 543 za2ct(ji,jj) + za*(cangvg*v_oce2(ji,jj)+zsang*u_oce2(ji,jj)) 592 544 zcca = z0+za*cangvg 593 545 zccb = zcorl2(ji,jj)+za*zsang 594 546 v_ice(ji,jj) = (zr-zccb*zu_ice2)/(zcca+epsd)*zmask 595 ! v_ice(ji,jj) = MAX( -1.0 , MIN( v_ice(ji,jj), 1.0 ) )596 547 597 548 END DO … … 600 551 CALL lbc_lnk( v_ice(:,:), 'V', -1. ) 601 552 602 553 ELSE 603 554 DO jj = k_j1+1, k_jpj-1 604 555 DO ji = 2, jpim1 605 !606 ! (ji,jj)607 ! |608 ! |609 ! (ji-1,jj) | (ji,jj)610 ! ---------611 ! | |612 ! | (ji,jj) |------(ji,jj)613 ! | |614 ! ---------615 !(ji-1,jj-1) (ji,jj-1)616 !617 556 zmask = (1.0-max(rzero,sign(rone,-zmass2(ji,jj))))*tmv(ji,jj) 618 557 zsang = sign(1.0,fcor(ji,jj))*sangvg … … 623 562 & /(e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj) 624 563 625 za = rhoco*sqrt((zu_ice2-u_oce2(ji,jj))**2+(v_ice(ji,jj)-v_oce2(ji,jj))**2)*(1.0-zfrld2(ji,jj)) 626 zr = z0*v_ice(ji,jj) + zf2(ji,jj) + za2ct(ji,jj) + za*(cangvg*v_oce2(ji,jj)+zsang*u_oce2(ji,jj)) 564 za = rhoco*SQRT((zu_ice2-u_oce2(ji,jj))**2 + & 565 (v_ice(ji,jj)-v_oce2(ji,jj))**2)*(1.0-zfrld2(ji,jj)) 566 zr = z0*v_ice(ji,jj) + zf2(ji,jj) + & 567 za2ct(ji,jj) + za*(cangvg*v_oce2(ji,jj)+zsang*u_oce2(ji,jj)) 627 568 zcca = z0+za*cangvg 628 569 zccb = zcorl2(ji,jj)+za*zsang 629 570 v_ice(ji,jj) = (zr-zccb*zu_ice2)/(zcca+epsd)*zmask 630 ! v_ice(ji,jj) = MAX( -1.0 , MIN( v_ice(ji,jj), 1.0 ) )631 571 632 572 END DO … … 637 577 DO jj = k_j1+1, k_jpj-1 638 578 DO ji = 2, jpim1 639 ! 640 ! (ji,jj) 641 ! | 642 ! | 643 ! (ji-1,jj) | (ji,jj) 644 ! --------- 645 ! | | 646 ! | (ji,jj) |------(ji,jj) 647 ! | | 648 ! --------- 649 ! (ji-1,jj-1) (ji,jj-1) 650 ! 651 652 zmask = (1.0-max(rzero,sign(rone,-zmass1(ji,jj))))*tmu(ji,jj) 653 zsang = sign(1.0,fcor(ji,jj))*sangvg 579 zmask = (1.0-MAX(rzero,SIGN(rone,-zmass1(ji,jj))))*tmu(ji,jj) 580 zsang = SIGN(1.0,fcor(ji,jj))*sangvg 654 581 z0 = zmass1(ji,jj)/dtevp 655 582 ! SB modif because ocean has no slip boundary condition … … 658 585 & /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj) 659 586 660 za = rhoco*sqrt((u_ice(ji,jj)-u_oce1(ji,jj))**2+(zv_ice1-v_oce1(ji,jj))**2)*(1.0-zfrld1(ji,jj)) 661 zr = z0*u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + za*(cangvg*u_oce1(ji,jj)-zsang*v_oce1(ji,jj)) 587 za = rhoco*SQRT((u_ice(ji,jj)-u_oce1(ji,jj))**2 + & 588 (zv_ice1-v_oce1(ji,jj))**2)*(1.0-zfrld1(ji,jj)) 589 zr = z0*u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + & 590 za*(cangvg*u_oce1(ji,jj)-zsang*v_oce1(ji,jj)) 662 591 zcca = z0+za*cangvg 663 592 zccb = zcorl1(ji,jj)+za*zsang 664 593 u_ice(ji,jj) = (zr+zccb*zv_ice1)/(zcca+epsd)*zmask 665 ! u_ice(ji,jj) = MAX( -1.0 , MIN( u_ice(ji,jj), 1.0 ) ) 666 667 END DO 668 END DO 594 END DO ! ji 595 END DO ! jj 669 596 670 597 CALL lbc_lnk( u_ice(:,:), 'U', -1. ) 671 598 672 ENDIF 673 674 !--- Convergence test. 675 DO jj = k_j1+1 , k_jpj-1 676 zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ) , & 677 ABS( v_ice(:,jj) - zv_ice(:,jj) ) ) 678 END DO 679 zresm = MAXVAL( zresr( 1:jpi , k_j1+1:k_jpj-1 ) ) 680 681 !------------------------------------------- 682 ! Compute dynamical inputs of the itd model 683 !------------------------------------------- 684 ! Mean over all iterations 685 686 ! IF ( stress_mean_swi .EQ. 1 ) THEN 687 ! DO jj = k_j1+1, k_jpj-1 688 ! DO ji = 2, jpim1 689 ! divu_i(ji,jj) = divu_i(ji,jj) + zdd(ji,jj) / nevp 690 ! delta_i(ji,jj) = delta_i(ji,jj) + deltat (ji,jj) / nevp 691 ! shear_i(ji,jj) = shear_i(ji,jj) + zds (ji,jj) / nevp 692 ! END DO 693 ! END DO 694 ! ENDIF 599 ENDIF 600 601 !--- Convergence test. 602 DO jj = k_j1+1 , k_jpj-1 603 zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ) , & 604 ABS( v_ice(:,jj) - zv_ice(:,jj) ) ) 605 END DO 606 zresm = MAXVAL( zresr( 1:jpi , k_j1+1:k_jpj-1 ) ) 695 607 696 608 ! ! ==================== ! … … 698 610 ! ! ==================== ! 699 611 700 !------------------------- 701 ! Prevent high velocities 702 !------------------------- 612 ! 613 !------------------------------------------------------------------------------! 614 ! 4) Prevent ice velocities when the ice is thin 615 !------------------------------------------------------------------------------! 616 ! 703 617 ! If the ice thickness is below 1cm then ice velocity should equal the 704 ! ocean velocity 618 ! ocean velocity, 619 ! This prevents high velocity when ice is thin 705 620 DO jj = k_j1+1, k_jpj-1 706 621 DO ji = 2, jpim1 … … 713 628 END DO 714 629 END DO 715 630 ! 631 !------------------------------------------------------------------------------! 632 ! 5) Compute stress rain invariants and store strain tensor 633 !------------------------------------------------------------------------------! 634 ! 635 ! Compute delta, shear and div, inputs for mechanical redistribution 716 636 DO jj = k_j1+1, k_jpj-1 717 637 DO ji = 2, jpim1 718 ! Recompute stress tensor invariants719 638 !- zdd(:,:), zdt(:,:): divergence and tension at centre 720 ! of grid cells721 639 !- zds(:,:): shear on northeast corner of grid cells 722 723 640 zindb = MAX( 0.0, SIGN( 1.0, at_i(ji,jj) - 1.0e-6 ) ) 724 641 zdummy = zindb * vt_i(ji,jj) / MAX(at_i(ji,jj) , 1.0e-06 ) … … 782 699 END DO !ji 783 700 784 ! dynamical invariants 701 ! MV This should be removed as it is the same as previous lines 702 ! ! dynamical invariants 785 703 ! IF ( stress_mean_swi .EQ. 0 ) THEN 704 ! the following should not be necessary 786 705 DO jj = k_j1+1, k_jpj-1 787 706 DO ji = 2, jpim1 … … 797 716 CALL lbc_lnk( shear_i(:,:), 'F', 1. ) 798 717 799 IF(ln_ctl) THEN 800 WRITE(charout,FMT="('lim_rhg : res =',D23.16, ' iter =',I4)") zresm, jter 801 CALL prt_ctl_info(charout) 802 CALL prt_ctl(tab2d_1=u_ice, clinfo1=' lim_rhg : u_ice :', tab2d_2=v_ice, clinfo2=' v_ice :') 803 ENDIF 804 805 ! Stress tensor 718 719 ! Store the stress tensor for next time step 720 ! this accelerates convergence and improves stability 806 721 IF ( stress_mean_swi .EQ. 1 ) THEN 807 722 DO jj = k_j1+1, k_jpj-1 … … 814 729 ENDIF 815 730 816 ! Ice internal stresses731 ! Lateral boundary condition 817 732 CALL lbc_lnk( stress1_i(:,:), 'T', 1. ) 818 733 CALL lbc_lnk( stress2_i(:,:), 'T', 1. ) 819 734 CALL lbc_lnk( stress12_i(:,:), 'F', 1. ) 820 735 821 !------------------------ 822 ! Write charge ellipse 823 !------------------------ 824 736 ! 737 !------------------------------------------------------------------------------! 738 ! 6) Control prints of residual and charge ellipse 739 !------------------------------------------------------------------------------! 740 ! 741 ! print the residual for convergence 742 IF(ln_ctl) THEN 743 WRITE(charout,FMT="('lim_rhg : res =',D23.16, ' iter =',I4)") zresm, iter 744 CALL prt_ctl_info(charout) 745 CALL prt_ctl(tab2d_1=u_ice, clinfo1=' lim_rhg : u_ice :', tab2d_2=v_ice, clinfo2=' v_ice :') 746 ENDIF 747 748 ! print charge ellipse 749 ! This can be desactivated once the user is sure that the stress state 750 ! lie on the charge ellipse. See Bouillon et al. 08 for more details 825 751 IF(ln_ctl) THEN 826 752 CALL prt_ctl_info('lim_rhg : numit :',ivar1=numit) -
trunk/NEMO/LIM_SRC_3/limrst.F90
r825 r834 6 6 #if defined key_lim3 7 7 !!---------------------------------------------------------------------- 8 !! 'key_lim3' : LIMsea-ice model8 !! 'key_lim3' : LIM3 sea-ice model 9 9 !!---------------------------------------------------------------------- 10 10 !! lim_rst_write : write of the restart file … … 19 19 USE ice_oce ! ice variables 20 20 USE daymod 21 USE limicepoints22 21 !USE limvar 23 22 … … 65 64 llbon 66 65 INTEGER :: & 67 ji, jj, jk , jl, jm, index66 ji, jj, jk , jl 68 67 INTEGER :: & 69 68 inumwrs, it0, itime -
trunk/NEMO/LIM_SRC_3/limtab.F90
r825 r834 5 5 !!====================================================================== 6 6 #if defined key_lim3 7 !!---------------------------------------------------------------------- 8 !! 'key_lim3' LIM3 sea-ice model 7 9 !!---------------------------------------------------------------------- 8 10 !! tab_2d_1d : 2-D to 1-D … … 20 22 21 23 !!---------------------------------------------------------------------- 22 !! LIM 2.0, UCL-LOCEAN-IPSL (2005)24 !! LIM 3.0, UCL-ASTR-LOCEAN-IPSL (2008) 23 25 !! $Header: /home/opalod/NEMOCVSROOT/NEMO/LIM_SRC/limtab.F90,v 1.2 2005/03/27 18:34:42 opalod Exp $ 24 26 !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt -
trunk/NEMO/LIM_SRC_3/limthd.F90
r825 r834 6 6 #if defined key_lim3 7 7 !!---------------------------------------------------------------------- 8 !! 'key_lim3' : LIMsea-ice model8 !! 'key_lim3' LIM3 sea-ice model 9 9 !!---------------------------------------------------------------------- 10 10 !! lim_thd : thermodynamic of sea ice … … 22 22 USE dom_ice ! LIM sea-ice domain 23 23 USE iceini 24 USE limicepoints25 24 USE limthd_dif 26 25 USE limthd_dh … … 84 83 !! * Local variables 85 84 INTEGER :: ji, jj, jk, jl, & 86 zjid , zjjd , &87 85 zji , zjj, & ! dummy loop indices 88 86 nbpb , & ! nb of icy pts for thermo. cal. … … 95 93 REAL(wp) :: & 96 94 zinda , & ! switch for test. the val. of concen. 97 zindb, zindg ,& ! switches for test. the val of arg98 z a , zh, zthsnice, &95 zindb, & ! switches for test. the val of arg 96 zthsnice , & 99 97 zfric_u , & ! friction velocity 100 98 zfnsol , & ! total non solar heat … … 104 102 105 103 REAL(wp) :: & 106 ztmelts , & ! ice layers melting point107 zaaa , & ! numeric value for 2nd degree polynomial108 zbbb , & ! " " " " " "109 zccc , & ! " " " " " "110 zdiscrim , & ! " " " " " "111 zq , & ! sea ice temporary energy of melting112 zqold , & ! test variable for ice heat content113 104 zareamin 114 105 … … 224 215 225 216 ! still need to be updated : fdtcn !!!! 226 ! !-- Lead heat budget (part 1, next one is in limthd_dh , then limthd_lab217 ! !-- Lead heat budget (part 1, next one is in limthd_dh 227 218 ! !-- qldif -- (or qldif_1d in 1d routines) 228 219 zfontn = sprecip(ji,jj) * lfus ! energy of melting … … 510 501 REAL(wp) :: & 511 502 zdes, & ! snow heat content increment (dummy) 512 ztmelts, & ! layer melting point (dummy)513 503 zeps ! very small value (1.e-10) 514 504 … … 560 550 jl !: category number 561 551 562 REAL(wp), DIMENSION(jpl) :: & !: ! goes to trash563 zdummy564 565 552 REAL(wp) :: & !: ! goes to trash 566 553 meance, & !: mean conservation error … … 573 560 INTEGER :: & 574 561 ji,jj,jk, & !: loop indices 575 zji, zjj , zjiindex562 zji, zjj 576 563 !!--------------------------------------------------------------------- 577 564 … … 745 732 jl !: category number 746 733 747 REAL(wp), DIMENSION(jpl) :: & !: ! goes to trash748 zdummy749 750 734 REAL(wp) :: & !: ! goes to trash 751 735 meance, & !: mean conservation error 752 max_cons_err, & !: maximum tolerated conservation error 753 max_surf_err !: maximum tolerated surface error 736 max_cons_err !: maximum tolerated conservation error 754 737 755 738 INTEGER :: & … … 758 741 INTEGER :: & 759 742 ji,jj,jk, & !: loop indices 760 zji, zjj , zjiindex743 zji, zjj 761 744 762 745 !!--------------------------------------------------------------------- … … 974 957 ENDIF 975 958 976 uscomi = 1.0 / ( 1.0 - amax ) ! inverse of minimum lead fraction977 959 rcdsn = hakdif * rcdsn 978 960 rcdic = hakdif * rcdic -
trunk/NEMO/LIM_SRC_3/limthd_dh.F90
r825 r834 1 1 MODULE limthd_dh 2 2 #if defined key_lim3 3 !!---------------------------------------------------------------------- 4 !! 'key_lim3' LIM3 sea-ice model 5 !!---------------------------------------------------------------------- 3 6 !!====================================================================== 4 7 !! *** MODULE limthd_dh *** … … 19 22 USE ice 20 23 USE par_ice 21 USE limicepoints22 24 23 25 IMPLICIT NONE … … 36 38 37 39 !!---------------------------------------------------------------------- 38 !! LIM 3.0, UCL-ASTR-LOCEAN-IPSL (200 5)40 !! LIM 3.0, UCL-ASTR-LOCEAN-IPSL (2008) 39 41 !!---------------------------------------------------------------------- 40 42 … … 70 72 !! ** References : Bitz and Lipscomb, JGR 99 71 73 !! Fichefet T. and M. Maqueda 1997, J. Geophys. Res., 102(C6), 12609-12646 72 !! Vancoppenolle Fichefet and Bitz, GRL 2005 74 !! Vancoppenolle, Fichefet and Bitz, GRL 2005 75 !! Vancoppenolle et al., OM08 73 76 !! 74 77 !! ** History : 75 78 !! original code 01-04 (LIM) 76 79 !! original routine 77 !! (05-2003) Martin Vancoppenolle, Louvain-La-Neuve, Belgium 78 !! (06/07-2005) Martin Vancoppenolle for the 3D version 80 !! (05-2003) M. Vancoppenolle, Louvain-La-Neuve, Belgium 81 !! (06/07-2005) 3D version 82 !! (03-2008) Clean code 79 83 !! 80 84 !!------------------------------------------------------------------ … … 90 94 jk , & !: ice layer index 91 95 isnow , & !: switch for presence (1) or absence (0) of snow 92 index , & !: ice point index (limicepoints)93 96 zji , & !: 2D corresponding indices to ji 94 97 zjj , & … … 107 110 zhnfi , & 108 111 zihg , & 109 zdhgm , &110 112 zdhnm , & 111 113 zhnnew , & 112 zdfseqv , &113 114 zeps = 1.0e-13, & 114 115 zhisn , & … … 116 117 !: entrapment 117 118 zds , & !: increment of bottom ice salinity 118 zoldsinew , & !: dummy argument for error diagnosis119 119 zcoeff , & !: dummy argument for snowfall partitioning 120 120 !: over ice and leads 121 zjiindex , & !: zdhs_temp, &122 121 zsm_snowice, & !: snow-ice salinity 123 122 zswi1 , & !: switch for computation of bottom salinity … … 159 158 160 159 REAL(wp), DIMENSION(jpij,jkmax) :: & 161 zqt_i_lay , & !: total ice heat content 162 zqt_s_lay !: total snow heat content 160 zqt_i_lay !: total ice heat content 163 161 164 162 ! Heat conservation … … 181 179 WRITE(numout,*) 'lim_thd_dh : computation of vertical snow/ice accretion/ablation' 182 180 WRITE(numout,*) '~~~~~~~~~' 181 183 182 zfsalt_melt(:) = 0.0 184 183 ftotal_fin(:) = 0.0 … … 211 210 dsm_i_se_1d(:) = 0.0 212 211 dsm_i_si_1d(:) = 0.0 213 ! ! heat conservation 214 ! IF ( con_i ) THEN 215 ! IF (jiindex_1d .GT.0) WRITE(numout,*) ' zf_surf : ', z_f_surf(jiindex_1d) 216 ! IF (jiindex_1d .GT.0) WRITE(numout,*) ' qnsr_ice : ', qnsr_ice_1d(jiindex_1d) 217 ! IF (jiindex_1d .GT.0) WRITE(numout,*) ' qsr_ice : ', qsr_ice_1d(jiindex_1d) 218 ! IF (jiindex_1d .GT.0) WRITE(numout,*) ' fc_su : ', fc_su(jiindex_1d) 219 ! IF (jiindex_1d .GT.0) WRITE(numout,*) ' i0 : ', i0(jiindex_1d) 220 ! ENDIF 221 212 ! 222 213 !------------------------------------------------------------------------------! 223 214 ! 2) Computing layer thicknesses and snow and sea-ice enthalpies. ! 224 215 !------------------------------------------------------------------------------! 225 226 !----------------- 216 ! 227 217 ! Layer thickness 228 !-----------------229 218 DO ji = kideb,kiut 230 219 zh_i(ji) = ht_i_b(ji) / nlay_i … … 232 221 END DO 233 222 234 !-------------------------------------- 235 ! Snow energy of melting, heat content 236 !-------------------------------------- 223 ! Total enthalpy of the snow 237 224 zqt_s(:) = 0.0 238 225 DO jk = 1, nlay_s 239 226 DO ji = kideb,kiut 240 ! this line could be removed241 ! q_s_b(ji,jk) = rhosn * ( cpic * ( rtt - t_s_b(ji,jk) ) + lfus )242 ! heat conservation243 227 zqt_s(ji) = zqt_s(ji) + q_s_b(ji,jk) * ht_s_b(ji) / nlay_s 244 228 END DO 245 229 END DO 246 230 247 ! ! heat conservation 248 ! qt_s_in(:,jl) = zqt_s(:) 249 ! IF (jiindex_1d.GT.0) & 250 ! WRITE(numout,*) 'jl : ', jl, ' zqts : ', zqt_s(jiindex_1d) / & 251 ! rdt_ice 252 253 ! IF (jiindex_1d.GT.0) THEN 254 ! WRITE(numout,*) ' Vertical profile : ' 255 ! WRITE(numout,*) ' ht_s_b : ', ht_s_b(jiindex_1d) 256 ! WRITE(numout,*) ' qm0 : ', q_s_b(jiindex_1d,1:nlay_s) * ht_s_b(jiindex_1d) / & 257 ! rdt_ice 258 ! WRITE(numout,*) ' q_s_b : ', q_s_b(jiindex_1d,1) 259 ! WRITE(numout,*) ' t_s_b : ', t_s_b(jiindex_1d,1) 260 ! WRITE(numout,*) ' zthick0: ', ht_s_b(jiindex_1d) 261 ! ENDIF 262 263 !------------------------------------- 264 ! Ice energy of melting, heat content 265 !------------------------------------- 231 ! Total enthalpy of the ice 266 232 zqt_i(:) = 0.0 267 233 DO jk = 1, nlay_i 268 234 DO ji = kideb,kiut 269 ! Melting point in K270 ! ztmelts = - tmut * s_i_b(ji,jk) + rtt271 ! this line could be removed272 ! q_i_b(ji,jk) = rhoic * &273 ! ( cpic * ( ztmelts - t_i_b(ji,jk) ) &274 ! + lfus * ( 1.0 - (ztmelts-rtt) / &275 ! MIN( ( t_i_b(ji,jk) - rtt ) , - zeps ) ) &276 ! - rcp * ( ztmelts-rtt ) )277 235 zqt_i(ji) = zqt_i(ji) + q_i_b(ji,jk) * ht_i_b(ji) / nlay_i 278 236 zqt_i_lay(ji,jk) = q_i_b(ji,jk) * ht_i_b(ji) / nlay_i 279 237 END DO 280 238 END DO 281 ! IF (jiindex_1d.GT.0) &282 ! WRITE(numout,*) 'jl : ', jl, ' zqti : ', zqt_i(jiindex_1d) / &283 ! rdt_ice284 285 ! IF (jiindex_1d.GT.0) THEN286 ! WRITE(numout,*) ' --- Vertical profile : '287 ! WRITE(numout,*) ' ht_s_b : ', ht_s_b(jiindex_1d)288 ! WRITE(numout,*) ' ht_i_b : ', ht_i_b(jiindex_1d)289 ! WRITE(numout,*) ' q_i_b : ', q_i_b(jiindex_1d, 1:nlay_i)290 ! WRITE(numout,*) ' t_i_b : ', t_i_b(jiindex_1d, 1:nlay_i)291 ! WRITE(numout,*) ' qm0 : ', zqt_i_lay(jiindex_1d,1:nlay_i) / rdt_ice292 ! WRITE(numout,*) ' zthick0: ', ht_i_b(jiindex_1d) / nlay_i293 ! ENDIF294 295 239 ! 296 240 !------------------------------------------------------------------------------| … … 298 242 !------------------------------------------------------------------------------| 299 243 ! 300 !--------------------- 301 ! Snow precips / melt302 !--------------------- 244 !------------------------- 245 ! 3.1 Snow precips / melt 246 !------------------------- 303 247 ! Snow accumulation in one thermodynamic time step 304 248 ! snowfall is partitionned between leads and ice … … 306 250 ! but because of the winds, more snow falls on leads than on sea ice 307 251 ! and a greater fraction (1-at_i)^beta of the total mass of snow 308 ! (beta < 1) falls in leads 252 ! (beta < 1) falls in leads. 309 253 ! In reality, beta depends on wind speed, 310 254 ! and should decrease with increasing wind speed but here, it is 311 ! considered as a constant 312 ! an average value is 0.66 255 ! considered as a constant. an average value is 0.66 313 256 ! Martin Vancoppenolle, December 2006 314 257 … … 322 265 ! Melt of fallen snow 323 266 DO ji = kideb, kiut 324 325 267 ! tatm_ice is now in K 326 268 zqprec(ji) = rhosn * ( cpic * ( rtt - tatm_ice_1d(ji) ) + lfus ) … … 334 276 qt_s_in(ji,jl) = qt_s_in(ji,jl) + zqprec(ji) * zdh_s_pre(ji) 335 277 zqt_s(ji) = zqt_s(ji) + zqprec(ji) * zdh_s_pre(ji) 336 337 END DO 338 339 ! IF (jiindex_1d .GT. 0) THEN 340 ! WRITE(numout,*) ' zqfont_su : ', zqfont_su(jiindex_1d) 341 ! WRITE(numout,*) ' zqprec : ', zqprec(jiindex_1d)*( zdh_s_pre(jiindex_1d) + zdh_s_mel(jiindex_1d) ) / rdt_ice 342 ! WRITE(numout,*) ' tatm : ', tatm_ice_1d(jiindex_1d) 343 ! WRITE(numout,*) 'jl : ', jl, ' zqts : ', zqt_s(jiindex_1d) / & 344 ! rdt_ice 345 ! ENDIF 346 347 ! updated total snow heat content 278 END DO 279 280 ! Update total snow heat content 348 281 zqt_s(ji) = MAX ( zqt_s(ji) - zqfont_su(ji) , 0.0 ) 349 350 ! IF (jiindex_1D .GT. 0) &351 ! WRITE(numout,*) 'jl : ', jl, ' zqts : ', zqt_s(jiindex_1d) / &352 ! rdt_ice353 282 354 283 ! Snow melt due to surface heat imbalance … … 358 287 zqfont_su(ji) = MAX( 0.0 , - zh_s(ji) - zdeltah(ji,jk) ) * & 359 288 q_s_b(ji,jk) 360 ! IF (ji.eq.jiindex_1d) WRITE(numout,*) ' zqfont_su : ', &361 ! zqfont_su(jiindex_1d)362 289 zdeltah(ji,jk) = MAX( zdeltah(ji,jk) , - zh_s(ji) ) 363 290 zdh_s_mel(ji) = zdh_s_mel(ji) + zdeltah(ji,jk) !resulting melt of snow 364 291 END DO 365 292 END DO 366 367 ! !+++++368 ! IF (jiindex_1d .GT. 0 ) THEN369 ! WRITE(numout,*) ' dhs_pre :', zdh_s_pre(jiindex_1d)370 ! WRITE(numout,*) ' dhs_mel :', zdh_s_mel(jiindex_1d)371 ! WRITE(numout,*) ' zqfont_su : ', zqfont_su(jiindex_1d)372 ! ENDIF373 ! !+++++374 293 375 294 ! Apply snow melt to snow depth … … 389 308 END DO ! ji 390 309 391 ! !+++++ 392 ! IF (jiindex_1d .GT. 1) WRITE(numout,*) ' dhs_tot :', dh_s_tot(jiindex_1d) 393 ! IF (jiindex_1d .GT. 1) WRITE(numout,*) ' zqfont_su :', zqfont_su(jiindex_1d) 394 ! !+++++ 395 396 !---------------------- 397 ! Surface ice ablation 398 !---------------------- 310 !-------------------------- 311 ! 3.2 Surface ice ablation 312 !-------------------------- 399 313 DO ji = kideb, kiut 400 ! IF (ji.eq.jiindex_1d) WRITE(numout,*) ' zqfont_su : ', &401 ! zqfont_su(jiindex_1d)402 314 dh_i_surf(ji) = 0.0 403 ! Heat conservation test315 ! For heat conservation test 404 316 z_f_surf(ji) = zqfont_su(ji) / rdt_ice ! heat conservation test 405 317 zdq_i(ji) = 0.0 … … 408 320 DO jk = 1, nlay_i 409 321 DO ji = kideb, kiut 322 ! melt of layer jk 410 323 zdeltah(ji,jk) = - zqfont_su(ji) / q_i_b(ji,jk) 324 ! recompute heat available 411 325 zqfont_su(ji) = MAX( 0.0 , - zh_i(ji) - zdeltah(ji,jk) ) * & 412 326 q_i_b(ji,jk) 327 ! melt of layer jk cannot be higher than its thickness 413 328 zdeltah(ji,jk) = MAX( zdeltah(ji,jk) , - zh_i(ji) ) 329 ! update surface melt 414 330 dh_i_surf(ji) = dh_i_surf(ji) + zdeltah(ji,jk) 331 ! for energy conservation 415 332 zdq_i(ji) = zdq_i(ji) + zdeltah(ji,jk) * & 416 333 q_i_b(ji,jk) / rdt_ice 417 ! Contribution to salt flux -> to change334 ! contribution to ice-ocean salt flux 418 335 zji = MOD( npb(ji) - 1, jpi ) + 1 419 336 zjj = ( npb(ji) - 1 ) / jpi + 1 420 337 zfsalt_melt(ji) = zfsalt_melt(ji) + & 421 ! ( sss_io(zji,zjj) - s_i_b(ji,jk) ) * &422 338 ( sss_io(zji,zjj) - sm_i_b(ji) ) * & 423 339 a_i_b(ji) * & 424 340 MIN( zdeltah(ji,jk) , 0.0 ) * rhoic / rdt_ice 425 !+++++++++++426 ! IF ( (zji.eq.jiindex) .AND. (zjj.eq.jjindex) ) THEN427 ! WRITE(numout,*) ' surf abl '428 ! WRITE(numout,*) ' zfsalt_melt : ', zfsalt_melt(ji) / ( sss_io(zji,zjj)+epsi16)429 ! WRITE(numout,*) ' sss_io : ', sss_io(zji,zjj)430 ! ENDIF431 !+++++++++++432 341 END DO ! ji 433 342 END DO ! jk … … 446 355 ENDIF 447 356 448 IF ( z_f_surf(ji) + zdq_i(ji) .GE. 1.0e-3 ) THEN!.GE. 1.0e-3 ) .AND. & 449 ! ( ( ht_i_b(ji) + dh_i_bott(ji) ) .GE. 1.0e-6 ) ) THEN 357 IF ( z_f_surf(ji) + zdq_i(ji) .GE. 1.0e-3 ) THEN! 450 358 WRITE(numout,*) ' ALERTE heat loss for surface melt ' 451 359 WRITE(numout,*) ' zji, zjj, jl :', zji, zjj, jl … … 471 379 ENDIF ! con_i 472 380 473 !------------------ 474 ! Snow sublimation 475 !------------------ 476 ! !++++++ 477 ! zqt_dummy(:) = 0.0 478 ! DO jk = 1, nlay_s 479 ! DO ji = kideb,kiut 480 ! q_s_b(ji,jk) = rhosn * ( cpic * ( rtt - t_s_b(ji,jk) ) + lfus ) 481 ! ! heat conservation 482 ! zqt_dummy(ji) = zqt_dummy(ji) + q_s_b(ji,jk) * old_ht_s_b(ji) / nlay_s 483 ! END DO 484 ! END DO 485 486 ! IF (jiindex_1d .GT. 0) & 487 ! WRITE(numout,*) 'jl : ', jl, ' Old zqts : ', zqt_dummy(jiindex_1d) / & 488 ! rdt_ice 489 !++++++ 490 ! !++++++ 491 ! zqt_dummy(:) = 0.0 492 ! DO jk = 1, nlay_s 493 ! DO ji = kideb,kiut 494 ! q_s_b(ji,jk) = rhosn * ( cpic * ( rtt - t_s_b(ji,jk) ) + lfus ) 495 ! ! heat conservation 496 ! zqt_dummy(ji) = zqt_dummy(ji) + q_s_b(ji,jk) * ht_s_b(ji) / nlay_s 497 ! END DO 498 ! END DO 499 ! IF (jiindex_1d .GT. 0) & 500 ! WRITE(numout,*) 'jl : ', jl, ' New zqts : ', zqt_dummy(jiindex_1d) / & 501 ! rdt_ice 502 ! !++++++ 381 !---------------------- 382 ! 3.3 Snow sublimation 383 !---------------------- 503 384 504 385 DO ji = kideb, kiut … … 510 391 ht_s_b(ji) = MAX( zzero , zdhcf ) 511 392 ! we recompute dh_s_tot 512 ! which may be different in case of total snow melt513 393 dh_s_tot(ji) = ht_s_b(ji) - zhsold(ji) 514 ! IF ( ht_s_b(ji) .LE. 0.0 ) THEN515 ! dh_s_tot(ji) = MAX( 0.0 , dh_s_tot(ji) )516 ! ENDIF517 518 394 qt_s_in(ji,jl) = qt_s_in(ji,jl) + zdh_s_sub(ji)*q_s_b(ji,1) 519 395 END DO !ji 520 396 521 ! IF ( jiindex_1d .GT. 0 ) THEN522 ! WRITE(numout,*) ' dh_s_sub : ', zdh_s_sub(jiindex_1d)523 ! WRITE(numout,*) ' dhs_tot : ', dh_s_tot(jiindex_1d)524 ! ENDIF525 526 !++++++527 397 zqt_dummy(:) = 0.0 528 398 DO jk = 1, nlay_s … … 533 403 END DO 534 404 END DO 535 ! IF (jiindex_1d .GT. 0) &536 ! WRITE(numout,*) 'jl : ', jl, ' Old zqts : ', zqt_dummy(jiindex_1d) / &537 ! rdt_ice538 ! !++++++539 540 ! IF (jiindex_1d .GT. 0) THEN541 ! WRITE(numout,*) ' t_s_b : ', t_s_b(jiindex_1d,1)542 ! WRITE(numout,*) ' q_s_b : ', q_s_b(jiindex_1d,1)543 ! ENDIF544 405 545 406 DO jk = 1, nlay_s !n … … 553 414 END DO 554 415 555 ! IF (jiindex_1d .GT. 0) THEN556 ! WRITE(numout,*) ' t_s_b : ', t_s_b(jiindex_1d,1)557 ! WRITE(numout,*) ' q_s_b : ', q_s_b(jiindex_1d,1)558 ! ENDIF559 560 416 ! 561 417 !------------------------------------------------------------------------------! … … 567 423 ! the inner conductive flux (fc_bo_i), from the open water heat flux 568 424 ! (qlbbqb) and the turbulent ocean flux (fbif). 569 ! fc_bo_i is positive downwards 570 ! fbif and qlbbq are positive to the ice 571 572 !--------------------------------------------- 573 ! Basal growth - salinity not varying in time 574 !--------------------------------------------- 575 IF ( num_sal .NE. 2 ) THEN 425 ! fc_bo_i is positive downwards. fbif and qlbbq are positive to the ice 426 427 !----------------------------------------------------- 428 ! 4.1 Basal growth - (a) salinity not varying in time 429 !----------------------------------------------------- 430 IF ( ( num_sal .NE. 2 ) .AND. ( num_sal .NE. 4 ) ) THEN 576 431 DO ji = kideb, kiut 577 432 IF ( ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) .LT. 0.0 ) THEN … … 590 445 dh_i_bott(ji) = - rdt_ice*( fc_bo_i(ji) + fbif_1d(ji) + & 591 446 qlbbq_1d(ji) ) / q_i_b(ji,nlay_i+1) 592 ! IF ( ji .EQ. jiindex_1d ) THEN593 ! WRITE(numout,*) ' s_i_new : ', s_i_new(ji)594 ! WRITE(numout,*) ' ztmelts : ', ztmelts595 ! WRITE(numout,*) ' t_bo : ', t_bo_b(ji)596 ! WRITE(numout,*) ' q_i_b : ', q_i_b(ji,nlay_i+1)597 ! WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(ji)598 ! ENDIF599 447 ENDIF ! heat budget 600 448 END DO ! ji 601 449 ENDIF ! num_sal 602 450 603 ! IF ( jiindex_1d .GT. 0 ) THEN 604 ! WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(jiindex_1d) 605 ! WRITE(numout,*) ' q_i_b : ', q_i_b(jiindex_1d,nlay_i+1) 606 ! WRITE(numout,*) ' fc_bo_i : ', fc_bo_i(jiindex_1d) 607 ! WRITE(numout,*) ' fbif_1d : ', fbif_1d(jiindex_1d) 608 ! WRITE(numout,*) ' qlbbq_1d : ', qlbbq_1d(jiindex_1d) 609 ! ENDIF 610 611 !----------------------------------------- 612 ! Basal growth - salinity varying in time 613 !----------------------------------------- 451 !------------------------------------------------- 452 ! 4.1 Basal growth - (b) salinity varying in time 453 !------------------------------------------------- 614 454 IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) THEN 615 455 ! the growth rate (dh_i_bott) is function of the new ice … … 617 457 ! salinity (snewice). snewice depends on dh_i_bott 618 458 ! it converges quickly, so, no problem 459 ! See Vancoppenolle et al., OM08 for more info on this 619 460 620 461 ! Initial value (tested 1D, can be anything between 1 and 20) … … 636 477 ( t_bo_b(ji) - rtt ) ) & 637 478 - rcp * ( ztmelts-rtt ) ) 638 ! !+++++++639 ! IF (ji.EQ.jiindex_1d) THEN640 ! WRITE(numout,*) ' ztmelts : ', ztmelts641 ! WRITE(numout,*) ' t_bo_b : ', t_bo_b(ji)642 ! WRITE(numout,*) ' q_i_b(ji,nlay_i+1) : ', q_i_b(ji,nlay_i+1)643 ! ENDIF644 ! !+++++++645 479 ! Bottom growth rate = - F*dt / q 646 480 dh_i_bott(ji) = - rdt_ice * ( fc_bo_i(ji) + fbif_1d(ji) & … … 650 484 ! zswi12 (1) if dh_i_bott/rdt .LT. 3.6e-7 and .GT. 2.0e-8 651 485 ! zswi1 (1) if dh_i_bott/rdt .LT. 2.0e-8 652 ! zgrr = MAX( dh_i_bott(ji) / rdt_ice, 0.0 )653 486 zgrr = MIN( 1.0e-3, MAX ( dh_i_bott(ji) / rdt_ice , zeps ) ) 654 487 zswi2 = MAX( zzero , SIGN( zone , zgrr - 3.6e-7 ) ) 655 488 zswi12 = MAX( zzero , SIGN( zone , zgrr - 2.0e-8 ) ) * ( 1.0 - zswi2 ) 656 489 zswi1 = 1. - zswi2 * zswi12 657 ! IF (ji.EQ.jiindex_1d) THEN658 ! WRITE(numout,*) ' zgrr : ', zgrr659 ! WRITE(numout,*) ' zswi2 : ', zswi2, ' zswi12 :', zswi12, ' zswi1 :', zswi1660 ! ENDIF661 490 zfracs = zswi1 * 0.12 + & 662 491 zswi12 * ( 0.8925 + 0.0568 * LOG( 100.0 * zgrr ) ) + & … … 673 502 IF ( ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) .LT. 0.0 ) THEN 674 503 ! New ice salinity must not exceed 15 psu 675 zoldsinew = s_i_new(ji)676 504 s_i_new(ji) = MIN( s_i_new(ji), s_i_max ) 677 505 ! Metling point in K … … 684 512 - rcp * ( ztmelts-rtt ) ) 685 513 ! Basal growth rate = - F*dt / q 686 ! sometimes, growth is very high because of very high bottom cond687 ! flux... this is dangerous688 514 dh_i_bott(ji) = - rdt_ice*( fc_bo_i(ji) + fbif_1d(ji) + & 689 515 qlbbq_1d(ji) ) / q_i_b(ji,nlay_i+1) 690 ! new lines 691 !----------------- 692 ! Salinity update 693 !----------------- 516 ! Salinity update 694 517 ! entrapment during bottom growth 695 518 dsm_i_se_1d(ji) = ( s_i_new(ji)*dh_i_bott(ji) + & … … 701 524 ENDIF ! num_sal 702 525 703 ! IF ( jiindex_1d .GT. 0 ) THEN 704 ! WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(jiindex_1d) 705 ! WRITE(numout,*) ' s_i_new : ', s_i_new(jiindex_1d) 706 ! ENDIF 707 708 !------------ 709 ! Basal melt 710 !------------ 711 !++++++ 526 !---------------- 527 ! 4.2 Basal melt 528 !---------------- 712 529 meance_dh = 0.0 713 530 numce_dh = 0 714 531 innermelt(:) = 0 715 !++++++716 532 717 533 DO ji = kideb, kiut … … 729 545 END DO 730 546 731 ! IF ( jiindex_1d .GT. 0 ) THEN732 ! WRITE(numout,*) ' zqfont_bo : ', zqfont_bo(jiindex_1d) / rdt_ice733 ! WRITE(numout,*) ' fc_bo_i : ', fc_bo_i(jiindex_1d)734 ! WRITE(numout,*) ' fbif : ', fbif_1d(jiindex_1d)735 ! WRITE(numout,*) ' qlbbq : ', qlbbq_1d(jiindex_1d)736 ! WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(jiindex_1d)737 ! ENDIF738 739 547 DO jk = nlay_i, 1, -1 740 548 DO ji = kideb, kiut 741 549 IF ( ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) .GE. 0.0 ) THEN 742 550 ztmelts = - tmut * s_i_b(ji,jk) + rtt 743 ! sometimes internal temperature gt melting point744 ! the whole layer is removed745 ! IF (ji.EQ.jiindex_1d) THEN746 ! WRITE(numout,*) ' jk : ', jk747 ! WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(ji)748 ! WRITE(numout,*) ' zqfont_bo : ', zqfont_bo(ji)749 ! WRITE(numout,*) ' t_i_b : ', t_i_b(ji,jk)750 ! WRITE(numout,*) ' q_i_b : ', q_i_b(ji,jk)751 ! ENDIF752 551 IF ( t_i_b(ji,jk) .GE. ztmelts ) THEN 753 552 zdeltah(ji,jk) = - zh_i(ji) 754 ! zqfont_bo(ji) = zqfont_bo(ji)755 553 dh_i_bott(ji) = dh_i_bott(ji) + zdeltah(ji,jk) 756 !++++++757 554 innermelt(ji) = 1 758 !++++++759 ! IF (ji.EQ.jiindex_1d) THEN760 ! WRITE(numout,*) ' Inner melt: ', innermelt(ji)761 ! WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(ji)762 ! WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(ji)763 ! ENDIF764 555 ELSE ! normal ablation 765 556 zdeltah(ji,jk) = - zqfont_bo(ji) / q_i_b(ji,jk) … … 774 565 zjj = ( npb(ji) - 1 ) / jpi + 1 775 566 zfsalt_melt(ji) = zfsalt_melt(ji) + & 776 ! ( sss_io(zji,zjj) - s_i_b(ji,jk) ) * &777 567 ( sss_io(zji,zjj) - sm_i_b(ji) ) * & 778 568 a_i_b(ji) * & 779 569 MIN( zdeltah(ji,jk) , 0.0 ) * rhoic / rdt_ice 780 ! IF (ji.EQ.jiindex_1d) THEN781 ! WRITE(numout,*) ' No inner melt '782 ! WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(ji)783 ! WRITE(numout,*) ' zqfont_bo : ', zqfont_bo(ji)784 ! WRITE(numout,*) ' t_i_b : ', t_i_b(ji,jk)785 ! WRITE(numout,*) ' q_i_b : ', q_i_b(ji,jk)786 ! WRITE(numout,*)787 ! ENDIF788 570 ENDIF 789 ! IF ( ji .EQ. jiindex_1d ) THEN790 ! WRITE(numout,*) ' basal melt '791 ! WRITE(numout,*) ' sss_io : ', sss_io(zji,zjj)792 ! WRITE(numout,*) ' zfsalt_melt : ', zfsalt_melt(ji)793 ! WRITE(numout,*) ' zfsalt_melt good units : ', zfsalt_melt(ji) / ( sss_io(zji,zjj) + epsi16 )794 ! ENDIF795 !+++++796 571 ENDIF 797 572 END DO ! ji 798 573 END DO ! jk 799 574 575 !------------------- 576 ! Conservation test 577 !------------------- 800 578 IF ( con_i ) THEN 801 579 DO ji = kideb, kiut 802 580 IF ( ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) .GE. 0.0 ) THEN 803 !-------------------804 ! Conservation test805 !-------------------806 581 IF ( ( zfbase(ji) + zdq_i(ji) ) .GE. 1.0e-3 ) THEN 807 582 numce_dh = numce_dh + 1 808 583 meance_dh = meance_dh + zfbase(ji) + zdq_i(ji) 809 584 ENDIF 810 IF ( zfbase(ji) + zdq_i(ji) .GE. 1.0e-3 ) THEN!.GE. 1.0e-3 ) .AND. & 811 ! ( ( ht_i_b(ji) + dh_i_bott(ji) ) .GE. 1.0e-6 ) ) THEN 585 IF ( zfbase(ji) + zdq_i(ji) .GE. 1.0e-3 ) THEN 812 586 WRITE(numout,*) ' ALERTE heat loss for basal melt ' 813 587 WRITE(numout,*) ' zji, zjj, jl :', zji, zjj, jl … … 834 608 ENDIF ! con_i 835 609 836 ! IF ( jiindex_1d .GT. 0 ) THEN837 ! WRITE(numout,*) ' zqfont_bo : ', zqfont_bo(jiindex_1d) / rdt_ice838 ! WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(jiindex_1d)839 ! ENDIF840 610 ! 841 611 !------------------------------------------------------------------------------! … … 843 613 !------------------------------------------------------------------------------! 844 614 ! 845 !------------------------------------------ 846 ! Excessive ablation in a 1-category model847 !------------------------------------------ 615 !---------------------------------------------- 616 ! 5.1 Excessive ablation in a 1-category model 617 !---------------------------------------------- 848 618 849 619 DO ji = kideb, kiut … … 872 642 END DO 873 643 874 ! IF (jiindex_1d .GT. 0 ) THEN 875 ! WRITE(numout,*) ' ht_i_b : ', ht_i_b(jiindex_1d) 876 ! WRITE(numout,*) ' ht_s_b : ', ht_s_b(jiindex_1d) 877 ! ENDIF 878 879 !-------------------------- 880 ! If too much ice melts 881 !-------------------------- 644 !----------------------------------- 645 ! 5.2 More than available ice melts 646 !----------------------------------- 882 647 ! then heat applied minus heat content at previous time step 883 648 ! should equal heat remaining … … 890 655 zhgnew(ji) = zihgnew * zhgnew(ji) ! ice thickness is put to 0 891 656 ! remaining heat 892 zfdt_final(ji) = ( 1.0 - zihgnew ) * ( zqfont_su(ji) + zqfont_bo(ji) ) ! & 893 ! + zihgnew * MAX ( zfdt_init(ji) - zqt_i(ji) - zqt_s(ji), 0.0 ) 894 895 ! !++++++++++ 896 ! IF (ji.EQ.jiindex_1d) THEN 897 ! WRITE(numout,*) ' zfdt_init : ', zfdt_init(ji) / rdt_ice 898 ! WRITE(numout,*) ' zqt_i : ', zqt_i(ji) / rdt_ice 899 ! WRITE(numout,*) ' 1. zqt_s : ', zqt_s(ji) / rdt_ice 900 ! WRITE(numout,*) ' zqt : ', ( zqt_i(ji) + zqt_s(ji) ) / rdt_ice 901 ! WRITE(numout,*) ' zhgnew : ', zhgnew(ji) 902 ! WRITE(numout,*) ' zihgnew : ', zihgnew 903 ! WRITE(numout,*) ' dh_i_surf : ', dh_i_surf(ji) 904 ! WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(ji) 905 ! WRITE(numout,*) ' ht_s_b : ', ht_s_b(ji) 906 ! WRITE(numout,*) ' 1. zfdt_final: ', zfdt_final(ji) / rdt_ice 907 ! WRITE(numout,*) ' zqfont_su : ', zqfont_su(jiindex_1d) / rdt_ice 908 ! WRITE(numout,*) ' zqfont_bo : ', zqfont_bo(jiindex_1d) / rdt_ice 909 ! ENDIF 910 ! !++++++++++ 657 zfdt_final(ji) = ( 1.0 - zihgnew ) * ( zqfont_su(ji) + zqfont_bo(ji) ) 911 658 912 659 ! If snow remains, energy is used to melt snow … … 923 670 zqt_s(ji) = zqt_s(ji) * ht_s_b(ji) 924 671 925 ! IF (ji.EQ.jiindex_1d) THEN926 ! WRITE(numout,*) ' 2. zfdt_final : ', zfdt_final(ji) / rdt_ice927 ! WRITE(numout,*) ' ht_s_b : ', ht_s_b(ji)928 ! WRITE(numout,*) ' zhgnew : ', zhgnew(ji)929 ! WRITE(numout,*) ' 2. zqt_s : ', zqt_s(ji) / rdt_ice930 ! WRITE(numout,*) ' zdhnm : ', zdhnm931 ! WRITE(numout,*)932 ! ENDIF933 934 672 ! Mass variations of ice and snow 935 673 !--------------------------------- … … 942 680 ! Remaining heat to the ocean 943 681 !--------------------------------- 944 ! check the switches here945 682 ! focea is in W.m-2 * dt 946 947 683 focea(ji) = - zfdt_final(ji) / rdt_ice 948 684 949 ! IF ( (zji.eq.jiindex) .AND. (zjj.eq.jjindex) ) THEN 950 ! WRITE(numout,*) ' focea : ', focea(ji) 951 ! ENDIF 952 953 END DO 954 955 ! IF (jiindex_1d .GT. 0 ) THEN 956 ! WRITE(numout,*) ' ht_i_b : ', ht_i_b(jiindex_1d) 957 ! WRITE(numout,*) ' ht_s_b : ', ht_s_b(jiindex_1d) 958 ! ENDIF 685 END DO 959 686 960 687 ftotal_fin (:) = zfdt_final(:) / rdt_ice 961 ! IF (jiindex_1d .GT. 1 ) WRITE(numout,*) ' ftotal_fin : ', ftotal_fin(jiindex_1d)962 688 963 689 !--------------------------- … … 986 712 focea(ji) * a_i_b(ji) * rdt_ice 987 713 988 ! ! astarojna989 714 zihic = 1.0 - MAX( zzero , SIGN( zone , -ht_i_b(ji) ) ) 990 715 ! equals 0 if ht_i = 0, 1 if ht_i gt 0 991 ! fstbif_1d est cumulatif merde992 716 fscbq_1d(ji) = a_i_b(ji) * fstbif_1d(ji) 993 ! here there is a bug994 717 qldif_1d(ji) = qldif_1d(ji) & 995 718 + fsup(ji) + ( 1.0 - zihgnew ) * focea(ji) * a_i_b(ji) & … … 1011 734 t_i_b(ji,jk) = zihgnew * t_i_b(ji,jk) + ( 1.0 - zihgnew ) * rtt 1012 735 q_i_b(ji,jk) = zihgnew * q_i_b(ji,jk) 1013 ! IF (ji.eq.jiindex_1d) THEN1014 ! WRITE(numout,*) ' jk : ', jk1015 ! WRITE(numout,*) ' q_i_b : ', q_i_b(ji,jk)1016 ! WRITE(numout,*) ' t_i_b : ', t_i_b(jiindex_1d, 1:nlay_i)1017 ! WRITE(numout,*) ' zihgnew : ', zihgnew, ' zhgnew : ', &1018 ! zhgnew(ji)1019 ! ENDIF1020 736 END DO 1021 737 END DO ! ji 1022 738 1023 ! IF (jiindex_1d.GT.0) THEN1024 ! WRITE(numout,*) ' --- Vertical profile : '1025 ! WRITE(numout,*) ' q_i_b : ', q_i_b(jiindex_1d, 1:nlay_i)1026 ! ENDIF1027 1028 739 DO ji = kideb, kiut 1029 740 ht_i_b(ji) = zhgnew(ji) 1030 741 END DO ! ji 1031 1032 ! IF (jiindex_1d .GT. 0 ) THEN1033 ! WRITE(numout,*) ' ht_i_b : ', ht_i_b(jiindex_1d)1034 ! WRITE(numout,*) ' ht_s_b : ', ht_s_b(jiindex_1d)1035 ! ENDIF1036 742 ! 1037 743 !------------------------------------------------------------------------------| … … 1065 771 zji = MOD( npb(ji) - 1, jpi ) + 1 1066 772 zjj = ( npb(ji) - 1 ) / jpi + 1 1067 ! zsm_snowice = MIN ( s_i_max, ( rhoic - rhosn ) / rhoic * &1068 ! sss_io(zji,zjj) )1069 773 1070 774 zsm_snowice = ( rhoic - rhosn ) / rhoic * & … … 1115 819 END DO !ji 1116 820 1117 ! IF (jiindex_1d .GT. 1 ) THEN1118 ! WRITE(numout,*) ' ht_i_b : ', ht_i_b(jiindex_1d)1119 ! WRITE(numout,*) ' ht_s_b : ', ht_s_b(jiindex_1d)1120 ! ENDIF1121 1122 821 END SUBROUTINE lim_thd_dh 1123 822 #else -
trunk/NEMO/LIM_SRC_3/limthd_dif.F90
r825 r834 1 1 MODULE limthd_dif 2 2 #if defined key_lim3 3 !!---------------------------------------------------------------------- 4 !! 'key_lim3' LIM3 sea-ice model 5 !!---------------------------------------------------------------------- 3 6 !!====================================================================== 4 7 !! *** MODULE limthd_dif *** … … 14 17 USE thd_ice 15 18 USE iceini 16 USE limicepoints17 19 USE limistate 18 20 USE in_out_manager … … 34 36 35 37 !!---------------------------------------------------------------------- 36 !! LIM 3.0, UCL-ASTR-LOCEAN-IPSL (200 5)38 !! LIM 3.0, UCL-ASTR-LOCEAN-IPSL (2008) 37 39 !! (c) UCL-ASTR and Martin Vancoppenolle 38 40 !!---------------------------------------------------------------------- … … 88 90 !! ** History : 89 91 !! (02-2003) Martin Vancoppenolle, Louvain-la-Neuve, Belgium 90 !! (06-2005) Martin Vancoppenolle still!!! for the 3d version...92 !! (06-2005) Martin Vancoppenolle, 3d version 91 93 !! (11-2006) Vectorized by Xavier Fettweis (UCL-ASTR) 92 94 !! (04-2007) Energy conservation tested by M. Vancoppenolle … … 106 108 layer, & ! vertical dummy loop index 107 109 nconv, & ! number of iterations in iterative procedure 108 ! nconvmax = 50, &! maximum number of iterations for temperature computations (50)109 index, & !110 110 minnumeqmin, & ! 111 111 maxnumeqmax … … 158 158 zfsw , & !solar radiation absorbed at the surface 159 159 zf , & ! surface flux function 160 dzf , & ! derivative of the surface flux function 161 zksn !effective snow conductivity 160 dzf ! derivative of the surface flux function 162 161 163 162 REAL(wp) :: & ! constant values … … 168 167 zbeta = 0.117, & !: for thermal conductivity (could be 0.13) 169 168 zraext_s = 1.0e08, & !: extinction coefficient of radiation in the snow 170 ! zraext_i = 1.5, & !extinction coefficient of radiation in the ice171 ! zraext_i = 1.0, & !extinction coefficient of radiation in the ice172 ! MY 0.87, FY 1.2 (Grenfell, Perovich and Light173 ! 07)174 ! zerrmax = 1.0e-4 , & !: error criterion for convergence of the iterative process175 zfc_int , & !: conductive heat flux at the interface176 169 zkimin = 0.10 , & !: minimum ice thermal conductivity 177 170 zht_smin = 1.0e-4 !: minimum snow depth 178 171 179 172 REAL(wp) :: & ! local variables 180 zheshth, & ! = zhe(ji) / thth181 173 ztmelt_i, & ! ice melting temperature 182 zerritmax, & ! current maximal error on temperature 183 zexp 174 zerritmax ! current maximal error on temperature 184 175 185 176 REAL(wp), DIMENSION(jpij) :: & 186 177 zerrit, & ! current error on temperature 187 178 zdifcase, & ! case of the equation resolution (1->4) 188 zghe, & ! correction factor of the thermal conductivity189 zhe, & ! effective thickness for compu. of equ. thermal conductivity190 179 zftrice, & ! solar radiation transmitted through the ice 191 zihic, zihe, zhsu 192 193 CHARACTER (len=50) :: charout 180 zihic, zhsu 181 194 182 !!-- End of declarations 195 183 !!---------------------------------------------------------------------------------------------- 196 184 197 ! new namelist parameters198 ! maxer_i_thd = maxer_i_thd199 ! thcon_i_swi = thcon_i_swi200 ! kappa_i = kappa_i201 ! nconv_i_thd = nconv_i_thd202 203 185 IF(lwp) WRITE(numout,*)'lim_thd_dif : Heat diffusion in sea ice for cat :', jl 204 186 205 thcon_i_swi = 1 ! Pringle et al. (2007)206 187 ! 207 188 !------------------------------------------------------------------------------! … … 339 320 zradtr_i(ji,nlay_i) 340 321 END DO 322 ! +++++ 323 341 324 DO layer = 1, nlay_i 342 325 DO ji = kideb , kiut … … 345 328 END DO 346 329 347 ! +++++348 330 349 331 ! … … 379 361 DO WHILE ((zerritmax > maxer_i_thd).AND.(nconv < nconv_i_thd)) 380 362 381 ! ! +++++382 ! DO ji = kideb , kiut383 ! zji = MOD( npb(ji) - 1, jpi ) + 1384 ! zjj = ( npb(ji) - 1 ) / jpi + 1385 ! IF ( (zji.EQ.jiindex) .AND. (zjj.EQ.jjindex) ) THEN386 ! WRITE(numout,*) ' Iteration : ', nconv387 ! WRITE(numout,*) ' t_su : ', t_su_b(ji)388 ! WRITE(numout,*) ' t_s : ', t_s_b(ji,1)389 ! WRITE(numout,*) ' t_i : ', t_i_b(ji,:)390 ! WRITE(numout,*) ' difcase : ', zdifcase(ji)391 ! ENDIF392 ! END DO393 ! ! +++++394 395 363 nconv = nconv+1 396 364 … … 410 378 411 379 IF ( thcon_i_swi .EQ. 1 ) THEN 412 ! Pringle et al formula included, should be tested380 ! Pringle et al formula included, 413 381 ! 2.11 + 0.09 S/T - 0.011.T 414 382 DO ji = kideb , kiut … … 489 457 zkappa_i(ji,0) = ztcond_i(ji,0)/MAX(zeps,zh_i(ji)) 490 458 zkappa_i(ji,nlay_i) = ztcond_i(ji,nlay_i) / MAX(zeps,zh_i(ji)) 491 492 459 !-- Interface 493 460 zkappa_s(ji,nlay_s) = 2.0*rcdsn*ztcond_i(ji,0)/MAX(zeps, & … … 616 583 numeqmin(ji) = 1 617 584 numeqmax(ji) = nlay_i + nlay_s + 1 585 618 586 !!surface equation 619 587 ztrid(ji,1,1) = 0.0 620 588 ztrid(ji,1,2) = dzf(ji) - zg1s*zkappa_s(ji,0) 621 589 ztrid(ji,1,3) = zg1s*zkappa_s(ji,0) 622 ! bug again ???623 ! zindterm(ji,1) = dzf(ji)*ztsuold(ji) - zf(ji)624 590 zindterm(ji,1) = dzf(ji)*t_su_b(ji) - zf(ji) 625 591 … … 853 819 END DO 854 820 821 !-------------------------! 822 ! Heat conservation ! 823 !-------------------------! 855 824 IF ( con_i ) THEN 856 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++857 ! ++++ heat conservation ++++ '858 ! Internal fluxes ... goes to trash!!!859 825 860 826 DO ji = kideb, kiut … … 884 850 zji = MOD( npb(ji) - 1, jpi ) + 1 885 851 zjj = ( npb(ji) - 1 ) / jpi + 1 886 ! IF ( ( zji.EQ.jiindex ) .AND. ( zjj.EQ.jjindex) ) THEN887 ! WRITE(numout,*) ' layer : ', layer888 ! WRITE(numout,*) ' fc_i : ', fc_i(ji,layer)889 ! WRITE(numout,*) ' zkappa_i : ', zkappa_i(ji,layer)890 ! WRITE(numout,*) ' ztcond_i : ', ztcond_i(ji,layer)891 ! WRITE(numout,*) ' zspeche_i: ', zspeche_i(ji,layer)892 ! WRITE(numout,*) ' ti2 : ', t_i_b(ji,layer+1)893 ! WRITE(numout,*) ' ti1 : ', t_i_b(ji,layer)894 ! ENDIF895 852 END DO 896 853 END DO … … 903 860 904 861 ENDIF 905 906 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++907 862 908 863 END SUBROUTINE lim_thd_dif -
trunk/NEMO/LIM_SRC_3/limthd_ent.F90
r825 r834 1 1 MODULE limthd_ent 2 2 #if defined key_lim3 3 !!---------------------------------------------------------------------- 4 !! 'key_lim3' LIM3 sea-ice model 5 !!---------------------------------------------------------------------- 3 6 !!====================================================================== 4 7 !! *** MODULE limthd_ent *** … … 7 10 !! after vertical growth/decay 8 11 !!====================================================================== 9 10 !!---------------------------------------------------------------------- 11 !! lim_thd_ent : ice temperature redistribution 12 !! lim_thd_ent : ice redistribution of enthalpy 12 13 13 14 !! * Modules used … … 22 23 USE limistate 23 24 USE ice 24 ! USE limthd25 25 USE limvar 26 26 USE par_ice 27 USE limicepoints28 27 29 28 IMPLICIT NONE … … 42 41 43 42 !!---------------------------------------------------------------------- 44 !! LIM 3.0, UCL-ASTR-LOCEAN-IPSL (200 5)43 !! LIM 3.0, UCL-ASTR-LOCEAN-IPSL (2008) 45 44 !!---------------------------------------------------------------------- 46 45 CONTAINS … … 59 58 !! ** Method : linear conservative remapping 60 59 !! 61 !! ** Steps : 1) Switches62 !! 2) S now redistribution63 !! 3) Ice enthalpyredistribution64 !! 4) Ice salinity65 !! 5) I nverse temperature computation from enthalpy60 !! ** Steps : 1) Grid 61 !! 2) Switches 62 !! 3) Snow redistribution 63 !! 4) Ice enthalpy redistribution 64 !! 5) Ice salinity, recover temperature 66 65 !! 67 66 !! ** Arguments … … 76 75 !! (07-2005) Martin for 3d adapatation 77 76 !! (11-2006) Vectorized by Xavier Fettweis (ASTR) 77 !! (03-2008) Energy conservation and clean code 78 78 !! * Arguments 79 79 … … 84 84 85 85 INTEGER :: & 86 ji,jk ,jk_a, & ! dummy loop indices87 zji, zjj , & !dummy indices86 ji,jk , & ! dummy loop indices 87 zji, zjj , & ! dummy indices 88 88 ntop0 , & ! old layer top index 89 89 nbot1 , & ! new layer bottom index … … 140 140 ! Energy conservation 141 141 REAL(wp), DIMENSION(jpij) :: & 142 zqti_in, zqts_in, zqt_in, & 143 zqti_fin, zqts_fin, zqt_fin 142 zqti_in, zqts_in, & 143 zqti_fin, zqts_fin 144 145 !------------------------------------------------------------------------------| 144 146 145 147 zeps = 1.0d-20 146 148 zeps6 = 1.0d-06 147 148 IF(lwp) THEN149 WRITE(numout,*)150 WRITE(numout,*) 'lim_thd_ent : Redistribution of energy in the ice '151 WRITE(numout,*) '~~~~~~~~~~~'152 WRITE(numout,*) ' kideb : ', kideb153 WRITE(numout,*) ' kiut : ', kiut154 ENDIF155 !156 !------------------------------------------------------------------------------|157 ! 1) Number of layers / Thickness of the layers |158 !------------------------------------------------------------------------------|159 !160 ! Number of layers depends on the ice thickness161 ! 1 layer if ice thinner than 30 cm162 ! 2 layers if ice thinner than 60 cm163 ! 4 layers otherwise164 165 nlays0 = nlay_s166 nlayi0 = nlay_i167 168 DO ji = kideb, kiut169 zh_i(ji) = old_ht_i_b(ji) / nlay_i170 zh_s(ji) = old_ht_s_b(ji) / nlay_s171 ENDDO172 173 149 zthick0(:,:) = 0.0 174 150 zm0(:,:) = 0.0 … … 178 154 z_i(:,:) = 0.0 179 155 z_s(:,:) = 0.0 156 157 ! 158 !------------------------------------------------------------------------------| 159 ! 1) Grid | 160 !------------------------------------------------------------------------------| 161 ! 162 nlays0 = nlay_s 163 nlayi0 = nlay_i 164 165 DO ji = kideb, kiut 166 zh_i(ji) = old_ht_i_b(ji) / nlay_i 167 zh_s(ji) = old_ht_s_b(ji) / nlay_s 168 ENDDO 169 180 170 ! 181 171 !------------------------------------------------------------------------------| … … 294 284 ENDDO 295 285 296 ! IF (jiindex_1d .GT. 0 ) THEN297 ! WRITE(numout,*) ' icsuind : ', icsuind(jiindex_1d)298 ! WRITE(numout,*) ' icsuswi : ', icsuswi(jiindex_1d)299 ! WRITE(numout,*) ' icboind : ', icboind(jiindex_1d)300 ! WRITE(numout,*) ' icboswi : ', icboswi(jiindex_1d)301 ! WRITE(numout,*) ' snicind : ', snicind(jiindex_1d)302 ! WRITE(numout,*) ' snicswi : ', snicswi(jiindex_1d)303 ! ENDIF304 286 ! 305 287 !------------------------------------------------------------------------------| … … 361 343 ENDDO 362 344 363 ! IF (jiindex_1d .GT. 0 ) THEN364 ! WRITE(numout,*) ' zqts_in : ', zqts_in(jiindex_1d) / rdt_ice365 ! WRITE(numout,*) ' q_s_b : ', q_s_b(jiindex_1d,1)366 ! WRITE(numout,*) ' t_s_b : ', t_s_b(jiindex_1d,1)367 ! WRITE(numout,*) ' zthick0: ', zthick0(jiindex_1d,1)368 ! ENDIF369 370 ! WRITE(numout,*) ' TEST : '371 ! WRITE(numout,*) ' maxnbot0 : ', maxnbot0372 ! WRITE(numout,*) ' kideb : ', kideb373 ! WRITE(numout,*) ' kiut : ', kiut374 345 DO jk = 2, maxnbot0 375 346 DO ji = kideb, kiut … … 381 352 zswitch = 1.0 - MAX (0.0, SIGN ( 1.0, zeps - ht_s_b(ji) ) ) 382 353 zqts_in(ji) = zqts_in(ji) + ( 1. - snswi(ji) ) * qm0(ji,jk) * zswitch 383 ! IF (ji.EQ.jiindex_1d) THEN384 ! WRITE(numout,*) ' limsum ', limsum385 ! WRITE(numout,*) ' snswi : ', snswi(ji)386 ! WRITE(numout,*) ' snind : ', snind(ji)387 ! WRITE(numout,*) ' t_s_b : ', t_s_b(ji,1)388 ! WRITE(numout,*) ' q_s_b : ', rhosn * ( cpic * ( rtt - &389 ! t_s_b(ji,limsum) ) + lfus )390 ! WRITE(numout,*) ' zthick0: ', zthick0(jiindex_1d,jk)391 ! ENDIF392 354 END DO ! jk 393 355 ENDDO ! ji 394 395 ! Heat conservation396 ! IF (jiindex_1d .GT. 0 ) THEN397 ! WRITE(numout,*) ' zqts_in : ', zqts_in(jiindex_1d) / rdt_ice398 ! WRITE(numout,*) ' q_s_b : ', q_s_b(jiindex_1d,1)399 ! WRITE(numout,*) ' t_s_b : ', t_s_b(jiindex_1d,1)400 ! ENDIF401 356 402 357 !------------------------------------------------ … … 423 378 424 379 !------------------ 425 ! NEW SNOW PROFILE380 ! new snow profile 426 381 !------------------ 427 382 … … 460 415 END DO 461 416 ENDDO 462 463 ! IF (jiindex_1d .GT. 0) THEN464 ! WRITE(numout,*) ' qm0 : ', qm0(jiindex_1d,0:maxnbot0) / rdt_ice465 ! WRITE(numout,*) ' maxnbot0 : ', maxnbot0466 ! WRITE(numout,*) ' zm0 : ', zm0(jiindex_1d,0:maxnbot0)467 ! WRITE(numout,*) ' z_s : ', z_s(jiindex_1d,0:1)468 ! WRITE(numout,*) ' zhl0 : ', zhl0(jiindex_1d,0:maxnbot0)469 ! WRITE(numout,*) ' zhl0 : ', zhl0(jiindex_1d,0:maxnbot0)470 ! ENDIF471 417 472 418 !---------------- … … 480 426 q_s_b(ji,layer1) = q_s_b(ji,layer1) + zrl01(layer1,layer0)*qm0(ji,layer0) & 481 427 * MAX(0.0,SIGN(1.0,nbot0(ji)-layer0+zeps)) 482 ! IF (ji.EQ.jiindex_1d) THEN483 ! WRITE(numout,*) ' layer0, layer1 ', layer0, layer1484 ! WRITE(numout,*) ' q_s_b : ', q_s_b(ji,layer1) / rdt_ice485 ! WRITE(numout,*) ' nbot0 : ', nbot0(ji)486 ! WRITE(numout,*)487 ! ENDIF488 428 END DO 489 429 END DO … … 516 456 ENDIF 517 457 518 ! WRITE(numout,*) ' ht_s_b : ', ht_s_b(jiindex_1d)519 ! WRITE(numout,*) ' dh_s_tot : ', dh_s_tot(jiindex_1d)520 ! WRITE(numout,*) ' dh_snowi : ', dh_snowice(jiindex_1d)521 ! WRITE(numout,*) ' zqts_fin : ', zqts_fin(jiindex_1d) / rdt_ice522 ! WRITE(numout,*) ' difference : ', ( zqts_fin(jiindex_1d) - &523 ! zqts_in(jiindex_1d) ) / rdt_ice524 525 ! IF (jiindex_1d .GT. 1) THEN526 ! WRITE(numout,*) ' Vertical profile : '527 ! WRITE(numout,*) ' ht_s_b : ', ht_s_b(jiindex_1d)528 ! WRITE(numout,*) ' nbot0 : ', nbot0(jiindex_1d)529 ! WRITE(numout,*) ' qm0 : ', qm0(jiindex_1d,1:nbot0(jiindex_1d)) / &530 ! rdt_ice531 ! WRITE(numout,*) ' zthick0: ', zthick0(jiindex_1d,1:nbot0(jiindex_1d))532 ! WRITE(numout,*) ' zm0 : ', zm0(jiindex_1d,0:nbot0(jiindex_1d))533 ! ENDIF534 535 458 !--------------------- 536 459 ! Recover heat content … … 541 464 END DO !ji 542 465 ENDDO !jk 543 ! IF (jiindex_1d .GT. 1) THEN544 ! WRITE(numout,*) ' q_s_b : ', q_s_b(jiindex_1d,1)545 ! WRITE(numout,*) ' t_s_b : ', t_s_b(jiindex_1d,1)546 ! ENDIF547 466 548 467 !--------------------- 549 468 ! Recover temperature 550 469 !--------------------- 551 ! is it really necessary ???? since we don't use snow temperature anymore552 ! it is now since, lim_thd_glohec uses t_s_b and not q_s_b!!!553 ! it should not be necessary anymore554 555 470 zfac1 = 1. / ( rhosn * cpic ) 556 471 zfac2 = lfus / cpic … … 563 478 END DO 564 479 ENDDO 565 ! IF (jiindex_1d .GT. 1) THEN566 ! WRITE(numout,*) ' q_s_b : ', q_s_b(jiindex_1d,1)567 ! WRITE(numout,*) ' t_s_b : ', t_s_b(jiindex_1d,1)568 ! ENDIF569 480 ! 570 481 !------------------------------------------------------------------------------| … … 587 498 ( 1 - icsuind(ji) ) * icsuswi(ji) + snicswi(ji) , & 588 499 nlay_i + 2 ) ) 589 ! IF (ji.EQ.jiindex_1D) THEN590 ! WRITE(numout,*) ' nbot0 : ', nbot0(ji)591 ! WRITE(numout,*) ' nlayi0 : ', nlayi0592 ! WRITE(numout,*) ' icboind: ', icboind(ji)593 ! WRITE(numout,*) ' icsuind: ', icsuind(ji)594 ! WRITE(numout,*) ' icsuswi: ', icsuswi(ji)595 ! WRITE(numout,*) ' snicswi: ', snicswi(ji)596 ! WRITE(numout,*) ' jkmax : ', jkmax597 ! WRITE(numout,*) ' term : ', &598 ! nlayi0+(1-icboind(ji))+(1-icsuind(ji))*icsuswi(ji)+snicswi(ji)599 ! ENDIF600 500 ! maximum reference number of the bottommost layer over all domain 601 501 maxnbot0 = MAX( maxnbot0 , nbot0(ji) ) … … 612 512 ! the ice layer number goes from 1 to nlay_i 613 513 ! limsum is the real ice layer number corresponding to present jk 614 limsum = ( (icsuswi(ji)*(icsuind(ji)+jk-1) + (1-icsuswi(ji))*jk))*(1-snicswi(ji)) + (jk-1)*snicswi(ji) 514 limsum = ( (icsuswi(ji)*(icsuind(ji)+jk-1) + & 515 (1-icsuswi(ji))*jk))*(1-snicswi(ji)) + (jk-1)*snicswi(ji) 615 516 zm0(ji,jk)= icsuswi(ji)*dh_i_surf(ji) + snicswi(ji)*dh_snowice(ji) & 616 517 + limsum * zh_i(ji) … … 647 548 MIN((t_i_b(ji,limsum)-rtt),-zeps) ) - rcp*(ztmelts-rtt) ) & 648 549 * zthick0(ji,jk) 649 650 ! IF (ji .EQ. jiindex_1d ) THEN651 ! WRITE(numout,*) 'jk : ', jk, ' t_i_b : ', t_i_b(ji,limsum), &652 ! ' s_i_b : ', s_i_b(ji,limsum)653 ! WRITE(numout,*) ' q_i_b : ', &654 ! rhoic*(cpic*(ztmelts-t_i_b(ji,limsum)) + lfus *( &655 ! 1.0-(ztmelts-rtt)/ &656 ! MIN((t_i_b(ji,limsum)-rtt),-zeps) ) - &657 ! rcp*(ztmelts-rtt) )658 ! WRITE(numout,*) ' zthick0: ', zthick0(ji,jk)659 ! WRITE(numout,*) ' limsum : ', limsum660 ! WRITE(numout,*) ' icsuswi: ', icsuswi(ji), ' icsuind : ', &661 ! icsuind(ji)662 ! WRITE(numout,*) ' snicswi: ', snicswi(ji)663 ! WRITE(numout,*) ' FAC1 : ', snicswi(ji)*(jk-1)664 ! WRITE(numout,*) ' FAC2 : ', icsuswi(ji)*(jk-1+icsuind(ji))665 ! WRITE(numout,*) ' FAC3 : ', (1-icsuswi(ji))*(1-snicswi(ji))*jk666 ! ENDIF667 550 END DO 668 551 ENDDO 669 ! IF (jiindex_1d .GT.0) WRITE(numout,*) ' qm0 : ', qm0(jiindex_1d,1:nbot0(jiindex_1d)) / &670 ! rdt_ice671 672 ! ! Heat conservation673 ! WRITE(numout,*) ' zqts_in : ', zqts_in(jiindex_1d) / rdt_ice674 552 675 553 !---------------------------- … … 696 574 ! Snow ice layer heat content 697 575 !----------------------------- 698 699 576 DO ji = kideb, kiut 700 577 ! energy of the flooding seawater … … 704 581 qt_i_in(ji,jl) = qt_i_in(ji,jl) + zqsnic 705 582 706 ! The heat of the flooding seawater is removed from the lead707 ! default line708 ! qldif_1d(ji) = qldif_1d(ji) - zqsnic * a_i_b(ji)709 ! LIM 2 line710 ! qldif_1d(ji) = qldif_1d(ji) + zdrmh * ( 1.0 - alphs * ( rhosn/rhoic ) ) * xlic * ( 1.0 - frld_1d(ji) )711 ! debug line : test a plus (should also be positive for the ice)712 583 qldif_1d(ji) = qldif_1d(ji) + zqsnic * a_i_b(ji) 713 ! qldif_1d(ji) = qldif_1d(ji) + rhosn * cpic *( t_bo_b(ji)-t_s_b(ji,1) ) * &714 ! a_i_b(ji)*dh_snowice(ji)715 !++++++716 ! IF ( ji.EQ.jiindex_1D ) THEN717 ! WRITE(numout,*) ' zqsnic : ', zqsnic / rdt_ice718 ! WRITE(numout,*) ' zqsnow : ', zqsnow(ji) / rdt_ice719 ! ENDIF720 ! !++++++721 584 722 585 ! enthalpy of the newly formed snow-ice layer … … 725 588 qm0(ji,1) = snicswi(ji) * zqsnic + ( 1 - snicswi(ji) ) * qm0(ji,1) 726 589 727 ! !++++++728 ! IF (ji.EQ.jiindex_1D) THEN729 ! WRITE(numout,*) ' snicswi : ', snicswi(ji)730 ! WRITE(numout,*) ' zqsnic : ', zqsnic / rdt_ice731 ! WRITE(numout,*) ' qldif_1d : ', qldif_1d(ji)732 ! WRITE(numout,*) ' dhsnowice :', dh_snowice(ji)733 ! WRITE(numout,*) ' t_bo_b :', t_bo_b(ji)734 ! WRITE(numout,*) ' t_s_b :', t_s_b(ji,1)735 ! WRITE(numout,*) ' limsum :', limsum736 ! WRITE(numout,*) ' a_i_b :', a_i_b(ji)737 ! ENDIF738 ! !++++++739 740 590 ENDDO ! ji 741 742 ! IF (jiindex_1D.GT.1) THEN743 ! WRITE(numout,*) ' Vertical profile : '744 ! WRITE(numout,*) ' ht_s_b : ', ht_s_b(jiindex_1d)745 ! WRITE(numout,*) ' ht_i_b : ', ht_i_b(jiindex_1d)746 ! WRITE(numout,*) ' nbot0 : ', nbot0(jiindex_1d)747 ! WRITE(numout,*) ' qm0 : ', qm0(jiindex_1d,1:nbot0(jiindex_1d)) / &748 ! rdt_ice749 ! WRITE(numout,*) ' zthick0: ', zthick0(jiindex_1d,1:nbot0(jiindex_1d))750 ! WRITE(numout,*) ' zm0 : ', zm0(jiindex_1d,0:nbot0(jiindex_1d))751 ! ENDIF752 591 753 592 DO jk = ntop0, maxnbot0 … … 757 596 * MAX( 0.0 , SIGN(1.0,ht_i_b(ji)-zeps6+zeps) ) & 758 597 * MAX( 0.0 , SIGN( 1. , nbot0(ji) - jk + zeps ) ) 759 ! IF (ji.EQ.jiindex_1D) THEN 760 ! WRITE(numout,*) ' jk : ', jk 761 ! WRITE(numout,*) ' qm0 : ', qm0(ji,jk) 762 ! WRITE(numout,*) ' nbot0:', nbot0(ji) 763 ! WRITE(numout,*) ' fac 1 : ', MAX( 0.0 , SIGN( 1.0, nbot0(ji) & 764 ! - jk + zeps ) ) 765 ! ENDIF 766 END DO 767 ENDDO 768 ! IF (jiindex_1D .GT. 0 ) THEN 769 ! WRITE(numout,*) ' qm0 : ', qm0(jiindex_1d,1:jkmax) / rdt_ice 770 ! WRITE(numout,*) ' zqti_in : ', zqti_in(jiindex_1d) / rdt_ice 771 ! ENDIF 772 ! Heat conservation 773 ! WRITE(numout,*) ' zm0 : ', zm0(jiindex_1d,:) 774 ! WRITE(numout,*) ' zthick0 : ', zthick0(jiindex_1d,:) 775 ! WRITE(numout,*) ' qm0 : ', qm0(jiindex_1d,0:jkmax) / rdt_ice 598 END DO 599 ENDDO 776 600 777 601 !------------- … … 803 627 ENDDO 804 628 805 ! What is the difference with zthick0 ???806 629 !--thicknesses of the layers 807 630 DO layer0 = ntop0, maxnbot0 … … 810 633 END DO 811 634 ENDDO 812 813 814 ! Use less ??815 ! DO layer1 = ntop1, nbot1816 ! DO ji = kideb, kiut817 ! IF (ji.eq.jiindex_1d) THEN818 ! WRITE(numout,*) ' jk : ', layer1819 ! WRITE(numout,*) ' q_i_b : ', q_i_b(ji,layer1)820 ! ENDIF821 ! q_i_b(ji,layer1) = q_i_b(ji,layer1) * (1.0 - MAX(0.0,SIGN(1.0,ht_i_b(ji)-zeps6+zeps)))822 ! IF (ji.eq.jiindex_1d) THEN823 ! WRITE(numout,*) ' ht_i_b ', ht_i_b(ji)824 ! WRITE(numout,*) ' factor2 ', ht_i_b(ji)-zeps6+zeps825 ! WRITE(numout,*) ' factor .... ', (1.0 - &826 ! MAX(0.0,SIGN(1.0,ht_i_b(ji)-zeps6+zeps)))827 ! WRITE(numout,*) ' q_i_b : ', q_i_b(ji,layer1)828 ! ENDIF829 ! END DO830 ! ENDDO831 635 832 636 !------------------------ … … 848 652 ENDDO 849 653 850 851 852 ! WRITE(numout,*) ' s_i : ', s_i_b(jiindex_1d,1:nlay_i)853 854 654 !------------------------- 855 655 ! Heat conservation check … … 861 661 END DO 862 662 END DO 863 ! IF ( jiindex_1d .GT. 0 ) &864 ! WRITE(numout,*) ' zqti_fin : ', zqti_fin(jiindex_1d) / rdt_ice865 663 ! 866 664 DO ji = kideb, kiut … … 882 680 ENDIF 883 681 END DO 884 !885 ! WRITE(numout,*) ' old_ht_i : ', old_ht_i_b(jiindex_1d)886 ! WRITE(numout,*) ' ht_i_b : ', ht_i_b(jiindex_1d)887 ! WRITE(numout,*) ' icsuswi : ', icsuswi(jiindex_1d)888 ! WRITE(numout,*) ' icboswi : ', icboswi(jiindex_1d)889 ! WRITE(numout,*) ' snicswi : ', snicswi(jiindex_1d)890 ! WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(jiindex_1d)891 ! WRITE(numout,*) ' dh_i_surf : ', dh_i_surf(jiindex_1d)892 ! WRITE(numout,*) ' dh_snowice : ', dh_snowice(jiindex_1d)893 ! WRITE(numout,*) ' zqti_fin : ', zqti_fin(jiindex_1d) / rdt_ice894 ! WRITE(numout,*) ' difference : ', ( zqti_fin(jiindex_1d) - &895 ! zqti_in(jiindex_1d) ) / rdt_ice896 682 897 683 !---------------------- … … 904 690 ENDDO !jk 905 691 906 !++++++++++++++++907 692 ! Heat conservation 908 693 zqti_fin(:) = 0.0 … … 912 697 END DO 913 698 END DO 914 ! IF (jiindex_1d .GT. 0 ) & 915 ! WRITE(numout,*) ' zqti_fin : ', zqti_fin(jiindex_1d) / rdt_ice916 ! WRITE(numout,*) ' q_i_b : ', q_i_b(jiindex_1d,1:nlay_i)917 !++++++++++++++++ 918 !------------------------------------------------- 919 ! Update salinity, bottom and snow ice entrapment 920 ! -------------------------------------------------699 700 ! 701 !------------------------------------------------------------------------------| 702 ! 5) Update salinity and recover temperature | 703 !------------------------------------------------------------------------------| 704 ! 705 ! Update salinity (basal entrapment, snow ice formation) 921 706 DO ji = kideb, kiut 922 707 sm_i_b(ji) = sm_i_b(ji) & … … 924 709 END DO !ji 925 710 926 !---------------------927 711 ! Recover temperature 928 !---------------------929 712 DO jk = 1, nlay_i 930 713 … … 944 727 END DO !jk 945 728 946 ! Fin de la routine Enth947 729 END SUBROUTINE lim_thd_ent 948 730 … … 957 739 #endif 958 740 END MODULE limthd_ent 959 -
trunk/NEMO/LIM_SRC_3/limthd_lac.F90
r825 r834 1 1 MODULE limthd_lac 2 2 #if defined key_lim3 3 !!---------------------------------------------------------------------- 4 !! 'key_lim3' LIM3 sea-ice model 5 !!---------------------------------------------------------------------- 3 6 !!====================================================================== 4 7 !! *** MODULE limthd_lac *** … … 20 23 USE iceini 21 24 USE limtab 22 USE limicepoints23 25 USE taumod 24 26 USE blk_oce … … 43 45 44 46 !!---------------------------------------------------------------------- 45 !! LIM 2.0, UCL-LOCEAN-IPSL (2005)47 !! LIM 3.0, UCL-ASTR-LOCEAN-IPSL (2008) 46 48 !! $Header: /home/opalod/NEMOCVSROOT/NEMO/LIM_SRC/limthd_lac.F90,v 1.5 2005/03/27 18:34:42 opalod Exp $ 47 49 !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt … … 79 81 !! 80 82 !! History : 81 !! 1.0 ! 01-04 (T. Fichefet, M. A. Morales Maqueda, H. Goosse)82 !! ! original code83 !! 2.0 ! 02-08 (C. Ethe, G. Madec) F90, mpp84 83 !! 3.0 ! 12-05 (M. Vancoppenolle) Thorough rewrite of the routine 85 84 !! Salinity variations in sea ice, 86 85 !! Multi-layer code 87 !! 3.1 ! 01-06 (M. Vancoppenolle) I ce thickness distribution86 !! 3.1 ! 01-06 (M. Vancoppenolle) ITD 88 87 !! 3.2 ! 04-07 (M. Vancoppenolle) Mass and energy conservation tested 89 88 !!------------------------------------------------------------------------ … … 185 184 186 185 REAL(wp) :: & 187 zde , & ! :increment of energy in category jl 188 zde_old , & ! :dummy variable for energy conservation 189 zde_new , & ! :dummy variable for energy conservation 190 zde_nice , & ! :dummy variable for energy conservation 191 zde_diff ! :dummy variable for energy conservation 186 zde ! :increment of energy in category jl 192 187 193 188 CHARACTER (len = 15) :: fieldid … … 199 194 vt_i_init(:,:) = 0.0 200 195 vt_s_init(:,:) = 0.0 201 IF(lwp) THEN 202 WRITE(numout,*) 'lim_thd_lac : Creating new ice' 203 WRITE(numout,*) '~~~~~~~~~~~' 196 zeps6 = 1.0e-6 197 198 !------------------------------------------------------------------------------! 199 ! 1) Conservation check and changes in each ice category 200 !------------------------------------------------------------------------------! 201 IF ( con_i ) THEN 202 CALL lim_column_sum (jpl, v_i, vt_i_init) 203 CALL lim_column_sum (jpl, v_s, vt_s_init) 204 CALL lim_column_sum_energy (jpl, nlay_i, e_i, et_i_init) 205 CALL lim_column_sum (jpl, e_s(:,:,1,:) , et_s_init) 204 206 ENDIF 205 zeps6 = 1.0e-6206 207 !++++++++++208 WRITE(numout,*) ' Step 0 '209 WRITE(numout,*) ' v_i : ', v_i(jiindex,jjindex,1:jpl)210 WRITE(numout,*) ' a_i : ', a_i(jiindex,jjindex,1:jpl)211 DO jk = 1, nlay_i212 WRITE(numout,*) ' e_i : ', jk, e_i(jiindex,jjindex,jk,1:jpl)213 END DO214 !++++++++++215 216 !------------------------------------------------------------------------------!217 ! 1) Conservation check and changes in each ice category218 !------------------------------------------------------------------------------!219 CALL lim_column_sum (jpl, v_i, vt_i_init)220 CALL lim_column_sum (jpl, v_s, vt_s_init)221 CALL lim_column_sum_energy (jpl, nlay_i, e_i, et_i_init)222 CALL lim_column_sum (jpl, e_s(:,:,1,:) , et_s_init)223 207 224 208 !------------------------------------------------------------------------------| … … 497 481 ! 5) Compute thickness, salinity, enthalpy, age, area and volume of new ice 498 482 !------------------------------------------------------------------------------! 483 499 484 !---------------------- 500 485 ! Thickness of new ice … … 511 496 zjj = ( npac(ji) - 1 ) / jpi + 1 512 497 WRITE(numout,*) ' collection thickness <= 0 ', zh_newice(ji), ji, zji, zjj 513 ! WRITE(numout,*) ' LATITUDE ', gphit(zji,zjj), ' LONGITUDE ', &514 ! glamt(zji,zjj)515 ! WRITE(numout,*) ' zh_newice ', zh_newice(ji)516 ! WRITE(numout,*) ' ji ', ji517 ! WRITE(numout,*) ' a_i ', a_i(zji,zjj,1:jpl)518 ! WRITE(numout,*) ' v_i ', v_i(zji,zjj,1:jpl)519 498 ENDIF 520 499 END DO … … 554 533 / ( t_bo_b(ji) - rtt ) ) & 555 534 - rcp * ( ztmelts-rtt ) ) 556 !+++++++++++557 IF ( ji .EQ. jiindex_1D ) THEN558 WRITE(numout,*) ' ze_newice : ', ze_newice(jiindex_1d)559 ENDIF560 !+++++++++++561 535 ze_newice(ji) = MAX( ze_newice(ji) , 0.0 ) + & 562 536 MAX( 0.0 , SIGN( 1.0 , - ze_newice(ji) ) ) & 563 537 * rhoic * lfus 564 !+++++++++++565 IF ( ji .EQ. jiindex_1D ) THEN566 WRITE(numout,*) ' ze_newice : ', ze_newice(jiindex_1d)567 ENDIF568 !+++++++++++569 538 END DO ! ji 570 !+++++++++++571 IF ( jiindex_1d .GT. 0) WRITE(numout,*) ' ze_newice : ', ze_newice(jiindex_1d)572 !+++++++++++573 574 539 !---------------- 575 540 ! Age of new ice … … 585 550 zqbgow(ji) = qldif_1d(ji) - qcmif_1d(ji) !<0 586 551 END DO ! ji 587 !+++++++++++588 IF (jiindex_1d .GT. 0 ) THEN589 WRITE(numout,*) ' qldif : ', qldif_1d(jiindex_1d)590 WRITE(numout,*) ' qcmif : ', qcmif_1d(jiindex_1d)591 WRITE(numout,*) ' zqbgow : ', zqbgow(jiindex_1d)592 ENDIF593 !+++++++++++594 552 595 553 !------------------- … … 605 563 zv_newice(ji) = ( 1.0 - zfrazb ) * zv_newice(ji) 606 564 END DO 607 !+++++++++++608 IF (jiindex_1d .GT. 0) THEN609 WRITE(numout,*) ' zv_newice : ', zv_newice(jiindex_1d)610 ENDIF611 !+++++++++++612 565 613 566 !------------------------------------ … … 663 616 za_newice(ji) = za_newice(ji) - zda_res(ji) 664 617 zv_newice(ji) = zv_newice(ji) - zdv_res(ji) 665 IF ( ji .EQ. jiindex_1D ) THEN666 WRITE(numout,*) ' zv_newice : ', zv_newice(ji)667 WRITE(numout,*) ' zdv_res : ', zdv_res (ji)668 ENDIF669 618 ELSE 670 619 zda_res(ji) = 0.0 … … 678 627 zat_i_ac(:) = 0.0 679 628 680 WRITE(numout,*) ' zv_i_ac : ', zv_i_ac(jiindex, 1:jpl)681 629 DO jl = 1, jpl 682 630 DO ji = 1, nbpac … … 691 639 END DO ! ji 692 640 END DO ! jl 693 WRITE(numout,*) ' zv_i_ac : ', zv_i_ac(jiindex, 1:jpl)694 641 695 ! !++++++++++++++++696 DO ji = 1, nbpac697 ! !+++++698 IF (zat_i_ac(ji).gt.1.0) THEN699 zji = MOD( npac(ji) - 1, jpi ) + 1700 zjj = ( npac(ji) - 1 ) / jpi + 1701 WRITE(numout,*) ' *** ERROR MESSAGE *** '702 WRITE(numout,*) ' something MUST be wrong '703 WRITE(numout,*) ' at_i : ', zat_i_ac(ji)704 WRITE(numout,*) ' is wrong '705 WRITE(numout,*) ' point : ', zji, zjj706 ENDIF707 END DO ! ji708 !++++++++++++++++709 710 642 !---------------------------------- 711 643 ! Heat content - lateral accretion … … 849 781 zindb = 1.0 - MAX( 0.0 , SIGN( 1.0 , - & 850 782 za_i_ac(ji,jl) ) ) ! 0 if no ice and 1 if yes 851 !zo_i_ac(ji,jl) = zv_old(ji,jl) * zo_i_ac(ji,jl) / &852 !MAX(zv_i_ac(ji,jl),zeps) * zindb853 783 zoa_i_ac(ji,jl) = za_old(ji,jl) * zoa_i_ac(ji,jl) / & 854 784 MAX( za_i_ac(ji,jl) , zeps ) * zindb … … 940 870 END DO 941 871 942 !++++943 WRITE(numout,*) 'lim_thd_lac : Salt flux diagnostic '944 WRITE(numout,*) '~~~~~~~~~~~~'945 WRITE(numout,*) ' *** Salt fluxes at bottom interface ***'946 WRITE(numout,*) ' fseqv : ', fseqv(jiindex,jjindex)947 WRITE(numout,*)948 !++++949 950 872 !------------------------------------------------------------------------------| 951 873 ! 10) Conservation check and changes in each ice category 952 874 !------------------------------------------------------------------------------| 953 875 876 IF ( con_i ) THEN 954 877 CALL lim_column_sum (jpl, v_i, vt_i_final) 955 878 fieldid = 'v_i, limthd_lac' … … 972 895 WRITE(numout,*) ' et_i_init : ', et_i_init(jiindex,jjindex) 973 896 WRITE(numout,*) ' et_i_final: ', et_i_final(jiindex,jjindex) 897 898 ENDIF 974 899 975 900 END SUBROUTINE lim_thd_lac -
trunk/NEMO/LIM_SRC_3/limthd_sal.F90
r825 r834 1 1 MODULE limthd_sal 2 2 #if defined key_lim3 3 !!---------------------------------------------------------------------- 4 !! 'key_lim3' LIM3 sea-ice model 5 !!---------------------------------------------------------------------- 3 6 !!====================================================================== 4 7 !! *** MODULE limthd_sal *** … … 20 23 USE limvar 21 24 USE par_ice 22 USE limicepoints23 25 24 26 IMPLICIT NONE … … 60 62 !! 61 63 !! ** History : 64 !! 65 !! "Je ne suis pour l'instant qu'a 80% de ma condition, mais c'est 66 !! les 30% qui restent qui seront les plus difficiles" 67 !! E. Mpenza 68 !! 62 69 !!------------------------------------------------------------------- 63 70 !! History : … … 65 72 !! 3.0 ! 05-12 Routine rewritten for the 3-D version 66 73 !!--------------------------------------------------------------------- 74 !! 67 75 !! * Local variables 68 76 INTEGER, INTENT(in) :: & … … 74 82 75 83 REAL(wp) :: & 76 ! ds_i_gd, & !: salt loss by gravity drainage77 ! ds_i_seg, & !: salt gain by segregation of seawater at the bottom78 ! ds_i_flush, & !: summer salt loss by flushing79 zargtemp, & !: used to compute Scharzacher profile80 84 zsold, & !: old salinity 81 ! ds_i_snice, & !: increase of salinity due to snow ice82 85 zeps=1.0e-06 , & !: very small 83 86 iflush , & !: flushing (1) or not (0) … … 87 90 i_ice_switch , & !: ice thickness above a certain treshold or not 88 91 ztmelts , & !: freezing point of sea ice 89 zs_snowice , & !: salinity of forming snow ice90 92 zaaa , & !: dummy factor 91 93 zbbb , & !: dummy factor … … 95 97 REAL(wp), DIMENSION(jpij) :: & 96 98 ze_init , & !initial total enthalpy 97 ze_end , & !final total enthalpy98 99 zhiold , & 99 100 zsiold … … 157 158 END DO ! ji 158 159 END DO ! jk 159 160 ! IF (jiindex_1d .GT. 0 ) THEN161 ! WRITE(numout,*) ' s_i_b : ', s_i_b(jiindex_1d,1:nlay_i)162 ! WRITE(numout,*) ' t_i_b : ', t_i_b(jiindex_1d,1:nlay_i)163 ! WRITE(numout,*) ' q_i_b : ', q_i_b(jiindex_1d,1:nlay_i)164 ! ENDIF165 160 166 161 DO ji = kideb, kiut … … 215 210 CALL lim_var_salprof1d(kideb,kiut) 216 211 217 !--------------------------------------218 ! Energy of melting for each ice layer219 !--------------------------------------220 ! q=q(S,T) J.m-3221 ! DO jk = 1, nlay_i222 ! DO ji = kideb, kiut223 ! ! New heat content after desalination224 ! ztmelts = - tmut * s_i_b(ji,jk) + rtt ! Tfi in K225 ! q_i_b(ji,jk) = rhoic * &226 ! ( cpic * ( ztmelts-t_i_b(ji,jk) ) &227 ! + lfus * ( 1.0 - (ztmelts-rtt) / &228 ! MIN( ( t_i_b(ji,jk) - rtt ) , -zeps ) ) &229 ! - rcp * ( ztmelts-rtt ) )230 ! END DO ! ji231 ! END DO232 233 !--------------------234 ! Total heat content235 !--------------------236 ! ze_end(:) = 0.0237 ! DO jk = 1, nlay_i238 ! DO ji = kideb, kiut239 ! ze_end(ji) = ze_end(ji) + q_i_b(ji,jk) * ht_i_b(ji) / nlay_i240 ! END DO241 ! END DO242 243 ! IF (jiindex_1d .GT. 0 ) THEN244 ! WRITE(numout,*) ' ze_init : ', ze_init(jiindex_1d) / rdt_ice245 ! WRITE(numout,*) ' ze_end : ', ze_end (jiindex_1d) / rdt_ice246 ! WRITE(numout,*) ' s_i_b : ', s_i_b(jiindex_1d,1:nlay_i)247 ! WRITE(numout,*) ' t_i_b : ', t_i_b(jiindex_1d,1:nlay_i)248 ! WRITE(numout,*) ' q_i_b : ', q_i_b(jiindex_1d,1:nlay_i)249 ! ENDIF250 251 212 !---------------------------- 252 213 ! Heat flux - brine drainage 253 214 !---------------------------- 254 255 ! WRITE(numout,*) ' jiindex, jjindex : ', jiindex, jjindex256 ! WRITE(numout,*) ' jiindex_1d : ', jiindex_1d257 ! WRITE(numout,*) ' kideb, kiut: ', kideb, kiut258 215 259 216 DO ji = kideb, kiut … … 265 222 iaccrbo = MAX( 0.0 , SIGN ( 1.0 , dh_i_bott(ji) ) ) 266 223 267 ! IF ( ( iflush .EQ. 1 ) .OR. ( iaccrbo .EQ. 1 ) ) THEN268 ! fhbri_1d(ji) = fhbri_1d(ji) + ( ze_end(ji) - ze_init(ji) ) * &269 ! a_i_b(ji) / at_i_b(ji) / rdt_ice270 ! ze_end GT ze_init... ?271 ! ENDIF272 224 fhbri_1d(ji) = 0.0 273 274 !++++++ 275 ! zji = MOD( npb(ji) - 1, jpi ) + 1 276 ! zjj = ( npb(ji) - 1 ) / jpi + 1 277 ! fhbricat(zji,zjj,jl) = ( ze_end(ji) - ze_init(ji) ) / & 278 ! rdt_ice 279 !++++++ 280 ! ftotal_fin(ji) = ftotal_fin(ji) + fhbricat(zji,zjj,jl) 281 ! IF ( ji .EQ. jiindex_1d) WRITE(numout,*) ' fhbri : ', & 282 ! fhbricat(zji,zjj,jl) * rdt_ice 283 ! IF ( ji .EQ. jiindex_1d ) & 284 ! WRITE(numout,*) ' ftotal_fin : ', ftotal_fin(jiindex_1d) 285 ! !+++++ 286 287 END DO ! ji 288 289 ! IF (jiindex_1d .GT. 0 ) THEN 290 ! WRITE(numout,*) ' Category : ', jl 291 ! WRITE(numout,*) ' jiindex_1d : ', jiindex_1d 292 ! WRITE(numout,*) ' fhbricat : ', fhbricat(jiindex,jjindex,jl) 293 ! WRITE(numout,*) ' fhbri_1d : ', fhbri_1d(jiindex_1d) 294 ! WRITE(numout,*) ' ze_end : ', ze_end(jiindex_1d) 295 ! WRITE(numout,*) ' ze_init : ', ze_init(jiindex_1d) 296 ! WRITE(numout,*) ' a_i_b : ', a_i_b(jiindex_1d) 297 ! WRITE(numout,*) ' at_i_b : ', at_i_b(jiindex_1d) 298 ! WRITE(numout,*) ' rdt_ice : ', rdt_ice 299 ! ENDIF 225 END DO ! ji 300 226 301 227 !---------------------------- … … 377 303 ENDIF 378 304 379 !!--multiyear sea-ice salinity380 !!-- ca va pas alors on vire bordel de djeuille381 ! if (ht.gt.360) then382 ! do jk = 1, nlay_i383 ! s_i_b(ji,jk) = 1.58 + 0.18*ht_i_b(ji)384 ! end do385 ! endif386 305 DO jk = 1, nlay_i 387 306 s_i_b(ji,jk) = sm_i_b(ji) … … 415 334 END DO ! ji 416 335 ENDIF 417 418 !+++++419 ! IF ( (zji.EQ.jiindex) .AND. (zjj.EQ.jjindex) ) THEN420 ! WRITE(numout,*) ' *** FS_EQV *** '421 ! WRITE(numout,*) ' limthd_sal '422 ! WRITE(numout,*) ' fseqv ', fseqv_1d(ji)423 ! WRITE(numout,*) ' sss_io ', sss_io(zji,zjj)424 ! WRITE(numout,*) ' s_i_new ', s_i_new(ji)425 ! WRITE(numout,*) ' dh_i_bott ', dh_i_bott(ji)426 ! WRITE(numout,*) ' a_i_b ', a_i_b(ji)427 ! WRITE(numout,*) ' *** FS_BRI *** '428 ! WRITE(numout,*) ' fsbri ', fsbri_1d(ji)429 ! WRITE(numout,*) ' ht_i ', ht_i_b(ji)430 ! WRITE(numout,*) ' dsm_i_gd ', dsm_i_gd_1d(ji)431 ! WRITE(numout,*) ' dsm_i_fl ', dsm_i_fl_1d(ji)432 ! ENDIF433 434 !+++++435 336 436 337 !-- End of salinity computations -
trunk/NEMO/LIM_SRC_3/limtrp.F90
r825 r834 6 6 #if defined key_lim3 7 7 !!---------------------------------------------------------------------- 8 !! 'key_lim3' : LIMsea-ice model8 !! 'key_lim3' LIM3 sea-ice model 9 9 !!---------------------------------------------------------------------- 10 10 !! lim_trp : advection/diffusion process of sea ice … … 26 26 USE lib_mpp 27 27 USE par_ice 28 USE limicepoints29 28 30 29 IMPLICIT NONE … … 50 49 # include "vectopt_loop_substitute.h90" 51 50 !!---------------------------------------------------------------------- 52 !! LIM 2.0, UCL-LOCEAN-IPSL (2005)51 !! LIM 3.0, UCL-ASTR-LOCEAN-IPSL (2008) 53 52 !! $Header: /home/opalod/NEMOCVSROOT/NEMO/LIM_SRC/limtrp.F90,v 1.5 2005/03/27 18:34:42 opalod Exp $ 54 53 !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt … … 75 74 !!--------------------------------------------------------------------- 76 75 !! * Local Variables 77 INTEGER :: ji, jj, jk, jl, index,layer, & ! dummy loop indices76 INTEGER :: ji, jj, jk, jl, layer, & ! dummy loop indices 78 77 initad ! number of sub-timestep for the advection 79 78 INTEGER :: ji_maxu, ji_maxv, jj_maxu, jj_maxv … … 81 80 REAL(wp) :: & 82 81 zindb , & 83 zacrith, & 84 zindsn , & 85 zindic , & 86 zusvosn, & 87 zusvoic, & 88 zignm , & 89 zindhe , & 90 zvbord , & 91 zcfl , & 92 zusnit , & 93 zrtt, ztsn, ztic1, ztic2, zsal, zage, & 94 zq, zbigval, ze1, ze2, ze, & 82 zindsn , & 83 zindic , & 84 zusvosn, & 85 zusvoic, & 86 zvbord , & 87 zcfl , & 88 zusnit , & 89 zrtt, zsal, zage, & 90 zbigval, ze, & 95 91 zmaxu, zmaxv 96 92 … … 101 97 REAL(wp), DIMENSION(jpi,jpj,jpl):: & ! temporary workspace 102 98 zs0ice, zs0sn, zs0a , & 103 zs0c0 , zs0c1 , zs0c2 ,&99 zs0c0 , & 104 100 zs0sm , zs0oi 105 101 … … 554 550 END DO ! jl 555 551 556 !-----------------------------------------------------------------------------!557 ! 6) Correct age558 !-----------------------------------------------------------------------------!559 560 ! DO jl = 1, jpl561 ! DO jj = 1, jpj562 ! DO ji = 1, jpi563 ! IF (old_a_i(ji,jj,jl).eq.0.00) THEN564 ! o_i(ji,jj,jl) = MAX( MIN( rdt_ice*float(numit)/86400.0, o_i(ji,jj,jl) ), 0.0 )565 ! ENDIF566 ! END DO567 ! END DO568 ! END DO569 570 571 552 ENDIF 572 553 -
trunk/NEMO/LIM_SRC_3/limupdate.F90
r825 r834 9 9 !! from advection and ice thermodynamics 10 10 !!====================================================================== 11 #if defined key_lim3 11 12 !!---------------------------------------------------------------------- 12 !! 'key_lim3' : LIMsea-ice model13 !! 'key_lim3' LIM3 sea-ice model 13 14 !!---------------------------------------------------------------------- 14 15 !! lim_update : computes update of sea-ice global variables … … 42 43 USE thd_ice ! LIM thermodynamic sea-ice variables 43 44 USE par_ice 44 USE limicepoints45 45 USE limitd_th 46 46 USE limvar … … 80 80 INTEGER :: & 81 81 ji, jj, & ! geographical indices 82 jk, jl, jm, & ! layer, category and type indices 83 index, num_ex_cat, zji, zjj ! dummy loop indices 82 jk, jl, jm ! layer, category and type indices 84 83 INTEGER :: & 85 84 jbnd1, jbnd2 … … 96 95 rzero = 0.e0 , & 97 96 rone = 1.e0 , & 98 z_slope_s , & ! slope of the salinity profile99 97 zhimax ! maximum thickness tolerated for advection of 100 98 ! in an ice-free cell 101 99 REAL(wp) :: & ! dummy switches and arguments 102 zind he, zconc, zindb, zindsn, zindic, zacrith, zignm, &103 zrtt, z tsn, zindg, zthsnice, za, zh, zdvres, zviold,ztic,&104 zbigvalue, zvsold, z vold, z_dv_ex, z_da_ex, zamax, zdv,&105 z_prescr_hi, zat_i_old, z_a_i, zalpha, zs_inf, zsmin, zsmax,&100 zindb, zindsn, zindic, zacrith, & 101 zrtt, zindg, zh, zdvres, zviold, & 102 zbigvalue, zvsold, z_da_ex, zamax, & 103 z_prescr_hi, zat_i_old, & 106 104 ztmelts, ze_s 107 108 REAL(wp), DIMENSION(nlay_i) :: &109 z_de, zs_zero110 105 111 106 REAL(wp), DIMENSION(jpl) :: z_da_i, z_dv_i … … 113 108 LOGICAL, DIMENSION(jpi,jpj,jpl) :: & 114 109 internal_melt 110 115 111 INTEGER :: & 116 112 ind_im, layer ! indices for internal melt … … 734 730 ( 1.0 - zindsn ) * 0.0 735 731 736 ! ztsn = MIN( MAX( zrtt, t_s(ji,jj,1,jl) ) , rtt )737 ! t_s(ji,jj,1,jl) = zindsn*ztsn + (1.0 - zindsn)*t_bo(ji,jj)738 739 732 END DO ! ji 740 733 END DO ! jj … … 864 857 865 858 ! 1) Count the number of existing categories 866 ! num_ex_cat = 0867 859 DO jl = 1, jpl 868 860 zindb = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi03 ) ) … … 1005 997 1006 998 END SUBROUTINE lim_update 999 #else 1000 !!---------------------------------------------------------------------- 1001 !! Default option Empty Module No sea-ice model 1002 !!---------------------------------------------------------------------- 1003 CONTAINS 1004 SUBROUTINE lim_update ! Empty routine 1005 END SUBROUTINE lim_update 1006 1007 #endif 1008 1007 1009 END MODULE limupdate -
trunk/NEMO/LIM_SRC_3/limvar.F90
r825 r834 1 1 MODULE limvar 2 2 #if defined key_lim3 3 !!---------------------------------------------------------------------- 4 !! 'key_lim3' LIM3 sea-ice model 5 !!---------------------------------------------------------------------- 3 6 !!====================================================================== 4 7 !! *** MODULE limvar *** … … 38 41 USE ice_oce ! ice variables 39 42 USE thd_ice 40 USE limicepoints41 43 USE in_out_manager 42 44 USE ice … … 62 64 63 65 !!---------------------------------------------------------------------- 64 !! LIM -@ 4.0, UCL-ASTR (2006)66 !! LIM3.0, UCL-ASTR-LOCEAN-IPSL (2008) 65 67 !! (c) UCL-ASTR and Martin Vancoppenolle 66 68 !!---------------------------------------------------------------------- … … 91 93 !! ** History : 92 94 !! (01-2006) Martin Vancoppenolle, UCL-ASTR 93 !! Take it easy man94 !! Life is just a simple game, between ups / and downs \ :@)95 95 !! 96 96 !! note : you could add an argument when you need only at_i, vt_i … … 103 103 jj, & ! spatial dummy loop index 104 104 jk, & ! vertical layering dummy loop index 105 jl, & ! ice category dummy loop index 106 jm, & ! ice types dummy loop index 107 index 105 jl ! ice category dummy loop index 108 106 109 107 REAL :: zeps, epsi16, zinda, epsi06 … … 116 114 epsi16 = 1.0e-16 117 115 epsi06 = 1.0e-6 118 119 IF(lwp) THEN120 WRITE(numout,*)121 WRITE(numout,*) 'lim_var_agg : Aggregate ice variables '122 WRITE(numout,*) '~~~~~~~~~~~'123 ENDIF124 116 125 117 !------------------ … … 228 220 jj, & ! spatial dummy loop index 229 221 jk, & ! vertical layering dummy loop index 230 jl, & ! ice category dummy loop index 231 jm, & ! ice types dummy loop index 232 index 222 jl ! ice category dummy loop index 233 223 234 224 REAL :: zq_i, zaaa, zbbb, zccc, zdiscrim, & … … 242 232 !!-- End of declarations 243 233 !!------------------------------------------------------------------------------ 244 ! IF (lwp) THEN245 ! WRITE(numout,*)246 ! WRITE(numout,*) 'lim_var_glo2eqv : transform global variables into equivalent ice variables '247 ! WRITE(numout,*) '~~~~~~~~~~~~~~~'248 ! ENDIF249 234 250 235 !------------------------------------------------------- … … 377 362 !! 378 363 !!------------------------------------------------------------------ 379 !! * Arguments380 381 !! * Local variables382 INTEGER :: ji, & ! spatial dummy loop index383 jj, & ! spatial dummy loop index384 jk, & ! vertical layering dummy loop index385 jl, & ! ice category dummy loop index386 jm, & ! ice types dummy loop index387 index388 389 REAL :: zq_i, zaaa, zbbb, zccc, zdiscrim, &390 ztmelts391 392 REAL :: zeps393 394 !Ice specific heat, massive latent heat395 zeps = 1.0e-10396 397 !!-- End of declarations398 !!------------------------------------------------------------------------------399 400 ! IF (lwp) THEN401 ! WRITE(numout,*)402 ! WRITE(numout,*) 'lim_var_eqv2glo : transform equivalent variables into global ice variables '403 ! WRITE(numout,*) '~~~~~~~~~~~~~~~'404 ! ENDIF405 364 406 365 v_i(:,:,:) = ht_i(:,:,:) * a_i(:,:,:) … … 408 367 smv_i(:,:,:) = sm_i(:,:,:) * v_i(:,:,:) 409 368 oa_i (:,:,:) = o_i (:,:,:) * a_i(:,:,:) 410 411 ! DO jl = 1, jpl412 ! DO jk = 1, nlay_i413 ! DO jj = 1, jpj414 ! DO ji = 1, jpi415 ! ! heat content per unit volume416 ! e_i(ji,jj,jk,jl) = zidto * rhoic * &417 ! ( cpic * ( ztmelts - t_i(ji,jj,jk,jl) ) &418 ! + lfus * ( 1.0 - (ztmelts-rtt) / MIN((t_i(ji,jj,jk,jl)-rtt),-zeps) ) &419 ! - rcp * ( ztmelts - rtt ) &420 ! )421 422 ! ! Correct dimensions to avoid big values423 ! e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / unit_fac424 425 ! ! Mutliply by ice volume, and divide by number of layers to get heat content in 10^9 J426 ! e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * &427 ! area(ji,jj) * a_i(ji,jj,jl) * ht_i(ji,jj,jl) / &428 ! FLOAT(nlay_i)429 430 ! END DO431 ! END DO432 ! END DO433 ! END DO434 369 435 370 END SUBROUTINE lim_var_eqv2glo … … 787 722 788 723 #endif 789 724 END MODULE limvar -
trunk/NEMO/LIM_SRC_3/limwri.F90
r825 r834 1 1 MODULE limwri 2 #if defined key_lim3 3 !!---------------------------------------------------------------------- 4 !! 'key_lim3' LIM3 sea-ice model 5 !!---------------------------------------------------------------------- 2 6 !!====================================================================== 3 7 !! *** MODULE limwri *** … … 8 12 !! $Header: /home/opalod/NEMOCVSROOT/NEMO/LIM_SRC/limwri.F90,v 1.4 2005/03/27 18:34:42 opalod Exp $ 9 13 !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 10 !!----------------------------------------------------------------------11 #if defined key_lim312 !!----------------------------------------------------------------------13 !! 'key_lim3' LIM sea-ice model14 14 !!---------------------------------------------------------------------- 15 15 !! lim_wri : write of the diagnostics variables in ouput file … … 29 29 USE iceini 30 30 USE lbclnk 31 USE limicepoints32 31 USE par_ice 33 32 USE limvar … … 98 97 99 98 REAL(wp), DIMENSION(jpi,jpj) :: & 100 zfield , zfielda99 zfield 101 100 102 101 REAL(wp), DIMENSION(jpi,jpj,jpl) :: & … … 107 106 CHARACTER(len = 40) :: & 108 107 clhstnam, clop, & 109 clhstnama , clopa108 clhstnama 110 109 111 110 INTEGER , SAVE :: & 112 111 nice, nhorid, ndim, niter, ndepid 113 112 INTEGER , SAVE :: & 114 nicea, nhorida, n dima, nitera, ndepida, ndimitd113 nicea, nhorida, nitera, ndimitd 115 114 INTEGER , DIMENSION( jpij ) , SAVE :: & 116 115 ndex51 -
trunk/NEMO/LIM_SRC_3/thd_ice.F90
r825 r834 17 17 PRIVATE 18 18 19 !!--------------------------- 19 20 !! * Share Module variables 21 !!--------------------------- 20 22 REAL(wp) , PUBLIC :: & !!! ** ice-thermo namelist (namicethd) ** 21 hmelt = -0.15 , & !: maximum melting at the bottom 22 hicmin = 0.2 , & !: ice th. corr. to max. ener. in brine pocket 23 hmelt = -0.15 , & !: maximum melting at the bottom; active only for 24 !: one category 25 hicmin = 0.2 , & !: (REMOVE) 23 26 hiclim = 0.05 , & !: minimum ice thickness 24 27 amax = 0.999 , & !: maximum lead fraction 25 sbeta = 1.0 , & !: numerical scheme for diffusion in ice 26 parlat = 0.0 , & !: percent. of energy used for lateral ablation27 hakspl = 0.5 , & !: slope of distr. for Hakkinen-Mellro's lat. melt28 hibspl = 0.5 , & !: slope of distribution for Hibler's lat. melt29 exld = 2.0 , & !: exponent for leads-closure rate30 hakdif = 1.0 , & !: coefficient for diffusions of ice and snow31 thth = 0.2 , & !: thick. for comp. of eq. thermal conduct28 sbeta = 1.0 , & !: numerical scheme for diffusion in ice (REMOVE) 29 parlat = 0.0 , & !: (REMOVE) 30 hakspl = 0.5 , & !: (REMOVE) 31 hibspl = 0.5 , & !: (REMOVE) 32 exld = 2.0 , & !: (REMOVE) 33 hakdif = 1.0 , & !: (REMOVE) 34 thth = 0.2 , & !: (REMOVE) 32 35 hnzst = 0.1 , & !: thick. of the surf. layer in temp. comp. 33 36 parsub = 1.0 , & !: switch for snow sublimation or not … … 41 44 hiccrit = (/0.3,0.3/) !: ice th. for lateral accretion in the NH (SH) (m) 42 45 43 REAL(wp) , PUBLIC :: & !: 44 uscomi !: inverse of minimum lead fraction 46 !!----------------------------- 47 !! * Share 1D Module variables 48 !!----------------------------- 49 !: In ice thermodynamics, to spare memory, the vectors are folded 50 !: from 1D to 2D vectors. The following variables, with ending _1d (or _b) 51 !: are the variables corresponding to 2d vectors 45 52 46 53 INTEGER , PUBLIC, DIMENSION(jpij) :: & !: 47 54 npb , & !: number of points where computations has to be done 48 npac !: correspondance between the points 55 npac !: correspondance between the points (lateral accretion) 49 56 50 57 REAL(wp), PUBLIC, DIMENSION(jpij) :: & !: … … 123 130 q_s_b !: Snow enthalpy per unit volume 124 131 125 ! goes to trash, just for conservation checks 126 REAL(wp), PUBLIC, DIMENSION(jpij,jpl) :: & !: ! goes to trash 132 ! Clean the following ... 133 ! These variables are coded for conservation checks 134 REAL(wp), PUBLIC, DIMENSION(jpij,jpl) :: & ! 127 135 qt_i_in , & !: ice energy summed over categories (initial) 128 136 qt_i_fin , & !: ice energy summed over categories (final)
Note: See TracChangeset
for help on using the changeset viewer.