Changeset 5048 for branches/2015/dev_r5044_CNRS_LIM3CLEAN
- Timestamp:
- 2015-02-02T11:28:50+01:00 (9 years ago)
- Location:
- branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM
- Files:
-
- 13 edited
-
CONFIG/SHARED/field_def.xml (modified) (1 diff)
-
CONFIG/SHARED/namelist_ice_lim3_ref (modified) (2 diffs)
-
NEMO/LIM_SRC_3/ice.F90 (modified) (5 diffs)
-
NEMO/LIM_SRC_3/iceini.F90 (modified) (2 diffs)
-
NEMO/LIM_SRC_3/limitd_th.F90 (modified) (1 diff)
-
NEMO/LIM_SRC_3/limthd.F90 (modified) (6 diffs)
-
NEMO/LIM_SRC_3/limthd_dh.F90 (modified) (1 diff)
-
NEMO/LIM_SRC_3/limthd_dif.F90 (modified) (30 diffs)
-
NEMO/LIM_SRC_3/limtrp.F90 (modified) (1 diff)
-
NEMO/LIM_SRC_3/limwri.F90 (modified) (2 diffs)
-
NEMO/LIM_SRC_3/par_ice.F90 (modified) (1 diff)
-
NEMO/LIM_SRC_3/thd_ice.F90 (modified) (4 diffs)
-
NEMO/OPA_SRC/SBC/sbcice_lim.F90 (modified) (7 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/CONFIG/SHARED/field_def.xml
r4996 r5048 297 297 <field id="vfxsub" long_name="snw sublimation" unit="m/day" /> 298 298 <field id="vfxspr" long_name="snw precipitation on ice" unit="m/day" /> 299 300 <field id="afxtot" long_name="area tendency (total)" unit="s-1" /> 301 <field id="afxdyn" long_name="area tendency (dynamics)" unit="s-1" /> 302 <field id="afxthd" long_name="area tendency (thermo)" unit="s-1" /> 299 303 300 304 <field id="hfxsum" long_name="heat fluxes causing surface ice melt" unit="W/m2" /> -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/CONFIG/SHARED/namelist_ice_lim3_ref
r5046 r5048 2 2 !! NEMO/LIM3 : 1 - dynamics/advection/thermo (namicerun) 3 3 !! namelists 2 - ice intialisation (namiceini) 4 !! 3 - ice dynamic 4 !! 3 - ice dynamics (namicedyn) 5 5 !! 4 - ice advection (namicetrp) 6 !! 5 - thermodynamic 6 !! 5 - thermodynamics (namicethd) 7 7 !! 6 - ice salinity (namicesal) 8 8 !! 7 - mechanical redistribution of ice (namiceitdme) … … 17 17 cn_icerst_out = "restart_ice" ! suffix of ice restart name (output) 18 18 ln_limdyn = .true. ! ice dynamics (T) or thermodynamics only (F) 19 amax = 0.999 ! maximum ice concentration20 ln_nicep = .false. ! Ice points output for debug (yes or no)21 ln_limdiahsb = .false. 19 amax = 0.999 ! maximum tolerated ice concentration 20 ln_nicep = .false. ! ice points output for debug (yes or no) 21 ln_limdiahsb = .false. ! check the heat and salt budgets (T) or not (F) 22 22 ln_limdiaout = .true. ! output the heat and salt budgets (T) or not (F) 23 23 / 24 24 !----------------------------------------------------------------------- 25 &namiceini ! ice initiali sation25 &namiceini ! ice initialization 26 26 !----------------------------------------------------------------------- 27 27 ln_limini = .false. ! activate ice initialization (T) or not (F) 28 28 thres_sst = 0.0 ! threshold water temperature for initial sea ice 29 hts_ini_n = 0.3 ! initial snow thickness in the north30 hts_ini_s = 0.3 ! " " south31 hti_ini_n = 3.0 ! initial ice thickness in the north32 hti_ini_s = 1.0 ! " " south33 ati_ini_n = 0.9 ! initial ice concentration in the north34 ati_ini_s = 0.9 ! " " south35 smi_ini_n = 6.3 ! initial ice salinity in the north36 smi_ini_s = 6.3 ! " " south37 tmi_ini_n = 270. ! initial ice/snw temp in the north38 tmi_ini_s = 270. ! initial ice/snw temp in the south29 hts_ini_n = 0.3 ! initial snow thickness (North) 30 hts_ini_s = 0.3 ! " " (South) 31 hti_ini_n = 3.0 ! initial ice thickness (North) 32 hti_ini_s = 1.0 ! " " (South) 33 ati_ini_n = 0.9 ! initial ice concentration (North) 34 ati_ini_s = 0.9 ! " " (South) 35 smi_ini_n = 6.3 ! initial ice salinity (North) 36 smi_ini_s = 6.3 ! " " (South) 37 tmi_ini_n = 270. ! initial ice/snw temp (North) 38 tmi_ini_s = 270. ! initial ice/snw temp (South) 39 39 / 40 40 !----------------------------------------------------------------------- 41 &namicedyn ! ice dynamic 41 &namicedyn ! ice dynamics 42 42 !----------------------------------------------------------------------- 43 cw = 5.0e-03 ! drag coefficient for oceanic stress44 pstar = 2.0e+04 ! 1st bulk-rheology parameter 43 cw = 5.0e-03 ! ice-ocean drag coefficient 44 pstar = 2.0e+04 ! 1st bulk-rheology parameter (N/m) 45 45 c_rhg = 20.0 ! 2nd bulk-rhelogy parameter 46 creepl = 1.0e-12 ! creep limit 46 creepl = 1.0e-12 ! creep limit (s-1) 47 47 ecc = 2.0 ! eccentricity of the elliptical yield curve 48 ahi0 = 350.e0 ! horizontal eddy diffusivity coefficient for sea-ice [m2/s] !depends on resolution49 nevp = 120 ! number of iterations for subcycling in EVP50 relast = 0.333 ! ratio of elastic timescale over ice time step (1/3 if nevp=120 ;1/9 if nevp=300)51 hminrhg = 0.001 ! ice volume (a*h in m) below which ice velocity equal ocean velocity48 ahi0 = 350.e0 ! horizontal eddy diffusivity coefficient for sea-ice (m2/s), depends on resolution 49 nevp = 120 ! number of EVP subcycles 50 relast = 0.333 ! ratio of elastic timescale to ice time step (1/3 if nevp=120, 1/9 if nevp=300) 51 hminrhg = 0.001 ! ice volume (a*h in m) below which ice velocity equals ocean velocity 52 52 / 53 53 !----------------------------------------------------------------------- 54 &namicethd ! ice thermodynamic 54 &namicethd ! ice thermodynamics 55 55 !----------------------------------------------------------------------- 56 hiccrit = 0.1 ! ice thickness for lateral accretion 57 ! caution 1.0, 1.0 best value to be used!!! (gilles G.) ???? 58 fraz_swi = 0 ! use of frazil ice collection thickness in function of wind (1.0) or not (0.0) 56 hiccrit = 0.1 ! thickness for new ice formation in open water 57 fraz_swi = 0 ! use frazil ice collection thickness as a function of wind (1.0) or not (0.0) 59 58 maxfrazb = 0.0 ! maximum portion of frazil ice collecting at the ice bottom 60 vfrazb = 0.41 66667! thresold drift speed for frazil ice collecting at the ice bottom59 vfrazb = 0.417 ! thresold drift speed for frazil ice collecting at the ice bottom 61 60 Cfrazb = 5.0 ! squeezing coefficient for frazil ice collecting at the ice bottom 62 61 hiclim = 0.10 ! minimum thickness of the first thickness category, must be smaller than hiccrit 63 62 betas = 0.6 ! exponent in lead-ice fractionation of snow precipitation 0.66 64 63 ! betas = 1 -> equipartition, betas < 1 -> more on leads 65 kappa_i = 1.0 ! extinction radiation parameter in sea ice ( 1.0)64 kappa_i = 1.0 ! extinction radiation parameter in sea ice (m-1) 66 65 nconv_i_thd = 50 ! maximal number of iterations for heat diffusion computation 67 maxer_i_thd = 0.0001 ! maximal error in temperature for heat diffusion computation 68 thcon_i_swi = 1 ! switch for computation of thermal conductivity in the ice 69 ! (0) Untersteiner (1964), (1) Pringle et al. (2007) 66 maxer_i_thd = 0.0001 ! temperature error tolerance after heat diffusion 67 thcon_i_swi = 1 ! sea ice thermal conductivity, Untersteiner_64 (0) or Pringle_et_al_07 (1) 68 nn_monocat = 0 ! virtual ITD mono-category parameterizations (1) or not (0); 1 ok for jpl=1 only 69 ! 2=replace ridging by piling, 3=activate G(he) only, 4=activate lateral melting only) --- those are temporary test options 70 70 / 71 71 !----------------------------------------------------------------------- 72 72 &namicesal ! ice salinity 73 73 !----------------------------------------------------------------------- 74 num_sal = 2 ! salinity option: 1 -> S = bulk_sal 75 ! 2 -> S = S(z,t) with a simple parameterization 76 ! 3 -> S = S(z) profile of Scwharzacher [1959] 77 ! 4 -> S = S(h) Cox and Weeks [1974] 74 num_sal = 2 ! S=bulk_sal (1), S(z,t) - (2), S(z) Schwarzacher_59 (3) 78 75 bulk_sal = 4.0 ! if 1 is used, it represents the ice salinity 79 sal_G = 5.00 ! restoring salinity for GD80 time_G = 1.728e+6 ! restoring time for GD81 sal_F = 2.00 ! restoring salinity for flushing82 time_F = 8.640e+5 ! restoring time for flushing83 s_i_max = 20.0 ! Maximum salinity84 s_i_min = 0.1 ! Minimum tolerated ice salinity76 sal_G = 5.00 ! restoring salinity, gravity drainage (g/kg) 77 time_G = 1.728e+6 ! restoring time scale, gravity drainage (s) 78 sal_F = 2.00 ! restoring salinity, flushing (g/kg) 79 time_F = 8.640e+5 ! restoring time, flushing (s) 80 s_i_max = 20.0 ! maximum tolerated ice salinity (g/kg) 81 s_i_min = 0.1 ! minimum tolerated ice salinity (g/kg) 85 82 / 86 83 !----------------------------------------------------------------------- 87 84 &namiceitdme ! parameters for mechanical redistribution of ice 88 85 !----------------------------------------------------------------------- 89 ridge_scheme_swi = 0 ! which ridging scheme using (1=Rothrock,else=Hibler79)90 Cs = 0.50 ! shearing energy contributionto ridging91 Cf = 17.0 ! ratio of ridging work to PEchange in ridging92 fsnowrdg = 0.5 ! snow fraction that survives in ridging93 fsnowrft = 0.5 ! snow fraction that survives in rafting94 Gstar = 0.15 ! fractional area of thin ice being ridged 95 astar = 0.05 ! equivalent of gstar (0.05 for TH75 and 0.03 for weaker ice)96 Hstar = 100.0 ! parameter determining the maximum thickness of ridged ice97 raft_swi = 1 ! rafting or not98 hparmeter = 0.75 ! threshold thickness for rafting or not86 ridge_scheme_swi = 0 ! ice strength parameteriztaion, Hibler_79 (0), Rothrock_75 (1) 87 Cs = 0.50 ! ratio of shearing energy contributing to ridging 88 Cf = 17.0 ! ratio of ridging work to potential energy change in ridging 89 fsnowrdg = 0.5 ! snow volume fraction that survives in ridging 90 fsnowrft = 0.5 ! snow volume fraction that survives in rafting 91 Gstar = 0.15 ! fractional area of thin ice being ridged (partfun_swi = 0) 92 astar = 0.05 ! measures fractional area of thin ice being ridged (partfun_swi = 1) 93 Hstar = 100.0 ! determines the maximum thickness of ridged ice (m) 94 raft_swi = 1 ! rafting activated (1) or not (0) 95 hparmeter = 0.75 ! threshold thickness for rafting (m) 99 96 Craft = 5.0 ! coefficient used in the rafting function 100 ridge_por = 0.3 ! initial porosity of the ridged ice (typically 0.30)101 partfun_swi = 1 ! participation function linear, TH75 (0) or exponential Letal07 (1)102 brinstren_swi = 0 ! (1) use brine volume to diminish ice strength97 ridge_por = 0.3 ! porosity of newly ridged ice 98 partfun_swi = 1 ! participation function: linear, Thorndike et al 75 (0) or exponential, Lipscomb 07 (1) 99 brinstren_swi = 0 ! ice strength is a function of brine volume fraction (1) or not (0) -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/ice.F90
r5047 r5048 165 165 REAL(wp), PUBLIC :: r1_rdtice !: = 1. / rdt_ice 166 166 167 ! !!** ice-dynamic namelist (namicedyn) ** 167 ! Clem je sais pas ou les mettre 168 INTEGER , PUBLIC :: nn_hibnd !: thickness category boundaries: tanh function (1), or h^(-alpha) function (2) 169 REAL(wp), PUBLIC :: rn_hibnd !: mean thickness of the domain (used to compute category boundaries, nn_hibnd = 2 only) 170 171 ! !!** ice-dynamics namelist (namicedyn) ** 168 172 INTEGER , PUBLIC :: nevp !: number of iterations for subcycling 169 173 REAL(wp), PUBLIC :: cw !: drag coefficient for oceanic stress … … 191 195 ! ! 3 - salinity profile, constant in time 192 196 INTEGER , PUBLIC :: thcon_i_swi !: thermal conductivity: =0 Untersteiner (1964) ; =1 Pringle et al (2007) 197 INTEGER , PUBLIC :: nn_monocat !: virtual ITD mono-category parameterizations (1) or not (0) 193 198 194 199 ! !!** ice-mechanical redistribution namelist (namiceitdme) … … 260 265 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: wfx_res !: residual component of wfx_ice [kg/m2] 261 266 267 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: afx_tot !: ice concentration tendency (total) [s-1] 268 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: afx_thd !: ice concentration tendency (thermodynamics) [s-1] 269 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: afx_dyn !: ice concentration tendency (dynamics) [s-1] 270 262 271 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: sfx_bog !: salt flux due to ice growth/melt [PSU/m2/s] 263 272 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: sfx_bom !: salt flux due to ice growth/melt [PSU/m2/s] … … 290 299 291 300 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ftr_ice !: transmitted solar radiation under ice 301 302 ! lateral melting (mono-category) 303 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: dh_i_melt !: net melting (m) 292 304 293 305 ! temporary arrays for dummy version of the code … … 436 448 & wfx_bog(jpi,jpj) , wfx_dyn(jpi,jpj) , wfx_bom(jpi,jpj) , wfx_sum(jpi,jpj) , & 437 449 & wfx_res(jpi,jpj) , wfx_sni(jpi,jpj) , wfx_opw(jpi,jpj) , wfx_spr(jpi,jpj) , qlead (jpi,jpj) , & 438 & fhtur (jpi,jpj) , ftr_ice(jpi,jpj,jpl) , & 450 & afx_tot(jpi,jpj) , afx_thd(jpi,jpj), afx_dyn(jpi,jpj) , & 451 & fhtur (jpi,jpj) , ftr_ice(jpi,jpj,jpl) , dh_i_melt(jpi,jpj,jpl) , & 439 452 & sfx_res(jpi,jpj) , sfx_bri(jpi,jpj) , sfx_dyn(jpi,jpj) , & 440 453 & sfx_bog(jpi,jpj) , sfx_bom(jpi,jpj) , sfx_sum(jpi,jpj) , sfx_sni(jpi,jpj) , sfx_opw(jpi,jpj) , & -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/iceini.F90
r5047 r5048 172 172 !! limistate (only) and is changed to 99 m in ice_init 173 173 !!------------------------------------------------------------------ 174 INTEGER :: jl ! dummy loop index 175 REAL(wp) :: zc1, zc2, zc3, zx1 ! local scalars 174 INTEGER :: jl ! dummy loop index 175 REAL(wp) :: zc1, zc2, zc3, zx1 ! local scalars 176 REAL(wp) :: zhmax, znum, zden, zalpha ! 176 177 !!------------------------------------------------------------------ 177 178 … … 189 190 !- Thickness categories boundaries 190 191 !---------------------------------- 192 193 ! Clem - je sais pas encore dans quelle namelist les mettre, ca depend des chgts liés à bdy 194 nn_hibnd = 2 ! thickness category boundaries: tanh function (1) h^(-alpha) (2) 195 rn_himean = 2.5 ! mean thickness of the domain (used to compute category boundaries, nn_hibnd = 2 only) 196 191 197 hi_max(:) = 0._wp 192 198 193 zc1 = 3._wp / REAL( jpl, wp ) 194 zc2 = 10._wp * zc1 195 zc3 = 3._wp 196 197 DO jl = 1, jpl 198 zx1 = REAL( jl-1, wp ) / REAL( jpl, wp ) 199 hi_max(jl) = hi_max(jl-1) + zc1 + zc2 * (1._wp + TANH( zc3 * (zx1 - 1._wp ) ) ) 200 END DO 199 SELECT CASE ( nn_hbnd ) 200 !---------------------- 201 CASE (1) ! tanh function (CICE) 202 !---------------------- 203 204 zc1 = 3._wp / REAL( jpl, wp ) 205 zc2 = 10._wp * zc1 206 zc3 = 3._wp 207 208 DO jl = 1, jpl 209 zx1 = REAL( jl-1, wp ) / REAL( jpl, wp ) 210 hi_max(jl) = hi_max(jl-1) + zc1 + zc2 * (1._wp + TANH( zc3 * (zx1 - 1._wp ) ) ) 211 END DO 212 213 !---------------------- 214 CASE (2) ! h^(-alpha) function 215 !---------------------- 216 217 zalpha = 0.05 ! exponent of the transform function 218 219 zhmax = 3.*rn_himean 220 221 DO jl = 1, jpl 222 znum = jpl * ( zhmax+1 )**zalpha 223 zden = ( jpl - jl ) * ( zhmax+1 )**zalpha + jl 224 hi_max(jl) = ( znum / zden )**(1./zalpha) - 1 225 END DO 226 227 END SELECT 228 229 ! mv- would be nice if we could put the instruction hi_max(jpl) = 99 here 201 230 202 231 IF(lwp) WRITE(numout,*) ' Thickness category boundaries ' -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90
r4990 r5048 93 93 CALL lim_thd_lac 94 94 CALL lim_var_glo2eqv ! only for info 95 96 ! MV: Could put lateral melting here, would be better I think ??? 97 95 98 96 99 IF(ln_ctl) THEN ! Control print -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90
r5047 r5048 83 83 !! 84 84 INTEGER :: ji, jj, jk, jl ! dummy loop indices 85 INTEGER :: nbpb ! nb of icy pts for thermo. cal. 85 INTEGER :: nbpb ! nb of icy pts for vertical thermo calculations 86 INTEGER :: nbplm ! nb of icy pts for lateral melting calculations (mono-cat) 86 87 INTEGER :: ii, ij ! temporary dummy loop index 87 88 REAL(wp) :: zfric_umin = 0._wp ! lower bound for the friction velocity (cice value=5.e-04) … … 434 435 CALL tab_1d_2d( nbpb, hfx_res , npb, hfx_res_1d(1:nbpb) , jpi, jpj ) 435 436 CALL tab_1d_2d( nbpb, hfx_err_rem , npb, hfx_err_rem_1d(1:nbpb) , jpi, jpj ) 437 438 IF ( ( ( nn_monocat .EQ. 1 ) .OR. ( nn_monocat .EQ. 4 ) ) .AND. ( jpl == 1 ) ) THEN 439 CALL tab_1d_2d( nbpb, dh_i_melt(:,:,jl) , npb, dh_i_melt_1d(1:nbpb) , jpi, jpj ) 440 ENDIF 436 441 ! 437 442 CALL tab_1d_2d( nbpb, qns_ice(:,:,jl), npb, qns_ice_1d(1:nbpb) , jpi, jpj) … … 471 476 !---------------------------------- 472 477 CALL lim_var_eqv2glo 478 479 !---------------------------------- 480 ! 5.X) Lateral melting 481 !---------------------------------- 482 !!! declare dh_i_melt (ok), dh_i_melt_1d (ok), nbplm (ok), nplm (ok), zda(ok) 483 484 IF ( ( ( nn_monocat .EQ. 1 ) .OR. ( nn_monocat .EQ. 4 ) ) .AND. ( jpl == 1 ) ) THEN 485 486 WRITE(numout,*) ' Lateral melting ON ' 487 488 ! select points where lateral melting occurs 489 jl = 1 490 491 nbplm = 0 492 DO jj = 1, jpj 493 DO ji = 1, jpi 494 IF ( ( dh_i_melt(ji,jj,jl) .LT.-epsi10 ) .AND. & 495 & ( ht_i(ji,jj,jl) - dh_i_melt(ji,jj,jl) .GT. epsi10 ) .AND. & 496 & ( ht_i(ji,jj,jl) .GT. epsi10 ) ) THEN 497 nbplm = nbplm + 1 498 nplm(nbplm) = (jj - 1) * jpi + ji 499 ENDIF 500 END DO 501 END DO 502 503 IF( nbplm > 0 ) THEN ! If there is no net melting, do nothing 504 505 ! Move to 1D arrays 506 !------------------------- 507 508 CALL tab_2d_1d( nbplm, a_i_1d (1:nbplm), a_i(:,:,jl) , jpi, jpj, nplm(1:nbplm) ) 509 CALL tab_2d_1d( nbplm, ht_i_1d (1:nbplm), ht_i(:,:,jl) , jpi, jpj, nplm(1:nbplm) ) 510 CALL tab_2d_1d( nbplm, dh_i_melt_1d(1:nbplm), dh_i_melt(:,:,jl) , jpi, jpj, nplm(1:nbplm) ) 511 512 ! Compute lateral melting (dA = A/2h dh ) 513 DO ji = 1, nbplm 514 zda = a_i_1d(ji) * dh_i_melt_1d(ji) / ( 2._wp * ht_i_1d(ji) ) 515 a_i_1d(ji) = a_i_1d(ji) + zda 516 END DO 517 518 ! Move back to 2D arrays 519 !------------------------- 520 CALL tab_1d_2d( nbplm, a_i (:,:,jl) , nplm, a_i_1d (1:nbplm) , jpi, jpj ) 521 at_i(:,:) = a_i(:,:,jl) 522 523 ENDIF 524 525 ENDIF 473 526 474 527 !-------------------------------------------- … … 563 616 !!------------------------------------------------------------------- 564 617 INTEGER :: ios ! Local integer output status for namelist read 565 NAMELIST/namicethd/ hiccrit, fraz_swi, maxfrazb, vfrazb, Cfrazb, &618 NAMELIST/namicethd/ hiccrit, fraz_swi, maxfrazb, vfrazb, Cfrazb, & 566 619 & hiclim, parsub, betas, & 567 & kappa_i, nconv_i_thd, maxer_i_thd, thcon_i_swi 620 & kappa_i, nconv_i_thd, maxer_i_thd, thcon_i_swi, & 621 & nn_monocat 568 622 !!------------------------------------------------------------------- 569 623 ! … … 582 636 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicethd in configuration namelist', lwp ) 583 637 IF(lwm) WRITE ( numoni, namicethd ) 638 ! 639 IF ( ( jpl > 1 ) .AND. ( nn_monocat == 1 ) ) THEN 640 nn_monocat = 0 641 WRITE(numout, *) ' nn_monocat must be 0 in multi-category case ' 642 ENDIF 584 643 585 644 IF( lk_cpl .AND. parsub /= 0.0 ) CALL ctl_stop( 'In coupled mode, use parsub = 0. or send dqla' ) … … 597 656 WRITE(numout,*)' switch for snow sublimation (=1) or not (=0) parsub = ', parsub 598 657 WRITE(numout,*)' coefficient for ice-lead partition of snowfall betas = ', betas 599 WRITE(numout,*)' extinction radiation parameter in sea ice (1.0)kappa_i = ', kappa_i658 WRITE(numout,*)' extinction radiation parameter in sea ice kappa_i = ', kappa_i 600 659 WRITE(numout,*)' maximal n. of iter. for heat diffusion computation nconv_i_thd = ', nconv_i_thd 601 660 WRITE(numout,*)' maximal err. on T for heat diffusion computation maxer_i_thd = ', maxer_i_thd 602 661 WRITE(numout,*)' switch for comp. of thermal conductivity in the ice thcon_i_swi = ', thcon_i_swi 603 662 WRITE(numout,*)' check heat conservation in the ice/snow con_i = ', con_i 663 WRITE(numout,*)' virtual ITD mono-category parameterizations (1) or not nn_monocat = ', nn_monocat 604 664 ENDIF 605 665 ! -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90
r5047 r5048 688 688 END DO 689 689 END DO 690 691 ! 692 !------------------------------------------------------------- 693 ! Net melting (input for thin ice melting, mono-category case) 694 !------------------------------------------------------------- 695 ! 696 IF ( ( ( nn_monocat .EQ. 1 ) .OR. ( nn_monocat .EQ. 4 ) ) .AND. ( jpl == 1 ) ) THEN 697 DO ji = kideb, kiut 698 dh_i_melt_1d(ji) = MIN( dh_i_surf(ji) + dh_i_bott(ji) + dh_snowice(ji), 0._wp ) 699 END DO 700 ENDIF 690 701 691 702 CALL wrk_dealloc( jpij, zh_s, zqprec, zq_su, zq_bo, zf_tt, zq_rema ) -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90
r5047 r5048 100 100 INTEGER :: nconv ! number of iterations in iterative procedure 101 101 INTEGER :: minnumeqmin, maxnumeqmax 102 102 103 INTEGER, POINTER, DIMENSION(:) :: numeqmin ! reference number of top equation 103 104 INTEGER, POINTER, DIMENSION(:) :: numeqmax ! reference number of bottom equation 104 105 INTEGER, POINTER, DIMENSION(:) :: isnow ! switch for presence (1) or absence (0) of snow 106 105 107 REAL(wp) :: zg1s = 2._wp ! for the tridiagonal system 106 108 REAL(wp) :: zg1 = 2._wp ! … … 113 115 REAL(wp) :: zerritmax ! current maximal error on temperature 114 116 REAL(wp) :: zhsu 115 REAL(wp), POINTER, DIMENSION(:) :: ztfs ! ice melting point 116 REAL(wp), POINTER, DIMENSION(:) :: ztsub ! old surface temperature (before the iterative procedure ) 117 REAL(wp), POINTER, DIMENSION(:) :: ztsubit ! surface temperature at previous iteration 118 REAL(wp), POINTER, DIMENSION(:) :: zh_i ! ice layer thickness 119 REAL(wp), POINTER, DIMENSION(:) :: zh_s ! snow layer thickness 120 REAL(wp), POINTER, DIMENSION(:) :: zfsw ! solar radiation absorbed at the surface 121 REAL(wp), POINTER, DIMENSION(:) :: zf ! surface flux function 122 REAL(wp), POINTER, DIMENSION(:) :: dzf ! derivative of the surface flux function 123 REAL(wp), POINTER, DIMENSION(:) :: zerrit ! current error on temperature 124 REAL(wp), POINTER, DIMENSION(:) :: zdifcase ! case of the equation resolution (1->4) 125 REAL(wp), POINTER, DIMENSION(:) :: zftrice ! solar radiation transmitted through the ice 126 REAL(wp), POINTER, DIMENSION(:) :: zihic 127 REAL(wp), POINTER, DIMENSION(:,:) :: ztcond_i ! Ice thermal conductivity 128 REAL(wp), POINTER, DIMENSION(:,:) :: zradtr_i ! Radiation transmitted through the ice 129 REAL(wp), POINTER, DIMENSION(:,:) :: zradab_i ! Radiation absorbed in the ice 130 REAL(wp), POINTER, DIMENSION(:,:) :: zkappa_i ! Kappa factor in the ice 131 REAL(wp), POINTER, DIMENSION(:,:) :: ztib ! Old temperature in the ice 132 REAL(wp), POINTER, DIMENSION(:,:) :: zeta_i ! Eta factor in the ice 133 REAL(wp), POINTER, DIMENSION(:,:) :: ztitemp ! Temporary temperature in the ice to check the convergence 134 REAL(wp), POINTER, DIMENSION(:,:) :: zspeche_i ! Ice specific heat 135 REAL(wp), POINTER, DIMENSION(:,:) :: z_i ! Vertical cotes of the layers in the ice 136 REAL(wp), POINTER, DIMENSION(:,:) :: zradtr_s ! Radiation transmited through the snow 137 REAL(wp), POINTER, DIMENSION(:,:) :: zradab_s ! Radiation absorbed in the snow 138 REAL(wp), POINTER, DIMENSION(:,:) :: zkappa_s ! Kappa factor in the snow 139 REAL(wp), POINTER, DIMENSION(:,:) :: zeta_s ! Eta factor in the snow 140 REAL(wp), POINTER, DIMENSION(:,:) :: ztstemp ! Temporary temperature in the snow to check the convergence 141 REAL(wp), POINTER, DIMENSION(:,:) :: ztsb ! Temporary temperature in the snow 142 REAL(wp), POINTER, DIMENSION(:,:) :: z_s ! Vertical cotes of the layers in the snow 143 REAL(wp), POINTER, DIMENSION(:,:) :: zswiterm ! Independent term 144 REAL(wp), POINTER, DIMENSION(:,:) :: zswitbis ! temporary independent term 145 REAL(wp), POINTER, DIMENSION(:,:) :: zdiagbis 146 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrid ! tridiagonal system terms 117 118 REAL(wp), POINTER, DIMENSION(:) :: ztfs ! ice melting point 119 REAL(wp), POINTER, DIMENSION(:) :: ztsub ! old surface temperature (before the iterative procedure ) 120 REAL(wp), POINTER, DIMENSION(:) :: ztsubit ! surface temperature at previous iteration 121 REAL(wp), POINTER, DIMENSION(:) :: zh_i ! ice layer thickness 122 REAL(wp), POINTER, DIMENSION(:) :: zh_s ! snow layer thickness 123 REAL(wp), POINTER, DIMENSION(:) :: zfsw ! solar radiation absorbed at the surface 124 REAL(wp), POINTER, DIMENSION(:) :: zf ! surface flux function 125 REAL(wp), POINTER, DIMENSION(:) :: dzf ! derivative of the surface flux function 126 REAL(wp), POINTER, DIMENSION(:) :: zerrit ! current error on temperature 127 REAL(wp), POINTER, DIMENSION(:) :: zdifcase ! case of the equation resolution (1->4) 128 REAL(wp), POINTER, DIMENSION(:) :: zftrice ! solar radiation transmitted through the ice 129 REAL(wp), POINTER, DIMENSION(:) :: zihic 130 131 REAL(wp), POINTER, DIMENSION(:,:) :: ztcond_i ! Ice thermal conductivity 132 REAL(wp), POINTER, DIMENSION(:,:) :: zradtr_i ! Radiation transmitted through the ice 133 REAL(wp), POINTER, DIMENSION(:,:) :: zradab_i ! Radiation absorbed in the ice 134 REAL(wp), POINTER, DIMENSION(:,:) :: zkappa_i ! Kappa factor in the ice 135 REAL(wp), POINTER, DIMENSION(:,:) :: ztib ! Old temperature in the ice 136 REAL(wp), POINTER, DIMENSION(:,:) :: zeta_i ! Eta factor in the ice 137 REAL(wp), POINTER, DIMENSION(:,:) :: ztitemp ! Temporary temperature in the ice to check the convergence 138 REAL(wp), POINTER, DIMENSION(:,:) :: zspeche_i ! Ice specific heat 139 REAL(wp), POINTER, DIMENSION(:,:) :: z_i ! Vertical cotes of the layers in the ice 140 REAL(wp), POINTER, DIMENSION(:,:) :: zradtr_s ! Radiation transmited through the snow 141 REAL(wp), POINTER, DIMENSION(:,:) :: zradab_s ! Radiation absorbed in the snow 142 REAL(wp), POINTER, DIMENSION(:,:) :: zkappa_s ! Kappa factor in the snow 143 REAL(wp), POINTER, DIMENSION(:,:) :: zeta_s ! Eta factor in the snow 144 REAL(wp), POINTER, DIMENSION(:,:) :: ztstemp ! Temporary temperature in the snow to check the convergence 145 REAL(wp), POINTER, DIMENSION(:,:) :: ztsb ! Temporary temperature in the snow 146 REAL(wp), POINTER, DIMENSION(:,:) :: z_s ! Vertical cotes of the layers in the snow 147 REAL(wp), POINTER, DIMENSION(:,:) :: zindterm ! 'Ind'ependent term 148 REAL(wp), POINTER, DIMENSION(:,:) :: zindtbis ! Temporary 'ind'ependent term 149 REAL(wp), POINTER, DIMENSION(:,:) :: zdiagbis ! Temporary 'dia'gonal term 150 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrid ! Tridiagonal system terms 151 147 152 ! diag errors on heat 148 REAL(wp), POINTER, DIMENSION(:) :: zdq, zq_ini, zhfx_err 153 REAL(wp), POINTER, DIMENSION(:) :: zdq, zq_ini, zhfx_err 154 155 ! Mono-category 156 REAL(wp) :: zepsilon ! determines thres. above which computation of G(h) is done 157 REAL(wp) :: zratio_s ! dummy factor 158 REAL(wp) :: zratio_i ! dummy factor 159 REAL(wp) :: zh_thres ! thickness thres. for G(h) computation 160 REAL(wp) :: zhe ! dummy factor 161 REAL(wp) :: zswitch ! dummy switch 162 REAL(wp) :: zkimean ! mean sea ice thermal conductivity 163 REAL(wp) :: zfac ! dummy factor 164 REAL(wp) :: zihe ! dummy factor 165 REAL(wp) :: zheshth ! dummy factor 166 167 REAL(wp), POINTER, DIMENSION(:) :: zghe ! G(he), th. conduct enhancement factor, mono-cat 168 149 169 !!------------------------------------------------------------------ 150 170 ! 151 171 CALL wrk_alloc( jpij, numeqmin, numeqmax, isnow ) 152 172 CALL wrk_alloc( jpij, ztfs, ztsub, ztsubit, zh_i, zh_s, zfsw ) 153 CALL wrk_alloc( jpij, zf, dzf, zerrit, zdifcase, zftrice, zihic )173 CALL wrk_alloc( jpij, zf, dzf, zerrit, zdifcase, zftrice, zihic, zghe ) 154 174 CALL wrk_alloc( jpij, nlay_i+1, ztcond_i, zradtr_i, zradab_i, zkappa_i, ztib, zeta_i, ztitemp, z_i, zspeche_i, kjstart=0) 155 175 CALL wrk_alloc( jpij, nlay_s+1, zradtr_s, zradab_s, zkappa_s, ztsb, zeta_s, ztstemp, z_s, kjstart=0) 156 CALL wrk_alloc( jpij, nlay_i+3, z switerm, zswitbis, zdiagbis )176 CALL wrk_alloc( jpij, nlay_i+3, zindterm, zindtbis, zdiagbis ) 157 177 CALL wrk_alloc( jpij, nlay_i+3, 3, ztrid ) 158 178 … … 200 220 ! 201 221 !------------------------------------------------------------------------------| 202 ! 2) Radiation s|222 ! 2) Radiation | 203 223 !------------------------------------------------------------------------------| 204 224 ! … … 213 233 ! zftrice = io.qsr_ice is below the surface 214 234 ! ftr_ice = io.qsr_ice.exp(-k(h_i)) transmitted below the ice 235 ! fr1_i0_1d = i0 for a thin ice cover, fr1_i0_2d = i0 for a thick ice cover 215 236 zhsu = 0.1_wp ! threshold for the computation of i0 216 !fr1_i0_1d = i0 for a thin ice surface217 !fr1_i0_2d = i0 for a thick ice surface218 237 DO ji = kideb , kiut 219 238 ! switches … … 339 358 END DO 340 359 ENDIF 341 ! 342 !------------------------------------------------------------------------------| 343 ! 5) kappa factors | 344 !------------------------------------------------------------------------------| 345 ! 360 361 ! 362 !------------------------------------------------------------------------------| 363 ! 6) G(he) - enhancement of thermal conductivity in mono-category case | 364 !------------------------------------------------------------------------------| 365 ! 366 ! Computation of effective thermal conductivity G(h) 367 ! Used in mono-category case only to simulate an ITD implicitly 368 ! Fichefet and Morales Maqueda, JGR 1997 369 370 zghe(:) = 0._wp 371 372 SELECT CASE ( nn_monocat ) 373 374 CASE (0,2,4) 375 376 zghe(kideb:kiut) = 1._wp 377 378 CASE (1,3) ! LIM3 379 380 zepsilon = 0.1 381 zh_thres = EXP( 1._wp ) * zepsilon / 2. 382 383 DO ji = kideb, kiut 384 385 ! Mean sea ice thermal conductivity 386 zkimean = SUM( ztcond_i(ji,0:nlay_i) ) / REAL(nlay_i+1,wp) 387 388 ! Effective thickness he (zhe) 389 zfac = 1._wp / ( rcdsn + zkimean ) 390 zratio_s = rcdsn * zfac 391 zratio_i = zkimean * zfac 392 zhe = zratio_s * ht_i_1d(ji) + zratio_i * ht_s_1d(ji) 393 394 ! G(he) 395 zswitch = MAX( 0._wp , SIGN( 1._wp , zhe - zh_thres ) ) ! =0 if zhe < zh_thres, if > 396 zghe(ji) = ( 1.0 - zswitch ) + zswitch * ( 0.5 + 0.5 * LOG( 2.*zhe / zepsilon ) ) 397 398 ! Impose G(he) < 2. 399 zghe(ji) = MIN( zghe(ji), 2.0 ) 400 401 END DO 402 403 END SELECT 404 405 ! 406 !------------------------------------------------------------------------------| 407 ! 7) kappa factors | 408 !------------------------------------------------------------------------------| 409 ! 410 !--- Snow 346 411 DO ji = kideb, kiut 347 348 !-- Snow kappa factors 349 zkappa_s(ji,0) = rcdsn / MAX(epsi10,zh_s(ji)) 350 zkappa_s(ji,nlay_s) = rcdsn / MAX(epsi10,zh_s(ji)) 412 zfac = 1. / MAX( epsi10 , zh_s(ji) ) 413 zkappa_s(ji,0) = zghe(ji) * rcdsn * zfac 414 zkappa_s(ji,nlay_s) = zghe(ji) * rcdsn * zfac 351 415 END DO 352 416 353 417 DO jk = 1, nlay_s-1 354 418 DO ji = kideb , kiut 355 zkappa_s(ji,jk) = 2.0 * rcdsn / &356 MAX(epsi10,2.0*zh_s(ji))357 358 END DO 359 419 zkappa_s(ji,jk) = zghe(ji) * 2.0 * rcdsn / MAX( epsi10, 2.0*zh_s(ji) ) 420 END DO 421 END DO 422 423 !--- Ice 360 424 DO jk = 1, nlay_i-1 361 425 DO ji = kideb , kiut 362 !-- Ice kappa factors 363 zkappa_i(ji,jk) = 2.0*ztcond_i(ji,jk)/ & 364 MAX(epsi10,2.0*zh_i(ji)) 365 END DO 366 END DO 367 368 DO ji = kideb , kiut 369 zkappa_i(ji,0) = ztcond_i(ji,0)/MAX(epsi10,zh_i(ji)) 370 zkappa_i(ji,nlay_i) = ztcond_i(ji,nlay_i) / MAX(epsi10,zh_i(ji)) 371 !-- Interface 372 zkappa_s(ji,nlay_s) = 2.0*rcdsn*ztcond_i(ji,0)/MAX(epsi10, & 373 (ztcond_i(ji,0)*zh_s(ji) + rcdsn*zh_i(ji))) 374 zkappa_i(ji,0) = zkappa_s(ji,nlay_s)*REAL( isnow(ji) ) & 375 + zkappa_i(ji,0)*REAL( 1 - isnow(ji) ) 376 END DO 377 ! 378 !------------------------------------------------------------------------------| 379 ! 6) Sea ice specific heat, eta factors | 426 zkappa_i(ji,jk) = zghe(ji) * 2.0 * ztcond_i(ji,jk) / MAX( epsi10 , 2.0*zh_i(ji) ) 427 END DO 428 END DO 429 430 !--- Snow-ice interface 431 DO ji = kideb , kiut 432 zfac = 1./ MAX( epsi10 , zh_i(ji) ) 433 zkappa_i(ji,0) = zghe(ji) * ztcond_i(ji,0) * zfac 434 zkappa_i(ji,nlay_i) = zghe(ji) * ztcond_i(ji,nlay_i) * zfac 435 zkappa_s(ji,nlay_s) = zghe(ji) * zghe(ji) * 2.0 * rcdsn * ztcond_i(ji,0) / & 436 & MAX( epsi10, ( zghe(ji) * ztcond_i(ji,0) * zh_s(ji) + zghe(ji) * rcdsn * zh_i(ji) ) ) 437 zkappa_i(ji,0) = zkappa_s(ji,nlay_s)*REAL( isnow(ji), wp ) + zkappa_i(ji,0)*REAL( 1 - isnow(ji), wp ) 438 END DO 439 440 ! 441 !------------------------------------------------------------------------------| 442 ! 8) Sea ice specific heat, eta factors | 380 443 !------------------------------------------------------------------------------| 381 444 ! … … 383 446 DO ji = kideb , kiut 384 447 ztitemp(ji,jk) = t_i_1d(ji,jk) 385 zspeche_i(ji,jk) = cpic + zgamma*s_i_1d(ji,jk)/ & 386 MAX((t_i_1d(ji,jk)-rtt)*(ztib(ji,jk)-rtt),epsi10) 387 zeta_i(ji,jk) = rdt_ice / MAX(rhoic*zspeche_i(ji,jk)*zh_i(ji), & 388 epsi10) 448 zspeche_i(ji,jk) = cpic + zgamma*s_i_1d(ji,jk)/ MAX( (t_i_1d(ji,jk)-rtt)*(ztib(ji,jk)-rtt) , epsi10 ) 449 zeta_i(ji,jk) = rdt_ice / MAX( rhoic*zspeche_i(ji,jk)*zh_i(ji), epsi10 ) 389 450 END DO 390 451 END DO … … 396 457 END DO 397 458 END DO 398 ! 399 !------------------------------------------------------------------------------| 400 ! 7) surface flux computation | 459 460 ! 461 !------------------------------------------------------------------------------| 462 ! 9) surface flux computation | 401 463 !------------------------------------------------------------------------------| 402 464 ! … … 418 480 ! 419 481 !------------------------------------------------------------------------------| 420 ! 8) tridiagonal system terms|482 ! 10) tridiagonal system terms | 421 483 !------------------------------------------------------------------------------| 422 484 ! … … 433 495 ztrid(ji,numeq,2) = 0. 434 496 ztrid(ji,numeq,3) = 0. 435 z switerm(ji,numeq)= 0.436 z switbis(ji,numeq)= 0.497 zindterm(ji,numeq)= 0. 498 zindtbis(ji,numeq)= 0. 437 499 zdiagbis(ji,numeq)= 0. 438 500 ENDDO … … 446 508 zkappa_i(ji,jk)) 447 509 ztrid(ji,numeq,3) = - zeta_i(ji,jk)*zkappa_i(ji,jk) 448 z switerm(ji,numeq) = ztib(ji,jk) + zeta_i(ji,jk)* &510 zindterm(ji,numeq) = ztib(ji,jk) + zeta_i(ji,jk)* & 449 511 zradab_i(ji,jk) 450 512 END DO … … 458 520 + zkappa_i(ji,nlay_i-1) ) 459 521 ztrid(ji,numeq,3) = 0.0 460 z switerm(ji,numeq) = ztib(ji,nlay_i) + zeta_i(ji,nlay_i)* &522 zindterm(ji,numeq) = ztib(ji,nlay_i) + zeta_i(ji,nlay_i)* & 461 523 ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i)*zg1 & 462 524 * t_bo_1d(ji) ) … … 478 540 zkappa_s(ji,jk) ) 479 541 ztrid(ji,numeq,3) = - zeta_s(ji,jk)*zkappa_s(ji,jk) 480 z switerm(ji,numeq) = ztsb(ji,jk) + zeta_s(ji,jk)* &542 zindterm(ji,numeq) = ztsb(ji,jk) + zeta_s(ji,jk)* & 481 543 zradab_s(ji,jk) 482 544 END DO … … 485 547 IF ( nlay_i.eq.1 ) THEN 486 548 ztrid(ji,nlay_s+2,3) = 0.0 487 z switerm(ji,nlay_s+2) = zswiterm(ji,nlay_s+2) + zkappa_i(ji,1)* &549 zindterm(ji,nlay_s+2) = zindterm(ji,nlay_s+2) + zkappa_i(ji,1)* & 488 550 t_bo_1d(ji) 489 551 ENDIF … … 502 564 ztrid(ji,1,2) = dzf(ji) - zg1s*zkappa_s(ji,0) 503 565 ztrid(ji,1,3) = zg1s*zkappa_s(ji,0) 504 z switerm(ji,1) = dzf(ji)*t_su_1d(ji) - zf(ji)566 zindterm(ji,1) = dzf(ji)*t_su_1d(ji) - zf(ji) 505 567 506 568 !!first layer of snow equation … … 508 570 ztrid(ji,2,2) = 1.0 + zeta_s(ji,1)*(zkappa_s(ji,1) + zkappa_s(ji,0)*zg1s) 509 571 ztrid(ji,2,3) = - zeta_s(ji,1)* zkappa_s(ji,1) 510 z switerm(ji,2) = ztsb(ji,1) + zeta_s(ji,1)*zradab_s(ji,1)572 zindterm(ji,2) = ztsb(ji,1) + zeta_s(ji,1)*zradab_s(ji,1) 511 573 512 574 ELSE … … 525 587 zkappa_s(ji,0) * zg1s ) 526 588 ztrid(ji,2,3) = - zeta_s(ji,1)*zkappa_s(ji,1) 527 z switerm(ji,2) = ztsb(ji,1) + zeta_s(ji,1) * &589 zindterm(ji,2) = ztsb(ji,1) + zeta_s(ji,1) * & 528 590 ( zradab_s(ji,1) + & 529 591 zkappa_s(ji,0) * zg1s * t_su_1d(ji) ) … … 549 611 ztrid(ji,numeqmin(ji),2) = dzf(ji) - zkappa_i(ji,0)*zg1 550 612 ztrid(ji,numeqmin(ji),3) = zkappa_i(ji,0)*zg1 551 z switerm(ji,numeqmin(ji)) = dzf(ji)*t_su_1d(ji) - zf(ji)613 zindterm(ji,numeqmin(ji)) = dzf(ji)*t_su_1d(ji) - zf(ji) 552 614 553 615 !!first layer of ice equation … … 556 618 + zkappa_i(ji,0) * zg1 ) 557 619 ztrid(ji,numeqmin(ji)+1,3) = - zeta_i(ji,1)*zkappa_i(ji,1) 558 z switerm(ji,numeqmin(ji)+1)= ztib(ji,1) + zeta_i(ji,1)*zradab_i(ji,1)620 zindterm(ji,numeqmin(ji)+1)= ztib(ji,1) + zeta_i(ji,1)*zradab_i(ji,1) 559 621 560 622 !!case of only one layer in the ice (surface & ice equations are altered) … … 569 631 ztrid(ji,numeqmin(ji)+1,3) = 0.0 570 632 571 z switerm(ji,numeqmin(ji)+1) = ztib(ji,1) + zeta_i(ji,1)* &633 zindterm(ji,numeqmin(ji)+1) = ztib(ji,1) + zeta_i(ji,1)* & 572 634 ( zradab_i(ji,1) + zkappa_i(ji,1)*t_bo_1d(ji) ) 573 635 ENDIF … … 589 651 zg1) 590 652 ztrid(ji,numeqmin(ji),3) = - zeta_i(ji,1) * zkappa_i(ji,1) 591 z switerm(ji,numeqmin(ji)) = ztib(ji,1) + zeta_i(ji,1)*( zradab_i(ji,1) + &653 zindterm(ji,numeqmin(ji)) = ztib(ji,1) + zeta_i(ji,1)*( zradab_i(ji,1) + & 592 654 zkappa_i(ji,0) * zg1 * t_su_1d(ji) ) 593 655 … … 598 660 zkappa_i(ji,1)) 599 661 ztrid(ji,numeqmin(ji),3) = 0.0 600 z switerm(ji,numeqmin(ji)) = ztib(ji,1) + zeta_i(ji,1)* &662 zindterm(ji,numeqmin(ji)) = ztib(ji,1) + zeta_i(ji,1)* & 601 663 (zradab_i(ji,1) + zkappa_i(ji,1)*t_bo_1d(ji)) & 602 664 + t_su_1d(ji)*zeta_i(ji,1)*zkappa_i(ji,0)*2.0 … … 610 672 ! 611 673 !------------------------------------------------------------------------------| 612 ! 9) tridiagonal system solving|674 ! 11) tridiagonal system solving | 613 675 !------------------------------------------------------------------------------| 614 676 ! … … 622 684 623 685 DO ji = kideb , kiut 624 z switbis(ji,numeqmin(ji)) = zswiterm(ji,numeqmin(ji))686 zindtbis(ji,numeqmin(ji)) = zindterm(ji,numeqmin(ji)) 625 687 zdiagbis(ji,numeqmin(ji)) = ztrid(ji,numeqmin(ji),2) 626 688 minnumeqmin = MIN(numeqmin(ji),minnumeqmin) … … 633 695 zdiagbis(ji,numeq) = ztrid(ji,numeq,2) - ztrid(ji,numeq,1)* & 634 696 ztrid(ji,numeq-1,3)/zdiagbis(ji,numeq-1) 635 z switbis(ji,numeq) = zswiterm(ji,numeq) - ztrid(ji,numeq,1)* &636 z switbis(ji,numeq-1)/zdiagbis(ji,numeq-1)697 zindtbis(ji,numeq) = zindterm(ji,numeq) - ztrid(ji,numeq,1)* & 698 zindtbis(ji,numeq-1)/zdiagbis(ji,numeq-1) 637 699 END DO 638 700 END DO … … 640 702 DO ji = kideb , kiut 641 703 ! ice temperatures 642 t_i_1d(ji,nlay_i) = z switbis(ji,numeqmax(ji))/zdiagbis(ji,numeqmax(ji))704 t_i_1d(ji,nlay_i) = zindtbis(ji,numeqmax(ji))/zdiagbis(ji,numeqmax(ji)) 643 705 END DO 644 706 … … 646 708 DO ji = kideb , kiut 647 709 jk = numeq - nlay_s - 1 648 t_i_1d(ji,jk) = (z switbis(ji,numeq) - ztrid(ji,numeq,3)* &710 t_i_1d(ji,jk) = (zindtbis(ji,numeq) - ztrid(ji,numeq,3)* & 649 711 t_i_1d(ji,jk+1))/zdiagbis(ji,numeq) 650 712 END DO … … 654 716 ! snow temperatures 655 717 IF (ht_s_1d(ji).GT.0._wp) & 656 t_s_1d(ji,nlay_s) = (z switbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) &718 t_s_1d(ji,nlay_s) = (zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) & 657 719 * t_i_1d(ji,1))/zdiagbis(ji,nlay_s+1) & 658 720 * MAX(0.0,SIGN(1.0,ht_s_1d(ji))) … … 662 724 ztsubit(ji) = t_su_1d(ji) 663 725 IF( t_su_1d(ji) < ztfs(ji) ) & 664 t_su_1d(ji) = ( z switbis(ji,numeqmin(ji)) - ztrid(ji,numeqmin(ji),3)* ( REAL( isnow(ji) )*t_s_1d(ji,1) &726 t_su_1d(ji) = ( zindtbis(ji,numeqmin(ji)) - ztrid(ji,numeqmin(ji),3)* ( REAL( isnow(ji) )*t_s_1d(ji,1) & 665 727 & + REAL( 1 - isnow(ji) )*t_i_1d(ji,1) ) ) / zdiagbis(ji,numeqmin(ji)) 666 728 END DO 667 729 ! 668 730 !-------------------------------------------------------------------------- 669 ! 1 0) Has the scheme converged ?, end of the iterative procedure |731 ! 12) Has the scheme converged ?, end of the iterative procedure | 670 732 !-------------------------------------------------------------------------- 671 733 ! … … 709 771 ! 710 772 !-------------------------------------------------------------------------! 711 ! 1 1) Fluxes at the interfaces !773 ! 13) Fluxes at the interfaces ! 712 774 !-------------------------------------------------------------------------! 713 775 DO ji = kideb, kiut … … 771 833 CALL wrk_dealloc( jpij, numeqmin, numeqmax, isnow ) 772 834 CALL wrk_dealloc( jpij, ztfs, ztsub, ztsubit, zh_i, zh_s, zfsw ) 773 CALL wrk_dealloc( jpij, zf, dzf, zerrit, zdifcase, zftrice, zihic )835 CALL wrk_dealloc( jpij, zf, dzf, zerrit, zdifcase, zftrice, zihic, zghe ) 774 836 CALL wrk_dealloc( jpij, nlay_i+1, ztcond_i, zradtr_i, zradab_i, zkappa_i, & 775 837 & ztib, zeta_i, ztitemp, z_i, zspeche_i, kjstart = 0 ) 776 838 CALL wrk_dealloc( jpij, nlay_s+1, zradtr_s, zradab_s, zkappa_s, ztsb, zeta_s, ztstemp, z_s, kjstart = 0 ) 777 CALL wrk_dealloc( jpij, nlay_i+3, z switerm, zswitbis, zdiagbis )839 CALL wrk_dealloc( jpij, nlay_i+3, zindterm, zindtbis, zdiagbis ) 778 840 CALL wrk_dealloc( jpij, nlay_i+3, 3, ztrid ) 779 841 CALL wrk_dealloc( jpij, zdq, zq_ini, zhfx_err ) -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90
r4990 r5048 420 420 END DO 421 421 ! ------------------------------------------------- 422 423 !-------------------------------------- 424 ! Impose a_i < amax in mono-category 425 !-------------------------------------- 426 ! 427 IF ( ( nn_monocat == 2 ) .AND. ( jpl == 1 ) ) THEN ! simple conservative piling, comparable with LIM2 428 DO jj = 1, jpj 429 DO ji = 1, jpi 430 a_i(ji,jj,1) = MIN( a_i(ji,jj,1), amax ) 431 END DO 432 END DO 433 ENDIF 422 434 423 435 ! --- diags --- -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90
r4990 r5048 200 200 201 201 ztmp = rday / rhoic 202 CALL iom_put( "vfxres" , wfx_res * ztmp ) ! daily prod./melting due to limupdate 203 CALL iom_put( "vfxopw" , wfx_opw * ztmp ) ! daily lateral thermodynamic ice production 204 CALL iom_put( "vfxsni" , wfx_sni * ztmp ) ! daily snowice ice production 205 CALL iom_put( "vfxbog" , wfx_bog * ztmp ) ! daily bottom thermodynamic ice production 206 CALL iom_put( "vfxdyn" , wfx_dyn * ztmp ) ! daily dynamic ice production (rid/raft) 207 CALL iom_put( "vfxsum" , wfx_sum * ztmp ) ! surface melt 208 CALL iom_put( "vfxbom" , wfx_bom * ztmp ) ! bottom melt 209 CALL iom_put( "vfxice" , wfx_ice * ztmp ) ! total ice growth/melt 210 CALL iom_put( "vfxsnw" , wfx_snw * ztmp ) ! total snw growth/melt 211 CALL iom_put( "vfxsub" , wfx_sub * ztmp ) ! sublimation (snow) 212 CALL iom_put( "vfxspr" , wfx_spr * ztmp ) ! precip (snow) 202 CALL iom_put( "vfxres" , wfx_res * ztmp ) ! daily prod./melting due to limupdate 203 CALL iom_put( "vfxopw" , wfx_opw * ztmp ) ! daily lateral thermodynamic ice production 204 CALL iom_put( "vfxsni" , wfx_sni * ztmp ) ! daily snowice ice production 205 CALL iom_put( "vfxbog" , wfx_bog * ztmp ) ! daily bottom thermodynamic ice production 206 CALL iom_put( "vfxdyn" , wfx_dyn * ztmp ) ! daily dynamic ice production (rid/raft) 207 CALL iom_put( "vfxsum" , wfx_sum * ztmp ) ! surface melt 208 CALL iom_put( "vfxbom" , wfx_bom * ztmp ) ! bottom melt 209 CALL iom_put( "vfxice" , wfx_ice * ztmp ) ! total ice growth/melt 210 CALL iom_put( "vfxsnw" , wfx_snw * ztmp ) ! total snw growth/melt 211 CALL iom_put( "vfxsub" , wfx_sub * ztmp ) ! sublimation (snow) 212 CALL iom_put( "vfxspr" , wfx_spr * ztmp ) ! precip (snow) 213 214 CALL iom_put( "afxtot" , afx_tot ) ! concentration tendency (total) 215 CALL iom_put( "afxdyn" , afx_dyn ) ! concentration tendency (dynamics) 216 CALL iom_put( "afxthd" , afx_thd ) ! concentration tendency (thermo) 213 217 214 218 CALL iom_put ('hfxthd', hfx_thd(:,:) ) ! … … 216 220 CALL iom_put ('hfxres', hfx_res(:,:) ) ! 217 221 CALL iom_put ('hfxout', hfx_out(:,:) ) ! 218 CALL iom_put ('hfxin' , hfx_in(:,:) ) !222 CALL iom_put ('hfxin' , hfx_in(:,:) ) ! 219 223 CALL iom_put ('hfxsnw', hfx_snw(:,:) ) ! 220 224 CALL iom_put ('hfxsub', hfx_sub(:,:) ) ! -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/par_ice.F90
r4873 r5048 12 12 13 13 ! !!! ice thermodynamics 14 INTEGER, PUBLIC, PARAMETER :: nlay_i = 5!: number of ice layers14 INTEGER, PUBLIC, PARAMETER :: nlay_i = 2 !: number of ice layers 15 15 INTEGER, PUBLIC, PARAMETER :: nlay_s = 1 !: number of snow layers 16 16 -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90
r5047 r5048 35 35 !: are the variables corresponding to 2d vectors 36 36 37 INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: npb !: number of points where computations has to be done 38 INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: npac !: correspondance between points (lateral accretion) 37 INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: npb !: address vector for 1d vertical thermo computations 38 INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: nplm !: address vector for mono-category lateral melting 39 INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: npac !: address vector for new ice formation 39 40 40 41 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: qlead_1d … … 109 110 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: dh_i_bott !: Ice bottom accretion/ablation [m] 110 111 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: dh_snowice !: Snow ice formation [m of ice] 112 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: dh_i_melt_1d !: Net melting [m], for mono-cat lateral melting 111 113 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: sm_i_1d !: Ice bulk salinity [ppt] 112 114 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: s_i_new !: Salinity of new ice at the bottom … … 138 140 !!---------------------------------------------------------------------! 139 141 140 ALLOCATE( npb (jpij) , np ac (jpij),&142 ALLOCATE( npb (jpij) , nplm (jpij) , npac (jpij), & 141 143 ! ! 142 144 & qlead_1d (jpij) , ftr_ice_1d (jpij) , & … … 165 167 & ht_s_1d (jpij) , fc_su (jpij) , fc_bo_i (jpij) , & 166 168 & dh_s_tot (jpij) , dh_i_surf(jpij) , dh_i_bott(jpij) , & 167 & dh_snowice(jpij) , sm_i_1d (jpij) , s_i_new (jpij) , & 169 & dh_snowice(jpij) , dh_i_melt_1d(jpij) , & 170 & sm_i_1d (jpij) , s_i_new (jpij) , & 168 171 & t_s_1d(jpij,nlay_s), & 169 172 & t_i_1d(jpij,nlay_i+1), s_i_1d(jpij,nlay_i+1) , & -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90
r4990 r5048 191 191 CALL sbc_cpl_ice_tau( utau_ice , vtau_ice ) 192 192 193 ! MV -> seb194 ! CALL sbc_cpl_ice_flx( p_frld=ato_i, palbi=zalb_ice, psst=sst_m, pist=t_su )195 196 ! IF( nn_limflx == 2 ) CALL ice_lim_flx( t_su, zalb_ice, qns_ice, qsr_ice , &197 ! & dqns_ice, qla_ice, dqla_ice, nn_limflx )198 ! ! Latent heat flux is forced to 0 in coupled :199 ! ! it is included in qns (non-solar heat flux)200 ! qla_ice (:,:,:) = 0._wp201 ! dqla_ice (:,:,:) = 0._wp202 ! END MV -> seb203 !204 193 END SELECT 205 194 … … 236 225 wfx_spr(:,:) = 0._wp ; 237 226 227 afx_tot(:,:) = at_i(:,:) ; afx_dyn(:,:) = 0._wp 228 afx_thd(:,:) = 0._wp 229 238 230 hfx_in (:,:) = 0._wp ; hfx_out(:,:) = 0._wp 239 231 hfx_thd(:,:) = 0._wp ; … … 252 244 ! ---------------------------------------------- 253 245 IF( .NOT. lk_c1d ) THEN 246 247 IF ( ln_limdyn ) afx_dyn(:,:) = at_i(:,:) 248 254 249 CALL lim_dyn( kt ) ! Ice dynamics ( rheology/dynamics ) 255 250 CALL lim_trp( kt ) ! Ice transport ( Advection/diffusion ) 256 251 CALL lim_var_glo2eqv ! equivalent variables, requested for rafting 257 252 IF( ln_nicep ) CALL lim_prt_state( kt, jiindx, jjindx,-1, ' - ice dyn & trp - ' ) ! control print 258 253 IF( nn_monocat /= 2 ) CALL lim_itd_me ! Mechanical redistribution ! (ridging/rafting) 259 254 CALL lim_var_agg( 1 ) 260 255 #if defined key_bdy … … 278 273 oa_i_b (:,:,:) = oa_i (:,:,:) 279 274 smv_i_b(:,:,:) = smv_i(:,:,:) 275 276 IF ( ln_limdyn ) afx_dyn(:,:) = ( at_i(:,:) - afx_dyn(:,:) ) * r1_rdtice 277 afx_thd(:,:) = at_i(:,:) 280 278 281 279 ! ---------------------------------------------- 282 ! ice thermodynamic 280 ! ice thermodynamics 283 281 ! ---------------------------------------------- 284 282 CALL lim_var_glo2eqv ! equivalent variables … … 288 286 phicif(:,:) = vt_i(:,:) 289 287 290 ! MV -> seb291 288 SELECT CASE( kblk ) 292 289 CASE ( jp_cpl ) … … 299 296 dqla_ice (:,:,:) = 0._wp 300 297 END SELECT 301 ! END MV -> seb302 298 ! 303 299 CALL lim_var_bv ! bulk brine volume (diag) … … 320 316 ! ! Diagnostics and outputs 321 317 IF (ln_limdiaout) CALL lim_diahsb 318 319 afx_thd(:,:) = ( at_i(:,:) - afx_thd(:,:) ) * r1_rdtice 320 afx_tot(:,:) = ( at_i(:,:) - afx_tot(:,:) ) * r1_rdtice 322 321 323 322 CALL lim_wri( 1 ) ! Ice outputs
Note: See TracChangeset
for help on using the changeset viewer.
