Changeset 8517 for branches/2017
- Timestamp:
- 2017-09-12T20:46:13+02:00 (7 years ago)
- Location:
- branches/2017/dev_r8183_ICEMODEL/NEMOGCM
- Files:
-
- 17 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/CONFIG/ORCA2_SAS_LIM3/EXP00/namelist_ice_cfg
r8321 r8517 1 1 !!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 2 !! LIM3 configuration namelist: Overwrites SHARED/namelist_ice_lim3_ref 3 !! 1 - Generic parameters (namicerun) 4 !! 2 - Diagnostics (namicediag) 5 !! 3 - Ice initialization (namiceini) 6 !! 4 - Ice discretization (namiceitd) 7 !! 5 - Ice dynamics and transport (namicedyn) 8 !! 6 - Ice thermodynamics (namicethd) 9 !! 7 - Ice salinity (namicesal) 10 !! 8 - Ice mechanical redistribution (namiceitdme) 11 !! 9 - Ice/snow albedos (namicealb) 12 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 2 !! ESIM configuration namelist: Overwrites SHARED/namelist_ice_lim3_ref 3 !! 1 - Generic parameters (namice_run) 4 !! 2 - Ice thickness discretization (namice_itd) 5 !! 3 - Ice dynamics (namice_dyn) 6 !! 4 - Ice ridging/rafting (namice_rdgrft) 7 !! 5 - Ice rheology (namice_rhg) 8 !! 6 - Ice advection (namice_adv) 9 !! 7 - Ice thermodynamics (namice_thd) 10 !! 8 - Ice salinity (namice_sal) 11 !! 9 - Ice melt ponds (namice_mp) 12 !! 10 - Ice initialization (namice_ini) 13 !! 11 - Ice/snow albedos (namice_alb) 14 !! 12 - Ice diagnostics (namice_dia) 15 !!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 16 ! 13 17 !------------------------------------------------------------------------------ 14 &namice run ! Generic parameters18 &namice_run ! Generic parameters 15 19 !------------------------------------------------------------------------------ 16 20 / 17 21 !------------------------------------------------------------------------------ 18 &namice diag ! Diagnostics22 &namice_itd ! Ice discretization 19 23 !------------------------------------------------------------------------------ 20 24 / 21 25 !------------------------------------------------------------------------------ 22 &namice ini ! Ice initialization26 &namice_dyn ! Ice dynamics 23 27 !------------------------------------------------------------------------------ 24 28 / 25 29 !------------------------------------------------------------------------------ 26 &namice itd ! Ice discretization30 &namice_rdgrft ! Ice ridging/rafting 27 31 !------------------------------------------------------------------------------ 28 32 / 29 33 !------------------------------------------------------------------------------ 30 &namice dyn ! Ice dynamics and transport34 &namice_rhg ! Ice rheology 31 35 !------------------------------------------------------------------------------ 32 36 / 33 37 !------------------------------------------------------------------------------ 34 &namice thd ! Ice thermodynamics38 &namice_adv ! Ice advection 35 39 !------------------------------------------------------------------------------ 36 40 / 37 41 !------------------------------------------------------------------------------ 38 &namice sal ! Ice salinity42 &namice_thd ! Ice thermodynamics 39 43 !------------------------------------------------------------------------------ 40 44 / 41 45 !------------------------------------------------------------------------------ 42 &namice itdme ! Ice mechanical redistribution (ridging and rafting)46 &namice_sal ! Ice salinity 43 47 !------------------------------------------------------------------------------ 44 48 / 45 !----------------------------------------------------------------------- 46 &namice alb ! albedo parameters47 !----------------------------------------------------------------------- 49 !------------------------------------------------------------------------------ 50 &namicemp ! Melt ponds 51 !------------------------------------------------------------------------------ 48 52 / 53 !------------------------------------------------------------------------------ 54 &namice_ini ! Ice initialization 55 !------------------------------------------------------------------------------ 56 / 57 !------------------------------------------------------------------------------ 58 &namice_alb ! albedo parameters 59 !------------------------------------------------------------------------------ 60 / 61 !------------------------------------------------------------------------------ 62 &namice_dia ! Diagnostics 63 !------------------------------------------------------------------------------ 64 / -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/CONFIG/SHARED/namelist_ice_lim3_ref
r8516 r8517 18 18 &namice_run ! Generic parameters 19 19 !------------------------------------------------------------------------------ 20 jpl = 5! number of ice categories21 nlay_i = 2! number of ice layers22 nlay_s = 1! number of snow layers (only 1 is working)23 nn_monocat = 0! virtual ITD mono-category parameterizations (1-4 => jpl = 1 only) or not (0)20 jpl = 5 ! number of ice categories 21 nlay_i = 2 ! number of ice layers 22 nlay_s = 1 ! number of snow layers (only 1 is working) 23 nn_monocat = 0 ! virtual ITD mono-category parameterizations (1-4 => jpl = 1 only) or not (0) 24 24 ! 2: simple piling instead of ridging --- temporary option 25 25 ! 3: activate G(he) only --- temporary option 26 26 ! 4: activate extra lateral melting only --- temporary option 27 ln_icedyn = .true. ! ice dynamics (T) or not (F) 28 ln_icethd = .true. ! ice thermo (T) or not (F) 27 29 rn_amax_n = 0.997 ! maximum tolerated ice concentration NH 28 30 rn_amax_s = 0.997 ! maximum tolerated ice concentration SH … … 41 43 &namice_dyn ! Ice dynamics 42 44 !------------------------------------------------------------------------------ 43 ln_icedyn = .true. ! ice dynamics (T) or not (F) => DO NOT TOUCH UNLESS U KNOW WHAT U DO 44 nn_icedyn = 2 ! switch for ice dynamics 45 ! 2: total 46 ! 1: advection only (no diffusion, no ridging/rafting) 47 ! 0: advection only (as 1 but with prescribed velocity, bypass rheology) 48 rn_uice = 0.00001 ! prescribed ice u-velocity 49 rn_vice = -0.00001 ! prescribed ice v-velocity 45 ln_dynFULL = .true. ! dyn.: full ice dynamics (rheology + advection + ridging/rafting + correction) 46 ln_dynRHGADV = .false. ! dyn.: no ridge/raft & no corrections (rheology + advection) 47 ln_dynADV = .false. ! dyn.: only advection w prescribed vel.(rn_uvice + advection) 48 rn_uice = 0.00001 ! prescribed ice u-velocity 49 rn_vice = 0. ! prescribed ice v-velocity 50 50 rn_ishlat = 2. ! free slip (0) ; partial slip (0-2) ; no slip (2) ; strong slip (>2) 51 51 rn_cio = 5.0e-03 ! ice-ocean drag coefficient (-) … … 53 53 rn_gamma = 0.15 ! fraction of ocean depth that ice must reach to initiate landfast 54 54 ! recommended range: [0.1 ; 0.25] 55 rn_icebfr = 10. ! maximum bottom stress per unit area of contact (N/m2)55 rn_icebfr = 10. ! maximum bottom stress per unit area of contact [N/m2] 56 56 ! a very large value ensures ice velocity=0 even with a small contact area 57 57 ! recommended range: ?? (should be greater than atm-ice stress => >0.1 N/m2) 58 rn_lfrelax = 1.e-5 ! relaxation time scale to reach static friction (s-1)58 rn_lfrelax = 1.e-5 ! relaxation time scale to reach static friction [s-1] 59 59 / 60 60 !------------------------------------------------------------------------------ … … 63 63 ! -- ice_rdgrft_strength -- ! 64 64 ln_str_H79 = .true. ! ice strength param.: Hibler_79 => P = pstar*<h>*exp(-c_rhg*A) 65 rn_pstar = 2.0e+04 ! ice strength thickness parameter (N/m2)65 rn_pstar = 2.0e+04 ! ice strength thickness parameter [N/m2] 66 66 rn_crhg = 20.0 ! ice strength conc. parameter (-) 67 67 ln_str_R75 = .false. ! ice strength param.: Rothrock_75 => P = Cf*coeff*integral(wr.h^2) … … 74 74 ln_partf_exp = .true. ! Exponential ridging participation function (Lipscomb, 2007) 75 75 rn_astar = 0.03 ! exponential measure of ridging ice fraction [set to 0.05 if hstar=100] 76 ln_ridging = .true. ! ridging activated (T) or not (F) => DO NOT TOUCH UNLESS U KNOW WHAT U DO77 rn_hstar = 25.0 ! determines the maximum thickness of ridged ice (m)(Hibler, 1980)76 ln_ridging = .true. ! ridging activated (T) or not (F) 77 rn_hstar = 25.0 ! determines the maximum thickness of ridged ice [m] (Hibler, 1980) 78 78 rn_porordg = 0.3 ! porosity of newly ridged ice (Lepparanta et al., 1995) 79 79 rn_fsnwrdg = 0.5 ! snow volume fraction that survives in ridging 80 80 rn_fpndrdg = 1.0 ! pond fraction that survives in ridging (small a priori) 81 ln_rafting = .true. ! rafting activated (T) or not (F) => DO NOT TOUCH UNLESS U KNOW WHAT U DO82 rn_hraft = 0.75 ! threshold thickness for rafting (m)81 ln_rafting = .true. ! rafting activated (T) or not (F) 82 rn_hraft = 0.75 ! threshold thickness for rafting [m] 83 83 rn_craft = 5.0 ! squeezing coefficient used in the rafting function 84 84 rn_fsnwrft = 0.5 ! snow volume fraction that survives in rafting … … 89 89 !------------------------------------------------------------------------------ 90 90 ln_rhg_EVP = .true. ! EVP rheology 91 rn_creepl = 1.0e-12 ! creep limit (s-1)91 rn_creepl = 1.0e-12 ! creep limit [1/s] 92 92 rn_ecc = 2.0 ! eccentricity of the elliptical yield curve 93 93 nn_nevp = 120 ! number of EVP subcycles … … 105 105 &namice_thd ! Ice thermodynamics 106 106 !------------------------------------------------------------------------------ 107 ln_icethd = .true. ! ice thermo (T) or not (F) => DO NOT TOUCH UNLESS U KNOW WHAT U DO108 107 ! -- icethd_dif -- ! 109 rn_kappa_i = 1.0 ! radiation attenuation coefficient in sea ice (m-1)108 rn_kappa_i = 1.0 ! radiation attenuation coefficient in sea ice [1/m] 110 109 ln_cndi_U64 = .false. ! sea ice thermal conductivity: k = k0 + beta.S/T (Untersteiner, 1964) 111 110 ln_cndi_P07 = .true. ! sea ice thermal conductivity: k = k0 + beta1.S/T - beta2.T (Pringle et al., 2007) … … 114 113 ! Obs: 0.1-0.5 (Lecomte et al, JAMES 2013) 115 114 ! -- icethd_dh -- ! 116 ln_icedH = .true. ! activate ice thickness change from growing/melting (T) or not (F) => DO NOT TOUCH UNLESS U KNOW WHAT U DO115 ln_icedH = .true. ! activate ice thickness change from growing/melting (T) or not (F) 117 116 rn_blow_s = 0.66 ! mesure of snow blowing into the leads 118 117 ! = 1 => no snow blowing, < 1 => some snow blowing 119 118 ! -- icethd_da -- ! 120 ln_icedA = .true. ! activate lateral melting param. (T) or not (F) => DO NOT TOUCH UNLESS U KNOW WHAT U DO119 ln_icedA = .true. ! activate lateral melting param. (T) or not (F) 121 120 rn_beta = 1.0 ! coef. beta for lateral melting param. Recommended range=[0.8-1.2] 122 121 ! => decrease = more melt and melt peaks toward higher concentration (A~0.5 for beta=1 ; A~0.8 for beta=0.2) … … 127 126 ! 10 vs 8m = -20% melting 128 127 ! -- icethd_lac -- ! 129 ln_icedO = .true. ! activate ice growth in open-water (T) or not (F) => DO NOT TOUCH UNLESS U KNOW WHAT U DO128 ln_icedO = .true. ! activate ice growth in open-water (T) or not (F) 130 129 rn_hinew = 0.1 ! thickness for new ice formation in open water (m), must be larger than rn_hnewice 131 130 ln_frazil = .false. ! Frazil ice parameterization (ice collection as a function of wind) … … 148 147 !------------------------------------------------------------------------------ 149 148 ! -- icethd_sal -- ! 150 ln_icedS = .true. ! activate gravity drainage and flushing (T) or not (F) => DO NOT TOUCH UNLESS U KNOW WHAT U DO149 ln_icedS = .true. ! activate gravity drainage and flushing (T) or not (F) 151 150 nn_icesal = 2 ! ice salinity option 152 151 ! 1: constant ice salinity (S=rn_icesal) -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/CONFIG/SHARED/namelist_ref
r8321 r8517 429 429 / 430 430 !----------------------------------------------------------------------- 431 &namsbc_alb ! albedo parameters432 !-----------------------------------------------------------------------433 nn_ice_alb = 1 ! parameterization of ice/snow albedo434 ! 0: Shine & Henderson-Sellers (JGR 1985), giving clear-sky albedo435 ! 1: "home made" based on Brandt et al. (JClim 2005) and Grenfell & Perovich (JGR 2004),436 ! giving cloud-sky albedo437 rn_alb_sdry = 0.85 ! dry snow albedo : 0.80 (nn_ice_alb = 0); 0.85 (nn_ice_alb = 1); obs 0.85-0.87 (cloud-sky)438 rn_alb_smlt = 0.75 ! melting snow albedo : 0.65 ( '' ) ; 0.75 ( '' ) ; obs 0.72-0.82 ( '' )439 rn_alb_idry = 0.60 ! dry ice albedo : 0.72 ( '' ) ; 0.60 ( '' ) ; obs 0.54-0.65 ( '' )440 rn_alb_imlt = 0.50 ! bare puddled ice albedo : 0.53 ( '' ) ; 0.50 ( '' ) ; obs 0.49-0.58 ( '' )441 rn_alb_dpnd = 0.27 ! ponded ice albedo : 0.25 ( '' ) ; 0.27 ( '' ) ; obs 0.10-0.30 ( '' )442 /443 !-----------------------------------------------------------------------444 431 &namsbc_wave ! External fields from wave model (ln_wave=T) 445 432 !----------------------------------------------------------------------- -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/ice.F90
r8516 r8517 157 157 INTEGER , PUBLIC :: nlay_i !: number of ice layers 158 158 INTEGER , PUBLIC :: nlay_s !: number of snow layers 159 LOGICAL , PUBLIC :: ln_icedyn !: flag for ice dynamics (T) or not (F) 160 LOGICAL , PUBLIC :: ln_icethd !: flag for ice thermo (T) or not (F) 159 161 REAL(wp) , PUBLIC :: rn_amax_n !: maximum ice concentration Northern hemisphere 160 162 REAL(wp) , PUBLIC :: rn_amax_s !: maximum ice concentration Southern hemisphere … … 168 170 169 171 ! !!** ice-dynamics namelist (namice_dyn) ** 170 LOGICAL , PUBLIC :: ln_icedyn !: flag for ice dynamics (T) or not (F) 171 INTEGER , PUBLIC :: nn_icedyn !: flag for ice dynamics 172 REAL(wp), PUBLIC :: rn_uice !: prescribed u-vel (case nn_icedyn=0) 173 REAL(wp), PUBLIC :: rn_vice !: prescribed v-vel (case nn_icedyn=0) 172 LOGICAL, PUBLIC :: ln_dynFULL !: dyn.: full ice dynamics (rheology + advection + ridging/rafting + correction) 173 LOGICAL, PUBLIC :: ln_dynRHGADV !: dyn.: no ridge/raft & no corrections (rheology + advection) 174 LOGICAL, PUBLIC :: ln_dynADV !: dyn.: only advection w prescribed vel.(rn_uvice + advection) 174 175 REAL(wp), PUBLIC :: rn_ishlat !: lateral boundary condition for sea-ice 175 176 REAL(wp), PUBLIC :: rn_cio !: drag coefficient for oceanic stress 176 177 LOGICAL , PUBLIC :: ln_landfast !: landfast ice parameterization (T or F) 177 REAL(wp), PUBLIC :: rn_gamma !: fraction of ocean depth that ice must reach to initiate landfast ice178 REAL(wp), PUBLIC :: rn_icebfr !: maximum bottom stress per unit area of contact (landfast ice)179 REAL(wp), PUBLIC :: rn_lfrelax !: relaxation time scale (s-1) to reach static friction (landfast ice)178 REAL(wp), PUBLIC :: rn_gamma !: fraction of ocean depth that ice must reach to initiate landfast ice 179 REAL(wp), PUBLIC :: rn_icebfr !: maximum bottom stress per unit area of contact (landfast ice) 180 REAL(wp), PUBLIC :: rn_lfrelax !: relaxation time scale (s-1) to reach static friction (landfast ice) 180 181 ! 181 182 ! !!** ice-rdige/raft namelist (namice_rdgrft) ** … … 199 200 ! 200 201 ! !!** ice-thermodynamics namelist (namice_thd) ** 201 LOGICAL , PUBLIC :: ln_icethd !: flag for ice thermo (T) or not (F)202 202 ! -- icethd_dif -- ! 203 203 REAL(wp), PUBLIC :: rn_kappa_i !: coef. for the extinction of radiation Grenfell et al. (2006) [1/m] … … 266 266 REAL(wp), PUBLIC :: r1_nlay_s !: 1 / nlay_s 267 267 REAL(wp), PUBLIC :: rswitch !: switch for the presence of ice (1) or not (0) 268 REAL(wp), PUBLIC :: rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft !: conservation diagnostics 268 269 REAL(wp), PUBLIC, PARAMETER :: epsi06 = 1.e-06_wp !: small number 269 270 REAL(wp), PUBLIC, PARAMETER :: epsi10 = 1.e-10_wp !: small number 270 271 REAL(wp), PUBLIC, PARAMETER :: epsi20 = 1.e-20_wp !: small number 271 ! 272 LOGICAL , PUBLIC :: l_piling !: =T simple conservative piling, comparable with LIM2 272 273 273 274 274 ! !!** define arrays -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceadv.F90
r8516 r8517 37 37 PUBLIC ice_adv_init ! called by icestp 38 38 39 INTEGER :: nice_adv ! choice of the type of advection scheme 40 ! ! associated indices: 41 INTEGER, PARAMETER :: np_advPRA = 1 ! Prather scheme 42 INTEGER, PARAMETER :: np_advUMx = 2 ! Ultimate-Macho scheme 43 39 44 !! * Substitution 40 45 # include "vectopt_loop_substitute.h90" … … 61 66 !!---------------------------------------------------------------------- 62 67 INTEGER, INTENT(in) :: kt ! number of iteration 63 !64 INTEGER :: ji, jj, jk, jl, jt ! dummy loop indices65 INTEGER :: initad ! number of sub-timestep for the advection66 REAL(wp) :: zcfl , zusnit ! - -67 CHARACTER(len=80) :: cltmp68 !69 REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b70 REAL(wp) :: zdv71 REAL(wp), DIMENSION(jpi,jpj) :: zatold, zeiold, zesold, zsmvold72 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zhimax, zviold, zvsold73 68 !!--------------------------------------------------------------------- 74 69 IF( nn_timing == 1 ) CALL timing_start('iceadv') … … 80 75 ENDIF 81 76 82 CALL ice_var_agg( 1 ) ! integrated values + ato_i83 84 77 ! conservation test 85 IF( ln_icediachk ) CALL ice_cons_hsm(0, 'iceadv', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 86 87 ! store old values for diag 88 zviold (:,:,:) = v_i(:,:,:) 89 zvsold (:,:,:) = v_s(:,:,:) 90 zsmvold(:,:) = SUM( smv_i(:,:,:), dim=3 ) 91 zeiold (:,:) = et_i(:,:) 92 zesold (:,:) = et_s(:,:) 93 94 ! Thickness correction init. 95 zatold(:,:) = at_i 96 WHERE( a_i(:,:,:) >= epsi20 ) 97 ht_i(:,:,:) = v_i(:,:,:) / a_i(:,:,:) 98 ht_s(:,:,:) = v_s(:,:,:) / a_i(:,:,:) 99 ELSEWHERE 100 ht_i(:,:,:) = 0._wp 101 ht_s(:,:,:) = 0._wp 102 END WHERE 103 104 ! Record max of the surrounding ice thicknesses for correction in case advection creates ice too thick 105 zhimax(:,:,:) = ht_i(:,:,:) + ht_s(:,:,:) 106 DO jl = 1, jpl 107 DO jj = 2, jpjm1 108 DO ji = 2, jpim1 109 !!gm use of MAXVAL here is very probably less efficient than expending the 9 values 110 zhimax(ji,jj,jl) = MAX( epsi20, MAXVAL( ht_i(ji-1:ji+1,jj-1:jj+1,jl) + ht_s(ji-1:ji+1,jj-1:jj+1,jl) ) ) 111 END DO 112 END DO 113 END DO 114 CALL lbc_lnk( zhimax(:,:,:), 'T', 1. ) 115 78 IF( ln_icediachk ) CALL ice_cons_hsm(0, 'iceadv', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) 79 116 80 !---------- 117 81 ! Advection 118 82 !---------- 119 IF( ln_adv_UMx ) THEN !-- ULTIMATE-MACHO scheme 83 SELECT CASE( nice_adv ) 84 ! !-----------------------! 85 CASE( np_advUMx ) ! ULTIMATE-MACHO scheme ! 86 ! !-----------------------! 120 87 CALL ice_adv_umx( kt, u_ice, v_ice, & 121 88 & ato_i, v_i, v_s, smv_i, oa_i, a_i, a_ip, v_ip, e_s, e_i ) 122 123 ELSEIF( ln_adv_Pra ) THEN !-- PRATHER scheme 89 ! !-----------------------! 90 CASE( np_advPRA ) ! PRATHER scheme ! 91 ! !-----------------------! 124 92 CALL ice_adv_prather( kt, u_ice, v_ice, & 125 93 & ato_i, v_i, v_s, smv_i, oa_i, a_i, a_ip, v_ip, e_s, e_i ) 126 127 ENDIF 94 END SELECT 128 95 129 ! total ice fraction130 at_i(:,:) = a_i(:,:,1)131 DO jl = 2, jpl132 at_i(:,:) = at_i(:,:) + a_i(:,:,jl)133 END DO134 135 96 !------------ 136 97 ! diagnostics 137 98 !------------ 138 DO jj = 1, jpj 139 DO ji = 1, jpi 140 diag_trp_ei (ji,jj) = ( SUM( e_i (ji,jj,1:nlay_i,:) ) - zeiold(ji,jj) ) * r1_rdtice 141 diag_trp_es (ji,jj) = ( SUM( e_s (ji,jj,1:nlay_s,:) ) - zesold(ji,jj) ) * r1_rdtice 142 diag_trp_smv(ji,jj) = ( SUM( smv_i(ji,jj,:) ) - zsmvold(ji,jj) ) * r1_rdtice 143 diag_trp_vi (ji,jj) = SUM( v_i(ji,jj,:) - zviold(ji,jj,:) ) * r1_rdtice 144 diag_trp_vs (ji,jj) = SUM( v_s(ji,jj,:) - zvsold(ji,jj,:) ) * r1_rdtice 145 END DO 146 END DO 99 diag_trp_ei (:,:) = SUM(SUM( e_i (:,:,1:nlay_i,:) - e_i_b (:,:,1:nlay_i,:), dim=4 ), dim=3 ) * r1_rdtice 100 diag_trp_es (:,:) = SUM(SUM( e_s (:,:,1:nlay_s,:) - e_s_b (:,:,1:nlay_s,:), dim=4 ), dim=3 ) * r1_rdtice 101 diag_trp_smv(:,:) = SUM( smv_i(:,:,:) - smv_i_b(:,:,:) , dim=3 ) * r1_rdtice 102 diag_trp_vi (:,:) = SUM( v_i (:,:,:) - v_i_b (:,:,:) , dim=3 ) * r1_rdtice 103 diag_trp_vs (:,:) = SUM( v_s (:,:,:) - v_s_b (:,:,:) , dim=3 ) * r1_rdtice 147 104 IF( iom_use('icetrp') ) CALL iom_put( "icetrp" , diag_trp_vi * rday ) ! ice volume transport 148 105 IF( iom_use('snwtrp') ) CALL iom_put( "snwtrp" , diag_trp_vs * rday ) ! snw volume transport … … 150 107 IF( iom_use('deitrp') ) CALL iom_put( "deitrp" , diag_trp_ei ) ! advected ice enthalpy (W/m2) 151 108 IF( iom_use('destrp') ) CALL iom_put( "destrp" , diag_trp_es ) ! advected snw enthalpy (W/m2) 152 153 !-------------------------------------- 154 ! Thickness correction in case too high 155 !-------------------------------------- 156 IF( nn_icedyn == 2 ) THEN 157 ! 158 CALL ice_var_zapsmall !-- zap small areas 159 ! 160 DO jl = 1, jpl 161 DO jj = 1, jpj 162 DO ji = 1, jpi 163 IF ( v_i(ji,jj,jl) > 0._wp ) THEN !-- bound to zhimax 164 ! 165 ht_i (ji,jj,jl) = v_i (ji,jj,jl) / a_i(ji,jj,jl) 166 ht_s (ji,jj,jl) = v_s (ji,jj,jl) / a_i(ji,jj,jl) 167 zdv = v_i(ji,jj,jl) + v_s(ji,jj,jl) - zviold(ji,jj,jl) - zvsold(ji,jj,jl) 168 ! 169 IF ( ( zdv > 0.0 .AND. (ht_i(ji,jj,jl)+ht_s(ji,jj,jl)) > zhimax(ji,jj,jl) .AND. zatold(ji,jj) < 0.80 ) .OR. & 170 & ( zdv <= 0.0 .AND. (ht_i(ji,jj,jl)+ht_s(ji,jj,jl)) > zhimax(ji,jj,jl) ) ) THEN 171 a_i (ji,jj,jl) = ( v_i(ji,jj,jl) + v_s(ji,jj,jl) ) / zhimax(ji,jj,jl) 172 ht_i(ji,jj,jl) = v_i(ji,jj,jl) / a_i(ji,jj,jl) 173 ENDIF 174 ! 175 ENDIF 176 END DO 177 END DO 178 END DO 179 180 WHERE( ht_i(:,:,jpl) > hi_max(jpl) ) !-- bound ht_i to hi_max (99 m) 181 ht_i(:,:,jpl) = hi_max(jpl) 182 a_i (:,:,jpl) = v_i(:,:,jpl) / hi_max(jpl) 183 END WHERE 184 185 IF ( nn_pnd_scheme > 0 ) THEN !-- correct pond fraction to avoid a_ip > a_i 186 WHERE( a_ip(:,:,:) > a_i(:,:,:) ) a_ip(:,:,:) = a_i(:,:,:) 187 ENDIF 188 ! 189 ENDIF 190 191 !------------------------------------------------------------ 192 ! Impose a_i < amax if no ridging/rafting or in mono-category 193 !------------------------------------------------------------ 194 IF( l_piling ) THEN !-- simple conservative piling, comparable with 1-cat models 195 at_i(:,:) = SUM( a_i(:,:,:), dim=3 ) 196 DO jl = 1, jpl 197 WHERE( at_i(:,:) > epsi20 ) 198 a_i(:,:,jl) = a_i(:,:,jl) * ( 1._wp + MIN( rn_amax_2d(:,:) - at_i(:,:) , 0._wp ) / at_i(:,:) ) 199 END WHERE 200 END DO 201 ENDIF 202 203 ! agglomerate variables 204 vt_i(:,:) = SUM( v_i(:,:,:), dim=3 ) 205 vt_s(:,:) = SUM( v_s(:,:,:), dim=3 ) 206 at_i(:,:) = SUM( a_i(:,:,:), dim=3 ) 207 208 ! MV MP 2016 (remove once we get rid of a_i_frac and ht_i) 209 IF ( nn_pnd_scheme > 0 ) THEN 210 at_ip(:,:) = SUM( a_ip(:,:,:), dim = 3 ) 211 vt_ip(:,:) = SUM( v_ip(:,:,:), dim = 3 ) 212 ENDIF 213 ! END MP 2016 214 215 ! open water = 1 if at_i=0 216 WHERE( at_i == 0._wp ) ato_i = 1._wp 217 218 ! conservation test 219 IF( ln_icediachk ) CALL ice_cons_hsm(1, 'iceadv', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 220 221 ! -------------- 222 ! control prints 223 ! -------------- 224 IF( ln_icectl ) CALL ice_prt( kt, iiceprt, jiceprt,-1, ' - ice dyn & trp - ' ) 225 ! 226 IF( nn_timing == 1 ) CALL timing_stop('iceadv') 109 110 ! controls 111 IF( ln_icediachk ) CALL ice_cons_hsm(1, 'iceadv', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 112 IF( ln_icectl ) CALL ice_prt (kt, iiceprt, jiceprt,-1, ' - ice dyn & trp - ') ! prints 113 IF( nn_timing == 1 ) CALL timing_stop ('iceadv') ! timing 227 114 ! 228 115 END SUBROUTINE ice_adv … … 241 128 !! ** input : Namelist namice_adv 242 129 !!------------------------------------------------------------------- 243 INTEGER :: ios ! Local integer output status for namelist read130 INTEGER :: ios, ioptio ! Local integer output status for namelist read 244 131 !! 245 132 NAMELIST/namice_adv/ ln_adv_Pra, ln_adv_UMx, nn_UMx … … 260 147 WRITE(numout,*) '~~~~~~~~~~~~' 261 148 WRITE(numout,*) ' Namelist namice_adv' 262 WRITE(numout,*) ' advection scheme for ice transport (limtrp)' 263 WRITE(numout,*) ' type of advection scheme (Prather) ln_adv_Pra = ', ln_adv_Pra 264 WRITE(numout,*) ' type of advection scheme (Ulimate-Macho) ln_adv_UMx = ', ln_adv_UMx 265 WRITE(numout,*) ' order of the Ultimate-Macho scheme nn_UMx = ', nn_UMx 149 WRITE(numout,*) ' type of advection scheme (Prather) ln_adv_Pra = ', ln_adv_Pra 150 WRITE(numout,*) ' type of advection scheme (Ulimate-Macho) ln_adv_UMx = ', ln_adv_UMx 151 WRITE(numout,*) ' order of the Ultimate-Macho scheme nn_UMx = ', nn_UMx 266 152 ENDIF 267 153 ! 268 IF ( ( ln_adv_Pra .AND. ln_adv_UMx ) .OR. ( .NOT.ln_adv_Pra .AND. .NOT.ln_adv_UMx ) ) THEN 269 CALL ctl_stop( 'ice_adv_init: choose one and only one ice advection scheme (ln_adv_Pra or ln_adv_UMx)' ) 270 ENDIF 154 ! !== set the choice of ice advection ==! 155 ioptio = 0 156 IF( ln_adv_Pra ) THEN ; ioptio = ioptio + 1 ; nice_adv = np_advPRA ; ENDIF 157 IF( ln_adv_UMx ) THEN ; ioptio = ioptio + 1 ; nice_adv = np_advUMx ; ENDIF 158 IF( ioptio /= 1 ) CALL ctl_stop( 'ice_adv_init: choose one and only one ice advection scheme (ln_adv_Pra or ln_adv_UMx)' ) 271 159 ! 272 160 END SUBROUTINE ice_adv_init -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icealb.F90
r8514 r8517 321 321 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namice_alb in reference namelist', lwp ) 322 322 323 REWIND( numnam_ice_cfg ) ! Namelist nam sbc_alb in configuration namelist : Albedo parameters323 REWIND( numnam_ice_cfg ) ! Namelist namice_alb in configuration namelist : Albedo parameters 324 324 READ ( numnam_ice_cfg, namice_alb, IOSTAT = ios, ERR = 902 ) 325 325 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namice_alb in configuration namelist', lwp ) -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icecor.F90
r8514 r8517 53 53 ! 54 54 INTEGER :: ji, jj, jk, jl ! dummy loop indices 55 REAL(wp) :: zsal, z vi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b, zzc55 REAL(wp) :: zsal, zzc 56 56 REAL(wp), DIMENSION(jpi,jpj) :: zafx ! concentration trends diag 57 57 !!---------------------------------------------------------------------- 58 IF( nn_timing == 1 ) CALL timing_start('icecor') 58 ! controls 59 IF( nn_timing == 1 ) CALL timing_start('icecor') ! timing 60 IF( ln_icediachk ) CALL ice_cons_hsm(0, 'icecor', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 59 61 ! 60 62 IF( kt == nit000 .AND. lwp .AND. kn == 2 ) THEN … … 63 65 WRITE(numout,*) '~~~~~~~' 64 66 ENDIF 65 ! !--- conservation test66 IF( ln_icediachk ) CALL ice_cons_hsm(0, 'icecor', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)67 !68 67 ! 69 68 IF( kn == 2 ) THEN … … 175 174 END SELECT 176 175 ! 177 ! !--- conservation test 178 IF( ln_icediachk ) CALL ice_cons_hsm(1, 'icecor', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 179 ! 180 ! !--- control prints 181 IF( ln_ctl ) CALL ice_prt3D( 'icecor' ) 182 IF( ln_icectl .AND. kn == 2 ) CALL ice_prt( kt, iiceprt, jiceprt, 2, ' - Final state - ' ) 183 ! 184 IF( nn_timing == 1 ) CALL timing_stop('icecor') 176 ! controls 177 IF( ln_icediachk ) CALL ice_cons_hsm(1, 'icecor', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 178 IF( ln_ctl ) CALL ice_prt3D ('icecor') ! prints 179 IF( ln_icectl .AND. kn == 2 ) CALL ice_prt( kt, iiceprt, jiceprt, 2, ' - Final state - ' ) ! prints 180 IF( nn_timing == 1 ) CALL timing_stop ('icecor') ! timing 185 181 ! 186 182 END SUBROUTINE ice_cor -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icectl.F90
r8514 r8517 47 47 CONTAINS 48 48 49 SUBROUTINE ice_cons_hsm( icount, cd_routine, zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)49 SUBROUTINE ice_cons_hsm( icount, cd_routine, pdiag_v, pdiag_s, pdiag_t, pdiag_fv, pdiag_fs, pdiag_ft ) 50 50 !!---------------------------------------------------------------------- 51 51 !! *** ROUTINE ice_cons_hsm *** … … 56 56 !! ** Method : This is an online diagnostics which can be activated with ln_icediachk=true 57 57 !! It prints in ocean.output if there is a violation of conservation at each time-step 58 !! The thresholds (zv_sill, zs_sill, z h_sill) which determine violations are set to58 !! The thresholds (zv_sill, zs_sill, zt_sill) which determine violations are set to 59 59 !! a minimum of 1 mm of ice (over the ice area) that is lost/gained spuriously during 100 years. 60 60 !! For salt and heat thresholds, ice is considered to have a salinity of 10 … … 63 63 INTEGER , INTENT(in) :: icount ! called at: =0 the begining of the routine, =1 the end 64 64 CHARACTER(len=*), INTENT(in) :: cd_routine ! name of the routine 65 REAL(wp) , INTENT(inout) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b ! ????66 !! 67 REAL(wp) :: zvi, zsmv, zei, zfs, zfw,zft68 REAL(wp) ::zvmin, zamin, zamax69 REAL(wp) ::zvtrp, zetrp70 REAL(wp) :: zarea, zv_sill, zs_sill, zh_sill71 REAL(wp), PARAMETER ::zconv = 1.e-9 ! convert W to GW and kg to Mt65 REAL(wp) , INTENT(inout) :: pdiag_v, pdiag_s, pdiag_t, pdiag_fv, pdiag_fs, pdiag_ft 66 !! 67 REAL(wp) :: zv, zs, zt, zfs, zfv, zft 68 REAL(wp) :: zvmin, zamin, zamax 69 REAL(wp) :: zvtrp, zetrp 70 REAL(wp) :: zarea, zv_sill, zs_sill, zt_sill 71 REAL(wp), PARAMETER :: zconv = 1.e-9 ! convert W to GW and kg to Mt 72 72 !!---------------------------------------------------------------------- 73 73 ! 74 74 IF( icount == 0 ) THEN 75 ! ! water flux 76 pdiag_fv = glob_sum( -( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + & 77 & wfx_opw(:,:) + wfx_res(:,:) + wfx_dyn(:,:) + wfx_lam(:,:) + wfx_ice_sub(:,:) + & 78 & wfx_snw_sni(:,:) + wfx_snw_sum(:,:) + wfx_snw_dyn(:,:) + wfx_snw_sub(:,:) + wfx_spr(:,:) & 79 & ) * e1e2t(:,:) ) * zconv 80 ! 75 81 ! ! salt flux 76 zfs_b = glob_sum( ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + & 77 & sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) + sfx_sub(:,:) + sfx_lam(:,:) & 78 & ) * e1e2t(:,:) ) * zconv 79 ! 80 ! ! water flux 81 zfw_b = glob_sum( -( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + & 82 & wfx_opw(:,:) + wfx_res(:,:) + wfx_dyn(:,:) + wfx_lam(:,:) + wfx_ice_sub(:,:) + & 83 & wfx_snw_sni(:,:) + wfx_snw_sum(:,:) + wfx_snw_dyn(:,:) + wfx_snw_sub(:,:) + wfx_spr(:,:) & 84 & ) * e1e2t(:,:) ) * zconv 82 pdiag_fs = glob_sum( ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + & 83 & sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) + sfx_sub(:,:) + sfx_lam(:,:) & 84 & ) * e1e2t(:,:) ) * zconv 85 85 ! 86 86 ! ! heat flux 87 zft_b= glob_sum( ( hfx_sum(:,:) + hfx_bom(:,:) + hfx_bog(:,:) + hfx_dif(:,:) + hfx_opw(:,:) + hfx_snw(:,:) &88 & - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) - hfx_sub(:,:) - hfx_spr(:,:) &89 & ) * e1e2t(:,:) ) * zconv90 91 zvi_b= glob_sum( SUM( v_i * rhoic + v_s * rhosn, dim=3 ) * e1e2t * zconv )92 93 zsmv_b= glob_sum( SUM( smv_i * rhoic , dim=3 ) * e1e2t * zconv )94 95 zei_b= glob_sum( ( SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 ) &96 & + SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 ) ) * e1e2t ) * zconv87 pdiag_ft = glob_sum( ( hfx_sum(:,:) + hfx_bom(:,:) + hfx_bog(:,:) + hfx_dif(:,:) + hfx_opw(:,:) + hfx_snw(:,:) & 88 & - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) - hfx_sub(:,:) - hfx_spr(:,:) & 89 & ) * e1e2t(:,:) ) * zconv 90 91 pdiag_v = glob_sum( SUM( v_i * rhoic + v_s * rhosn, dim=3 ) * e1e2t * zconv ) 92 93 pdiag_s = glob_sum( SUM( smv_i * rhoic , dim=3 ) * e1e2t * zconv ) 94 95 pdiag_t = glob_sum( ( SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 ) & 96 & + SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 ) ) * e1e2t ) * zconv 97 97 98 98 ELSEIF( icount == 1 ) THEN 99 100 ! water flux 101 zfv = glob_sum( -( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + & 102 & wfx_opw(:,:) + wfx_res(:,:) + wfx_dyn(:,:) + wfx_lam(:,:) + wfx_ice_sub(:,:) + & 103 & wfx_snw_sni(:,:) + wfx_snw_sum(:,:) + wfx_snw_dyn(:,:) + wfx_snw_sub(:,:) + wfx_spr(:,:) & 104 & ) * e1e2t(:,:) ) * zconv - pdiag_fv 99 105 100 106 ! salt flux 101 107 zfs = glob_sum( ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + & 102 108 & sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) + sfx_sub(:,:) + sfx_lam(:,:) & 103 & ) * e1e2t(:,:) ) * zconv - zfs_b 104 105 ! water flux 106 zfw = glob_sum( -( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + & 107 & wfx_opw(:,:) + wfx_res(:,:) + wfx_dyn(:,:) + wfx_lam(:,:) + wfx_ice_sub(:,:) + & 108 & wfx_snw_sni(:,:) + wfx_snw_sum(:,:) + wfx_snw_dyn(:,:) + wfx_snw_sub(:,:) + wfx_spr(:,:) & 109 & ) * e1e2t(:,:) ) * zconv - zfw_b 109 & ) * e1e2t(:,:) ) * zconv - pdiag_fs 110 110 111 111 ! heat flux 112 112 zft = glob_sum( ( hfx_sum(:,:) + hfx_bom(:,:) + hfx_bog(:,:) + hfx_dif(:,:) + hfx_opw(:,:) + hfx_snw(:,:) & 113 113 & - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) - hfx_sub(:,:) - hfx_spr(:,:) & 114 & ) * e1e2t(:,:) ) * zconv - zft_b114 & ) * e1e2t(:,:) ) * zconv - pdiag_ft 115 115 116 116 ! outputs 117 zv i= ( ( glob_sum( SUM( v_i * rhoic + v_s * rhosn, dim=3 ) * e1e2t ) * zconv &118 & - zvi_b ) * r1_rdtice - zfw) * rday119 120 zs mv= ( ( glob_sum( SUM( smv_i * rhoic , dim=3 ) * e1e2t ) * zconv &121 & - zsmv_b) * r1_rdtice + zfs ) * rday122 123 z ei= ( glob_sum( ( SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 ) &117 zv = ( ( glob_sum( SUM( v_i * rhoic + v_s * rhosn, dim=3 ) * e1e2t ) * zconv & 118 & - pdiag_v ) * r1_rdtice - zfv ) * rday 119 120 zs = ( ( glob_sum( SUM( smv_i * rhoic , dim=3 ) * e1e2t ) * zconv & 121 & - pdiag_s ) * r1_rdtice + zfs ) * rday 122 123 zt = ( glob_sum( ( SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 ) & 124 124 & + SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 ) ) * e1e2t ) * zconv & 125 & - zei_b) * r1_rdtice + zft125 & - pdiag_t ) * r1_rdtice + zft 126 126 127 127 ! zvtrp and zetrp must be close to 0 if the advection scheme is conservative … … 137 137 zv_sill = zarea * 2.5e-5 138 138 zs_sill = zarea * 25.e-5 139 z h_sill = zarea * 10.e-5139 zt_sill = zarea * 10.e-5 140 140 141 141 IF(lwp) THEN 142 IF ( ABS( zv i ) > zv_sill ) WRITE(numout,*) 'violation volume [Mt/day] (',cd_routine,') = ',zvi143 IF ( ABS( zs mv ) > zs_sill ) WRITE(numout,*) 'violation saline [psu*Mt/day] (',cd_routine,') = ',zsmv144 IF ( ABS( z ei ) > zh_sill ) WRITE(numout,*) 'violation enthalpy [GW] (',cd_routine,') = ',zei142 IF ( ABS( zv ) > zv_sill ) WRITE(numout,*) 'violation volume [Mt/day] (',cd_routine,') = ',zv 143 IF ( ABS( zs ) > zs_sill ) WRITE(numout,*) 'violation saline [psu*Mt/day] (',cd_routine,') = ',zs 144 IF ( ABS( zt ) > zt_sill ) WRITE(numout,*) 'violation enthalpy [GW] (',cd_routine,') = ',zt 145 145 IF ( ABS(zvtrp ) > zv_sill .AND. cd_routine == 'iceadv' ) THEN 146 WRITE(numout,*) 'violation vtrp [Mt/day] (',cd_routine,') = ',zvtrp 147 WRITE(numout,*) 'violation etrp [GW] (',cd_routine,') = ',zetrp 148 ENDIF 149 IF ( zvmin < -epsi10 ) WRITE(numout,*) 'violation v_i<0 [m] (',cd_routine,') = ',zvmin 150 IF ( zamax > MAX( rn_amax_n, rn_amax_s ) + epsi10 .AND. & 151 & cd_routine /= 'iceadv' .AND. cd_routine /= 'icerdgrft' ) THEN 152 WRITE(numout,*) 'violation a_i>amax (',cd_routine,') = ',zamax 153 IF ( zamax > 1._wp ) WRITE(numout,*) 'violation a_i>1 (',cd_routine,') = ',zamax 154 ENDIF 155 IF ( zamin < -epsi10 ) WRITE(numout,*) 'violation a_i<0 (',cd_routine,') = ',zamin 146 WRITE(numout,*) 'violation vtrp [Mt/day] (',cd_routine,') = ',zvtrp 147 WRITE(numout,*) 'violation etrp [GW] (',cd_routine,') = ',zetrp 148 ENDIF 149 IF ( zvmin < -epsi10 ) WRITE(numout,*) 'violation v_i<0 [m] (',cd_routine,') = ',zvmin 150 IF ( zamax > MAX(rn_amax_n,rn_amax_s) + epsi10 .AND. cd_routine /= 'iceadv' .AND. cd_routine /= 'icerdgrft' ) & 151 & WRITE(numout,*) 'violation a_i>amax (',cd_routine,') = ',zamax 152 IF ( zamin < -epsi10 ) WRITE(numout,*) 'violation a_i<0 (',cd_routine,') = ',zamin 156 153 ENDIF 157 154 ! … … 169 166 !! ** Method : This is an online diagnostics which can be activated with ln_icediachk=true 170 167 !! It prints in ocean.output if there is a violation of conservation at each time-step 171 !! The thresholds (zv_sill, zs_sill, z h_sill) which determine the violation are set to168 !! The thresholds (zv_sill, zs_sill, zt_sill) which determine the violation are set to 172 169 !! a minimum of 1 mm of ice (over the ice area) that is lost/gained spuriously during 100 years. 173 170 !! For salt and heat thresholds, ice is considered to have a salinity of 10 … … 176 173 CHARACTER(len=*), INTENT(in) :: cd_routine ! name of the routine 177 174 REAL(wp) :: zhfx, zsfx, zvfx 178 REAL(wp) :: zarea, zv_sill, zs_sill, z h_sill175 REAL(wp) :: zarea, zv_sill, zs_sill, zt_sill 179 176 REAL(wp), PARAMETER :: zconv = 1.e-9 ! convert W to GW and kg to Mt 180 177 !!---------------------------------------------------------------------- 178 179 ! water flux 180 zvfx = glob_sum( ( wfx_ice + wfx_snw + wfx_spr + wfx_sub + diag_vice + diag_vsnw ) * e1e2t ) * zconv * rday 181 182 ! salt flux 183 zsfx = glob_sum( ( sfx + diag_smvi ) * e1e2t ) * zconv * rday 181 184 182 185 ! heat flux … … 184 187 ! & - SUM( qevap_ice * a_i_b, dim=3 ) & !!clem: I think this line must be commented (but need check) 185 188 & ) * e1e2t ) * zconv 186 ! salt flux187 zsfx = glob_sum( ( sfx + diag_smvi ) * e1e2t ) * zconv * rday188 ! water flux189 zvfx = glob_sum( ( wfx_ice + wfx_snw + wfx_spr + wfx_sub + diag_vice + diag_vsnw ) * e1e2t ) * zconv * rday190 189 191 190 ! set threshold values and calculate the ice area (+epsi10 to set a threshold > 0 when there is no ice) … … 193 192 zv_sill = zarea * 2.5e-5 194 193 zs_sill = zarea * 25.e-5 195 z h_sill = zarea * 10.e-5194 zt_sill = zarea * 10.e-5 196 195 197 196 IF(lwp) THEN 198 IF( ABS( zvfx ) > zv_sill ) WRITE(numout,*) 'violation vfx [Mt/day] (',cd_routine,') = ',(zvfx)199 IF( ABS( zsfx ) > zs_sill ) WRITE(numout,*) 'violation sfx [psu*Mt/day] (',cd_routine,') = ',(zsfx)200 IF( ABS( zhfx ) > z h_sill ) WRITE(numout,*) 'violation hfx [GW] (',cd_routine,') = ',(zhfx)197 IF( ABS( zvfx ) > zv_sill ) WRITE(numout,*) 'violation vfx [Mt/day] (',cd_routine,') = ',zvfx 198 IF( ABS( zsfx ) > zs_sill ) WRITE(numout,*) 'violation sfx [psu*Mt/day] (',cd_routine,') = ',zsfx 199 IF( ABS( zhfx ) > zt_sill ) WRITE(numout,*) 'violation hfx [GW] (',cd_routine,') = ',zhfx 201 200 ENDIF 202 201 ! … … 215 214 INTEGER :: ialert_id ! number of the current alert 216 215 REAL(wp) :: ztmelts ! ice layer melting point 217 CHARACTER (len=30), DIMENSION(20) 218 INTEGER , DIMENSION(20) 216 CHARACTER (len=30), DIMENSION(20) :: cl_alname ! name of alert 217 INTEGER , DIMENSION(20) :: inb_alp ! number of alerts positive 219 218 !!------------------------------------------------------------------- 220 219 -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icedyn.F90
r8516 r8517 34 34 PUBLIC ice_dyn_init ! called by icestp.F90 35 35 36 INTEGER :: nice_dyn ! choice of the type of advection scheme36 INTEGER :: nice_dyn ! choice of the type of dynamics 37 37 ! ! associated indices: 38 INTEGER, PARAMETER :: np_dynNO = 0 ! no ice dynamics and ice advection 39 INTEGER, PARAMETER :: np_dynFULL = 1 ! full ice dynamics (rheology + advection + ridging/rafting + correction) 40 INTEGER, PARAMETER :: np_dyn = 2 ! no ridging/rafting (rheology + advection + correction) 41 INTEGER, PARAMETER :: np_dynPURE = 3 ! pure dynamics (rheology + advection) 42 38 INTEGER, PARAMETER :: np_dynFULL = 1 ! full ice dynamics (rheology + advection + ridging/rafting + correction) 39 INTEGER, PARAMETER :: np_dynRHGADV1 = 2 ! no ridging/rafting (rheology + advection + correction) 40 INTEGER, PARAMETER :: np_dynRHGADV2 = 3 ! pure dynamics (rheology + advection) 41 INTEGER, PARAMETER :: np_dynADV = 4 ! only advection w prescribed vel.(rn_uvice + advection) 42 ! 43 ! ** namelist (namdyn) ** 44 REAL(wp) :: rn_uice ! prescribed u-vel (case np_dynADV) 45 REAL(wp) :: rn_vice ! prescribed v-vel (case np_dynADV) 46 43 47 !! * Substitutions 44 48 # include "vectopt_loop_substitute.h90" … … 61 65 INTEGER, INTENT(in) :: kt ! ice time step 62 66 !! 63 INTEGER :: jl ! dummy loop indices 67 INTEGER :: ji, jj, jl ! dummy loop indices 68 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zhmax 64 69 !!-------------------------------------------------------------------- 65 70 ! … … 72 77 ENDIF 73 78 74 CALL ice_var_agg(1) ! -- aggregate ice categories75 79 ! 76 IF( ln_landfast ) THEN !-- Landfast ice parameterization: define max bottom friction80 IF( ln_landfast ) THEN !-- Landfast ice parameterization: define max bottom friction 77 81 tau_icebfr(:,:) = 0._wp 78 82 DO jl = 1, jpl … … 81 85 ENDIF 82 86 83 SELECT CASE( nice_dyn ) ! -- Set which dynamics is running 87 zhmax(:,:,:) = ht_i_b(:,:,:) !-- Record max of the surrounding 9-pts ice thick. (for CALL Hbig) 88 DO jl = 1, jpl 89 DO jj = 2, jpjm1 90 DO ji = 2, jpim1 91 !!gm use of MAXVAL here is very probably less efficient than expending the 9 values 92 zhmax(ji,jj,jl) = MAX( epsi20, MAXVAL( ht_i_b(ji-1:ji+1,jj-1:jj+1,jl) ) ) 93 END DO 94 END DO 95 END DO 96 CALL lbc_lnk( zhmax(:,:,:), 'T', 1. ) 97 ! 98 ! 99 SELECT CASE( nice_dyn ) !-- Set which dynamics is running 84 100 85 101 CASE ( np_dynFULL ) !== all dynamical processes ==! 86 CALL ice_rhg ( kt ) ! -- rheology 87 CALL ice_adv ( kt ) ! -- advection of ice 88 CALL ice_rdgrft( kt ) ! -- ridging/rafting 89 CALL ice_cor ( kt , 1 ) ! -- Corrections 90 91 CASE ( np_dyn ) !== pure dynamics only ==! (no ridging/rafting) (nono cat. case 2) 92 CALL ice_rhg ( kt ) ! -- rheology 93 CALL ice_adv ( kt ) ! -- advection of ice 94 CALL ice_cor ( kt , 1 ) ! -- Corrections 95 96 CASE ( np_dynPURE ) !== pure dynamics only ==! (nn_icedyn= 1 ) 97 CALL ice_rhg ( kt ) ! -- rheology 98 CALL ice_adv ( kt ) ! -- advection of ice 99 100 CASE ( np_dynNO ) !== prescribed ice velocities ==! (nn_icedyn= 0 ) 102 CALL ice_rhg ( kt ) ! -- rheology 103 CALL ice_adv ( kt ) ; CALL Hbig( zhmax ) ! -- advection of ice + correction on ice thickness 104 CALL ice_rdgrft( kt ) ! -- ridging/rafting 105 CALL ice_cor ( kt , 1 ) ! -- Corrections 106 107 CASE ( np_dynRHGADV1 ) !== no ridge/raft ==! (mono cat. case 2) 108 CALL ice_rhg ( kt ) ! -- rheology 109 CALL ice_adv ( kt ) ! -- advection of ice 110 CALL Hpiling ! -- simple pile-up (replaces ridging/rafting) 111 CALL ice_cor ( kt , 1 ) ! -- Corrections 112 113 CASE ( np_dynRHGADV2 ) !== no ridge/raft & no corrections ==! 114 CALL ice_rhg ( kt ) ! -- rheology 115 CALL ice_adv ( kt ) ! -- advection of ice 116 CALL Hpiling ! -- simple pile-up (replaces ridging/rafting) 117 118 CASE ( np_dynADV ) !== pure advection ==! (prescribed velocities) 101 119 u_ice(:,:) = rn_uice * umask(:,:,1) 102 120 v_ice(:,:) = rn_vice * vmask(:,:,1) 103 121 !!CALL RANDOM_NUMBER(u_ice(:,:)) 104 122 !!CALL RANDOM_NUMBER(v_ice(:,:)) 123 CALL ice_adv ( kt ) ! -- advection of ice 105 124 106 125 END SELECT … … 110 129 END SUBROUTINE ice_dyn 111 130 131 SUBROUTINE Hbig( phmax ) 132 !!------------------------------------------------------------------- 133 !! *** ROUTINE Hbig *** 134 !! 135 !! ** Purpose : Thickness correction in case advection scheme creates 136 !! abnormally tick ice 137 !! 138 !! ** Method : 1- check whether ice thickness resulting from advection is 139 !! larger than the surrounding 9-points before advection 140 !! and reduce it if a) divergence or b) convergence & at_i>0.8 141 !! 2- bound ice thickness with hi_max (99m) 142 !! 143 !! ** input : Max thickness of the surrounding 9-points 144 !!------------------------------------------------------------------- 145 REAL(wp), DIMENSION(:,:,:), INTENT(in) :: phmax ! max ice thick from surrounding 9-pts 146 ! 147 INTEGER :: ji, jj, jl ! dummy loop indices 148 REAL(wp) :: zh, zdv 149 !!------------------------------------------------------------------- 150 ! 151 CALL ice_var_zapsmall !-- zap small areas 152 ! 153 DO jl = 1, jpl 154 DO jj = 1, jpj 155 DO ji = 1, jpi 156 IF ( v_i(ji,jj,jl) > 0._wp ) THEN !-- bound to hmax 157 ! 158 zh = v_i (ji,jj,jl) / a_i(ji,jj,jl) 159 zdv = v_i(ji,jj,jl) - v_i_b(ji,jj,jl) 160 ! 161 IF ( ( zdv > 0.0 .AND. zh > phmax(ji,jj,jl) .AND. at_i_b(ji,jj) < 0.80 ) .OR. & 162 & ( zdv <= 0.0 .AND. zh > phmax(ji,jj,jl) ) ) THEN 163 a_i (ji,jj,jl) = v_i(ji,jj,jl) / MIN( phmax(ji,jj,jl), hi_max(jpl) ) !-- bound ht_i to hi_max (99 m) 164 ENDIF 165 ! 166 ENDIF 167 END DO 168 END DO 169 END DO 170 171 IF ( nn_pnd_scheme > 0 ) THEN !-- correct pond fraction to avoid a_ip > a_i 172 WHERE( a_ip(:,:,:) > a_i(:,:,:) ) a_ip(:,:,:) = a_i(:,:,:) 173 ENDIF 174 ! 175 END SUBROUTINE Hbig 176 177 SUBROUTINE Hpiling 178 !!------------------------------------------------------------------- 179 !! *** ROUTINE Hpiling *** 180 !! 181 !! ** Purpose : Simple conservative piling comparable with 1-cat models 182 !! 183 !! ** Method : pile-up ice when no ridging/rafting 184 !! 185 !! ** input : a_i 186 !!------------------------------------------------------------------- 187 INTEGER :: jl ! dummy loop indices 188 !!------------------------------------------------------------------- 189 ! 190 CALL ice_var_zapsmall !-- zap small areas 191 ! 192 at_i(:,:) = SUM( a_i(:,:,:), dim=3 ) 193 DO jl = 1, jpl 194 WHERE( at_i(:,:) > epsi20 ) 195 a_i(:,:,jl) = a_i(:,:,jl) * ( 1._wp + MIN( rn_amax_2d(:,:) - at_i(:,:) , 0._wp ) / at_i(:,:) ) 196 END WHERE 197 END DO 198 ! 199 END SUBROUTINE Hpiling 112 200 113 201 SUBROUTINE ice_dyn_init … … 123 211 !! ** input : Namelist namice_dyn 124 212 !!------------------------------------------------------------------- 125 INTEGER :: ios ! Local integer output status for namelist read126 !! 127 NAMELIST/namice_dyn/ ln_ icedyn , nn_icedyn, rn_uice , rn_vice ,&213 INTEGER :: ios, ioptio ! Local integer output status for namelist read 214 !! 215 NAMELIST/namice_dyn/ ln_dynFULL, ln_dynRHGADV, ln_dynADV, rn_uice, rn_vice, & 128 216 & rn_ishlat , rn_cio , & 129 217 & ln_landfast, rn_gamma , rn_icebfr, rn_lfrelax … … 144 232 WRITE(numout,*) '~~~~~~~~~~~~' 145 233 WRITE(numout,*) ' Namelist namice_dyn' 146 WRITE(numout,*) ' Ice dynamics (T) or not (F) ln_icedyn = ', ln_icedyn 147 WRITE(numout,*) ' associated switch nn_icedyn = ', nn_icedyn 148 WRITE(numout,*) ' =2 all processes (default option)' 149 WRITE(numout,*) ' =1 advection only (no ridging/rafting)' 150 WRITE(numout,*) ' =0 advection only with prescribed velocity given by ' 234 WRITE(numout,*) ' Full ice dynamics (rhg + adv + ridge/raft + corr) ln_dynFULL = ', ln_dynFULL 235 WRITE(numout,*) ' No ridge/raft & No cor (rhg + adv) ln_dynRHGADV = ', ln_dynRHGADV 236 WRITE(numout,*) ' Advection only (rn_uvice + adv) ln_dynADV = ', ln_dynADV 237 WRITE(numout,*) ' with prescribed velocity given by ' 151 238 WRITE(numout,*) ' a uniform field (u,v)_ice = (rn_uice,rn_vice) = (', rn_uice,',', rn_vice,')' 152 239 WRITE(numout,*) ' lateral boundary condition for sea ice dynamics rn_ishlat = ', rn_ishlat … … 158 245 ENDIF 159 246 ! !== set the choice of ice dynamics ==! 160 SELECT CASE( nn_icedyn ) 161 CASE( 2 ) 162 IF( nn_monocat /= 2 ) THEN !--- full dynamics (rheology + advection + ridging/rafting + correction) 163 nice_dyn = np_dynFULL 164 ELSE 165 nice_dyn = np_dyn !--- dynamics without ridging/rafting 166 ENDIF 167 CASE( 1 ) !--- dynamics without ridging/rafting and correction 168 nice_dyn = np_dynPURE 169 CASE( 0 ) !--- prescribed ice velocities (from namelist) 170 nice_dyn = np_dynNO 171 END SELECT 247 ioptio = 0 248 ! !--- full dynamics (rheology + advection + ridging/rafting + correction) 249 IF( ln_dynFULL .AND. nn_monocat /= 2 ) THEN ; ioptio = ioptio + 1 ; nice_dyn = np_dynFULL ; ENDIF 250 ! !--- dynamics without ridging/rafting (rheology + advection + correction) 251 IF( ln_dynFULL .AND. nn_monocat == 2 ) THEN ; ioptio = ioptio + 1 ; nice_dyn = np_dynRHGADV1 ; ENDIF 252 ! !--- dynamics without ridging/rafting and corr (rheology + advection) 253 IF( ln_dynRHGADV ) THEN ; ioptio = ioptio + 1 ; nice_dyn = np_dynRHGADV2 ; ENDIF 254 ! !--- advection only with prescribed ice velocities (from namelist) 255 IF( ln_dynADV ) THEN ; ioptio = ioptio + 1 ; nice_dyn = np_dynADV ; ENDIF 256 ! 257 IF( ioptio /= 1 ) CALL ctl_stop( 'ice_dyn_init: one and only one ice dynamics option has to be defined ' ) 172 258 ! 173 259 ! !--- Lateral boundary conditions … … 177 263 ELSEIF ( 2. < rn_ishlat ) THEN ; IF(lwp) WRITE(numout,*) ' ===>>> ice lateral strong-slip' 178 264 ENDIF 179 !180 265 ! !--- NO Landfast ice : set to zero once for all 181 266 IF( .NOT. ln_landfast ) tau_icebfr(:,:) = 0._wp 182 267 ! 183 ! !--- simple conservative piling, comparable with LIM2184 l_piling = nn_icedyn == 1 .OR. ( nn_monocat == 2 .AND. jpl == 1 )185 !186 268 END SUBROUTINE ice_dyn_init 187 269 -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceitd.F90
r8515 r8517 67 67 REAL(wp) :: zx3 68 68 REAL(wp) :: zslope ! used to compute local thermodynamic "speeds" 69 REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b ! conservation check70 69 71 70 INTEGER , DIMENSION(jpij) :: idxice2 ! compute remapping or not … … 81 80 IF( kt == nit000 .AND. lwp ) WRITE(numout,*) '-- ice_itd_rem: remapping ice thickness distribution' 82 81 83 IF( ln_icediachk ) CALL ice_cons_hsm(0, 'iceitd_rem', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)82 IF( ln_icediachk ) CALL ice_cons_hsm(0, 'iceitd_rem', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) 84 83 85 84 !----------------------------------------------------------------------------------------------- … … 311 310 ENDIF 312 311 ! 313 IF( ln_icediachk ) CALL ice_cons_hsm(1, 'iceitd_rem', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)312 IF( ln_icediachk ) CALL ice_cons_hsm(1, 'iceitd_rem', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) 314 313 ! 315 314 END SUBROUTINE ice_itd_rem -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerdgrft.F90
r8515 r8517 127 127 ! 128 128 INTEGER, PARAMETER :: nitermax = 20 129 !130 REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b131 129 !!----------------------------------------------------------------------------- 132 IF( nn_timing == 1 ) CALL timing_start('icerdgrft') 130 ! controls 131 IF( nn_timing == 1 ) CALL timing_start('icerdgrft') ! timing 132 IF( ln_icediachk ) CALL ice_cons_hsm(0, 'icerdgrft', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 133 133 134 134 IF( kt == nit000 ) THEN … … 140 140 ! 141 141 ENDIF 142 ! ! conservation test143 IF( ln_icediachk ) CALL ice_cons_hsm(0, 'icerdgrft', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)144 142 145 143 !-----------------------------------------------------------------------------! … … 284 282 CALL ice_var_agg( 1 ) 285 283 286 !-----------------------------------------------------------------------------! 287 ! control prints 288 !-----------------------------------------------------------------------------! 289 ! ! conservation test 290 IF( ln_icediachk ) CALL ice_cons_hsm(1, 'icerdgrft', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 291 292 ! ! control prints 293 IF( ln_ctl ) CALL ice_prt3D( 'icerdgrft' ) 294 ! 295 IF( nn_timing == 1 ) CALL timing_stop('icerdgrft') 284 ! controls 285 IF( ln_icediachk ) CALL ice_cons_hsm(1, 'icerdgrft', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 286 IF( ln_ctl ) CALL ice_prt3D ('icerdgrft') ! prints 287 IF( nn_timing == 1 ) CALL timing_stop ('icerdgrft') ! timing 296 288 ! 297 289 END SUBROUTINE ice_rdgrft … … 795 787 ELSEIF( ln_str_H79 ) THEN ! Ice strength => Hibler (1979) method ! 796 788 ! !--------------------------------------------------! 797 strength(:,:) = rn_pstar * vt_i(:,:) * EXP( - rn_crhg * ( 1._wp - at_i(:,:) ) ) * tmask(:,:,1)789 strength(:,:) = rn_pstar * SUM( v_i(:,:,:), dim=3 ) * EXP( -rn_crhg * ( 1._wp - SUM( a_i(:,:,:), dim=3 ) ) ) 798 790 ! 799 791 ismooth = 1 -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerhg.F90
r8516 r8517 32 32 PUBLIC ice_rhg ! called by icestp.F90 33 33 PUBLIC ice_rhg_init ! called by icestp.F90 34 35 INTEGER :: nice_rhg ! choice of the type of rheology 36 ! ! associated indices: 37 INTEGER, PARAMETER :: np_rhgEVP = 1 ! EVP rheology 38 !! INTEGER, PARAMETER :: np_rhgEAP = 2 ! EAP rheology 34 39 35 40 !! * Substitutions … … 53 58 !!-------------------------------------------------------------------- 54 59 INTEGER, INTENT(in) :: kt ! ice time step 55 !!56 REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b57 60 !!-------------------------------------------------------------------- 58 ! 59 IF( nn_timing == 1 ) CALL timing_start('icerhg') 61 ! controls 62 IF( nn_timing == 1 ) CALL timing_start('icerhg') ! timing 63 IF( ln_icediachk ) CALL ice_cons_hsm(0, 'icerhg', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 60 64 ! 61 65 IF( kt == nit000 .AND. lwp ) THEN … … 64 68 WRITE(numout,*)'~~~~~~~' 65 69 ENDIF 66 ! ! -- conservation test67 IF( ln_icediachk ) CALL ice_cons_hsm(0, 'icerhg', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)68 70 69 ! ----------------------- 70 ! Rheology (ice dynamics) 71 ! ----------------------- 72 CALL ice_rhg_evp( kt, stress1_i, stress2_i, stress12_i, u_ice, v_ice, shear_i, divu_i, delta_i ) 71 ! -------- 72 ! Rheology 73 ! -------- 74 SELECT CASE( nice_rhg ) 75 ! !------------------------! 76 CASE( np_rhgEVP ) ! Elasto-Viscous-Plastic ! 77 ! !------------------------! 78 CALL ice_rhg_evp( kt, stress1_i, stress2_i, stress12_i, u_ice, v_ice, shear_i, divu_i, delta_i ) 79 80 END SELECT 73 81 ! 74 ! !- conservation test75 IF( ln_icediachk ) CALL ice_cons_hsm(1, 'icerhg', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)76 IF( ln_ctl ) CALL ice_prt3D ('icerhg') !- Controlprints77 IF( nn_timing == 1 ) CALL timing_stop ('icerhg') !-timing82 ! controls 83 IF( ln_icediachk ) CALL ice_cons_hsm(1, 'icerhg', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 84 IF( ln_ctl ) CALL ice_prt3D ('icerhg') ! prints 85 IF( nn_timing == 1 ) CALL timing_stop ('icerhg') ! timing 78 86 ! 79 87 END SUBROUTINE ice_rhg … … 92 100 !! ** input : Namelist namice_rhg 93 101 !!------------------------------------------------------------------- 94 INTEGER :: ios ! Local integer output status for namelist read102 INTEGER :: ios, ioptio ! Local integer output status for namelist read 95 103 !! 96 104 NAMELIST/namice_rhg/ ln_rhg_EVP, rn_creepl, rn_ecc , nn_nevp, rn_relast … … 118 126 ENDIF 119 127 ! 128 ! !== set the choice of ice advection ==! 129 ioptio = 0 130 IF( ln_rhg_EVP ) THEN ; ioptio = ioptio + 1 ; nice_rhg = np_rhgEVP ; ENDIF 131 !! IF( ln_rhg_EAP ) THEN ; ioptio = ioptio + 1 ; nice_rhg = np_rhgEAP ; ENDIF 132 IF( ioptio /= 1 ) CALL ctl_stop( 'ice_rhg_init: choose one and only one ice rheology' ) 133 134 ! ! allocate tke arrays 135 !!clem example IF( zdf_tke_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'zdf_tke_init : unable to allocate arrays' ) 136 ! 120 137 END SUBROUTINE ice_rhg_init 121 138 -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerhg_evp.F90
r8515 r8517 23 23 USE sbc_oce , ONLY : ln_ice_embd, nn_fsbc, ssh_m 24 24 USE sbc_ice , ONLY : utau_ice, vtau_ice, snwice_mass, snwice_mass_b 25 USE ice ! ice variables26 USE icerdgrft ! ice strength25 USE ice ! sea-ice: ice variables 26 USE icerdgrft ! sea-ice: ice strength 27 27 ! 28 28 USE lbclnk ! Lateral Boundary Condition / MPP link -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icestp.F90
r8516 r8517 258 258 CALL ice_itd_init ! ice thickness distribution initialization 259 259 ! 260 CALL ice_dyn_init ! set ice dynamics parameters261 !262 CALL ice_rdgrft_init! set ice ridging/rafting parameters263 !264 CALL ice_rhg_init ! set ice rheologyparameters265 !266 CALL ice_adv_init ! set ice advection parameters 267 !268 CALL ice_thd_init! set ice thermodynics parameters269 !270 CALL ice_thd_sal_init ! set ice salinity parameters271 260 IF( ln_icedyn ) THEN 261 CALL ice_dyn_init ! set ice dynamics parameters 262 CALL ice_rdgrft_init ! set ice ridging/rafting parameters 263 CALL ice_rhg_init ! set ice rheology parameters 264 CALL ice_adv_init ! set ice advection parameters 265 ENDIF 266 267 IF( ln_icethd ) THEN 268 CALL ice_thd_init ! set ice thermodynics parameters 269 CALL ice_thd_sal_init ! set ice salinity parameters 270 ENDIF 271 272 272 ! MV MP 2016 273 273 CALL lim_mp_init ! set melt ponds parameters … … 283 283 CALL ice_var_glo2eqv 284 284 ! 285 CALL ice_update_init ! ice surface boundary condition 286 ! 287 CALL ice_alb_init ! ice surface albedo 288 ! 289 CALL ice_dia_init ! initialization for diags 290 ! 291 fr_i (:,:) = at_i(:,:) ! initialisation of sea-ice fraction 292 tn_ice(:,:,:) = t_su(:,:,:) ! initialisation of surface temp for coupled simu 293 ! 285 CALL ice_update_init ! ice surface boundary condition 286 ! 287 CALL ice_alb_init ! ice surface albedo 288 ! 289 CALL ice_dia_init ! initialization for diags 290 ! 291 fr_i (:,:) = at_i(:,:) ! initialisation of sea-ice fraction 292 tn_ice(:,:,:) = t_su(:,:,:) ! initialisation of surface temp for coupled simu 293 ! 294 ! ! set max concentration in both hemispheres 294 295 WHERE( gphit(:,:) > 0._wp ) ; rn_amax_2d(:,:) = rn_amax_n ! NH 295 296 ELSEWHERE ; rn_amax_2d(:,:) = rn_amax_s ! SH … … 306 307 !! 307 308 !! ** Method : Read the namice_run namelist and check the parameter 308 !! values called at the first timestep (nit000)309 !! values called at the first timestep (nit000) 309 310 !! 310 311 !! ** input : Namelist namice_run 311 312 !!------------------------------------------------------------------- 312 313 INTEGER :: ios ! Local integer output status for namelist read 313 NAMELIST/namice_run/ jpl, nlay_i, nlay_s, nn_monocat, rn_amax_n, rn_amax_s, cn_icerst_in, cn_icerst_indir,&314 & cn_icerst_out, cn_icerst_outdir314 NAMELIST/namice_run/ jpl, nlay_i, nlay_s, nn_monocat, ln_icedyn, ln_icethd, rn_amax_n, rn_amax_s, & 315 & cn_icerst_in, cn_icerst_indir, cn_icerst_out, cn_icerst_outdir 315 316 !!------------------------------------------------------------------- 316 317 ! … … 333 334 WRITE(numout,*) ' number of snow layers nlay_s = ', nlay_s 334 335 WRITE(numout,*) ' virtual ITD mono-category param (1-4) or not (0) nn_monocat = ', nn_monocat 336 WRITE(numout,*) ' Ice dynamics (T) or not (F) ln_icedyn = ', ln_icedyn 337 WRITE(numout,*) ' Ice thermodynamics (T) or not (F) ln_icethd = ', ln_icethd 335 338 WRITE(numout,*) ' maximum ice concentration for NH = ', rn_amax_n 336 339 WRITE(numout,*) ' maximum ice concentration for SH = ', rn_amax_s 337 340 ENDIF 338 341 ! 339 ! 342 ! !--- check consistency 340 343 IF ( jpl > 1 .AND. nn_monocat == 1 ) THEN 341 344 nn_monocat = 0 … … 349 352 IF( ln_bdy .AND. ln_icediachk ) CALL ctl_warn('ice_run_init: online conservation check does not work with BDY') 350 353 ! 351 rdt_ice = REAL(nn_fsbc) * rdt 354 rdt_ice = REAL(nn_fsbc) * rdt !--- sea-ice timestep and inverse 352 355 r1_rdtice = 1._wp / rdt_ice 353 IF( lwp ) WRITE(numout,*) ' ice timestep rdt_ice 354 ! 355 r1_nlay_i = 1._wp / REAL( nlay_i, wp ) 356 IF( lwp ) WRITE(numout,*) ' ice timestep rdt_ice = ', rdt_ice 357 ! 358 r1_nlay_i = 1._wp / REAL( nlay_i, wp ) !--- inverse of nlay_i and nlay_s 356 359 r1_nlay_s = 1._wp / REAL( nlay_s, wp ) 357 360 ! -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd.F90
r8515 r8517 81 81 INTEGER :: ji, jj, jk, jl ! dummy loop indices 82 82 REAL(wp) :: zfric_u, zqld, zqfr, zqfr_neg 83 REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b ! conservation check84 83 REAL(wp), PARAMETER :: zfric_umin = 0._wp ! lower bound for the friction velocity (cice value=5.e-04) 85 84 REAL(wp), PARAMETER :: zch = 0.0057_wp ! heat transfer coefficient … … 87 86 ! 88 87 !!------------------------------------------------------------------- 89 90 IF( nn_timing == 1 ) CALL timing_start('icethd') 88 ! controls 89 IF( nn_timing == 1 ) CALL timing_start('icethd') ! timing 90 IF( ln_icediachk ) CALL ice_cons_hsm(0, 'icethd', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 91 91 92 92 IF( kt == nit000 .AND. lwp ) THEN … … 96 96 ENDIF 97 97 98 ! conservation test99 IF( ln_icediachk ) CALL ice_cons_hsm( 0, 'icethd', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b )100 101 98 CALL ice_var_glo2eqv 102 99 … … 257 254 oa_i(:,:,:) = o_i(:,:,:) * a_i(:,:,:) 258 255 259 IF( ln_icediachk ) CALL ice_cons_hsm( 1, 'icethd', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)256 IF( ln_icediachk ) CALL ice_cons_hsm(1, 'icethd', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) 260 257 ! 261 258 CALL ice_var_zapsmall ! --- remove very small ice concentration (<1e-10) --- ! … … 266 263 IF( ln_icedO ) CALL ice_thd_lac ! --- frazil ice growing in leads --- ! 267 264 ! 268 IF( ln_icectl ) CALL ice_prt( kt, iiceprt, jiceprt, 1, ' - ice thermodyn. - ' ) ! control print269 IF( ln_ ctl ) CALL ice_prt3D( 'icethd' ) ! Control print270 !271 IF( nn_timing == 1 ) CALL timing_stop('icethd')265 ! controls 266 IF( ln_icectl ) CALL ice_prt (kt, iiceprt, jiceprt, 1, ' - ice thermodyn. - ') ! prints 267 IF( ln_ctl ) CALL ice_prt3D ('icethd') ! prints 268 IF( nn_timing == 1 ) CALL timing_stop('icethd') ! timing 272 269 ! 273 270 END SUBROUTINE ice_thd … … 539 536 INTEGER :: ios ! Local integer output status for namelist read 540 537 !! 541 NAMELIST/namice_thd/ ln_icethd,rn_kappa_i, ln_cndi_U64, ln_cndi_P07, ln_dqns_i, rn_cnd_s, &538 NAMELIST/namice_thd/ rn_kappa_i, ln_cndi_U64, ln_cndi_P07, ln_dqns_i, rn_cnd_s, & 542 539 & ln_icedH, rn_blow_s, & 543 540 & ln_icedA, rn_beta, rn_dmin, & … … 560 557 WRITE(numout,*) '~~~~~~~~~~~~~' 561 558 WRITE(numout,*) ' Namelist namice_thd' 562 WRITE(numout,*) ' Ice thermodynamics (T) or not (F) ln_icethd = ', ln_icethd563 559 WRITE(numout,*) ' -- icethd_dif --' 564 560 WRITE(numout,*) ' extinction radiation parameter in sea ice rn_kappa_i = ', rn_kappa_i -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_lac.F90
r8514 r8517 83 83 84 84 REAL(wp) :: zv_newfra 85 REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b ! conservation check86 85 87 86 INTEGER , DIMENSION(jpij) :: jcat ! indexes of categories where new ice grows … … 113 112 !!-----------------------------------------------------------------------! 114 113 115 IF( ln_icediachk ) CALL ice_cons_hsm(0, 'icethd_lac', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)114 IF( ln_icediachk ) CALL ice_cons_hsm(0, 'icethd_lac', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) 116 115 117 116 CALL ice_var_agg(1) … … 479 478 ENDIF ! nidx > 0 480 479 ! 481 IF( ln_icediachk ) CALL ice_cons_hsm(1, 'icethd_lac', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)480 IF( ln_icediachk ) CALL ice_cons_hsm(1, 'icethd_lac', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) 482 481 ! 483 482 END SUBROUTINE ice_thd_lac -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icevar.F90
r8514 r8517 104 104 ! END MP 2016 105 105 106 DO jj = 1, jpj ! open water fraction 107 DO ji = 1, jpi 108 ato_i(ji,jj) = MAX( 1._wp - at_i(ji,jj), 0._wp ) 109 END DO 110 END DO 111 !!gm I think this should do the work : 112 ! ato_i(:,:) = MAX( 1._wp - at_i(:,:), 0._wp ) 113 !!gm end 106 ato_i(:,:) = 1._wp - at_i(:,:) ! open water fraction 114 107 115 108 IF( kn > 1 ) THEN … … 543 536 544 537 ! to be sure that at_i is the sum of a_i(jl) 545 at_i (:,:) = a_i(:,:,1) 546 vt_i (:,:) = v_i(:,:,1) 547 DO jl = 2, jpl 548 at_i(:,:) = at_i(:,:) + a_i(:,:,jl) 549 vt_i(:,:) = vt_i(:,:) + v_i(:,:,jl) 550 END DO 551 552 ! open water = 1 if at_i=0 (no re-calculation of ato_i here) 553 DO jj = 1, jpj 554 DO ji = 1, jpi 555 rswitch = MAX( 0._wp , SIGN( 1._wp, - at_i(ji,jj) ) ) 556 ato_i(ji,jj) = rswitch + (1._wp - rswitch ) * ato_i(ji,jj) 557 END DO 558 END DO 538 at_i (:,:) = SUM( a_i(:,:,:), dim=3 ) 539 vt_i (:,:) = SUM( v_i(:,:,:), dim=3 ) 540 541 ! open water = 1 if at_i=0 542 WHERE( at_i(:,:) == 0._wp ) ato_i(:,:) = 1._wp 559 543 ! 560 544 END SUBROUTINE ice_var_zapsmall
Note: See TracChangeset
for help on using the changeset viewer.