Changeset 5067 for branches/2015
- Timestamp:
- 2015-02-06T19:12:57+01:00 (10 years ago)
- Location:
- branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO
- Files:
-
- 29 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/dom_ice.F90
r5064 r5067 5 5 !!====================================================================== 6 6 !! History : 3.0 ! 2003-08 (M. Vancoppenolle) LIM-3 original code 7 !! 4.0! 2011-02 (G. Madec) dynamical allocation7 !! 3.5 ! 2011-02 (G. Madec) dynamical allocation 8 8 !!---------------------------------------------------------------------- 9 9 USE in_out_manager ! I/O manager … … 20 20 INTEGER, PUBLIC :: njeq , njeqm1 !: j-index of the equator if it is inside the domain 21 21 22 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: fcor !: coriolis coefficient 23 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tms, tmi !: temperature mask, mask for stress 24 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tmu, tmv !: mask at u and v velocity points 25 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tmf !: mask at f-point 26 22 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: fcor !: coriolis coefficient 27 23 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: wght !: weight of the 4 neighbours to compute averages 28 24 … … 41 37 !!------------------------------------------------------------------- 42 38 ! 43 ALLOCATE( fcor(jpi,jpj) , & 44 & tms (jpi,jpj) , tmi (jpi,jpj) , & 45 & tmu (jpi,jpj) , tmv (jpi,jpj) , & 46 & tmf (jpi,jpj) , & 47 & wght(jpi,jpj,2,2) , STAT = dom_ice_alloc ) 39 ALLOCATE( fcor(jpi,jpj), wght(jpi,jpj,2,2), STAT = dom_ice_alloc ) 48 40 ! 49 41 IF( dom_ice_alloc /= 0 ) CALL ctl_warn( 'dom_ice_alloc: failed to allocate arrays.' ) -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/ice.F90
r5064 r5067 165 165 166 166 ! !!** ice-thickness distribution namelist (namiceitd) ** 167 INTEGER , PUBLIC :: nn_itdshp !: categories distribution following: tanh function (1), or h^(-alpha) function (2) 168 REAL(wp), PUBLIC :: rn_itmean !: mean thickness of the domain (used to compute the distribution, nn_itdshp = 2 only) 167 INTEGER , PUBLIC :: jpl !: number of ice categories 168 INTEGER , PUBLIC :: nlay_i !: number of ice layers 169 INTEGER , PUBLIC :: nlay_s !: number of snow layers 170 INTEGER , PUBLIC :: nn_catbnd !: categories distribution following: tanh function (1), or h^(-alpha) function (2) 171 REAL(wp), PUBLIC :: rn_himean !: mean thickness of the domain (used to compute the distribution, nn_itdshp = 2 only) 169 172 170 173 ! !!** ice-dynamics namelist (namicedyn) ** 171 INTEGER , PUBLIC :: nevp !: number of iterations for subcycling 172 REAL(wp), PUBLIC :: cw !: drag coefficient for oceanic stress 173 REAL(wp), PUBLIC :: pstar !: determines ice strength (N/M), Hibler JPO79 174 REAL(wp), PUBLIC :: c_rhg !: determines changes in ice strength 175 REAL(wp), PUBLIC :: creepl !: creep limit : has to be under 1.0e-9 176 REAL(wp), PUBLIC :: ecc !: eccentricity of the elliptical yield curve 177 REAL(wp), PUBLIC :: ahi0 !: sea-ice hor. eddy diffusivity coeff. (m2/s) 178 REAL(wp), PUBLIC :: relast !: ratio => telast/rdt_ice (1/3 or 1/9 depending on nb of subcycling nevp) 174 INTEGER , PUBLIC :: nn_icestr !: ice strength parameterization (0=Hibler79 1=Rothrock75) 175 INTEGER , PUBLIC :: nn_icestr_bvf !: use brine volume to diminish ice strength 176 INTEGER , PUBLIC :: nn_nevp !: number of iterations for subcycling 177 REAL(wp), PUBLIC :: rn_cio !: drag coefficient for oceanic stress 178 REAL(wp), PUBLIC :: rn_pstar !: determines ice strength (N/M), Hibler JPO79 179 REAL(wp), PUBLIC :: rn_crhg !: determines changes in ice strength 180 REAL(wp), PUBLIC :: rn_creepl !: creep limit : has to be under 1.0e-9 181 REAL(wp), PUBLIC :: rn_ecc !: eccentricity of the elliptical yield curve 182 REAL(wp), PUBLIC :: rn_ahi0_ref !: sea-ice hor. eddy diffusivity coeff. (m2/s) 183 REAL(wp), PUBLIC :: rn_relast !: ratio => telast/rdt_ice (1/3 or 1/9 depending on nb of subcycling nevp) 179 184 180 185 ! !!** ice-salinity namelist (namicesal) ** 181 REAL(wp), PUBLIC :: s_i_max!: maximum ice salinity [PSU]182 REAL(wp), PUBLIC :: s_i_min!: minimum ice salinity [PSU]183 REAL(wp), PUBLIC :: sal_G!: restoring salinity for gravity drainage [PSU]184 REAL(wp), PUBLIC :: sal_F!: restoring salinity for flushing [PSU]185 REAL(wp), PUBLIC :: time_G!: restoring time constant for gravity drainage (= 20 days) [s]186 REAL(wp), PUBLIC :: time_F!: restoring time constant for gravity drainage (= 10 days) [s]187 REAL(wp), PUBLIC :: bulk_sal!: bulk salinity (ppt) in case of constant salinity186 REAL(wp), PUBLIC :: rn_simax !: maximum ice salinity [PSU] 187 REAL(wp), PUBLIC :: rn_simin !: minimum ice salinity [PSU] 188 REAL(wp), PUBLIC :: rn_sal_gd !: restoring salinity for gravity drainage [PSU] 189 REAL(wp), PUBLIC :: rn_sal_fl !: restoring salinity for flushing [PSU] 190 REAL(wp), PUBLIC :: rn_time_gd !: restoring time constant for gravity drainage (= 20 days) [s] 191 REAL(wp), PUBLIC :: rn_time_fl !: restoring time constant for gravity drainage (= 10 days) [s] 192 REAL(wp), PUBLIC :: rn_icesal !: bulk salinity (ppt) in case of constant salinity 188 193 189 194 ! !!** ice-salinity namelist (namicesal) ** 190 INTEGER , PUBLIC :: n um_sal!: salinity configuration used in the model195 INTEGER , PUBLIC :: nn_icesal !: salinity configuration used in the model 191 196 ! ! 1 - constant salinity in both space and time 192 197 ! ! 2 - prognostic salinity (s(z,t)) 193 198 ! ! 3 - salinity profile, constant in time 194 INTEGER , PUBLIC :: thcon_i_swi!: thermal conductivity: =0 Untersteiner (1964) ; =1 Pringle et al (2007)199 INTEGER , PUBLIC :: nn_ice_thcon !: thermal conductivity: =0 Untersteiner (1964) ; =1 Pringle et al (2007) 195 200 INTEGER , PUBLIC :: nn_monocat !: virtual ITD mono-category parameterizations (1) or not (0) 196 201 197 202 ! !!** ice-mechanical redistribution namelist (namiceitdme) 198 REAL(wp), PUBLIC :: Cs!: fraction of shearing energy contributing to ridging199 REAL(wp), PUBLIC :: Cf !: ratio of ridging work to PE loss200 REAL(wp), PUBLIC :: fsnowrdg!: fractional snow loss to the ocean during ridging201 REAL(wp), PUBLIC :: fsnowrft!: fractional snow loss to the ocean during ridging202 REAL(wp), PUBLIC :: Gstar!: fractional area of young ice contributing to ridging203 REAL(wp), PUBLIC :: astar!: equivalent of G* for an exponential participation function204 REAL(wp), PUBLIC :: Hstar!: thickness that determines the maximal thickness of ridged ice205 REAL(wp), PUBLIC :: hparmeter!: threshold thickness (m) for rafting / ridging206 REAL(wp), PUBLIC :: Craft!: coefficient for smoothness of the hyperbolic tangent in rafting207 REAL(wp), PUBLIC :: r idge_por!: initial porosity of ridges (0.3 regular value)208 REAL(wp), PUBLIC :: betas!: coef. for partitioning of snowfall between leads and sea ice209 REAL(wp), PUBLIC :: kappa_i!: coef. for the extinction of radiation Grenfell et al. (2006) [1/m]210 REAL(wp), PUBLIC :: n conv_i_thd!: maximal number of iterations for heat diffusion211 REAL(wp), PUBLIC :: maxer_i_thd!: maximal tolerated error (C) for heat diffusion203 REAL(wp), PUBLIC :: rn_cs !: fraction of shearing energy contributing to ridging 204 REAL(wp), PUBLIC :: rn_pe_rdg !: ridging work divided by pot. energy change in ridging, nn_ice str = 1 205 REAL(wp), PUBLIC :: rn_fsnowrdg !: fractional snow loss to the ocean during ridging 206 REAL(wp), PUBLIC :: rn_fsnowrft !: fractional snow loss to the ocean during ridging 207 REAL(wp), PUBLIC :: rn_gstar !: fractional area of young ice contributing to ridging 208 REAL(wp), PUBLIC :: rn_astar !: equivalent of G* for an exponential participation function 209 REAL(wp), PUBLIC :: rn_hstar !: thickness that determines the maximal thickness of ridged ice 210 REAL(wp), PUBLIC :: rn_hraft !: threshold thickness (m) for rafting / ridging 211 REAL(wp), PUBLIC :: rn_craft !: coefficient for smoothness of the hyperbolic tangent in rafting 212 REAL(wp), PUBLIC :: rn_por_rdg !: initial porosity of ridges (0.3 regular value) 213 REAL(wp), PUBLIC :: rn_betas !: coef. for partitioning of snowfall between leads and sea ice 214 REAL(wp), PUBLIC :: rn_kappa_i !: coef. for the extinction of radiation Grenfell et al. (2006) [1/m] 215 REAL(wp), PUBLIC :: nn_conv_dif !: maximal number of iterations for heat diffusion 216 REAL(wp), PUBLIC :: rn_terr_dif !: maximal tolerated error (C) for heat diffusion 212 217 213 218 ! !!** ice-mechanical redistribution namelist (namiceitdme) 214 INTEGER , PUBLIC :: ridge_scheme_swi !: scheme used for ice ridging 215 INTEGER , PUBLIC :: raft_swi !: rafting of ice or not 216 INTEGER , PUBLIC :: partfun_swi !: participation function: =0 Thorndike et al. (1975), =1 Lipscomb et al. (2007) 217 INTEGER , PUBLIC :: brinstren_swi !: use brine volume to diminish ice strength 218 219 REAL(wp), PUBLIC :: usecc2 !: = 1.0 / ( ecc * ecc ) 220 REAL(wp), PUBLIC :: rhoco !: = rau0 * cw 219 INTEGER , PUBLIC :: nn_rafting !: rafting of ice or not 220 INTEGER , PUBLIC :: nn_partfun !: participation function: =0 Thorndike et al. (1975), =1 Lipscomb et al. (2007) 221 222 REAL(wp), PUBLIC :: usecc2 !: = 1.0 / ( rn_ecc * rn_ecc ) 223 REAL(wp), PUBLIC :: rhoco !: = rau0 * cio 221 224 222 225 ! !!** switch for presence of ice or not … … 364 367 !!-------------------------------------------------------------------------- 365 368 ! !!: ** Namelist namicerun read in sbc_lim_init ** 366 INTEGER , PUBLIC :: jpl !: number of ice categories367 INTEGER , PUBLIC :: nlay_i !: number of ice layers368 INTEGER , PUBLIC :: nlay_s !: number of snow layers369 369 CHARACTER(len=32) , PUBLIC :: cn_icerst_in !: suffix of ice restart name (input) 370 370 CHARACTER(len=32) , PUBLIC :: cn_icerst_out !: suffix of ice restart name (output) 371 371 LOGICAL , PUBLIC :: ln_limdyn !: flag for ice dynamics (T) or not (F) 372 372 LOGICAL , PUBLIC :: ln_nicep !: flag for sea-ice points output (T) or not (F) 373 REAL(wp) , PUBLIC :: amax!: maximum ice concentration373 REAL(wp) , PUBLIC :: rn_amax !: maximum ice concentration 374 374 ! 375 375 !!-------------------------------------------------------------------------- -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limadv.F90
r5064 r5067 85 85 zs2new = MIN( 2.0 * zslpmax - 0.3334 * ABS( zs1new ), & 86 86 & MAX( ABS( zs1new ) - zslpmax, psxx(ji,jj) ) ) 87 zin0 = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tm s(ji,jj) ! Case of empty boxes & Apply mask87 zin0 = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1) ! Case of empty boxes & Apply mask 88 88 89 89 ps0 (ji,jj) = zslpmax … … 270 270 zs2new = MIN( ( 2.0 * zslpmax - 0.3334 * ABS( zs1new ) ), & 271 271 & MAX( ABS( zs1new )-zslpmax, psyy(ji,jj) ) ) 272 zin0 = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tm s(ji,jj) ! Case of empty boxes & Apply mask272 zin0 = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1) ! Case of empty boxes & Apply mask 273 273 ! 274 274 ps0 (ji,jj) = zslpmax -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limcons.F90
r5064 r5067 173 173 IF( icount == 0 ) THEN 174 174 175 zvi_b = glob_sum( SUM( v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 ) * e12t(:,:) * tm s(:,:) )176 zsmv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * e12t(:,:) * tm s(:,:) )177 zei_b = glob_sum( SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 ) * e12t(:,:) * tm s* zconv + &178 & SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 ) * e12t(:,:) * tm s* zconv )175 zvi_b = glob_sum( SUM( v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 ) * e12t(:,:) * tmask(:,:,1) ) 176 zsmv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * e12t(:,:) * tmask(:,:,1) ) 177 zei_b = glob_sum( SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 ) * e12t(:,:) * tmask(:,:,1) * zconv + & 178 & SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 ) * e12t(:,:) * tmask(:,:,1) * zconv ) 179 179 zfw_b = glob_sum( - ( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) + & 180 180 & wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) + wfx_sub(:,:) + wfx_spr(:,:) & 181 & ) * e12t(:,:) * tm s(:,:) )181 & ) * e12t(:,:) * tmask(:,:,1) ) 182 182 zfs_b = glob_sum( ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + & 183 183 & sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) & 184 & ) * e12t(:,:) * tm s(:,:) )184 & ) * e12t(:,:) * tmask(:,:,1) ) 185 185 zft_b = glob_sum( ( hfx_sum(:,:) + hfx_bom(:,:) + hfx_bog(:,:) + hfx_dif(:,:) + hfx_opw(:,:) + hfx_snw(:,:) & 186 186 & - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) - hfx_sub(:,:) - hfx_spr(:,:) & 187 & ) * e12t(:,:) * tm s(:,:) * zconv )187 & ) * e12t(:,:) * tmask(:,:,1) * zconv ) 188 188 189 189 ELSEIF( icount == 1 ) THEN … … 191 191 zfs = glob_sum( ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + & 192 192 & sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) & 193 & ) * e12t(:,:) * tm s(:,:) ) - zfs_b193 & ) * e12t(:,:) * tmask(:,:,1) ) - zfs_b 194 194 zfw = glob_sum( - ( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) + & 195 195 & wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) + wfx_sub(:,:) + wfx_spr(:,:) & 196 & ) * e12t(:,:) * tm s(:,:) ) - zfw_b196 & ) * e12t(:,:) * tmask(:,:,1) ) - zfw_b 197 197 zft = glob_sum( ( hfx_sum(:,:) + hfx_bom(:,:) + hfx_bog(:,:) + hfx_dif(:,:) + hfx_opw(:,:) + hfx_snw(:,:) & 198 198 & - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) - hfx_sub(:,:) - hfx_spr(:,:) & 199 & ) * e12t(:,:) * tm s(:,:) * zconv ) - zft_b199 & ) * e12t(:,:) * tmask(:,:,1) * zconv ) - zft_b 200 200 201 zvi = ( glob_sum( SUM( v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 ) * e12t(:,:) * tms(:,:) ) - zvi_b ) * r1_rdtice - zfw 202 zsmv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * e12t(:,:) * tms(:,:) ) - zsmv_b ) * r1_rdtice + ( zfs / rhoic ) 203 zei = glob_sum( SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 ) * e12t(:,:) * tms * zconv + & 204 & SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 ) * e12t(:,:) * tms * zconv ) * r1_rdtice - zei_b * r1_rdtice + zft 201 zvi = ( glob_sum( SUM( v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 ) & 202 & * e12t(:,:) * tmask(:,:,1) ) - zvi_b ) * r1_rdtice - zfw 203 zsmv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * e12t(:,:) * tmask(:,:,1) ) - zsmv_b ) * r1_rdtice + ( zfs / rhoic ) 204 zei = glob_sum( SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 ) * e12t(:,:) * tmask(:,:,1) * zconv + & 205 & SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 ) * e12t(:,:) * tmask(:,:,1) * zconv ) & 206 & * r1_rdtice - zei_b * r1_rdtice + zft 205 207 206 208 zvmin = glob_min(v_i) … … 213 215 IF ( ABS( zei ) > 1.e-2 ) WRITE(numout,*) 'violation enthalpy [GW] (',cd_routine,') = ',(zei) 214 216 IF ( zvmin < -epsi10 ) WRITE(numout,*) 'violation v_i<0 [m] (',cd_routine,') = ',(zvmin) 215 IF( cd_routine /= 'limtrp' .AND. cd_routine /= 'limitd_me' .AND. zamax > amax+epsi10 ) THEN217 IF( cd_routine /= 'limtrp' .AND. cd_routine /= 'limitd_me' .AND. zamax > rn_amax+epsi10 ) THEN 216 218 WRITE(numout,*) 'violation a_i>amax (',cd_routine,') = ',zamax 217 219 ENDIF -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limctl.F90
r5064 r5067 124 124 DO jj = 1, jpj 125 125 DO ji = 1, jpi 126 IF( tm s(ji,jj) <= 0._wp .AND. at_i(ji,jj) > 0._wp ) THEN126 IF( tmask(ji,jj,1) <= 0._wp .AND. at_i(ji,jj) > 0._wp ) THEN 127 127 !CALL lim_prt( kt, ji, jj, 1, ' ALERTE 6 : Ice on continents ' ) 128 !WRITE(numout,*) ' masks s, u, v : ', tm s(ji,jj), tmu(ji,jj), tmv(ji,jj)128 !WRITE(numout,*) ' masks s, u, v : ', tmask(ji,jj,1), umask(ji,jj,1), vmask(ji,jj,1) 129 129 !WRITE(numout,*) ' sst : ', sst_m(ji,jj) 130 130 !WRITE(numout,*) ' sss : ', sss_m(ji,jj) … … 293 293 WRITE(numout,*) ' ~~~~~~~~~~~~~~ ' 294 294 WRITE(numout,*) ' Simple state ' 295 WRITE(numout,*) ' masks s,u,v : ', tm s(ji,jj), tmu(ji,jj), tmv(ji,jj)295 WRITE(numout,*) ' masks s,u,v : ', tmask(ji,jj,1), umask(ji,jj,1), vmask(ji,jj,1) 296 296 WRITE(numout,*) ' lat - long : ', gphit(ji,jj), glamt(ji,jj) 297 297 WRITE(numout,*) ' Time step : ', numit -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limdiahsb.F90
r5064 r5067 71 71 72 72 ! 1/area 73 z1_area = 1._wp / MAX( glob_sum( e12t(:,:) * tm s(:,:) ), epsi06 )74 75 rswitch = MAX( 0._wp , SIGN( 1._wp , glob_sum( e12t(:,:) * tm s(:,:) ) - epsi06 ) )73 z1_area = 1._wp / MAX( glob_sum( e12t(:,:) * tmask(:,:,1) ), epsi06 ) 74 75 rswitch = MAX( 0._wp , SIGN( 1._wp , glob_sum( e12t(:,:) * tmask(:,:,1) ) - epsi06 ) ) 76 76 ! ----------------------- ! 77 77 ! 1 - Content variations ! 78 78 ! ----------------------- ! 79 zbg_ivo = glob_sum( vt_i(:,:) * e12t(:,:) * tm s(:,:) ) ! volume ice80 zbg_svo = glob_sum( vt_s(:,:) * e12t(:,:) * tm s(:,:) ) ! volume snow81 zbg_are = glob_sum( at_i(:,:) * e12t(:,:) * tm s(:,:) ) ! area82 zbg_sal = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * e12t(:,:) * tm s(:,:) ) ! mean salt content83 zbg_tem = glob_sum( ( tm_i(:,:) - rtt ) * vt_i(:,:) * e12t(:,:) * tm s(:,:) ) ! mean temp content84 85 !zbg_ihc = glob_sum( et_i(:,:) * e12t(:,:) * tm s(:,:) ) / MAX( zbg_ivo,epsi06 ) ! ice heat content86 !zbg_shc = glob_sum( et_s(:,:) * e12t(:,:) * tm s(:,:) ) / MAX( zbg_svo,epsi06 ) ! snow heat content79 zbg_ivo = glob_sum( vt_i(:,:) * e12t(:,:) * tmask(:,:,1) ) ! volume ice 80 zbg_svo = glob_sum( vt_s(:,:) * e12t(:,:) * tmask(:,:,1) ) ! volume snow 81 zbg_are = glob_sum( at_i(:,:) * e12t(:,:) * tmask(:,:,1) ) ! area 82 zbg_sal = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * e12t(:,:) * tmask(:,:,1) ) ! mean salt content 83 zbg_tem = glob_sum( ( tm_i(:,:) - rtt ) * vt_i(:,:) * e12t(:,:) * tmask(:,:,1) ) ! mean temp content 84 85 !zbg_ihc = glob_sum( et_i(:,:) * e12t(:,:) * tmask(:,:,1) ) / MAX( zbg_ivo,epsi06 ) ! ice heat content 86 !zbg_shc = glob_sum( et_s(:,:) * e12t(:,:) * tmask(:,:,1) ) / MAX( zbg_svo,epsi06 ) ! snow heat content 87 87 88 88 ! Volume 89 89 ztmp = rswitch * z1_area * r1_rau0 * rday 90 zbg_vfx = ztmp * glob_sum( emp(:,:) * e12t(:,:) * tm s(:,:) )91 zbg_vfx_bog = ztmp * glob_sum( wfx_bog(:,:) * e12t(:,:) * tm s(:,:) )92 zbg_vfx_opw = ztmp * glob_sum( wfx_opw(:,:) * e12t(:,:) * tm s(:,:) )93 zbg_vfx_sni = ztmp * glob_sum( wfx_sni(:,:) * e12t(:,:) * tm s(:,:) )94 zbg_vfx_dyn = ztmp * glob_sum( wfx_dyn(:,:) * e12t(:,:) * tm s(:,:) )95 zbg_vfx_bom = ztmp * glob_sum( wfx_bom(:,:) * e12t(:,:) * tm s(:,:) )96 zbg_vfx_sum = ztmp * glob_sum( wfx_sum(:,:) * e12t(:,:) * tm s(:,:) )97 zbg_vfx_res = ztmp * glob_sum( wfx_res(:,:) * e12t(:,:) * tm s(:,:) )98 zbg_vfx_spr = ztmp * glob_sum( wfx_spr(:,:) * e12t(:,:) * tm s(:,:) )99 zbg_vfx_snw = ztmp * glob_sum( wfx_snw(:,:) * e12t(:,:) * tm s(:,:) )100 zbg_vfx_sub = ztmp * glob_sum( wfx_sub(:,:) * e12t(:,:) * tm s(:,:) )90 zbg_vfx = ztmp * glob_sum( emp(:,:) * e12t(:,:) * tmask(:,:,1) ) 91 zbg_vfx_bog = ztmp * glob_sum( wfx_bog(:,:) * e12t(:,:) * tmask(:,:,1) ) 92 zbg_vfx_opw = ztmp * glob_sum( wfx_opw(:,:) * e12t(:,:) * tmask(:,:,1) ) 93 zbg_vfx_sni = ztmp * glob_sum( wfx_sni(:,:) * e12t(:,:) * tmask(:,:,1) ) 94 zbg_vfx_dyn = ztmp * glob_sum( wfx_dyn(:,:) * e12t(:,:) * tmask(:,:,1) ) 95 zbg_vfx_bom = ztmp * glob_sum( wfx_bom(:,:) * e12t(:,:) * tmask(:,:,1) ) 96 zbg_vfx_sum = ztmp * glob_sum( wfx_sum(:,:) * e12t(:,:) * tmask(:,:,1) ) 97 zbg_vfx_res = ztmp * glob_sum( wfx_res(:,:) * e12t(:,:) * tmask(:,:,1) ) 98 zbg_vfx_spr = ztmp * glob_sum( wfx_spr(:,:) * e12t(:,:) * tmask(:,:,1) ) 99 zbg_vfx_snw = ztmp * glob_sum( wfx_snw(:,:) * e12t(:,:) * tmask(:,:,1) ) 100 zbg_vfx_sub = ztmp * glob_sum( wfx_sub(:,:) * e12t(:,:) * tmask(:,:,1) ) 101 101 102 102 ! Salt 103 zbg_sfx = ztmp * glob_sum( sfx(:,:) * e12t(:,:) * tm s(:,:) )104 zbg_sfx_bri = ztmp * glob_sum( sfx_bri(:,:) * e12t(:,:) * tm s(:,:) )105 zbg_sfx_res = ztmp * glob_sum( sfx_res(:,:) * e12t(:,:) * tm s(:,:) )106 zbg_sfx_dyn = ztmp * glob_sum( sfx_dyn(:,:) * e12t(:,:) * tm s(:,:) )107 108 zbg_sfx_bog = ztmp * glob_sum( sfx_bog(:,:) * e12t(:,:) * tm s(:,:) )109 zbg_sfx_opw = ztmp * glob_sum( sfx_opw(:,:) * e12t(:,:) * tm s(:,:) )110 zbg_sfx_sni = ztmp * glob_sum( sfx_sni(:,:) * e12t(:,:) * tm s(:,:) )111 zbg_sfx_bom = ztmp * glob_sum( sfx_bom(:,:) * e12t(:,:) * tm s(:,:) )112 zbg_sfx_sum = ztmp * glob_sum( sfx_sum(:,:) * e12t(:,:) * tm s(:,:) )103 zbg_sfx = ztmp * glob_sum( sfx(:,:) * e12t(:,:) * tmask(:,:,1) ) 104 zbg_sfx_bri = ztmp * glob_sum( sfx_bri(:,:) * e12t(:,:) * tmask(:,:,1) ) 105 zbg_sfx_res = ztmp * glob_sum( sfx_res(:,:) * e12t(:,:) * tmask(:,:,1) ) 106 zbg_sfx_dyn = ztmp * glob_sum( sfx_dyn(:,:) * e12t(:,:) * tmask(:,:,1) ) 107 108 zbg_sfx_bog = ztmp * glob_sum( sfx_bog(:,:) * e12t(:,:) * tmask(:,:,1) ) 109 zbg_sfx_opw = ztmp * glob_sum( sfx_opw(:,:) * e12t(:,:) * tmask(:,:,1) ) 110 zbg_sfx_sni = ztmp * glob_sum( sfx_sni(:,:) * e12t(:,:) * tmask(:,:,1) ) 111 zbg_sfx_bom = ztmp * glob_sum( sfx_bom(:,:) * e12t(:,:) * tmask(:,:,1) ) 112 zbg_sfx_sum = ztmp * glob_sum( sfx_sum(:,:) * e12t(:,:) * tmask(:,:,1) ) 113 113 114 114 ! Heat budget 115 115 zbg_ihc = glob_sum( et_i(:,:) * 1.e-20 ) ! ice heat content [1.e-20 J] 116 116 zbg_shc = glob_sum( et_s(:,:) * 1.e-20 ) ! snow heat content [1.e-20 J] 117 zbg_hfx_dhc = glob_sum( diag_heat_dhc(:,:) * e12t(:,:) * tm s(:,:) ) ! [in W]118 zbg_hfx_spr = glob_sum( hfx_spr(:,:) * e12t(:,:) * tm s(:,:) ) ! [in W]119 120 zbg_hfx_thd = glob_sum( hfx_thd(:,:) * e12t(:,:) * tm s(:,:) ) ! [in W]121 zbg_hfx_dyn = glob_sum( hfx_dyn(:,:) * e12t(:,:) * tm s(:,:) ) ! [in W]122 zbg_hfx_res = glob_sum( hfx_res(:,:) * e12t(:,:) * tm s(:,:) ) ! [in W]123 zbg_hfx_sub = glob_sum( hfx_sub(:,:) * e12t(:,:) * tm s(:,:) ) ! [in W]124 zbg_hfx_snw = glob_sum( hfx_snw(:,:) * e12t(:,:) * tm s(:,:) ) ! [in W]125 zbg_hfx_sum = glob_sum( hfx_sum(:,:) * e12t(:,:) * tm s(:,:) ) ! [in W]126 zbg_hfx_bom = glob_sum( hfx_bom(:,:) * e12t(:,:) * tm s(:,:) ) ! [in W]127 zbg_hfx_bog = glob_sum( hfx_bog(:,:) * e12t(:,:) * tm s(:,:) ) ! [in W]128 zbg_hfx_dif = glob_sum( hfx_dif(:,:) * e12t(:,:) * tm s(:,:) ) ! [in W]129 zbg_hfx_opw = glob_sum( hfx_opw(:,:) * e12t(:,:) * tm s(:,:) ) ! [in W]130 zbg_hfx_out = glob_sum( hfx_out(:,:) * e12t(:,:) * tm s(:,:) ) ! [in W]131 zbg_hfx_in = glob_sum( hfx_in(:,:) * e12t(:,:) * tm s(:,:) ) ! [in W]117 zbg_hfx_dhc = glob_sum( diag_heat_dhc(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W] 118 zbg_hfx_spr = glob_sum( hfx_spr(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W] 119 120 zbg_hfx_thd = glob_sum( hfx_thd(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W] 121 zbg_hfx_dyn = glob_sum( hfx_dyn(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W] 122 zbg_hfx_res = glob_sum( hfx_res(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W] 123 zbg_hfx_sub = glob_sum( hfx_sub(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W] 124 zbg_hfx_snw = glob_sum( hfx_snw(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W] 125 zbg_hfx_sum = glob_sum( hfx_sum(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W] 126 zbg_hfx_bom = glob_sum( hfx_bom(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W] 127 zbg_hfx_bog = glob_sum( hfx_bog(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W] 128 zbg_hfx_dif = glob_sum( hfx_dif(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W] 129 zbg_hfx_opw = glob_sum( hfx_opw(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W] 130 zbg_hfx_out = glob_sum( hfx_out(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W] 131 zbg_hfx_in = glob_sum( hfx_in(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W] 132 132 133 133 ! --------------------------------------------- ! 134 134 ! 2 - Trends due to forcing and ice growth/melt ! 135 135 ! --------------------------------------------- ! 136 z_frc_vol = r1_rau0 * glob_sum( - emp(:,:) * e12t(:,:) * tm s(:,:) ) ! volume fluxes137 z_frc_sal = r1_rau0 * glob_sum( sfx(:,:) * e12t(:,:) * tm s(:,:) ) ! salt fluxes136 z_frc_vol = r1_rau0 * glob_sum( - emp(:,:) * e12t(:,:) * tmask(:,:,1) ) ! volume fluxes 137 z_frc_sal = r1_rau0 * glob_sum( sfx(:,:) * e12t(:,:) * tmask(:,:,1) ) ! salt fluxes 138 138 z_bg_grme = glob_sum( - ( wfx_bog(:,:) + wfx_opw(:,:) + wfx_sni(:,:) + wfx_dyn(:,:) + & 139 & wfx_bom(:,:) + wfx_sum(:,:) + wfx_res(:,:) + wfx_snw(:,:) + wfx_sub(:,:) ) * e12t(:,:) * tm s(:,:) ) ! volume fluxes139 & wfx_bom(:,:) + wfx_sum(:,:) + wfx_res(:,:) + wfx_snw(:,:) + wfx_sub(:,:) ) * e12t(:,:) * tmask(:,:,1) ) ! volume fluxes 140 140 ! 141 141 frc_vol = frc_vol + z_frc_vol * rdt_ice -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90
r5064 r5067 85 85 IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limdyn', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 86 86 87 u_ice_b(:,:) = u_ice(:,:) * tmu(:,:)88 v_ice_b(:,:) = v_ice(:,:) * tmv(:,:)87 u_ice_b(:,:) = u_ice(:,:) * umask(:,:,1) 88 v_ice_b(:,:) = v_ice(:,:) * vmask(:,:,1) 89 89 90 90 ! Rheology (ice dynamics) … … 159 159 zv_io(:,:) = v_ice(:,:) - ssv_m(:,:) 160 160 ! frictional velocity at T-point 161 zcoef = 0.5_wp * cw161 zcoef = 0.5_wp * rn_cio 162 162 DO jj = 2, jpjm1 163 163 DO ji = fs_2, fs_jpim1 ! vector opt. 164 164 ust2s(ji,jj) = zcoef * ( zu_io(ji,jj) * zu_io(ji,jj) + zu_io(ji-1,jj) * zu_io(ji-1,jj) & 165 & + zv_io(ji,jj) * zv_io(ji,jj) + zv_io(ji,jj-1) * zv_io(ji,jj-1) ) * tms(ji,jj)165 & + zv_io(ji,jj) * zv_io(ji,jj) + zv_io(ji,jj-1) * zv_io(ji,jj-1) ) * tmask(ji,jj,1) 166 166 END DO 167 167 END DO … … 176 176 DO ji = fs_2, fs_jpim1 ! vector opt. 177 177 ust2s(ji,jj) = zcoef * SQRT( utau(ji,jj) * utau(ji,jj) + utau(ji-1,jj) * utau(ji-1,jj) & 178 & + vtau(ji,jj) * vtau(ji,jj) + vtau(ji,jj-1) * vtau(ji,jj-1) ) * tms(ji,jj)178 & + vtau(ji,jj) * vtau(ji,jj) + vtau(ji,jj-1) * vtau(ji,jj-1) ) * tmask(ji,jj,1) 179 179 END DO 180 180 END DO … … 243 243 !!------------------------------------------------------------------- 244 244 INTEGER :: ios ! Local integer output status for namelist read 245 NAMELIST/namicedyn/ cw, pstar, c_rhg, creepl, ecc, ahi0, nevp, relast245 NAMELIST/namicedyn/ nn_icestr, nn_icestr_bvf, rn_pstar, rn_crhg, rn_cio, rn_creepl, rn_ecc, nn_nevp, rn_relast, rn_ahi0_ref 246 246 INTEGER :: ji, jj 247 247 REAL(wp) :: za00, zd_max … … 261 261 WRITE(numout,*) 'lim_dyn_init : ice parameters for ice dynamics ' 262 262 WRITE(numout,*) '~~~~~~~~~~~~' 263 WRITE(numout,*) ' drag coefficient for oceanic stress cw = ', cw 264 WRITE(numout,*) ' first bulk-rheology parameter pstar = ', pstar 265 WRITE(numout,*) ' second bulk-rhelogy parameter c_rhg = ', c_rhg 266 WRITE(numout,*) ' creep limit creepl = ', creepl 267 WRITE(numout,*) ' eccentricity of the elliptical yield curve ecc = ', ecc 268 WRITE(numout,*) ' horizontal diffusivity coeff. (orca2 grid) ahi0 = ', ahi0 269 WRITE(numout,*) ' number of iterations for subcycling nevp = ', nevp 270 WRITE(numout,*) ' ratio of elastic timescale over ice time step relast = ', relast 263 WRITE(numout,*)' ice strength parameterization (0=Hibler 1=Rothrock) nn_icestr =', nn_icestr 264 WRITE(numout,*)' Including brine volume in ice strength comp. nn_icestr_bvf=', nn_icestr_bvf 265 WRITE(numout,*) ' drag coefficient for oceanic stress rn_cio = ', rn_cio 266 WRITE(numout,*) ' first bulk-rheology parameter rn_pstar = ', rn_pstar 267 WRITE(numout,*) ' second bulk-rhelogy parameter rn_crhg = ', rn_crhg 268 WRITE(numout,*) ' creep limit rn_creepl = ', rn_creepl 269 WRITE(numout,*) ' eccentricity of the elliptical yield curve rn_ecc = ', rn_ecc 270 WRITE(numout,*) ' number of iterations for subcycling nn_nevp = ', nn_nevp 271 WRITE(numout,*) ' ratio of elastic timescale over ice time step rn_relast = ', rn_relast 272 WRITE(numout,*) ' horizontal diffusivity coeff. (orca2 grid) rn_ahi0_ref = ', rn_ahi0_ref 271 273 ENDIF 272 274 ! 273 usecc2 = 1._wp / ( ecc *ecc )274 rhoco = rau0 * cw275 usecc2 = 1._wp / ( rn_ecc * rn_ecc ) 276 rhoco = rau0 * rn_cio 275 277 ! 276 278 ! Diffusion coefficients … … 278 280 IF( lk_mpp ) CALL mpp_max( zd_max ) ! max over the global domain 279 281 280 za00 = ahi0/ ( 1.e05_wp ) ! 1.e05 = 100km = max grid space at 60° latitude in orca2282 za00 = rn_ahi0_ref / ( 1.e05_wp ) ! 1.e05 = 100km = max grid space at 60° latitude in orca2 281 283 ! (60° = min latitude for ice cover) 282 284 DO jj = 1, jpj -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90
r5064 r5067 35 35 36 36 ! !!** init namelist (namiceini) ** 37 REAL(wp) :: thres_sst ! threshold water temperature for initial sea ice38 REAL(wp) :: hts_ini_n ! initial snow thickness in the north39 REAL(wp) :: hts_ini_s ! initial snow thickness in the south40 REAL(wp) :: hti_ini_n ! initial ice thickness in the north41 REAL(wp) :: hti_ini_s ! initial ice thickness in the south42 REAL(wp) :: ati_ini_n ! initial leads area in the north43 REAL(wp) :: ati_ini_s ! initial leads area in the south44 REAL(wp) :: smi_ini_n ! initial salinity45 REAL(wp) :: smi_ini_s ! initial salinity46 REAL(wp) :: tmi_ini_n ! initial temperature47 REAL(wp) :: tmi_ini_s ! initial temperature48 49 LOGICAL :: ln_ limini ! initialization or not37 REAL(wp) :: rn_thres_sst ! threshold water temperature for initial sea ice 38 REAL(wp) :: rn_hts_ini_n ! initial snow thickness in the north 39 REAL(wp) :: rn_hts_ini_s ! initial snow thickness in the south 40 REAL(wp) :: rn_hti_ini_n ! initial ice thickness in the north 41 REAL(wp) :: rn_hti_ini_s ! initial ice thickness in the south 42 REAL(wp) :: rn_ati_ini_n ! initial leads area in the north 43 REAL(wp) :: rn_ati_ini_s ! initial leads area in the south 44 REAL(wp) :: rn_smi_ini_n ! initial salinity 45 REAL(wp) :: rn_smi_ini_s ! initial salinity 46 REAL(wp) :: rn_tmi_ini_n ! initial temperature 47 REAL(wp) :: rn_tmi_ini_s ! initial temperature 48 49 LOGICAL :: ln_iceini ! initialization or not 50 50 !!---------------------------------------------------------------------- 51 51 !! LIM 3.0, UCL-LOCEAN-IPSL (2008) … … 112 112 ! surface temperature 113 113 DO jl = 1, jpl ! loop over categories 114 t_su (:,:,jl) = rtt * tm s(:,:)115 tn_ice(:,:,jl) = rtt * tm s(:,:)114 t_su (:,:,jl) = rtt * tmask(:,:,1) 115 tn_ice(:,:,jl) = rtt * tmask(:,:,1) 116 116 END DO 117 117 118 118 ! basal temperature (considered at freezing point) 119 t_bo(:,:) = ( eos_fzp( tsn(:,:,1,jp_sal) ) + rt0 ) * tm s(:,:)120 121 IF( ln_ limini ) THEN119 t_bo(:,:) = ( eos_fzp( tsn(:,:,1,jp_sal) ) + rt0 ) * tmask(:,:,1) 120 121 IF( ln_iceini ) THEN 122 122 123 123 !-------------------------------------------------------------------- … … 127 127 DO jj = 1, jpj ! ice if sst <= t-freez + ttest 128 128 DO ji = 1, jpi 129 IF( ( tsn(ji,jj,1,jp_tem) - ( t_bo(ji,jj) - rt0 ) ) * tm s(ji,jj) >=thres_sst ) THEN130 zswitch(ji,jj) = 0._wp * tm s(ji,jj) ! no ice129 IF( ( tsn(ji,jj,1,jp_tem) - ( t_bo(ji,jj) - rt0 ) ) * tmask(ji,jj,1) >= rn_thres_sst ) THEN 130 zswitch(ji,jj) = 0._wp * tmask(ji,jj,1) ! no ice 131 131 ELSE 132 zswitch(ji,jj) = 1._wp * tm s(ji,jj) ! ice132 zswitch(ji,jj) = 1._wp * tmask(ji,jj,1) ! ice 133 133 ENDIF 134 134 END DO … … 155 155 !----------------------------- 156 156 ! assign initial thickness, concentration, snow depth and salinity to an hemisphere-dependent array 157 zht_i_ini(1) = hti_ini_n ; zht_i_ini(2) =hti_ini_s ! ice thickness158 zht_s_ini(1) = hts_ini_n ; zht_s_ini(2) =hts_ini_s ! snow depth159 zat_i_ini(1) = ati_ini_n ; zat_i_ini(2) =ati_ini_s ! ice concentration160 zsm_i_ini(1) = smi_ini_n ; zsm_i_ini(2) =smi_ini_s ! bulk ice salinity161 ztm_i_ini(1) = tmi_ini_n ; ztm_i_ini(2) =tmi_ini_s ! temperature (ice and snow)157 zht_i_ini(1) = rn_hti_ini_n ; zht_i_ini(2) = rn_hti_ini_s ! ice thickness 158 zht_s_ini(1) = rn_hts_ini_n ; zht_s_ini(2) = rn_hts_ini_s ! snow depth 159 zat_i_ini(1) = rn_ati_ini_n ; zat_i_ini(2) = rn_ati_ini_s ! ice concentration 160 zsm_i_ini(1) = rn_smi_ini_n ; zsm_i_ini(2) = rn_smi_ini_s ! bulk ice salinity 161 ztm_i_ini(1) = rn_tmi_ini_n ; ztm_i_ini(2) = rn_tmi_ini_s ! temperature (ice and snow) 162 162 163 163 zvt_i_ini(:) = zht_i_ini(:) * zat_i_ini(:) ! ice volume … … 316 316 ht_i(ji,jj,jl) = zswitch(ji,jj) * zh_i_ini(jl,zhemis(ji,jj)) ! ice thickness 317 317 ht_s(ji,jj,jl) = ht_i(ji,jj,jl) * ( zht_s_ini( zhemis(ji,jj) ) / zht_i_ini( zhemis(ji,jj) ) ) ! snow depth 318 sm_i(ji,jj,jl) = zswitch(ji,jj) * zsm_i_ini(zhemis(ji,jj)) !+ ( 1._wp - zswitch(ji,jj) ) * s_i_min ! salinity318 sm_i(ji,jj,jl) = zswitch(ji,jj) * zsm_i_ini(zhemis(ji,jj)) !+ ( 1._wp - zswitch(ji,jj) ) * rn_simin ! salinity 319 319 o_i(ji,jj,jl) = zswitch(ji,jj) * 1._wp + ( 1._wp - zswitch(ji,jj) ) ! age 320 320 t_su(ji,jj,jl) = zswitch(ji,jj) * ztm_i_ini(zhemis(ji,jj)) + ( 1._wp - zswitch(ji,jj) ) * rtt ! surf temp … … 359 359 DO ji = 1, jpi 360 360 t_i(ji,jj,jk,jl) = zswitch(ji,jj) * ztm_i_ini(zhemis(ji,jj)) + ( 1._wp - zswitch(ji,jj) ) * rtt 361 s_i(ji,jj,jk,jl) = zswitch(ji,jj) * zsm_i_ini(zhemis(ji,jj)) !+ ( 1._wp - zswitch(ji,jj) ) * s_i_min361 s_i(ji,jj,jk,jl) = zswitch(ji,jj) * zsm_i_ini(zhemis(ji,jj)) !+ ( 1._wp - zswitch(ji,jj) ) * rn_simin 362 362 ztmelts = - tmut * s_i(ji,jj,jk,jl) + rtt !Melting temperature in K 363 363 … … 377 377 378 378 ELSE 379 ! if ln_ limini=false379 ! if ln_iceini=false 380 380 a_i (:,:,:) = 0._wp 381 381 v_i (:,:,:) = 0._wp … … 393 393 DO jl = 1, jpl 394 394 DO jk = 1, nlay_i 395 t_i(:,:,jk,jl) = rtt * tm s(:,:)395 t_i(:,:,jk,jl) = rtt * tmask(:,:,1) 396 396 END DO 397 397 DO jk = 1, nlay_s 398 t_s(:,:,jk,jl) = rtt * tm s(:,:)398 t_s(:,:,jk,jl) = rtt * tmask(:,:,1) 399 399 END DO 400 400 END DO 401 401 402 ENDIF ! ln_ limini402 ENDIF ! ln_iceini 403 403 404 404 at_i (:,:) = 0.0_wp … … 474 474 !! 8.5 ! 07-11 (M. Vancoppenolle) rewritten initialization 475 475 !!----------------------------------------------------------------------------- 476 NAMELIST/namiceini/ ln_ limini, thres_sst, hts_ini_n, hts_ini_s, hti_ini_n,hti_ini_s, &477 & ati_ini_n, ati_ini_s, smi_ini_n, smi_ini_s, tmi_ini_n,tmi_ini_s476 NAMELIST/namiceini/ ln_iceini, rn_thres_sst, rn_hts_ini_n, rn_hts_ini_s, rn_hti_ini_n, rn_hti_ini_s, & 477 & rn_ati_ini_n, rn_ati_ini_s, rn_smi_ini_n, rn_smi_ini_s, rn_tmi_ini_n, rn_tmi_ini_s 478 478 INTEGER :: ios ! Local integer output status for namelist read 479 479 !!----------------------------------------------------------------------------- … … 495 495 WRITE(numout,*) 'lim_istate_init : ice parameters inititialisation ' 496 496 WRITE(numout,*) '~~~~~~~~~~~~~~~' 497 WRITE(numout,*) ' initialization with ice (T) or not (F) ln_ limini = ', ln_limini498 WRITE(numout,*) ' threshold water temp. for initial sea-ice thres_sst = ',thres_sst499 WRITE(numout,*) ' initial snow thickness in the north hts_ini_n = ',hts_ini_n500 WRITE(numout,*) ' initial snow thickness in the south hts_ini_s = ',hts_ini_s501 WRITE(numout,*) ' initial ice thickness in the north hti_ini_n = ',hti_ini_n502 WRITE(numout,*) ' initial ice thickness in the south hti_ini_s = ',hti_ini_s503 WRITE(numout,*) ' initial ice concentr. in the north ati_ini_n = ',ati_ini_n504 WRITE(numout,*) ' initial ice concentr. in the north ati_ini_s = ',ati_ini_s505 WRITE(numout,*) ' initial ice salinity in the north smi_ini_n = ',smi_ini_n506 WRITE(numout,*) ' initial ice salinity in the south smi_ini_s = ',smi_ini_s507 WRITE(numout,*) ' initial ice/snw temp in the north tmi_ini_n = ',tmi_ini_n508 WRITE(numout,*) ' initial ice/snw temp in the south tmi_ini_s = ',tmi_ini_s497 WRITE(numout,*) ' initialization with ice (T) or not (F) ln_iceini = ', ln_iceini 498 WRITE(numout,*) ' threshold water temp. for initial sea-ice rn_thres_sst = ', rn_thres_sst 499 WRITE(numout,*) ' initial snow thickness in the north rn_hts_ini_n = ', rn_hts_ini_n 500 WRITE(numout,*) ' initial snow thickness in the south rn_hts_ini_s = ', rn_hts_ini_s 501 WRITE(numout,*) ' initial ice thickness in the north rn_hti_ini_n = ', rn_hti_ini_n 502 WRITE(numout,*) ' initial ice thickness in the south rn_hti_ini_s = ', rn_hti_ini_s 503 WRITE(numout,*) ' initial ice concentr. in the north rn_ati_ini_n = ', rn_ati_ini_n 504 WRITE(numout,*) ' initial ice concentr. in the north rn_ati_ini_s = ', rn_ati_ini_s 505 WRITE(numout,*) ' initial ice salinity in the north rn_smi_ini_n = ', rn_smi_ini_n 506 WRITE(numout,*) ' initial ice salinity in the south rn_smi_ini_s = ', rn_smi_ini_s 507 WRITE(numout,*) ' initial ice/snw temp in the north rn_tmi_ini_n = ', rn_tmi_ini_n 508 WRITE(numout,*) ' initial ice/snw temp in the south rn_tmi_ini_s = ', rn_tmi_ini_s 509 509 ENDIF 510 510 -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90
r5064 r5067 191 191 ! (thick, newly ridged ice). 192 192 193 closing_net(ji,jj) = Cs * 0.5 * ( Delta_i(ji,jj) - ABS( divu_i(ji,jj) ) ) - MIN( divu_i(ji,jj), 0._wp )193 closing_net(ji,jj) = rn_cs * 0.5 * ( Delta_i(ji,jj) - ABS( divu_i(ji,jj) ) ) - MIN( divu_i(ji,jj), 0._wp ) 194 194 195 195 ! 2.2 divu_adv … … 488 488 END DO !jl 489 489 490 zzc = Cf * Cp ! where Cp = (g/2)*(rhow-rhoi)*(rhoi/rhow) and Cfaccounts for frictional dissipation490 zzc = rn_pe_rdg * Cp ! where Cp = (g/2)*(rhow-rhoi)*(rhoi/rhow) and rn_pe_rdg accounts for frictional dissipation 491 491 strength(:,:) = zzc * strength(:,:) / aksum(:,:) 492 492 … … 498 498 ELSE ! kstrngth ne 1: Hibler (1979) form 499 499 ! 500 strength(:,:) = Pstar * vt_i(:,:) * EXP( - C_rhg * ( 1._wp - at_i(:,:) ) )500 strength(:,:) = rn_pstar * vt_i(:,:) * EXP( - rn_crhg * ( 1._wp - at_i(:,:) ) ) 501 501 ! 502 502 ksmooth = 1 … … 510 510 ! CAN BE REMOVED 511 511 ! 512 IF ( brinstren_swi== 1 ) THEN512 IF ( nn_icestr_bvf == 1 ) THEN 513 513 514 514 DO jj = 1, jpj … … 542 542 ! present 543 543 zworka(ji,jj) = 4.0 * strength(ji,jj) & 544 & + strength(ji-1,jj) * tm s(ji-1,jj) &545 & + strength(ji+1,jj) * tm s(ji+1,jj) &546 & + strength(ji,jj-1) * tm s(ji,jj-1) &547 & + strength(ji,jj+1) * tm s(ji,jj+1)548 549 zw1 = 4.0 + tm s(ji-1,jj) + tms(ji+1,jj) + tms(ji,jj-1) + tms(ji,jj+1)544 & + strength(ji-1,jj) * tmask(ji-1,jj,1) & 545 & + strength(ji+1,jj) * tmask(ji+1,jj,1) & 546 & + strength(ji,jj-1) * tmask(ji,jj-1,1) & 547 & + strength(ji,jj+1) * tmask(ji,jj+1,1) 548 549 zw1 = 4.0 + tmask(ji-1,jj,1) + tmask(ji+1,jj,1) + tmask(ji,jj-1,1) + tmask(ji,jj+1,1) 550 550 zworka(ji,jj) = zworka(ji,jj) / zw1 551 551 ELSE … … 619 619 CALL wrk_alloc( jpi,jpj,jpl+2, Gsum, kkstart = -1 ) 620 620 621 Gstari = 1.0/ Gstar622 astari = 1.0/ astar621 Gstari = 1.0/rn_gstar 622 astari = 1.0/rn_astar 623 623 aksum(:,:) = 0.0 624 624 athorn(:,:,:) = 0.0 … … 686 686 !----------------------------------------------------------------- 687 687 688 IF( partfun_swi== 0 ) THEN !--- Linear formulation (Thorndike et al., 1975)688 IF( nn_partfun == 0 ) THEN !--- Linear formulation (Thorndike et al., 1975) 689 689 DO jl = 0, jpl 690 690 DO jj = 1, jpj 691 691 DO ji = 1, jpi 692 IF( Gsum(ji,jj,jl) < Gstar) THEN692 IF( Gsum(ji,jj,jl) < rn_gstar) THEN 693 693 athorn(ji,jj,jl) = Gstari * (Gsum(ji,jj,jl)-Gsum(ji,jj,jl-1)) * & 694 694 (2.0 - (Gsum(ji,jj,jl-1)+Gsum(ji,jj,jl))*Gstari) 695 ELSEIF (Gsum(ji,jj,jl-1) < Gstar) THEN696 athorn(ji,jj,jl) = Gstari * ( Gstar-Gsum(ji,jj,jl-1)) * &697 (2.0 - (Gsum(ji,jj,jl-1)+ Gstar)*Gstari)695 ELSEIF (Gsum(ji,jj,jl-1) < rn_gstar) THEN 696 athorn(ji,jj,jl) = Gstari * (rn_gstar-Gsum(ji,jj,jl-1)) * & 697 (2.0 - (Gsum(ji,jj,jl-1)+rn_gstar)*Gstari) 698 698 ELSE 699 699 athorn(ji,jj,jl) = 0.0 … … 714 714 END DO 715 715 ! 716 ENDIF ! partfun_swi717 718 IF( raft_swi== 1 ) THEN ! Ridging and rafting ice participation functions716 ENDIF ! nn_partfun 717 718 IF( nn_rafting == 1 ) THEN ! Ridging and rafting ice participation functions 719 719 ! 720 720 DO jl = 1, jpl … … 723 723 IF ( athorn(ji,jj,jl) .GT. 0._wp ) THEN 724 724 !!gm TANH( -X ) = - TANH( X ) so can be computed only 1 time.... 725 aridge(ji,jj,jl) = ( TANH ( Craft * ( ht_i(ji,jj,jl) - hparmeter) ) + 1.0 ) * 0.5 * athorn(ji,jj,jl)726 araft (ji,jj,jl) = ( TANH ( - Craft * ( ht_i(ji,jj,jl) - hparmeter) ) + 1.0 ) * 0.5 * athorn(ji,jj,jl)725 aridge(ji,jj,jl) = ( TANH ( rn_craft * ( ht_i(ji,jj,jl) - rn_hraft ) ) + 1.0 ) * 0.5 * athorn(ji,jj,jl) 726 araft (ji,jj,jl) = ( TANH ( -rn_craft * ( ht_i(ji,jj,jl) - rn_hraft ) ) + 1.0 ) * 0.5 * athorn(ji,jj,jl) 727 727 IF ( araft(ji,jj,jl) < epsi06 ) araft(ji,jj,jl) = 0._wp 728 728 aridge(ji,jj,jl) = MAX( athorn(ji,jj,jl) - araft(ji,jj,jl), 0.0 ) … … 732 732 END DO ! jl 733 733 734 ELSE ! raft_swi= 0734 ELSE ! nn_rafting = 0 735 735 ! 736 736 DO jl = 1, jpl … … 740 740 ENDIF 741 741 742 IF ( raft_swi== 1 ) THEN742 IF ( nn_rafting == 1 ) THEN 743 743 744 744 IF( MAXVAL(aridge + araft - athorn(:,:,1:jpl)) .GT. epsi10 ) THEN … … 794 794 IF (a_i(ji,jj,jl) .GT. epsi10 .AND. athorn(ji,jj,jl) .GT. 0.0 ) THEN 795 795 hi = v_i(ji,jj,jl) / a_i(ji,jj,jl) 796 hrmean = MAX(SQRT( Hstar*hi), hi*krdgmin)796 hrmean = MAX(SQRT(rn_hstar*hi), hi*krdgmin) 797 797 hrmin(ji,jj,jl) = MIN(2.0*hi, 0.5*(hrmean + hi)) 798 798 hrmax(ji,jj,jl) = 2.0*hrmean - hrmin(ji,jj,jl) … … 1030 1030 !-------------------------------------------------------------------------- 1031 1031 vrdg1(ji,jj) = vicen_init(ji,jj,jl1) * afrac(ji,jj) 1032 vrdg2(ji,jj) = vrdg1(ji,jj) * ( 1. + r idge_por)1033 vsw (ji,jj) = vrdg1(ji,jj) * r idge_por1032 vrdg2(ji,jj) = vrdg1(ji,jj) * ( 1. + rn_por_rdg ) 1033 vsw (ji,jj) = vrdg1(ji,jj) * rn_por_rdg 1034 1034 1035 1035 vsrdg(ji,jj) = vsnwn_init(ji,jj,jl1) * afrac(ji,jj) … … 1061 1061 srdg2(ji,jj) = srdg1(ji,jj) + smsw(ji,jj) ! salt content of new ridge 1062 1062 1063 !srdg2(ji,jj) = MIN( s_i_max * vrdg2(ji,jj) , zsrdg2 ) ! impose a maximum salinity1063 !srdg2(ji,jj) = MIN( rn_simax * vrdg2(ji,jj) , zsrdg2 ) ! impose a maximum salinity 1064 1064 1065 1065 sfx_dyn(ji,jj) = sfx_dyn(ji,jj) - smsw(ji,jj) * rhoic * r1_rdtice … … 1090 1090 ! ij looping 1-icells 1091 1091 1092 msnow_mlt(ji,jj) = msnow_mlt(ji,jj) + rhosn*vsrdg(ji,jj)*(1.0- fsnowrdg) & ! rafting included1093 & + rhosn*vsrft(ji,jj)*(1.0- fsnowrft)1092 msnow_mlt(ji,jj) = msnow_mlt(ji,jj) + rhosn*vsrdg(ji,jj)*(1.0-rn_fsnowrdg) & ! rafting included 1093 & + rhosn*vsrft(ji,jj)*(1.0-rn_fsnowrft) 1094 1094 1095 1095 ! in J/m2 (same as e_s) 1096 esnow_mlt(ji,jj) = esnow_mlt(ji,jj) - esrdg(ji,jj)*(1.0- fsnowrdg) & !rafting included1097 & - esrft(ji,jj)*(1.0- fsnowrft)1096 esnow_mlt(ji,jj) = esnow_mlt(ji,jj) - esrdg(ji,jj)*(1.0-rn_fsnowrdg) & !rafting included 1097 & - esrft(ji,jj)*(1.0-rn_fsnowrft) 1098 1098 1099 1099 !----------------------------------------------------------------- … … 1206 1206 a_i (ji,jj ,jl2) = a_i (ji,jj ,jl2) + ardg2 (ji,jj) * farea 1207 1207 v_i (ji,jj ,jl2) = v_i (ji,jj ,jl2) + vrdg2 (ji,jj) * fvol(ji,jj) 1208 v_s (ji,jj ,jl2) = v_s (ji,jj ,jl2) + vsrdg (ji,jj) * fvol(ji,jj) * fsnowrdg1209 e_s (ji,jj,1,jl2) = e_s (ji,jj,1,jl2) + esrdg (ji,jj) * fvol(ji,jj) * fsnowrdg1208 v_s (ji,jj ,jl2) = v_s (ji,jj ,jl2) + vsrdg (ji,jj) * fvol(ji,jj) * rn_fsnowrdg 1209 e_s (ji,jj,1,jl2) = e_s (ji,jj,1,jl2) + esrdg (ji,jj) * fvol(ji,jj) * rn_fsnowrdg 1210 1210 smv_i(ji,jj ,jl2) = smv_i(ji,jj ,jl2) + srdg2 (ji,jj) * fvol(ji,jj) 1211 1211 oa_i (ji,jj ,jl2) = oa_i (ji,jj ,jl2) + oirdg2(ji,jj) * farea … … 1238 1238 a_i (ji,jj ,jl2) = a_i (ji,jj ,jl2) + arft2 (ji,jj) 1239 1239 v_i (ji,jj ,jl2) = v_i (ji,jj ,jl2) + virft (ji,jj) 1240 v_s (ji,jj ,jl2) = v_s (ji,jj ,jl2) + vsrft (ji,jj) * fsnowrft1241 e_s (ji,jj,1,jl2) = e_s (ji,jj,1,jl2) + esrft (ji,jj) * fsnowrft1240 v_s (ji,jj ,jl2) = v_s (ji,jj ,jl2) + vsrft (ji,jj) * rn_fsnowrft 1241 e_s (ji,jj,1,jl2) = e_s (ji,jj,1,jl2) + esrft (ji,jj) * rn_fsnowrft 1242 1242 smv_i(ji,jj ,jl2) = smv_i(ji,jj ,jl2) + smrft (ji,jj) 1243 1243 oa_i (ji,jj ,jl2) = oa_i (ji,jj ,jl2) + oirft2(ji,jj) … … 1331 1331 !!------------------------------------------------------------------- 1332 1332 INTEGER :: ios ! Local integer output status for namelist read 1333 NAMELIST/namiceitdme/ r idge_scheme_swi, Cs, Cf, fsnowrdg,fsnowrft, &1334 & Gstar, astar, Hstar, raft_swi, hparmeter, Craft, ridge_por, &1335 & partfun_swi, brinstren_swi1333 NAMELIST/namiceitdme/ rn_cs, rn_pe_rdg, rn_fsnowrdg, rn_fsnowrft, & 1334 & rn_gstar, rn_astar, rn_hstar, nn_rafting, rn_hraft, rn_craft, rn_por_rdg, & 1335 & nn_partfun 1336 1336 !!------------------------------------------------------------------- 1337 1337 ! … … 1349 1349 WRITE(numout,*)' lim_itd_me_init : ice parameters for mechanical ice redistribution ' 1350 1350 WRITE(numout,*)' ~~~~~~~~~~~~~~~' 1351 WRITE(numout,*)' Switch choosing the ice redistribution scheme ridge_scheme_swi', ridge_scheme_swi 1352 WRITE(numout,*)' Fraction of shear energy contributing to ridging Cs ', Cs 1353 WRITE(numout,*)' Ratio of ridging work to PotEner change in ridging Cf ', Cf 1354 WRITE(numout,*)' Fraction of snow volume conserved during ridging fsnowrdg ', fsnowrdg 1355 WRITE(numout,*)' Fraction of snow volume conserved during ridging fsnowrft ', fsnowrft 1356 WRITE(numout,*)' Fraction of total ice coverage contributing to ridging Gstar ', Gstar 1357 WRITE(numout,*)' Equivalent to G* for an exponential part function astar ', astar 1358 WRITE(numout,*)' Quantity playing a role in max ridged ice thickness Hstar ', Hstar 1359 WRITE(numout,*)' Rafting of ice sheets or not raft_swi ', raft_swi 1360 WRITE(numout,*)' Parmeter thickness (threshold between ridge-raft) hparmeter ', hparmeter 1361 WRITE(numout,*)' Rafting hyperbolic tangent coefficient Craft ', Craft 1362 WRITE(numout,*)' Initial porosity of ridges ridge_por ', ridge_por 1363 WRITE(numout,*)' Switch for part. function (0) linear (1) exponential partfun_swi ', partfun_swi 1364 WRITE(numout,*)' Switch for including brine volume in ice strength comp. brinstren_swi ', brinstren_swi 1351 WRITE(numout,*)' Fraction of shear energy contributing to ridging rn_cs ', rn_cs 1352 WRITE(numout,*)' Ratio of ridging work to PotEner change in ridging rn_pe_rdg ', rn_pe_rdg 1353 WRITE(numout,*)' Fraction of snow volume conserved during ridging rn_fsnowrdg ', rn_fsnowrdg 1354 WRITE(numout,*)' Fraction of snow volume conserved during ridging rn_fsnowrft ', rn_fsnowrft 1355 WRITE(numout,*)' Fraction of total ice coverage contributing to ridging rn_gstar ', rn_gstar 1356 WRITE(numout,*)' Equivalent to G* for an exponential part function rn_astar ', rn_astar 1357 WRITE(numout,*)' Quantity playing a role in max ridged ice thickness rn_hstar ', rn_hstar 1358 WRITE(numout,*)' Rafting of ice sheets or not nn_rafting ', nn_rafting 1359 WRITE(numout,*)' Parmeter thickness (threshold between ridge-raft) rn_hraft ', rn_hraft 1360 WRITE(numout,*)' Rafting hyperbolic tangent coefficient rn_craft ', rn_craft 1361 WRITE(numout,*)' Initial porosity of ridges rn_por_rdg ', rn_por_rdg 1362 WRITE(numout,*)' Switch for part. function (0) linear (1) exponential nn_partfun ', nn_partfun 1365 1363 ENDIF 1366 1364 ! -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90
r5064 r5067 361 361 ii = nind_i(ji) 362 362 ij = nind_j(ji) 363 IF ( a_i(ii,ij,1) > epsi10 .AND. ht_i(ii,ij,1) < hiclim) THEN364 a_i(ii,ij,1) = a_i(ii,ij,1) * ht_i(ii,ij,1) / hiclim365 ht_i(ii,ij,1) = hiclim363 IF ( a_i(ii,ij,1) > epsi10 .AND. ht_i(ii,ij,1) < rn_himin ) THEN 364 a_i(ii,ij,1) = a_i(ii,ij,1) * ht_i(ii,ij,1) / rn_himin 365 ht_i(ii,ij,1) = rn_himin 366 366 ENDIF 367 367 END DO !ji -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limmsh.F90
r5064 r5067 41 41 !! - Definition of some constants linked with the grid 42 42 !! - Definition of the metric coef. for the sea/ice 43 !! - Initialization of the ice masks (tmsk, umsk)44 43 !! 45 44 !! Reference : Deleersnijder et al. Ocean Modelling 100, 7-10 … … 103 102 !!gm end 104 103 105 ! !== ice masks ==!106 tms(:,:) = tmask(:,:,1) ! ice T-point : use surface tmask107 tmu(:,:) = umask(:,:,1) ! ice U-point : use surface umask (C-grid EVP)108 tmv(:,:) = vmask(:,:,1) ! ice V-point : use surface vmask (C-grid EVP)109 DO jj = 1, jpjm1 ! ice F-point : recompute fmask (due to nn_shlat)110 DO ji = 1 , jpim1 ! NO vector opt.111 tmf(ji,jj) = tms(ji,jj) * tms(ji+1,jj) * tms(ji,jj+1) * tms(ji+1,jj+1)112 END DO113 END DO114 CALL lbc_lnk( tmf(:,:), 'F', 1. ) ! lateral boundary conditions115 116 104 ! 117 105 END SUBROUTINE lim_msh -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90
r5064 r5067 141 141 REAL(wp), POINTER, DIMENSION(:,:) :: u_ice2, v_ice1 ! ice u/v component on V/U point 142 142 REAL(wp), POINTER, DIMENSION(:,:) :: zf1 , zf2 ! arrays for internal stresses 143 REAL(wp), POINTER, DIMENSION(:,:) :: zmask ! mask ocean grid points 143 144 144 145 REAL(wp), POINTER, DIMENSION(:,:) :: zdt ! tension at centre of grid cells … … 156 157 157 158 CALL wrk_alloc( jpi,jpj, zpresh, zfrld1, zmass1, zcorl1, za1ct , zpreshc, zfrld2, zmass2, zcorl2, za2ct ) 158 CALL wrk_alloc( jpi,jpj, u_oce2, u_ice2, v_oce1 , v_ice1 )159 CALL wrk_alloc( jpi,jpj, u_oce2, u_ice2, v_oce1 , v_ice1 , zmask ) 159 160 CALL wrk_alloc( jpi,jpj, zf1 , zu_ice, zf2 , zv_ice , zdt , zds ) 160 161 CALL wrk_alloc( jpi,jpj, zdt , zds , zs1 , zs2 , zs12 , zresr , zpice ) … … 187 188 188 189 #if defined key_lim3 189 CALL lim_itd_me_icestrength( ridge_scheme_swi) ! LIM-3: Ice strength on T-points190 CALL lim_itd_me_icestrength( nn_icestr ) ! LIM-3: Ice strength on T-points 190 191 #endif 191 192 … … 193 194 DO ji = 1 , jpi 194 195 #if defined key_lim3 195 zpresh(ji,jj) = tm s(ji,jj) * strength(ji,jj)196 zpresh(ji,jj) = tmask(ji,jj,1) * strength(ji,jj) 196 197 #endif 197 198 #if defined key_lim2 198 zpresh(ji,jj) = tm s(ji,jj) * pstar * vt_i(ji,jj) * EXP( -c_rhg * (1. - at_i(ji,jj) ) )199 #endif 200 ! tmi= 1 where there is ice or on land201 tmi(ji,jj) = 1._wp - ( 1._wp - MAX( 0._wp , SIGN ( 1._wp , vt_i(ji,jj) - zepsi ) ) ) * tms(ji,jj)199 zpresh(ji,jj) = tmask(ji,jj,1) * pstar * vt_i(ji,jj) * EXP( -c_rhg * (1. - at_i(ji,jj) ) ) 200 #endif 201 ! zmask = 1 where there is ice or on land 202 zmask(ji,jj) = 1._wp - ( 1._wp - MAX( 0._wp , SIGN ( 1._wp , vt_i(ji,jj) - zepsi ) ) ) * tmask(ji,jj,1) 202 203 END DO 203 204 END DO … … 207 208 DO jj = k_j1+1, k_jpj-1 208 209 DO ji = 2, jpim1 !RB caution no fs_ (ji+1,jj+1) 209 zstms = tm s(ji+1,jj+1) * wght(ji+1,jj+1,2,2) + tms(ji,jj+1) * wght(ji+1,jj+1,1,2) + &210 & tm s(ji+1,jj) * wght(ji+1,jj+1,2,1) + tms(ji,jj) * wght(ji+1,jj+1,1,1)210 zstms = tmask(ji+1,jj+1,1) * wght(ji+1,jj+1,2,2) + tmask(ji,jj+1,1) * wght(ji+1,jj+1,1,2) + & 211 & tmask(ji+1,jj,1) * wght(ji+1,jj+1,2,1) + tmask(ji,jj,1) * wght(ji+1,jj+1,1,1) 211 212 zpreshc(ji,jj) = ( zpresh(ji+1,jj+1) * wght(ji+1,jj+1,2,2) + zpresh(ji,jj+1) * wght(ji+1,jj+1,1,2) + & 212 213 & zpresh(ji+1,jj) * wght(ji+1,jj+1,2,1) + zpresh(ji,jj) * wght(ji+1,jj+1,1,1) & … … 250 251 DO ji = fs_2, fs_jpim1 251 252 252 zc1 = tm s(ji ,jj) * ( rhosn * vt_s(ji ,jj ) + rhoic * vt_i(ji ,jj ) )253 zc2 = tm s(ji+1,jj) * ( rhosn * vt_s(ji+1,jj ) + rhoic * vt_i(ji+1,jj ) )254 zc3 = tm s(ji ,jj+1) * ( rhosn * vt_s(ji ,jj+1) + rhoic * vt_i(ji ,jj+1) )255 256 zt11 = tm s(ji ,jj) * e1t(ji ,jj)257 zt12 = tm s(ji+1,jj) * e1t(ji+1,jj)258 zt21 = tm s(ji,jj) * e2t(ji,jj )259 zt22 = tm s(ji,jj+1) * e2t(ji,jj+1)253 zc1 = tmask(ji ,jj ,1) * ( rhosn * vt_s(ji ,jj ) + rhoic * vt_i(ji ,jj ) ) 254 zc2 = tmask(ji+1,jj ,1) * ( rhosn * vt_s(ji+1,jj ) + rhoic * vt_i(ji+1,jj ) ) 255 zc3 = tmask(ji ,jj+1,1) * ( rhosn * vt_s(ji ,jj+1) + rhoic * vt_i(ji ,jj+1) ) 256 257 zt11 = tmask(ji ,jj,1) * e1t(ji ,jj) 258 zt12 = tmask(ji+1,jj,1) * e1t(ji+1,jj) 259 zt21 = tmask(ji,jj ,1) * e2t(ji,jj ) 260 zt22 = tmask(ji,jj+1,1) * e2t(ji,jj+1) 260 261 261 262 ! Leads area. … … 274 275 v_oce1(ji,jj) = 0.5 * ( ( v_oce(ji ,jj) + v_oce(ji ,jj-1) ) * e1t(ji,jj) & 275 276 & + ( v_oce(ji+1,jj) + v_oce(ji+1,jj-1) ) * e1t(ji+1,jj) ) & 276 & / ( e1t(ji+1,jj) + e1t(ji,jj) ) * tmu(ji,jj)277 & / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1) 277 278 278 279 u_oce2(ji,jj) = 0.5 * ( ( u_oce(ji,jj ) + u_oce(ji-1,jj ) ) * e2t(ji,jj) & 279 280 & + ( u_oce(ji,jj+1) + u_oce(ji-1,jj+1) ) * e2t(ji,jj+1) ) & 280 & / ( e2t(ji,jj+1) + e2t(ji,jj) ) * tmv(ji,jj)281 & / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1) 281 282 282 283 ! Wind stress at U,V-point … … 305 306 ! 306 307 ! Time step for subcycling 307 dtevp = rdt_ice / n evp308 dtevp = rdt_ice / nn_nevp 308 309 #if defined key_lim3 309 dtotel = dtevp / ( 2._wp * r elast * rdt_ice )310 dtotel = dtevp / ( 2._wp * rn_relast * rdt_ice ) 310 311 #else 311 312 dtotel = dtevp / ( 2._wp * telast ) … … 314 315 z1_dtevp = 1._wp / dtevp 315 316 !-ecc2: square of yield ellipse eccenticrity (reminder: must become a namelist parameter) 316 ecc2 = ecc *ecc317 ecc2 = rn_ecc * rn_ecc 317 318 ecci = 1. / ecc2 318 319 … … 323 324 324 325 ! !----------------------! 325 DO jter = 1 , n evp! loop over jter !326 DO jter = 1 , nn_nevp ! loop over jter ! 326 327 ! !----------------------! 327 328 DO jj = k_j1, k_jpj-1 … … 331 332 332 333 DO jj = k_j1+1, k_jpj-1 333 DO ji = fs_2, fs_jpim1 !RB bug no vect opt due to tmi334 DO ji = fs_2, fs_jpim1 !RB bug no vect opt due to zmask 334 335 335 336 ! … … 363 364 zds(ji,jj) = ( ( u_ice(ji,jj+1) / e1u(ji,jj+1) - u_ice(ji,jj) / e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj) & 364 365 & + ( v_ice(ji+1,jj) / e2v(ji+1,jj) - v_ice(ji,jj) / e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) & 365 & ) * r1_e12f(ji,jj) * ( 2._wp - tmf(ji,jj) ) &366 & * tmi(ji,jj) * tmi(ji,jj+1) * tmi(ji+1,jj) * tmi(ji+1,jj+1)366 & ) * r1_e12f(ji,jj) * ( 2._wp - fmask(ji,jj,1) ) & 367 & * zmask(ji,jj) * zmask(ji,jj+1) * zmask(ji+1,jj) * zmask(ji+1,jj+1) 367 368 368 369 369 370 v_ice1(ji,jj) = 0.5_wp * ( ( v_ice(ji ,jj) + v_ice(ji ,jj-1) ) * e1t(ji+1,jj) & 370 371 & + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji ,jj) ) & 371 & / ( e1t(ji+1,jj) + e1t(ji,jj) ) * tmu(ji,jj)372 & / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1) 372 373 373 374 u_ice2(ji,jj) = 0.5_wp * ( ( u_ice(ji,jj ) + u_ice(ji-1,jj ) ) * e2t(ji,jj+1) & 374 375 & + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj ) ) & 375 & / ( e2t(ji,jj+1) + e2t(ji,jj) ) * tmv(ji,jj)376 & / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1) 376 377 END DO 377 378 END DO … … 387 388 388 389 delta = SQRT( divu_i(ji,jj)**2 + ( zdt(ji,jj)**2 + zdst**2 ) * usecc2 ) 389 delta_i(ji,jj) = delta + creepl390 delta_i(ji,jj) = delta + rn_creepl 390 391 391 392 !- Calculate Delta on corners … … 398 399 & ) * r1_e12f(ji,jj) 399 400 400 zddc = SQRT( zddc**2 + ( zdtc**2 + zds(ji,jj)**2 ) * usecc2 ) + creepl401 zddc = SQRT( zddc**2 + ( zdtc**2 + zds(ji,jj)**2 ) * usecc2 ) + rn_creepl 401 402 402 403 !-Calculate stress tensor components zs1 and zs2 at centre of grid cells (see section 3.5 of CICE user's guide). … … 438 439 DO jj = k_j1+1, k_jpj-1 439 440 DO ji = fs_2, fs_jpim1 440 rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass1(ji,jj) ) ) ) * tmu(ji,jj)441 rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass1(ji,jj) ) ) ) * umask(ji,jj,1) 441 442 z0 = zmass1(ji,jj) * z1_dtevp 442 443 … … 444 445 zv_ice1 = 0.5 * ( ( v_ice(ji ,jj) + v_ice(ji ,jj-1) ) * e1t(ji ,jj) & 445 446 & + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji+1,jj) ) & 446 & / ( e1t(ji+1,jj) + e1t(ji,jj) ) * tmu(ji,jj)447 & / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1) 447 448 za = rhoco * SQRT( ( u_ice(ji,jj) - u_oce(ji,jj) )**2 + & 448 449 & ( zv_ice1 - v_oce1(ji,jj) )**2 ) * ( 1.0 - zfrld1(ji,jj) ) … … 456 457 CALL lbc_lnk( u_ice(:,:), 'U', -1. ) 457 458 #if defined key_agrif && defined key_lim2 458 CALL agrif_rhg_lim2( jter, n evp, 'U' )459 CALL agrif_rhg_lim2( jter, nn_nevp, 'U' ) 459 460 #endif 460 461 #if defined key_bdy … … 465 466 DO ji = fs_2, fs_jpim1 466 467 467 rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass2(ji,jj) ) ) ) * tmv(ji,jj)468 rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass2(ji,jj) ) ) ) * vmask(ji,jj,1) 468 469 z0 = zmass2(ji,jj) * z1_dtevp 469 470 ! SB modif because ocean has no slip boundary condition 470 471 zu_ice2 = 0.5 * ( ( u_ice(ji,jj ) + u_ice(ji-1,jj ) ) * e2t(ji,jj) & 471 472 & + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj+1) ) & 472 & / ( e2t(ji,jj+1) + e2t(ji,jj) ) * tmv(ji,jj)473 & / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1) 473 474 za = rhoco * SQRT( ( zu_ice2 - u_oce2(ji,jj) )**2 + & 474 475 & ( v_ice(ji,jj) - v_oce(ji,jj))**2 ) * ( 1.0 - zfrld2(ji,jj) ) … … 482 483 CALL lbc_lnk( v_ice(:,:), 'V', -1. ) 483 484 #if defined key_agrif && defined key_lim2 484 CALL agrif_rhg_lim2( jter, n evp, 'V' )485 CALL agrif_rhg_lim2( jter, nn_nevp, 'V' ) 485 486 #endif 486 487 #if defined key_bdy … … 491 492 DO jj = k_j1+1, k_jpj-1 492 493 DO ji = fs_2, fs_jpim1 493 rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass2(ji,jj) ) ) ) * tmv(ji,jj)494 rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass2(ji,jj) ) ) ) * vmask(ji,jj,1) 494 495 z0 = zmass2(ji,jj) * z1_dtevp 495 496 ! SB modif because ocean has no slip boundary condition 496 497 zu_ice2 = 0.5 * ( ( u_ice(ji,jj ) + u_ice(ji-1,jj ) ) * e2t(ji,jj) & 497 498 & +( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj+1) ) & 498 & / ( e2t(ji,jj+1) + e2t(ji,jj) ) * tmv(ji,jj)499 & / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1) 499 500 500 501 za = rhoco * SQRT( ( zu_ice2 - u_oce2(ji,jj) )**2 + & … … 509 510 CALL lbc_lnk( v_ice(:,:), 'V', -1. ) 510 511 #if defined key_agrif && defined key_lim2 511 CALL agrif_rhg_lim2( jter, n evp, 'V' )512 CALL agrif_rhg_lim2( jter, nn_nevp, 'V' ) 512 513 #endif 513 514 #if defined key_bdy … … 517 518 DO jj = k_j1+1, k_jpj-1 518 519 DO ji = fs_2, fs_jpim1 519 rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass1(ji,jj) ) ) ) * tmu(ji,jj)520 rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass1(ji,jj) ) ) ) * umask(ji,jj,1) 520 521 z0 = zmass1(ji,jj) * z1_dtevp 521 522 zv_ice1 = 0.5 * ( ( v_ice(ji ,jj) + v_ice(ji ,jj-1) ) * e1t(ji,jj) & 522 523 & + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji+1,jj) ) & 523 & / ( e1t(ji+1,jj) + e1t(ji,jj) ) * tmu(ji,jj)524 & / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1) 524 525 525 526 za = rhoco * SQRT( ( u_ice(ji,jj) - u_oce(ji,jj) )**2 + & … … 534 535 CALL lbc_lnk( u_ice(:,:), 'U', -1. ) 535 536 #if defined key_agrif && defined key_lim2 536 CALL agrif_rhg_lim2( jter, n evp, 'U' )537 CALL agrif_rhg_lim2( jter, nn_nevp, 'U' ) 537 538 #endif 538 539 #if defined key_bdy … … 572 573 CALL lbc_lnk( v_ice(:,:), 'V', -1. ) 573 574 #if defined key_agrif && defined key_lim2 574 CALL agrif_rhg_lim2( n evp ,nevp, 'U' )575 CALL agrif_rhg_lim2( n evp ,nevp, 'V' )575 CALL agrif_rhg_lim2( nn_nevp , nn_nevp, 'U' ) 576 CALL agrif_rhg_lim2( nn_nevp , nn_nevp, 'V' ) 576 577 #endif 577 578 #if defined key_bdy … … 585 586 v_ice1(ji,jj) = 0.5_wp * ( ( v_ice(ji ,jj) + v_ice(ji, jj-1) ) * e1t(ji+1,jj) & 586 587 & + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji ,jj) ) & 587 & / ( e1t(ji+1,jj) + e1t(ji,jj) ) * tmu(ji,jj)588 & / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1) 588 589 589 590 u_ice2(ji,jj) = 0.5_wp * ( ( u_ice(ji,jj ) + u_ice(ji-1,jj ) ) * e2t(ji,jj+1) & 590 591 & + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj ) ) & 591 & / ( e2t(ji,jj+1) + e2t(ji,jj) ) * tmv(ji,jj)592 & / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1) 592 593 ENDIF 593 594 END DO … … 599 600 ! Recompute delta, shear and div, inputs for mechanical redistribution 600 601 DO jj = k_j1+1, k_jpj-1 601 DO ji = fs_2, jpim1 !RB bug no vect opt due to tmi602 DO ji = fs_2, jpim1 !RB bug no vect opt due to zmask 602 603 !- divu_i(:,:), zdt(:,:): divergence and tension at centre 603 604 !- zds(:,:): shear on northeast corner of grid cells … … 615 616 zds(ji,jj) = ( ( u_ice(ji,jj+1) / e1u(ji,jj+1) - u_ice(ji,jj) / e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj) & 616 617 & +( v_ice(ji+1,jj) / e2v(ji+1,jj) - v_ice(ji,jj) / e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) & 617 & ) * r1_e12f(ji,jj) * ( 2.0 - tmf(ji,jj) ) &618 & * tmi(ji,jj) * tmi(ji,jj+1) * tmi(ji+1,jj) * tmi(ji+1,jj+1)618 & ) * r1_e12f(ji,jj) * ( 2.0 - fmask(ji,jj,1) ) & 619 & * zmask(ji,jj) * zmask(ji,jj+1) * zmask(ji+1,jj) * zmask(ji+1,jj+1) 619 620 620 621 zdst = ( e2u(ji,jj) * v_ice1(ji,jj) - e2u(ji-1,jj ) * v_ice1(ji-1,jj ) & … … 622 623 623 624 delta = SQRT( divu_i(ji,jj)**2 + ( zdt(ji,jj)**2 + zdst**2 ) * usecc2 ) 624 delta_i(ji,jj) = delta + creepl625 delta_i(ji,jj) = delta + rn_creepl 625 626 626 627 ENDIF … … 690 691 ! 691 692 CALL wrk_dealloc( jpi,jpj, zpresh, zfrld1, zmass1, zcorl1, za1ct , zpreshc, zfrld2, zmass2, zcorl2, za2ct ) 692 CALL wrk_dealloc( jpi,jpj, u_oce2, u_ice2, v_oce1 , v_ice1 )693 CALL wrk_dealloc( jpi,jpj, u_oce2, u_ice2, v_oce1 , v_ice1 , zmask ) 693 694 CALL wrk_dealloc( jpi,jpj, zf1 , zu_ice, zf2 , zv_ice , zdt , zds ) 694 695 CALL wrk_dealloc( jpi,jpj, zdt , zds , zs1 , zs2 , zs12 , zresr , zpice ) -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limrst.F90
r5055 r5067 522 522 ! 523 523 ! clem: I do not understand why the following IF is needed 524 ! I suspect something inconsistent in the main code with option n um_sal=1525 IF( n um_sal == 1 ) THEN524 ! I suspect something inconsistent in the main code with option nn_icesal=1 525 IF( nn_icesal == 1 ) THEN 526 526 DO jl = 1, jpl 527 sm_i(:,:,jl) = bulk_sal527 sm_i(:,:,jl) = rn_icesal 528 528 DO jk = 1, nlay_i 529 s_i(:,:,jk,jl) = bulk_sal529 s_i(:,:,jk,jl) = rn_icesal 530 530 END DO 531 531 END DO -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90
r5064 r5067 26 26 USE phycst ! physical constants 27 27 USE dom_oce ! ocean domain 28 USE dom_ice, ONLY : tms29 28 USE ice ! LIM sea-ice variables 30 29 USE sbc_ice ! Surface boundary condition: sea-ice fields … … 170 169 zemp = emp(ji,jj) * pfrld(ji,jj) & ! evaporation over oceanic fraction 171 170 & - tprecip(ji,jj) * ( 1._wp - pfrld(ji,jj) ) & ! all precipitation reach the ocean 172 & + sprecip(ji,jj) * ( 1._wp - pfrld(ji,jj)** betas ) ! except solid precip intercepted by sea-ice171 & + sprecip(ji,jj) * ( 1._wp - pfrld(ji,jj)**rn_betas ) ! except solid precip intercepted by sea-ice 173 172 ENDIF 174 173 … … 197 196 snwice_mass_b(:,:) = snwice_mass(:,:) 198 197 ! new mass per unit area 199 snwice_mass (:,:) = tm s(:,:) * ( rhosn * vt_s(:,:) + rhoic * vt_i(:,:) )198 snwice_mass (:,:) = tmask(:,:,1) * ( rhosn * vt_s(:,:) + rhoic * vt_i(:,:) ) 200 199 ! time evolution of snow+ice mass 201 200 snwice_fmass (:,:) = ( snwice_mass(:,:) - snwice_mass_b(:,:) ) * r1_rdtice … … 346 345 ! ! embedded sea ice 347 346 IF( nn_ice_embd /= 0 ) THEN ! mass exchanges between ice and ocean (case 1 or 2) set the snow+ice mass 348 snwice_mass (:,:) = tm s(:,:) * ( rhosn * vt_s(:,:) + rhoic * vt_i(:,:) )347 snwice_mass (:,:) = tmask(:,:,1) * ( rhosn * vt_s(:,:) + rhoic * vt_i(:,:) ) 349 348 snwice_mass_b(:,:) = snwice_mass(:,:) 350 349 ELSE -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90
r5064 r5067 158 158 DO jj = 1, jpj 159 159 DO ji = 1, jpi 160 rswitch = tm s(ji,jj) * ( 1._wp - MAX( 0._wp , SIGN( 1._wp , - at_i(ji,jj) + epsi10 ) ) ) ! 0 if no ice160 rswitch = tmask(ji,jj,1) * ( 1._wp - MAX( 0._wp , SIGN( 1._wp , - at_i(ji,jj) + epsi10 ) ) ) ! 0 if no ice 161 161 ! 162 162 ! ! solar irradiance transmission at the mixed layer bottom and used in the lead heat budget … … 171 171 ! precip is included in qns but not in qns_ice 172 172 IF ( lk_cpl ) THEN 173 zqld = tm s(ji,jj) * rdt_ice * &173 zqld = tmask(ji,jj,1) * rdt_ice * & 174 174 & ( zqsr(ji,jj) * fraqsr_1lev(ji,jj) + zqns(ji,jj) & ! pfrld already included in coupled mode 175 & + ( pfrld(ji,jj)** betas - pfrld(ji,jj) ) * sprecip(ji,jj) * & ! heat content of precip175 & + ( pfrld(ji,jj)**rn_betas - pfrld(ji,jj) ) * sprecip(ji,jj) * & ! heat content of precip 176 176 & ( cpic * ( MIN( tatm_ice(ji,jj), rt0_snow ) - rtt ) - lfus ) & 177 177 & + ( 1._wp - pfrld(ji,jj) ) * ( tprecip(ji,jj) - sprecip(ji,jj) ) * rcp * ( tatm_ice(ji,jj) - rtt ) ) 178 178 ELSE 179 zqld = tm s(ji,jj) * rdt_ice * &179 zqld = tmask(ji,jj,1) * rdt_ice * & 180 180 & ( pfrld(ji,jj) * ( zqsr(ji,jj) * fraqsr_1lev(ji,jj) + zqns(ji,jj) ) & 181 & + ( pfrld(ji,jj)** betas - pfrld(ji,jj) ) * sprecip(ji,jj) * & ! heat content of precip181 & + ( pfrld(ji,jj)**rn_betas - pfrld(ji,jj) ) * sprecip(ji,jj) * & ! heat content of precip 182 182 & ( cpic * ( MIN( tatm_ice(ji,jj), rt0_snow ) - rtt ) - lfus ) & 183 183 & + ( 1._wp - pfrld(ji,jj) ) * ( tprecip(ji,jj) - sprecip(ji,jj) ) * rcp * ( tatm_ice(ji,jj) - rtt ) ) … … 185 185 186 186 !-- Energy needed to bring ocean surface layer until its freezing (<0, J.m-2) --- ! 187 zqfr = tm s(ji,jj) * rau0 * rcp * fse3t_m(ji,jj) * ( t_bo(ji,jj) - ( sst_m(ji,jj) + rt0 ) )187 zqfr = tmask(ji,jj,1) * rau0 * rcp * fse3t_m(ji,jj) * ( t_bo(ji,jj) - ( sst_m(ji,jj) + rt0 ) ) 188 188 189 189 !-- Energy Budget of the leads (J.m-2). Must be < 0 to form ice … … 229 229 & + pfrld(ji,jj) * qns(ji,jj) & 230 230 ! latent heat of precip (note that precip is included in qns but not in qns_ice) 231 & + ( pfrld(ji,jj)** betas - pfrld(ji,jj) ) * sprecip(ji,jj) &231 & + ( pfrld(ji,jj)**rn_betas - pfrld(ji,jj) ) * sprecip(ji,jj) & 232 232 & * ( cpic * ( MIN( tatm_ice(ji,jj), rt0_snow ) - rtt ) - lfus ) & 233 233 & + ( 1._wp - pfrld(ji,jj) ) * ( tprecip(ji,jj) - sprecip(ji,jj) ) * rcp * ( tatm_ice(ji,jj) - rtt ) & … … 651 651 !!------------------------------------------------------------------- 652 652 INTEGER :: ios ! Local integer output status for namelist read 653 NAMELIST/namicethd/ hiccrit, fraz_swi, maxfrazb, vfrazb,Cfrazb, &654 & hiclim, parsub,betas, &655 & kappa_i, nconv_i_thd, maxer_i_thd, thcon_i_swi, &653 NAMELIST/namicethd/ rn_hnewice, nn_frazil, rn_maxfrazb, rn_vfrazb, rn_Cfrazb, & 654 & rn_himin, parsub, rn_betas, & 655 & rn_kappa_i, nn_conv_dif, rn_terr_dif, nn_ice_thcon, & 656 656 & nn_monocat 657 657 !!------------------------------------------------------------------- … … 682 682 WRITE(numout,*) 683 683 WRITE(numout,*)' Namelist of ice parameters for ice thermodynamic computation ' 684 WRITE(numout,*)' ice thick. for lateral accretion hiccrit = ', hiccrit685 WRITE(numout,*)' Frazil ice thickness as a function of wind or not fraz_swi = ', fraz_swi686 WRITE(numout,*)' Maximum proportion of frazil ice collecting at bottom maxfrazb = ',maxfrazb687 WRITE(numout,*)' Thresold relative drift speed for collection of frazil vfrazb = ',vfrazb688 WRITE(numout,*)' Squeezing coefficient for collection of frazil Cfrazb = ',Cfrazb689 WRITE(numout,*)' minimum ice thickness hiclim = ', hiclim684 WRITE(numout,*)' ice thick. for lateral accretion rn_hnewice = ', rn_hnewice 685 WRITE(numout,*)' Frazil ice thickness as a function of wind or not nn_frazil = ', nn_frazil 686 WRITE(numout,*)' Maximum proportion of frazil ice collecting at bottom rn_maxfrazb = ', rn_maxfrazb 687 WRITE(numout,*)' Thresold relative drift speed for collection of frazil rn_vfrazb = ', rn_vfrazb 688 WRITE(numout,*)' Squeezing coefficient for collection of frazil rn_Cfrazb = ', rn_Cfrazb 689 WRITE(numout,*)' minimum ice thickness rn_himin = ', rn_himin 690 690 WRITE(numout,*)' numerical carac. of the scheme for diffusion in ice ' 691 691 WRITE(numout,*)' switch for snow sublimation (=1) or not (=0) parsub = ', parsub 692 WRITE(numout,*)' coefficient for ice-lead partition of snowfall betas = ',betas693 WRITE(numout,*)' extinction radiation parameter in sea ice kappa_i = ',kappa_i694 WRITE(numout,*)' maximal n. of iter. for heat diffusion computation n conv_i_thd = ', nconv_i_thd695 WRITE(numout,*)' maximal err. on T for heat diffusion computation maxer_i_thd = ', maxer_i_thd696 WRITE(numout,*)' switch for comp. of thermal conductivity in the ice thcon_i_swi = ', thcon_i_swi692 WRITE(numout,*)' coefficient for ice-lead partition of snowfall rn_betas = ', rn_betas 693 WRITE(numout,*)' extinction radiation parameter in sea ice rn_kappa_i = ', rn_kappa_i 694 WRITE(numout,*)' maximal n. of iter. for heat diffusion computation nn_conv_dif = ', nn_conv_dif 695 WRITE(numout,*)' maximal err. on T for heat diffusion computation rn_terr_dif = ', rn_terr_dif 696 WRITE(numout,*)' switch for comp. of thermal conductivity in the ice nn_ice_thcon = ', nn_ice_thcon 697 697 WRITE(numout,*)' check heat conservation in the ice/snow con_i = ', con_i 698 698 WRITE(numout,*)' virtual ITD mono-category parameterizations (1) or not nn_monocat = ', nn_monocat -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90
r5055 r5067 112 112 !!------------------------------------------------------------------ 113 113 114 ! Discriminate between varying salinity (n um_sal=2) and prescribed cases (other values)115 SELECT CASE( n um_sal ) ! varying salinity or not114 ! Discriminate between varying salinity (nn_icesal=2) and prescribed cases (other values) 115 SELECT CASE( nn_icesal ) ! varying salinity or not 116 116 CASE( 1, 3, 4 ) ; zswitch_sal = 0 ! prescribed salinity profile 117 117 CASE( 2 ) ; zswitch_sal = 1 ! varying salinity profile … … 227 227 !----------- 228 228 ! thickness change 229 zcoeff = ( 1._wp - ( 1._wp - at_i_1d(ji) )** betas ) / at_i_1d(ji)229 zcoeff = ( 1._wp - ( 1._wp - at_i_1d(ji) )**rn_betas ) / at_i_1d(ji) 230 230 zdh_s_pre(ji) = zcoeff * sprecip_1d(ji) * rdt_ice / rhosn 231 231 ! enthalpy of the precip (>0, J.m-3) (tatm_ice is now in K) … … 405 405 ! -> need for an iterative procedure, which converges quickly 406 406 407 IF ( n um_sal == 2 ) THEN407 IF ( nn_icesal == 2 ) THEN 408 408 num_iter_max = 5 409 409 ELSE … … 631 631 ! entrapment during snow ice formation 632 632 ! new salinity difference stored (to be used in limthd_ent.F90) 633 IF ( n um_sal == 2 ) THEN633 IF ( nn_icesal == 2 ) THEN 634 634 rswitch = MAX( 0._wp , SIGN( 1._wp , ht_i_1d(ji) - epsi10 ) ) 635 635 ! salinity dif due to snow-ice formation -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90
r5064 r5067 272 272 DO ji = kideb, kiut 273 273 ! ! radiation transmitted below the layer-th ice layer 274 zradtr_i(ji,jk) = zradtr_i(ji,0) * EXP( - kappa_i * ( MAX ( 0._wp , z_i(ji,jk) ) ) )274 zradtr_i(ji,jk) = zradtr_i(ji,0) * EXP( - rn_kappa_i * ( MAX ( 0._wp , z_i(ji,jk) ) ) ) 275 275 ! ! radiation absorbed by the layer-th ice layer 276 276 zradab_i(ji,jk) = zradtr_i(ji,jk-1) - zradtr_i(ji,jk) … … 309 309 zerritmax = 1000._wp ! maximal value of error on all points 310 310 311 DO WHILE ( zerritmax > maxer_i_thd .AND. nconv < nconv_i_thd)311 DO WHILE ( zerritmax > rn_terr_dif .AND. nconv < nn_conv_dif ) 312 312 ! 313 313 nconv = nconv + 1 … … 317 317 !------------------------------------------------------------------------------| 318 318 ! 319 IF( thcon_i_swi== 0 ) THEN ! Untersteiner (1964) formula319 IF( nn_ice_thcon == 0 ) THEN ! Untersteiner (1964) formula 320 320 DO ji = kideb , kiut 321 321 ztcond_i(ji,0) = rcdic + zbeta*s_i_1d(ji,1) / MIN(-epsi10,t_i_1d(ji,1)-rtt) … … 331 331 ENDIF 332 332 333 IF( thcon_i_swi== 1 ) THEN ! Pringle et al formula included: 2.11 + 0.09 S/T - 0.011.T333 IF( nn_ice_thcon == 1 ) THEN ! Pringle et al formula included: 2.11 + 0.09 S/T - 0.011.T 334 334 DO ji = kideb , kiut 335 335 ztcond_i(ji,0) = rcdic + 0.090_wp * s_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1)-rtt ) & … … 723 723 ! 724 724 ! check that nowhere it has started to melt 725 ! zerrit(ji) is a measure of error, it has to be under maxer_i_thd725 ! zerrit(ji) is a measure of error, it has to be under terr_dif 726 726 DO ji = kideb , kiut 727 727 t_su_1d(ji) = MAX( MIN( t_su_1d(ji) , rtt ) , 190._wp ) -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90
r5064 r5067 154 154 155 155 ! Default new ice thickness 156 hicol(:,:) = hiccrit157 158 IF( fraz_swi== 1 ) THEN156 hicol(:,:) = rn_hnewice 157 158 IF( nn_frazil == 1 ) THEN 159 159 160 160 !-------------------- … … 175 175 !------------- 176 176 ! C-grid wind stress components 177 ztaux = ( utau_ice(ji-1,jj ) * tmu(ji-1,jj) &178 & + utau_ice(ji ,jj ) * tmu(ji ,jj) ) * 0.5_wp179 ztauy = ( vtau_ice(ji ,jj-1) * tmv(ji ,jj-1) &180 & + vtau_ice(ji ,jj ) * tmv(ji ,jj) ) * 0.5_wp177 ztaux = ( utau_ice(ji-1,jj ) * umask(ji-1,jj ,1) & 178 & + utau_ice(ji ,jj ) * umask(ji ,jj ,1) ) * 0.5_wp 179 ztauy = ( vtau_ice(ji ,jj-1) * vmask(ji ,jj-1,1) & 180 & + vtau_ice(ji ,jj ) * vmask(ji ,jj ,1) ) * 0.5_wp 181 181 ! Square root of wind stress 182 182 ztenagm = SQRT( SQRT( ztaux**2 + ztauy**2 ) ) … … 194 194 ! C-grid ice velocity 195 195 rswitch = MAX( 0._wp, SIGN( 1._wp , at_i(ji,jj) ) ) 196 zvgx = rswitch * ( u_ice(ji-1,jj ) * tmu(ji-1,jj ) + u_ice(ji,jj) * tmu(ji,jj) ) * 0.5_wp197 zvgy = rswitch * ( v_ice(ji ,jj-1) * tmv(ji ,jj-1) + v_ice(ji,jj) * tmv(ji,jj) ) * 0.5_wp196 zvgx = rswitch * ( u_ice(ji-1,jj ) * umask(ji-1,jj ,1) + u_ice(ji,jj) * umask(ji,jj,1) ) * 0.5_wp 197 zvgy = rswitch * ( v_ice(ji ,jj-1) * vmask(ji ,jj-1,1) + v_ice(ji,jj) * vmask(ji,jj,1) ) * 0.5_wp 198 198 199 199 !----------------------------------- … … 319 319 !---------------------- 320 320 DO ji = 1, nbpac 321 zh_newice(ji) = hiccrit322 END DO 323 IF( fraz_swi== 1 ) zh_newice(1:nbpac) = hicol_1d(1:nbpac)321 zh_newice(ji) = rn_hnewice 322 END DO 323 IF( nn_frazil == 1 ) zh_newice(1:nbpac) = hicol_1d(1:nbpac) 324 324 325 325 !---------------------- 326 326 ! Salinity of new ice 327 327 !---------------------- 328 SELECT CASE ( n um_sal )328 SELECT CASE ( nn_icesal ) 329 329 CASE ( 1 ) ! Sice = constant 330 zs_newice(1:nbpac) = bulk_sal330 zs_newice(1:nbpac) = rn_icesal 331 331 CASE ( 2 ) ! Sice = F(z,t) [Vancoppenolle et al (2005)] 332 332 DO ji = 1, nbpac 333 333 ii = MOD( npac(ji) - 1 , jpi ) + 1 334 334 ij = ( npac(ji) - 1 ) / jpi + 1 335 zs_newice(ji) = MIN( 4.606 + 0.91 / zh_newice(ji) , s_i_max , 0.5 * sss_m(ii,ij) )335 zs_newice(ji) = MIN( 4.606 + 0.91 / zh_newice(ji) , rn_simax , 0.5 * sss_m(ii,ij) ) 336 336 END DO 337 337 CASE ( 3 ) ! Sice = F(z) [multiyear ice] … … 386 386 ! A fraction zfrazb of frazil ice is accreted at the ice bottom 387 387 rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp , - zat_i_1d(ji) ) ) 388 zfrazb = rswitch * ( TANH ( Cfrazb * ( zvrel_1d(ji) - vfrazb ) ) + 1.0 ) * 0.5 *maxfrazb388 zfrazb = rswitch * ( TANH ( rn_Cfrazb * ( zvrel_1d(ji) - rn_vfrazb ) ) + 1.0 ) * 0.5 * rn_maxfrazb 389 389 zv_frazb(ji) = zfrazb * zv_newice(ji) 390 390 zv_newice(ji) = ( 1.0 - zfrazb ) * zv_newice(ji) … … 408 408 ! we keep the excessive volume in memory and attribute it later to bottom accretion 409 409 DO ji = 1, nbpac 410 IF ( za_newice(ji) > ( amax - zat_i_1d(ji) ) ) THEN411 zda_res(ji) = za_newice(ji) - ( amax - zat_i_1d(ji) )410 IF ( za_newice(ji) > ( rn_amax - zat_i_1d(ji) ) ) THEN 411 zda_res(ji) = za_newice(ji) - ( rn_amax - zat_i_1d(ji) ) 412 412 zdv_res(ji) = zda_res (ji) * zh_newice(ji) 413 413 za_newice(ji) = za_newice(ji) - zda_res (ji) -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd_sal.F90
r5055 r5067 45 45 !! 46 46 !! ** Method : 3 possibilities 47 !! -> n um_sal = 1 -> Sice = cst [ice salinity constant in both time & space]48 !! -> n um_sal = 2 -> Sice = S(z,t) [Vancoppenolle et al. 2005]49 !! -> n um_sal = 3 -> Sice = S(z) [multiyear ice]47 !! -> nn_icesal = 1 -> Sice = cst [ice salinity constant in both time & space] 48 !! -> nn_icesal = 2 -> Sice = S(z,t) [Vancoppenolle et al. 2005] 49 !! -> nn_icesal = 3 -> Sice = S(z) [multiyear ice] 50 50 !!--------------------------------------------------------------------- 51 51 INTEGER, INTENT(in) :: kideb, kiut ! thickness category index … … 65 65 ! 1) Constant salinity, constant in time | 66 66 !------------------------------------------------------------------------------| 67 !!gm comment: if n um_sal = 1 s_i_new, s_i_1d and sm_i_1d can be set to bulk_sal one for all in the initialisation phase !!68 !!gm ===>>> simplification of almost all test on n um_sal value69 IF( n um_sal == 1 ) THEN70 s_i_1d (kideb:kiut,1:nlay_i) = bulk_sal71 sm_i_1d(kideb:kiut) = bulk_sal72 s_i_new(kideb:kiut) = bulk_sal67 !!gm comment: if nn_icesal = 1 s_i_new, s_i_1d and sm_i_1d can be set to rn_icesal one for all in the initialisation phase !! 68 !!gm ===>>> simplification of almost all test on nn_icesal value 69 IF( nn_icesal == 1 ) THEN 70 s_i_1d (kideb:kiut,1:nlay_i) = rn_icesal 71 sm_i_1d(kideb:kiut) = rn_icesal 72 s_i_new(kideb:kiut) = rn_icesal 73 73 ENDIF 74 74 … … 76 76 ! Module 2 : Constant salinity varying in time | 77 77 !------------------------------------------------------------------------------| 78 IF( n um_sal == 2 ) THEN78 IF( nn_icesal == 2 ) THEN 79 79 80 80 DO ji = kideb, kiut … … 89 89 !--------------------- 90 90 ! drainage by gravity drainage 91 dsm_i_gd_1d(ji) = - igravdr * MAX( sm_i_1d(ji) - sal_G , 0._wp ) / time_G* rdt_ice91 dsm_i_gd_1d(ji) = - igravdr * MAX( sm_i_1d(ji) - rn_sal_gd , 0._wp ) / rn_time_gd * rdt_ice 92 92 ! drainage by flushing 93 dsm_i_fl_1d(ji) = - iflush * MAX( sm_i_1d(ji) - sal_F , 0._wp ) / time_F* rdt_ice93 dsm_i_fl_1d(ji) = - iflush * MAX( sm_i_1d(ji) - rn_sal_fl , 0._wp ) / rn_time_fl * rdt_ice 94 94 95 95 !----------------- … … 115 115 ! Module 3 : Profile of salinity, constant in time | 116 116 !------------------------------------------------------------------------------| 117 IF( n um_sal == 3 ) CALL lim_var_salprof1d( kideb, kiut )117 IF( nn_icesal == 3 ) CALL lim_var_salprof1d( kideb, kiut ) 118 118 119 119 ! … … 133 133 !!------------------------------------------------------------------- 134 134 INTEGER :: ios ! Local integer output status for namelist read 135 NAMELIST/namicesal/ n um_sal, bulk_sal, sal_G, time_G, sal_F, time_F, &136 & s_i_max, s_i_min135 NAMELIST/namicesal/ nn_icesal, rn_icesal, rn_sal_gd, rn_time_gd, rn_sal_fl, rn_time_fl, & 136 & rn_simax, rn_simin 137 137 !!------------------------------------------------------------------- 138 138 ! … … 150 150 WRITE(numout,*) 'lim_thd_sal_init : Ice parameters for salinity ' 151 151 WRITE(numout,*) '~~~~~~~~~~~~~~~~' 152 WRITE(numout,*) ' switch for salinity n um_sal : ', num_sal153 WRITE(numout,*) ' bulk salinity value if n um_sal = 1 : ', bulk_sal154 WRITE(numout,*) ' restoring salinity for GD : ', sal_G155 WRITE(numout,*) ' restoring time for GD : ', time_G156 WRITE(numout,*) ' restoring salinity for flushing : ', sal_F157 WRITE(numout,*) ' restoring time for flushing : ', time_F158 WRITE(numout,*) ' Maximum tolerated ice salinity : ', s_i_max159 WRITE(numout,*) ' Minimum tolerated ice salinity : ', s_i_min152 WRITE(numout,*) ' switch for salinity nn_icesal : ', nn_icesal 153 WRITE(numout,*) ' bulk salinity value if nn_icesal = 1 : ', rn_icesal 154 WRITE(numout,*) ' restoring salinity for GD : ', rn_sal_gd 155 WRITE(numout,*) ' restoring time for GD : ', rn_time_gd 156 WRITE(numout,*) ' restoring salinity for flushing : ', rn_sal_fl 157 WRITE(numout,*) ' restoring time for flushing : ', rn_time_fl 158 WRITE(numout,*) ' Maximum tolerated ice salinity : ', rn_simax 159 WRITE(numout,*) ' Minimum tolerated ice salinity : ', rn_simin 160 160 ENDIF 161 161 ! -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90
r5064 r5067 400 400 DO jj = 1, jpj 401 401 DO ji = 1, jpi 402 a_i(ji,jj,1) = MIN( a_i(ji,jj,1), amax )402 a_i(ji,jj,1) = MIN( a_i(ji,jj,1), rn_amax ) 403 403 END DO 404 404 END DO -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90
r5064 r5067 92 92 DO jj = 1, jpj 93 93 DO ji = 1, jpi 94 IF( at_i(ji,jj) > amax .AND. a_i(ji,jj,jl) > 0._wp ) THEN95 a_i(ji,jj,jl) = a_i(ji,jj,jl) * ( 1._wp - ( 1._wp - amax / at_i(ji,jj) ) )94 IF( at_i(ji,jj) > rn_amax .AND. a_i(ji,jj,jl) > 0._wp ) THEN 95 a_i(ji,jj,jl) = a_i(ji,jj,jl) * ( 1._wp - ( 1._wp - rn_amax / at_i(ji,jj) ) ) 96 96 ht_i(ji,jj,jl) = v_i(ji,jj,jl) / a_i(ji,jj,jl) 97 97 ENDIF … … 118 118 ! Ice salinity bounds 119 119 !--------------------- 120 IF ( n um_sal == 2 ) THEN120 IF ( nn_icesal == 2 ) THEN 121 121 DO jl = 1, jpl 122 122 DO jj = 1, jpj … … 126 126 ! salinity stays in bounds 127 127 rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp, - v_i(ji,jj,jl) ) ) 128 smv_i(ji,jj,jl) = rswitch * MAX( MIN( s_i_max * v_i(ji,jj,jl), smv_i(ji,jj,jl) ), s_i_min * v_i(ji,jj,jl) )128 smv_i(ji,jj,jl) = rswitch * MAX( MIN( rn_simax * v_i(ji,jj,jl), smv_i(ji,jj,jl) ), rn_simin * v_i(ji,jj,jl) ) 129 129 ! associated salt flux 130 130 sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsal ) * rhoic * r1_rdtice -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90
r5064 r5067 82 82 83 83 !---------------------------------------------------------------------- 84 ! Constrain the thickness of the smallest category above hi clim84 ! Constrain the thickness of the smallest category above himin 85 85 !---------------------------------------------------------------------- 86 86 DO jj = 1, jpj 87 87 DO ji = 1, jpi 88 IF( v_i(ji,jj,1) > 0._wp .AND. ht_i(ji,jj,1) < hiclim) THEN89 zh = hiclim/ ht_i(ji,jj,1)88 IF( v_i(ji,jj,1) > 0._wp .AND. ht_i(ji,jj,1) < rn_himin ) THEN 89 zh = rn_himin / ht_i(ji,jj,1) 90 90 ht_s(ji,jj,1) = ht_s(ji,jj,1) * zh 91 91 ht_i(ji,jj,1) = ht_i(ji,jj,1) * zh … … 106 106 DO jj = 1, jpj 107 107 DO ji = 1, jpi 108 IF( at_i(ji,jj) > amax .AND. a_i(ji,jj,jl) > 0._wp ) THEN109 a_i(ji,jj,jl) = a_i(ji,jj,jl) * ( 1._wp - ( 1._wp - amax / at_i(ji,jj) ) )108 IF( at_i(ji,jj) > rn_amax .AND. a_i(ji,jj,jl) > 0._wp ) THEN 109 a_i(ji,jj,jl) = a_i(ji,jj,jl) * ( 1._wp - ( 1._wp - rn_amax / at_i(ji,jj) ) ) 110 110 ht_i(ji,jj,jl) = v_i(ji,jj,jl) / a_i(ji,jj,jl) 111 111 ENDIF … … 132 132 ! Ice salinity 133 133 !--------------------- 134 IF ( n um_sal == 2 ) THEN134 IF ( nn_icesal == 2 ) THEN 135 135 DO jl = 1, jpl 136 136 DO jj = 1, jpj … … 140 140 ! salinity stays in bounds 141 141 rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp, - v_i(ji,jj,jl) ) ) 142 smv_i(ji,jj,jl) = rswitch * MAX( MIN( s_i_max * v_i(ji,jj,jl), smv_i(ji,jj,jl) ), s_i_min * v_i(ji,jj,jl) ) !+ s_i_min * ( 1._wp - rswitch ) * v_i(ji,jj,jl)142 smv_i(ji,jj,jl) = rswitch * MAX( MIN( rn_simax * v_i(ji,jj,jl), smv_i(ji,jj,jl) ), rn_simin * v_i(ji,jj,jl) ) !+ rn_simin * ( 1._wp - rswitch ) * v_i(ji,jj,jl) 143 143 ! associated salt flux 144 144 sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsal ) * rhoic * r1_rdtice … … 167 167 CALL lbc_lnk( v_ice(:,:), 'V', -1. ) 168 168 !mask velocities 169 u_ice(:,:) = u_ice(:,:) * tmu(:,:)170 v_ice(:,:) = v_ice(:,:) * tmv(:,:)169 u_ice(:,:) = u_ice(:,:) * umask(:,:,1) 170 v_ice(:,:) = v_ice(:,:) * vmask(:,:,1) 171 171 172 172 ! for outputs -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90
r5064 r5067 1 1 2 MODULE limvar 2 3 !!====================================================================== … … 169 170 END DO 170 171 171 IF( n um_sal == 2 )THEN172 IF( nn_icesal == 2 )THEN 172 173 DO jl = 1, jpl 173 174 DO jj = 1, jpj … … 291 292 ! Vertically constant, constant in time 292 293 !--------------------------------------- 293 IF( n um_sal == 1 ) s_i(:,:,:,:) = bulk_sal294 IF( nn_icesal == 1 ) s_i(:,:,:,:) = rn_icesal 294 295 295 296 !----------------------------------- 296 297 ! Salinity profile, varying in time 297 298 !----------------------------------- 298 IF( n um_sal == 2 ) THEN299 IF( nn_icesal == 2 ) THEN 299 300 ! 300 301 DO jk = 1, nlay_i … … 344 345 END DO ! jl 345 346 ! 346 ENDIF ! n um_sal347 ENDIF ! nn_icesal 347 348 348 349 !------------------------------------------------------- … … 350 351 !------------------------------------------------------- 351 352 352 IF( n um_sal == 3 ) THEN ! Schwarzacher (1959) multiyear salinity profile (mean = 2.30)353 IF( nn_icesal == 3 ) THEN ! Schwarzacher (1959) multiyear salinity profile (mean = 2.30) 353 354 ! 354 355 sm_i(:,:,:) = 2.30_wp … … 362 363 END DO 363 364 ! 364 ENDIF ! n um_sal365 ENDIF ! nn_icesal 365 366 ! 366 367 CALL wrk_dealloc( jpi, jpj, jpl, z_slope_s, zalpha ) … … 451 452 ! Vertically constant, constant in time 452 453 !--------------------------------------- 453 IF( n um_sal == 1 ) s_i_1d(:,:) = bulk_sal454 IF( nn_icesal == 1 ) s_i_1d(:,:) = rn_icesal 454 455 455 456 !------------------------------------------------------ … … 457 458 !------------------------------------------------------ 458 459 459 IF( n um_sal == 2 ) THEN460 IF( nn_icesal == 2 ) THEN 460 461 ! 461 462 DO ji = kideb, kiut ! Slope of the linear profile zs_zero … … 495 496 !------------------------------------------------------- 496 497 497 IF( n um_sal == 3 ) THEN ! Schwarzacher (1959) multiyear salinity profile (mean = 2.30)498 IF( nn_icesal == 3 ) THEN ! Schwarzacher (1959) multiyear salinity profile (mean = 2.30) 498 499 ! 499 500 sm_i_1d(:) = 2.30_wp … … 570 571 571 572 ! ice salinity must stay in bounds 572 IF( n um_sal == 2 ) THEN573 smv_i(ji,jj,jl) = MAX( MIN( s_i_max * v_i(ji,jj,jl), smv_i(ji,jj,jl) ), s_i_min * v_i(ji,jj,jl) )573 IF( nn_icesal == 2 ) THEN 574 smv_i(ji,jj,jl) = MAX( MIN( rn_simax * v_i(ji,jj,jl), smv_i(ji,jj,jl) ), rn_simin * v_i(ji,jj,jl) ) 574 575 ENDIF 575 576 ! update exchanges with ocean -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90
r5055 r5067 106 106 DO jj = 2 , jpjm1 107 107 DO ji = 2 , jpim1 108 z2da(ji,jj) = ( u_ice(ji,jj) * tmu(ji,jj) + u_ice(ji-1,jj) * tmu(ji-1,jj) ) * 0.5_wp109 z2db(ji,jj) = ( v_ice(ji,jj) * tmv(ji,jj) + v_ice(ji,jj-1) * tmv(ji,jj-1) ) * 0.5_wp108 z2da(ji,jj) = ( u_ice(ji,jj) * umask(ji,jj,1) + u_ice(ji-1,jj) * umask(ji-1,jj,1) ) * 0.5_wp 109 z2db(ji,jj) = ( v_ice(ji,jj) * vmask(ji,jj,1) + v_ice(ji,jj-1) * vmask(ji,jj-1,1) ) * 0.5_wp 110 110 END DO 111 111 END DO -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limwri_dimg.h90
r4688 r5067 92 92 zinda = MAX( 0._wp , SIGN( 1._wp , ( 1.0 - frld(ji,jj) ) - 0.10 ) ) 93 93 zindb = zindh * zinda 94 ztmu = MAX( 0.5 * 1._wp , ( tmu(ji,jj) + tmu(ji+1,jj) + tmu(ji,jj+1) + tmu(ji+1,jj+1) ) )94 ztmu = MAX( 0.5 * 1._wp , ( umask(ji,jj,1) + umask(ji+1,jj,1) + umask(ji,jj+1,1) + umask(ji+1,jj+1,1) ) ) 95 95 zcmo(ji,jj,1) = ht_s (ji,jj,1) 96 96 zcmo(ji,jj,2) = ht_i (ji,jj,1) … … 99 99 zcmo(ji,jj,5) = sist (ji,jj) 100 100 zcmo(ji,jj,6) = fhtur (ji,jj) 101 zcmo(ji,jj,7) = zindb * ( u_ice(ji,jj ) * tmu(ji,jj ) + u_ice(ji+1,jj ) * tmu(ji+1,jj) &102 + u_ice(ji,jj+1) * tmu(ji,jj+1) + u_ice(ji+1,jj+1) * tmu(ji+1,jj+1) ) &101 zcmo(ji,jj,7) = zindb * ( u_ice(ji,jj ) * umask(ji,jj,1) + u_ice(ji+1,jj ) * umask(ji+1,jj,1) & 102 + u_ice(ji,jj+1) * umask(ji,jj+1,1) + u_ice(ji+1,jj+1) * umask(ji+1,jj+1,1) ) & 103 103 / ztmu 104 104 105 zcmo(ji,jj,8) = zindb * ( v_ice(ji,jj ) * tmu(ji,jj ) + v_ice(ji+1,jj ) * tmu(ji+1,jj) &106 + v_ice(ji,jj+1) * tmu(ji,jj+1) + v_ice(ji+1,jj+1) * tmu(ji+1,jj+1) ) &105 zcmo(ji,jj,8) = zindb * ( v_ice(ji,jj ) * umask(ji,jj,1) + v_ice(ji+1,jj ) * umask(ji+1,jj,1) & 106 + v_ice(ji,jj+1) * umask(ji,jj+1,1) + v_ice(ji+1,jj+1) * umask(ji+1,jj+1,1) ) & 107 107 / ztmu 108 108 zcmo(ji,jj,9) = sst_m(ji,jj) … … 135 135 zinda = MAX( 0._wp , SIGN( 1._wp , ( 1.0 - frld(ji,jj) ) - 0.10 ) ) 136 136 zindb = zindh * zinda 137 ztmu = MAX( 0.5 * 1._wp , ( tmu(ji,jj) + tmu(ji+1,jj) + tmu(ji,jj+1) + tmu(ji+1,jj+1) ) )137 ztmu = MAX( 0.5 * 1._wp , ( umask(ji,jj,1) + umask(ji+1,jj,1) + umask(ji,jj+1,1) + umask(ji+1,jj+1,1) ) ) 138 138 rcmoy(ji,jj,1) = ht_s (ji,jj,1) 139 139 rcmoy(ji,jj,2) = ht_i (ji,jj,1) … … 142 142 rcmoy(ji,jj,5) = sist (ji,jj) 143 143 rcmoy(ji,jj,6) = fhtur (ji,jj) 144 rcmoy(ji,jj,7) = zindb * ( u_ice(ji,jj ) * tmu(ji,jj ) + u_ice(ji+1,jj ) * tmu(ji+1,jj) &145 + u_ice(ji,jj+1) * tmu(ji,jj+1) + u_ice(ji+1,jj+1) * tmu(ji+1,jj+1) ) &144 rcmoy(ji,jj,7) = zindb * ( u_ice(ji,jj ) * umask(ji,jj,1) + u_ice(ji+1,jj ) * umask(ji+1,jj,1) & 145 + u_ice(ji,jj+1) * umask(ji,jj+1,1) + u_ice(ji+1,jj+1) * umask(ji+1,jj+1,1) ) & 146 146 / ztmu 147 147 148 rcmoy(ji,jj,8) = zindb * ( v_ice(ji,jj ) * tmu(ji,jj ) + v_ice(ji+1,jj ) * tmu(ji+1,jj) &149 + v_ice(ji,jj+1) * tmu(ji,jj+1) + v_ice(ji+1,jj+1) * tmu(ji+1,jj+1) ) &148 rcmoy(ji,jj,8) = zindb * ( v_ice(ji,jj ) * umask(ji,jj,1) + v_ice(ji+1,jj ) * umask(ji+1,jj,1) & 149 + v_ice(ji,jj+1) * umask(ji,jj+1,1) + v_ice(ji+1,jj+1) * umask(ji+1,jj+1,1) ) & 150 150 / ztmu 151 151 rcmoy(ji,jj,9) = sst_m(ji,jj) -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90
r5055 r5067 19 19 !!--------------------------- 20 20 ! !!! ** ice-thermo namelist (namicethd) ** 21 REAL(wp), PUBLIC :: hiclim!: minimum ice thickness21 REAL(wp), PUBLIC :: rn_himin !: minimum ice thickness 22 22 REAL(wp), PUBLIC :: parsub !: switch for snow sublimation or not 23 REAL(wp), PUBLIC :: maxfrazb!: maximum portion of frazil ice collecting at the ice bottom24 REAL(wp), PUBLIC :: vfrazb!: threshold drift speed for collection of bottom frazil ice25 REAL(wp), PUBLIC :: Cfrazb!: squeezing coefficient for collection of bottom frazil ice26 REAL(wp), PUBLIC :: hiccrit !: ice th. for lateral accretion in the NH (SH)(m)23 REAL(wp), PUBLIC :: rn_maxfrazb !: maximum portion of frazil ice collecting at the ice bottom 24 REAL(wp), PUBLIC :: rn_vfrazb !: threshold drift speed for collection of bottom frazil ice 25 REAL(wp), PUBLIC :: rn_Cfrazb !: squeezing coefficient for collection of bottom frazil ice 26 REAL(wp), PUBLIC :: rn_hnewice !: thickness for new ice formation (m) 27 27 28 INTEGER , PUBLIC :: fraz_swi !: use of frazil ice collection infunction of wind (1) or not (0)28 INTEGER , PUBLIC :: nn_frazil !: use of frazil ice collection as function of wind (1) or not (0) 29 29 30 30 !!----------------------------- -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90
r4990 r5067 162 162 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: gphit, gphiu !: latitude of t-, u-, v- and f-points (degre) 163 163 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: gphiv, gphif !: 164 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) :: e1t, e2t !: horizontal scale factorsat t-point (m)165 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) :: e1u, e2u !: horizontal scale factorsat u-point (m)166 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) :: e1v, e2v !: horizontal scale factorsat v-point (m)167 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) :: e1f, e2f !: horizontal scale factorsat f-point (m)164 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) :: e1t, e2t, r1_e1t, r1_e2t !: horizontal scale factors and inverse at t-point (m) 165 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) :: e1u, e2u, r1_e1u, r1_e2u !: horizontal scale factors and inverse at u-point (m) 166 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) :: e1v, e2v, r1_e1v, r1_e2v !: horizontal scale factors and inverse at v-point (m) 167 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) :: e1f, e2f, r1_e1f, r1_e2f !: horizontal scale factors and inverse at f-point (m) 168 168 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: e1e2t !: surface at t-point (m2) 169 169 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ff !: coriolis factor (2.*omega*sin(yphi) ) (s-1) … … 345 345 & tpol(jpiglo) , fpol(jpiglo) , STAT=ierr(2) ) 346 346 ! 347 ALLOCATE( glamt(jpi,jpj) , gphit(jpi,jpj) , e1t(jpi,jpj) , e2t(jpi,jpj) , & 348 & glamu(jpi,jpj) , gphiu(jpi,jpj) , e1u(jpi,jpj) , e2u(jpi,jpj) , & 349 & glamv(jpi,jpj) , gphiv(jpi,jpj) , e1v(jpi,jpj) , e2v(jpi,jpj) , e1e2t(jpi,jpj) , & 350 & glamf(jpi,jpj) , gphif(jpi,jpj) , e1f(jpi,jpj) , e2f(jpi,jpj) , ff (jpi,jpj) , STAT=ierr(3) ) 347 ALLOCATE( glamt(jpi,jpj) , gphit(jpi,jpj) , e1t(jpi,jpj) , e2t(jpi,jpj) , r1_e1t(jpi,jpj) , r1_e2t(jpi,jpj) , & 348 & glamu(jpi,jpj) , gphiu(jpi,jpj) , e1u(jpi,jpj) , e2u(jpi,jpj) , r1_e1u(jpi,jpj) , r1_e2u(jpi,jpj) , & 349 & glamv(jpi,jpj) , gphiv(jpi,jpj) , e1v(jpi,jpj) , e2v(jpi,jpj) , r1_e1v(jpi,jpj) , r1_e2v(jpi,jpj) , & 350 & glamf(jpi,jpj) , gphif(jpi,jpj) , e1f(jpi,jpj) , e2f(jpi,jpj) , r1_e1f(jpi,jpj) , r1_e2f(jpi,jpj) , & 351 & e1e2t(jpi,jpj) , ff (jpi,jpj) , STAT=ierr(3) ) 351 352 ! 352 353 ALLOCATE( gdep3w_0(jpi,jpj,jpk) , e3v_0(jpi,jpj,jpk) , e3f_0 (jpi,jpj,jpk) , & -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90
r4990 r5067 471 471 re2u_e1u(:,:) = e2u(:,:) / e1u(:,:) 472 472 re1v_e2v(:,:) = e1v(:,:) / e2v(:,:) 473 r1_e1t (:,:) = 1._wp / e1t(:,:) 474 r1_e1u (:,:) = 1._wp / e1u(:,:) 475 r1_e1v (:,:) = 1._wp / e1v(:,:) 476 r1_e1f (:,:) = 1._wp / e1f(:,:) 477 r1_e2t (:,:) = 1._wp / e2t(:,:) 478 r1_e2u (:,:) = 1._wp / e2u(:,:) 479 r1_e2v (:,:) = 1._wp / e2v(:,:) 480 r1_e2f (:,:) = 1._wp / e2f(:,:) 473 481 474 482 ! Control printing : Grid informations (if not restart) -
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90
r5059 r5067 236 236 CALL lim_thd( kt ) ! Ice thermodynamics 237 237 238 CALL lim_update2( kt ) ! Global variables update238 CALL lim_update2( kt ) ! Corrections 239 239 ! 240 240 CALL lim_sbc_flx( kt ) ! Update surface ocean mass, heat and salt fluxes … … 286 286 IF(lwm) CALL ctl_opn( numoni, 'output.namelist.ice', 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, 1 ) 287 287 288 CALL ice_run ! set some ice run parameters including jpl, nlay_i and nlay_s 288 CALL ice_run ! set some ice run parameters 289 ! 289 290 ! ! Allocate the ice arrays 290 291 ierr = ice_alloc () ! ice variables … … 303 304 & 'use more ocean levels or less ice/snow layers/categories.' ) 304 305 ! 306 CALL lim_itd_init ! ice thickness distribution initialization 305 307 ! 306 308 CALL lim_thd_init ! set ice thermodynics parameters … … 309 311 ! 310 312 CALL lim_msh ! ice mesh initialization 311 !312 CALL lim_itd_init ! ice thickness distribution initialization313 313 ! 314 314 CALL lim_itd_me_init ! ice thickness distribution initialization for mecanical deformation … … 352 352 INTEGER :: ios ! Local integer output status for namelist read 353 353 NAMELIST/namicerun/ jpl, nlay_i, nlay_s, cn_icerst_in, cn_icerst_out, & 354 & ln_limdyn, amax, ln_nicep, ln_limdiahsb, ln_limdiaout354 & ln_limdyn, rn_amax, ln_nicep, ln_limdiahsb, ln_limdiaout 355 355 !!------------------------------------------------------------------- 356 356 ! … … 369 369 WRITE(numout,*) 'ice_run : ice share parameters for dynamics/advection/thermo of sea-ice' 370 370 WRITE(numout,*) ' ~~~~~~' 371 WRITE(numout,*) ' number of ice categories = ', jpl372 WRITE(numout,*) ' number of ice layers = ', nlay_i373 WRITE(numout,*) ' number of snow layers = ', nlay_s371 WRITE(numout,*) ' number of ice categories = ', jpl 372 WRITE(numout,*) ' number of ice layers = ', nlay_i 373 WRITE(numout,*) ' number of snow layers = ', nlay_s 374 374 WRITE(numout,*) ' switch for ice dynamics (1) or not (0) ln_limdyn = ', ln_limdyn 375 WRITE(numout,*) ' maximum ice concentration = ', amax375 WRITE(numout,*) ' maximum ice concentration = ', rn_amax 376 376 WRITE(numout,*) ' Several ice points in the ice or not in ocean.output = ', ln_nicep 377 377 WRITE(numout,*) ' Diagnose heat/salt budget or not ln_limdiahsb = ', ln_limdiahsb … … 404 404 !!------------------------------------------------------------------- 405 405 INTEGER :: ios ! Local integer output status for namelist read 406 NAMELIST/namiceitd/ nn_ itdshp, rn_itmean406 NAMELIST/namiceitd/ nn_catbnd, rn_himean 407 407 ! 408 408 INTEGER :: jl ! dummy loop index … … 425 425 WRITE(numout,*) 'ice_itd : ice cat distribution' 426 426 WRITE(numout,*) ' ~~~~~~' 427 WRITE(numout,*) ' shape of ice categories distribution nn_ itdshp = ', nn_itdshp428 WRITE(numout,*) ' mean ice thickness in the domain (only active if nn_ itdshp=2) rn_itmean = ', rn_itmean427 WRITE(numout,*) ' shape of ice categories distribution nn_catbnd = ', nn_catbnd 428 WRITE(numout,*) ' mean ice thickness in the domain (only active if nn_catbnd=2) rn_himean = ', rn_himean 429 429 ENDIF 430 430 … … 438 438 hi_max(:) = 0._wp 439 439 440 SELECT CASE ( nn_ itdshp)440 SELECT CASE ( nn_catbnd ) 441 441 !---------------------- 442 442 CASE (1) ! tanh function (CICE) … … 456 456 zalpha = 0.05 ! exponent of the transform function 457 457 458 zhmax = 3.*rn_ itmean458 zhmax = 3.*rn_himean 459 459 460 460 DO jl = 1, jpl
Note: See TracChangeset
for help on using the changeset viewer.