Changeset 8515 for branches/2017/dev_r8183_ICEMODEL/NEMOGCM
- Timestamp:
- 2017-09-08T18:19:17+02:00 (7 years ago)
- Location:
- branches/2017/dev_r8183_ICEMODEL/NEMOGCM
- Files:
-
- 7 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/CONFIG/SHARED/namelist_ice_lim3_ref
r8514 r8515 1 1 !!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 2 !! LIM3namelist:2 !! ESIM namelist: 3 3 !! 1 - Generic parameters (namice_run) 4 4 !! 2 - Ice thickness discretization (namice_itd) … … 35 35 !------------------------------------------------------------------------------ 36 36 rn_himean = 2.0 ! expected domain-average ice thickness (m) 37 rn_himin = 0.1 ! minimum ice thickness (m) used in remapping 37 38 / 38 39 !------------------------------------------------------------------------------ … … 46 47 rn_perdg = 17.0 ! ridging work divided by pot. energy change in ridging 47 48 ! -- ice_rdgrft -- ! 48 rn_cs 49 rn_csrdg = 0.5 ! fraction of shearing energy contributing to ridging 49 50 ! -- ice_rdgrft_prep -- ! 50 51 ln_partf_lin = .false. ! Linear ridging participation function (Thorndike et al, 1975) … … 100 101 !------------------------------------------------------------------------------ 101 102 ln_icethd = .true. ! ice thermo (T) or not (F) => DO NOT TOUCH UNLESS U KNOW WHAT U DO 102 ! -- limthd_dif -- !103 ! -- icethd_dif -- ! 103 104 rn_kappa_i = 1.0 ! radiation attenuation coefficient in sea ice (m-1) 104 105 ln_cndi_U64 = .false. ! sea ice thermal conductivity: k = k0 + beta.S/T (Untersteiner, 1964) … … 107 108 rn_cnd_s = 0.31 ! thermal conductivity of the snow (0.31 W/m/K, Maykut and Untersteiner, 1971) 108 109 ! Obs: 0.1-0.5 (Lecomte et al, JAMES 2013) 109 ! -- limthd_dh -- !110 ! -- icethd_dh -- ! 110 111 ln_icedH = .true. ! activate ice thickness change from growing/melting (T) or not (F) => DO NOT TOUCH UNLESS U KNOW WHAT U DO 111 112 rn_blow_s = 0.66 ! mesure of snow blowing into the leads 112 113 ! = 1 => no snow blowing, < 1 => some snow blowing 113 ! -- limthd_da -- !114 ! -- icethd_da -- ! 114 115 ln_icedA = .true. ! activate lateral melting param. (T) or not (F) => DO NOT TOUCH UNLESS U KNOW WHAT U DO 115 116 rn_beta = 1.0 ! coef. beta for lateral melting param. Recommended range=[0.8-1.2] … … 120 121 ! => 6 vs 8m = +40% melting at the peak (A~0.5) 121 122 ! 10 vs 8m = -20% melting 122 ! -- limthd_lac -- !123 ! -- icethd_lac -- ! 123 124 ln_icedO = .true. ! activate ice growth in open-water (T) or not (F) => DO NOT TOUCH UNLESS U KNOW WHAT U DO 124 rn_hinew = 0.1 ! thickness for new ice formation in open water (m) 125 rn_hinew = 0.1 ! thickness for new ice formation in open water (m), must be larger than rn_hnewice 125 126 ln_frazil = .false. ! Frazil ice parameterization (ice collection as a function of wind) 126 127 rn_maxfraz = 1.0 ! maximum fraction of frazil ice collecting at the ice base 127 128 rn_vfraz = 0.417 ! thresold drift speed for frazil ice collecting at the ice bottom (m/s) 128 129 rn_Cfraz = 5.0 ! squeezing coefficient for frazil ice collecting at the ice bottom 129 ! -- limitd_th -- !130 rn_himin = 0.1 ! minimum ice thickness (m) used in remapping, must be smaller than rn_hnewice131 130 ! -- icestp -- ! 132 131 nn_iceflx = -1 ! Redistribute heat flux over ice categories … … 143 142 &namice_sal ! Ice salinity 144 143 !------------------------------------------------------------------------------ 145 ! -- limthd_sal -- !144 ! -- icethd_sal -- ! 146 145 ln_icedS = .true. ! activate gravity drainage and flushing (T) or not (F) => DO NOT TOUCH UNLESS U KNOW WHAT U DO 147 146 nn_icesal = 2 ! ice salinity option … … 171 170 &namice_ini ! Ice initialization 172 171 !------------------------------------------------------------------------------ 173 ! -- limistate -- !172 ! -- iceistate -- ! 174 173 ln_iceini = .true. ! activate ice initialization (T) or not (F) 175 174 ln_iceini_file = .false. ! netcdf file provided for initialization (T) or not (F) 176 rn_thres_sst = 2.0 ! max imum water temperature with initial ice (degC)175 rn_thres_sst = 2.0 ! max delta temp. above Tfreeze with initial ice = (sst - tfreeze) 177 176 rn_hts_ini_n = 0.3 ! initial real snow thickness (m), North 178 177 rn_hts_ini_s = 0.3 ! " " South -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceistate.F90
r8514 r8515 36 36 PRIVATE 37 37 38 PUBLIC ice_istate ! called by icestp.F90 38 PUBLIC ice_istate ! called by icestp.F90 39 PUBLIC ice_istate_init ! called by icestp.F90 39 40 40 41 INTEGER , PARAMETER :: jpfldi = 6 ! maximum number of files to read … … 102 103 103 104 IF(lwp) WRITE(numout,*) 104 IF(lwp) WRITE(numout,*) 'ice_istate 105 IF(lwp) WRITE(numout,*) '~~~~~~~~~~ 105 IF(lwp) WRITE(numout,*) 'ice_istate: sea-ice initialization ' 106 IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 106 107 107 108 !-------------------------------------------------------------------- … … 109 110 !-------------------------------------------------------------------- 110 111 ! 111 CALL ice_istate_init112 113 112 ! init surface temperature 114 113 DO jl = 1, jpl … … 154 153 !----------------------------- 155 154 ! assign initial thickness, concentration, snow depth and salinity to an hemisphere-dependent array 156 DO jj = 1, jpj 157 DO ji = 1, jpi 158 IF( ff_t(ji,jj) >= 0._wp ) THEN 159 zht_i_ini(ji,jj) = rn_hti_ini_n * zswitch(ji,jj) 160 zht_s_ini(ji,jj) = rn_hts_ini_n * zswitch(ji,jj) 161 zat_i_ini(ji,jj) = rn_ati_ini_n * zswitch(ji,jj) 162 zts_u_ini(ji,jj) = rn_tmi_ini_n * zswitch(ji,jj) 163 zsm_i_ini(ji,jj) = rn_smi_ini_n * zswitch(ji,jj) 164 ztm_i_ini(ji,jj) = rn_tmi_ini_n * zswitch(ji,jj) 165 ELSE 166 zht_i_ini(ji,jj) = rn_hti_ini_s * zswitch(ji,jj) 167 zht_s_ini(ji,jj) = rn_hts_ini_s * zswitch(ji,jj) 168 zat_i_ini(ji,jj) = rn_ati_ini_s * zswitch(ji,jj) 169 zts_u_ini(ji,jj) = rn_tmi_ini_s * zswitch(ji,jj) 170 zsm_i_ini(ji,jj) = rn_smi_ini_s * zswitch(ji,jj) 171 ztm_i_ini(ji,jj) = rn_tmi_ini_s * zswitch(ji,jj) 172 ENDIF 173 END DO 174 END DO 155 WHERE( ff_t(:,:) >= 0._wp ) 156 zht_i_ini(:,:) = rn_hti_ini_n * zswitch(:,:) 157 zht_s_ini(:,:) = rn_hts_ini_n * zswitch(:,:) 158 zat_i_ini(:,:) = rn_ati_ini_n * zswitch(:,:) 159 zts_u_ini(:,:) = rn_tmi_ini_n * zswitch(:,:) 160 zsm_i_ini(:,:) = rn_smi_ini_n * zswitch(:,:) 161 ztm_i_ini(:,:) = rn_tmi_ini_n * zswitch(:,:) 162 ELSEWHERE 163 zht_i_ini(:,:) = rn_hti_ini_s * zswitch(:,:) 164 zht_s_ini(:,:) = rn_hts_ini_s * zswitch(:,:) 165 zat_i_ini(:,:) = rn_ati_ini_s * zswitch(:,:) 166 zts_u_ini(:,:) = rn_tmi_ini_s * zswitch(:,:) 167 zsm_i_ini(:,:) = rn_smi_ini_s * zswitch(:,:) 168 ztm_i_ini(:,:) = rn_tmi_ini_s * zswitch(:,:) 169 END WHERE 175 170 ! 176 171 ENDIF ! ln_iceini_file … … 554 549 slf_i(jp_ati) = sn_ati ; slf_i(jp_tsu) = sn_tsu 555 550 slf_i(jp_tmi) = sn_tmi ; slf_i(jp_smi) = sn_smi 556 557 ! Define the initial parameters 558 ! ------------------------- 559 560 IF(lwp) THEN 551 ! 552 ! 553 IF(lwp) THEN ! control print 561 554 WRITE(numout,*) 562 555 WRITE(numout,*) 'ice_istate_init: ice parameters inititialisation ' 563 556 WRITE(numout,*) '~~~~~~~~~~~~~~~' 564 557 WRITE(numout,*) ' Namelist namice_ini' 565 WRITE(numout,*) ' initialization with ice (T) or not (F) ln_iceini = ', ln_iceini566 WRITE(numout,*) ' ice initialization from a netcdf file ln_iceini_file = ', ln_iceini_file567 WRITE(numout,*) ' threshold water temp. for initial sea-icern_thres_sst = ', rn_thres_sst568 WRITE(numout,*) ' initial snow thickness in the north rn_hts_ini_n = ', rn_hts_ini_n569 WRITE(numout,*) ' initial snow thickness in the south rn_hts_ini_s = ', rn_hts_ini_s570 WRITE(numout,*) ' initial ice thickness in the north rn_hti_ini_n = ', rn_hti_ini_n571 WRITE(numout,*) ' initial ice thickness in the south rn_hti_ini_s = ', rn_hti_ini_s572 WRITE(numout,*) ' initial ice concentr. in the north rn_ati_ini_n = ', rn_ati_ini_n573 WRITE(numout,*) ' initial ice concentr. in the north rn_ati_ini_s = ', rn_ati_ini_s574 WRITE(numout,*) ' initial ice salinity in the north rn_smi_ini_n = ', rn_smi_ini_n575 WRITE(numout,*) ' initial ice salinity in the south rn_smi_ini_s = ', rn_smi_ini_s576 WRITE(numout,*) ' initial ice/snw temp in the north rn_tmi_ini_n = ', rn_tmi_ini_n577 WRITE(numout,*) ' initial ice/snw temp in the south rn_tmi_ini_s = ', rn_tmi_ini_s558 WRITE(numout,*) ' initialization with ice (T) or not (F) ln_iceini = ', ln_iceini 559 WRITE(numout,*) ' ice initialization from a netcdf file ln_iceini_file = ', ln_iceini_file 560 WRITE(numout,*) ' max delta ocean temp. above Tfreeze with initial ice rn_thres_sst = ', rn_thres_sst 561 WRITE(numout,*) ' initial snow thickness in the north rn_hts_ini_n = ', rn_hts_ini_n 562 WRITE(numout,*) ' initial snow thickness in the south rn_hts_ini_s = ', rn_hts_ini_s 563 WRITE(numout,*) ' initial ice thickness in the north rn_hti_ini_n = ', rn_hti_ini_n 564 WRITE(numout,*) ' initial ice thickness in the south rn_hti_ini_s = ', rn_hti_ini_s 565 WRITE(numout,*) ' initial ice concentr. in the north rn_ati_ini_n = ', rn_ati_ini_n 566 WRITE(numout,*) ' initial ice concentr. in the north rn_ati_ini_s = ', rn_ati_ini_s 567 WRITE(numout,*) ' initial ice salinity in the north rn_smi_ini_n = ', rn_smi_ini_n 568 WRITE(numout,*) ' initial ice salinity in the south rn_smi_ini_s = ', rn_smi_ini_s 569 WRITE(numout,*) ' initial ice/snw temp in the north rn_tmi_ini_n = ', rn_tmi_ini_n 570 WRITE(numout,*) ' initial ice/snw temp in the south rn_tmi_ini_s = ', rn_tmi_ini_s 578 571 ENDIF 579 572 -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceitd.F90
r8514 r8515 650 650 REAL(wp) :: zhmax, znum, zden, zalpha ! - - 651 651 !! 652 NAMELIST/namice_itd/ rn_himean 652 NAMELIST/namice_itd/ rn_himean, rn_himin 653 653 !!------------------------------------------------------------------ 654 654 ! … … 668 668 WRITE(numout,*) ' Namelist namice_itd : ' 669 669 WRITE(numout,*) ' mean ice thickness in the domain rn_himean = ', rn_himean 670 WRITE(numout,*) ' minimum ice thickness rn_himin = ', rn_himin 670 671 ENDIF 671 672 ! -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerdgrft.F90
r8514 r8515 52 52 ! 53 53 ! ** namelist (namice_rdgrft) ** 54 REAL(wp) :: rn_cs 54 REAL(wp) :: rn_csrdg ! fraction of shearing energy contributing to ridging 55 55 LOGICAL :: ln_partf_lin ! participation function linear (Thorndike et al. (1975)) 56 56 REAL(wp) :: rn_gstar ! fractional area of young ice contributing to ridging … … 170 170 ! (thick, newly ridged ice). 171 171 172 closing_net(ji,jj) = rn_cs * 0.5_wp * ( delta_i(ji,jj) - ABS( divu_i(ji,jj) ) ) - MIN( divu_i(ji,jj), 0._wp )172 closing_net(ji,jj) = rn_csrdg * 0.5_wp * ( delta_i(ji,jj) - ABS( divu_i(ji,jj) ) ) - MIN( divu_i(ji,jj), 0._wp ) 173 173 174 174 ! 2.2 divu_adv … … 867 867 NAMELIST/namice_rdgrft/ ln_str_H79, rn_pstar, rn_crhg, & 868 868 & ln_str_R75, rn_perdg, & 869 & rn_cs 869 & rn_csrdg , & 870 870 & ln_partf_lin, rn_gstar, & 871 871 & ln_partf_exp, rn_astar, & … … 893 893 WRITE(numout,*) ' ice strength parameterization Rothrock (1975) ln_str_R75 = ', ln_str_R75 894 894 WRITE(numout,*) ' Ratio of ridging work to PotEner change in ridging rn_perdg = ', rn_perdg 895 WRITE(numout,*) ' Fraction of shear energy contributing to ridging rn_cs = ', rn_cs895 WRITE(numout,*) ' Fraction of shear energy contributing to ridging rn_csrdg = ', rn_csrdg 896 896 WRITE(numout,*) ' linear ridging participation function ln_partf_lin = ', ln_partf_lin 897 897 WRITE(numout,*) ' Fraction of ice coverage contributing to ridging rn_gstar = ', rn_gstar -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerhg_evp.F90
r8512 r8515 52 52 CONTAINS 53 53 54 SUBROUTINE ice_rhg_evp( kt, stress1_i, stress2_i, stress12_i, u_ice, v_ice, shear_i, divu_i,delta_i )54 SUBROUTINE ice_rhg_evp( kt, pstress1_i, pstress2_i, pstress12_i, pu_ice, pv_ice, pshear_i, pdivu_i, pdelta_i ) 55 55 !!------------------------------------------------------------------- 56 56 !! *** SUBROUTINE ice_rhg_evp *** … … 108 108 !!------------------------------------------------------------------- 109 109 INTEGER, INTENT(in) :: kt ! time step 110 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: stress1_i, stress2_i,stress12_i111 REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: u_ice, v_ice, shear_i, divu_i,delta_i110 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pstress1_i, pstress2_i, pstress12_i 111 REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: pu_ice, pv_ice, pshear_i, pdivu_i, pdelta_i 112 112 !! 113 113 INTEGER :: ji, jj ! dummy loop indices … … 247 247 248 248 ! Initialise stress tensor 249 zs1 (:,:) = stress1_i (:,:)250 zs2 (:,:) = stress2_i (:,:)251 zs12(:,:) = stress12_i(:,:)249 zs1 (:,:) = pstress1_i (:,:) 250 zs2 (:,:) = pstress2_i (:,:) 251 zs12(:,:) = pstress12_i(:,:) 252 252 253 253 ! Ice strength … … 340 340 IF(ln_ctl) THEN ! Convergence test 341 341 DO jj = 1, jpjm1 342 zu_ice(:,jj) = u_ice(:,jj) ! velocity at previous time step343 zv_ice(:,jj) = v_ice(:,jj)342 zu_ice(:,jj) = pu_ice(:,jj) ! velocity at previous time step 343 zv_ice(:,jj) = pv_ice(:,jj) 344 344 END DO 345 345 ENDIF … … 350 350 351 351 ! shear at F points 352 zds(ji,jj) = ( ( u_ice(ji,jj+1) * r1_e1u(ji,jj+1) -u_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj) &353 & + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) -v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) &352 zds(ji,jj) = ( ( pu_ice(ji,jj+1) * r1_e1u(ji,jj+1) - pu_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj) & 353 & + ( pv_ice(ji+1,jj) * r1_e2v(ji+1,jj) - pv_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) & 354 354 & ) * r1_e1e2f(ji,jj) * zfmask(ji,jj) 355 355 … … 367 367 368 368 ! divergence at T points 369 zdiv = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) *u_ice(ji-1,jj) &370 & + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) *v_ice(ji,jj-1) &369 zdiv = ( e2u(ji,jj) * pu_ice(ji,jj) - e2u(ji-1,jj) * pu_ice(ji-1,jj) & 370 & + e1v(ji,jj) * pv_ice(ji,jj) - e1v(ji,jj-1) * pv_ice(ji,jj-1) & 371 371 & ) * r1_e1e2t(ji,jj) 372 372 zdiv2 = zdiv * zdiv 373 373 374 374 ! tension at T points 375 zdt = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) -u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj) &376 & - ( v_ice(ji,jj) * r1_e1v(ji,jj) -v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj) &375 zdt = ( ( pu_ice(ji,jj) * r1_e2u(ji,jj) - pu_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj) & 376 & - ( pv_ice(ji,jj) * r1_e1v(ji,jj) - pv_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj) & 377 377 & ) * r1_e1e2t(ji,jj) 378 378 zdt2 = zdt * zdt … … 426 426 ! 427 427 ! !--- u_ice at V point 428 u_iceV(ji,jj) = 0.5_wp * ( ( u_ice(ji,jj ) +u_ice(ji-1,jj ) ) * e2t(ji,jj+1) &429 & + ( u_ice(ji,jj+1) +u_ice(ji-1,jj+1) ) * e2t(ji,jj ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1)428 u_iceV(ji,jj) = 0.5_wp * ( ( pu_ice(ji,jj ) + pu_ice(ji-1,jj ) ) * e2t(ji,jj+1) & 429 & + ( pu_ice(ji,jj+1) + pu_ice(ji-1,jj+1) ) * e2t(ji,jj ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1) 430 430 ! 431 431 ! !--- v_ice at U point 432 v_iceU(ji,jj) = 0.5_wp * ( ( v_ice(ji ,jj) +v_ice(ji ,jj-1) ) * e1t(ji+1,jj) &433 & + ( v_ice(ji+1,jj) +v_ice(ji+1,jj-1) ) * e1t(ji ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1)432 v_iceU(ji,jj) = 0.5_wp * ( ( pv_ice(ji ,jj) + pv_ice(ji ,jj-1) ) * e1t(ji+1,jj) & 433 & + ( pv_ice(ji+1,jj) + pv_ice(ji+1,jj-1) ) * e1t(ji ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1) 434 434 ! 435 435 END DO … … 444 444 DO ji = fs_2, fs_jpim1 445 445 ! !--- tau_io/(v_oce - v_ice) 446 zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice(ji,jj) - v_oce (ji,jj) ) &446 zTauO = zaV(ji,jj) * zrhoco * SQRT( ( pv_ice(ji,jj) - v_oce (ji,jj) ) * ( pv_ice(ji,jj) - v_oce (ji,jj) ) & 447 447 & + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) ) 448 448 ! !--- Ocean-to-Ice stress 449 ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) )449 ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - pv_ice(ji,jj) ) 450 450 ! 451 451 ! !--- tau_bottom/v_ice 452 zvel = MAX( zepsi, SQRT( v_ice(ji,jj) *v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) )452 zvel = MAX( zepsi, SQRT( pv_ice(ji,jj) * pv_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) ) 453 453 zTauB = - tau_icebfr(ji,jj) / zvel 454 454 ! 455 455 ! !--- Coriolis at V-points (energy conserving formulation) 456 456 zCory(ji,jj) = - 0.25_wp * r1_e2v(ji,jj) * & 457 & ( zmf(ji,jj ) * ( e2u(ji,jj ) * u_ice(ji,jj ) + e2u(ji-1,jj ) *u_ice(ji-1,jj ) ) &458 & + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) *u_ice(ji-1,jj+1) ) )457 & ( zmf(ji,jj ) * ( e2u(ji,jj ) * pu_ice(ji,jj ) + e2u(ji-1,jj ) * pu_ice(ji-1,jj ) ) & 458 & + zmf(ji,jj+1) * ( e2u(ji,jj+1) * pu_ice(ji,jj+1) + e2u(ji-1,jj+1) * pu_ice(ji-1,jj+1) ) ) 459 459 ! 460 460 ! !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io … … 465 465 ! 466 466 ! !--- ice velocity using implicit formulation (cf Madec doc & Bouillon 2009) 467 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * v_ice(ji,jj)& ! previous velocity468 & + zTauE + zTauO * v_ice(ji,jj)& ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)469 & ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )& ! m/dt + tau_io(only ice part) + landfast470 & + ( 1._wp - rswitch ) *v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0471 & ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )& ! v_ice = v_oce if mass < zmmin472 & ) * zmaskV(ji,jj)467 pv_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * pv_ice(ji,jj) & ! previous velocity 468 & + zTauE + zTauO * pv_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 469 & ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 470 & + ( 1._wp - rswitch ) * pv_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 471 & ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) ) & ! v_ice = v_oce if mass < zmmin 472 & ) * zmaskV(ji,jj) 473 473 ! 474 474 ! Bouillon 2013 475 !! v_ice(ji,jj) = ( zmV_t(ji,jj) * ( zbeta * v_ice(ji,jj) +v_ice_b(ji,jj) ) &475 !!pv_ice(ji,jj) = ( zmV_t(ji,jj) * ( zbeta * pv_ice(ji,jj) + pv_ice_b(ji,jj) ) & 476 476 !! & + zfV(ji,jj) + zCory(ji,jj) + zTauV_ia(ji,jj) + zTauO * v_oce(ji,jj) + zspgV(ji,jj) & 477 477 !! & ) / MAX( zmV_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchV(ji,jj) … … 479 479 END DO 480 480 END DO 481 CALL lbc_lnk( v_ice, 'V', -1. )481 CALL lbc_lnk( pv_ice, 'V', -1. ) 482 482 ! 483 483 #if defined key_agrif … … 491 491 492 492 ! tau_io/(u_oce - u_ice) 493 zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice(ji,jj) - u_oce (ji,jj) ) &493 zTauO = zaU(ji,jj) * zrhoco * SQRT( ( pu_ice(ji,jj) - u_oce (ji,jj) ) * ( pu_ice(ji,jj) - u_oce (ji,jj) ) & 494 494 & + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) ) 495 495 496 496 ! Ocean-to-Ice stress 497 ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) )497 ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - pu_ice(ji,jj) ) 498 498 499 499 ! tau_bottom/u_ice 500 zvel = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) *u_ice(ji,jj) ) )500 zvel = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + pu_ice(ji,jj) * pu_ice(ji,jj) ) ) 501 501 zTauB = - tau_icebfr(ji,jj) / zvel 502 502 503 503 ! Coriolis at U-points (energy conserving formulation) 504 504 zCorx(ji,jj) = 0.25_wp * r1_e1u(ji,jj) * & 505 & ( zmf(ji ,jj) * ( e1v(ji ,jj) * v_ice(ji ,jj) + e1v(ji ,jj-1) *v_ice(ji ,jj-1) ) &506 & + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) *v_ice(ji+1,jj-1) ) )505 & ( zmf(ji ,jj) * ( e1v(ji ,jj) * pv_ice(ji ,jj) + e1v(ji ,jj-1) * pv_ice(ji ,jj-1) ) & 506 & + zmf(ji+1,jj) * ( e1v(ji+1,jj) * pv_ice(ji+1,jj) + e1v(ji+1,jj-1) * pv_ice(ji+1,jj-1) ) ) 507 507 508 508 ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io … … 513 513 514 514 ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009) 515 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * u_ice(ji,jj)& ! previous velocity516 & + zTauE + zTauO * u_ice(ji,jj)& ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)517 & ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )& ! m/dt + tau_io(only ice part) + landfast518 & + ( 1._wp - rswitch ) *u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0519 & ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )& ! v_ice = v_oce if mass < zmmin520 & ) * zmaskU(ji,jj)515 pu_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * pu_ice(ji,jj) & ! previous velocity 516 & + zTauE + zTauO * pu_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 517 & ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 518 & + ( 1._wp - rswitch ) * pu_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 519 & ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) ) & ! v_ice = v_oce if mass < zmmin 520 & ) * zmaskU(ji,jj) 521 521 ! Bouillon 2013 522 !! u_ice(ji,jj) = ( zmU_t(ji,jj) * ( zbeta * u_ice(ji,jj) +u_ice_b(ji,jj) ) &522 !!pu_ice(ji,jj) = ( zmU_t(ji,jj) * ( zbeta * pu_ice(ji,jj) + pu_ice_b(ji,jj) ) & 523 523 !! & + zfU(ji,jj) + zCorx(ji,jj) + zTauU_ia(ji,jj) + zTauO * u_oce(ji,jj) + zspgU(ji,jj) & 524 524 !! & ) / MAX( zmU_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchU(ji,jj) 525 525 END DO 526 526 END DO 527 CALL lbc_lnk( u_ice, 'U', -1. )527 CALL lbc_lnk( pu_ice, 'U', -1. ) 528 528 ! 529 529 #if defined key_agrif … … 539 539 540 540 ! tau_io/(u_oce - u_ice) 541 zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice(ji,jj) - u_oce (ji,jj) ) &541 zTauO = zaU(ji,jj) * zrhoco * SQRT( ( pu_ice(ji,jj) - u_oce (ji,jj) ) * ( pu_ice(ji,jj) - u_oce (ji,jj) ) & 542 542 & + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) ) 543 543 544 544 ! Ocean-to-Ice stress 545 ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) )545 ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - pu_ice(ji,jj) ) 546 546 547 547 ! tau_bottom/u_ice 548 zvel = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) *u_ice(ji,jj) ) )548 zvel = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + pu_ice(ji,jj) * pu_ice(ji,jj) ) ) 549 549 zTauB = - tau_icebfr(ji,jj) / zvel 550 550 551 551 ! Coriolis at U-points (energy conserving formulation) 552 552 zCorx(ji,jj) = 0.25_wp * r1_e1u(ji,jj) * & 553 & ( zmf(ji ,jj) * ( e1v(ji ,jj) * v_ice(ji ,jj) + e1v(ji ,jj-1) *v_ice(ji ,jj-1) ) &554 & + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) *v_ice(ji+1,jj-1) ) )553 & ( zmf(ji ,jj) * ( e1v(ji ,jj) * pv_ice(ji ,jj) + e1v(ji ,jj-1) * pv_ice(ji ,jj-1) ) & 554 & + zmf(ji+1,jj) * ( e1v(ji+1,jj) * pv_ice(ji+1,jj) + e1v(ji+1,jj-1) * pv_ice(ji+1,jj-1) ) ) 555 555 556 556 ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io … … 561 561 562 562 ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009) 563 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * u_ice(ji,jj)& ! previous velocity564 & + zTauE + zTauO * u_ice(ji,jj)& ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)565 & ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )& ! m/dt + tau_io(only ice part) + landfast566 & + ( 1._wp - rswitch ) *u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0567 & ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )& ! v_ice = v_oce if mass < zmmin568 & ) * zmaskU(ji,jj)563 pu_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * pu_ice(ji,jj) & ! previous velocity 564 & + zTauE + zTauO * pu_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 565 & ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 566 & + ( 1._wp - rswitch ) * pu_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 567 & ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) ) & ! v_ice = v_oce if mass < zmmin 568 & ) * zmaskU(ji,jj) 569 569 ! Bouillon 2013 570 !! u_ice(ji,jj) = ( zmU_t(ji,jj) * ( zbeta * u_ice(ji,jj) +u_ice_b(ji,jj) ) &570 !!pu_ice(ji,jj) = ( zmU_t(ji,jj) * ( zbeta * pu_ice(ji,jj) + pu_ice_b(ji,jj) ) & 571 571 !! & + zfU(ji,jj) + zCorx(ji,jj) + zTauU_ia(ji,jj) + zTauO * u_oce(ji,jj) + zspgU(ji,jj) & 572 572 !! & ) / MAX( zmU_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchU(ji,jj) 573 573 END DO 574 574 END DO 575 CALL lbc_lnk( u_ice, 'U', -1. )575 CALL lbc_lnk( pu_ice, 'U', -1. ) 576 576 ! 577 577 #if defined key_agrif … … 585 585 586 586 ! tau_io/(v_oce - v_ice) 587 zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice(ji,jj) - v_oce (ji,jj) ) &587 zTauO = zaV(ji,jj) * zrhoco * SQRT( ( pv_ice(ji,jj) - v_oce (ji,jj) ) * ( pv_ice(ji,jj) - v_oce (ji,jj) ) & 588 588 & + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) ) 589 589 590 590 ! Ocean-to-Ice stress 591 ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) )591 ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - pv_ice(ji,jj) ) 592 592 593 593 ! tau_bottom/v_ice 594 zvel = MAX( zepsi, SQRT( v_ice(ji,jj) *v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) )594 zvel = MAX( zepsi, SQRT( pv_ice(ji,jj) * pv_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) ) 595 595 ztauB = - tau_icebfr(ji,jj) / zvel 596 596 597 597 ! Coriolis at V-points (energy conserving formulation) 598 598 zCory(ji,jj) = - 0.25_wp * r1_e2v(ji,jj) * & 599 & ( zmf(ji,jj ) * ( e2u(ji,jj ) * u_ice(ji,jj ) + e2u(ji-1,jj ) *u_ice(ji-1,jj ) ) &600 & + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) *u_ice(ji-1,jj+1) ) )599 & ( zmf(ji,jj ) * ( e2u(ji,jj ) * pu_ice(ji,jj ) + e2u(ji-1,jj ) * pu_ice(ji-1,jj ) ) & 600 & + zmf(ji,jj+1) * ( e2u(ji,jj+1) * pu_ice(ji,jj+1) + e2u(ji-1,jj+1) * pu_ice(ji-1,jj+1) ) ) 601 601 602 602 ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io … … 607 607 608 608 ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009) 609 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * v_ice(ji,jj)& ! previous velocity610 & + zTauE + zTauO * v_ice(ji,jj)& ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)611 & ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )& ! m/dt + tau_io(only ice part) + landfast612 & + ( 1._wp - rswitch ) *v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0613 & ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )& ! v_ice = v_oce if mass < zmmin614 & ) * zmaskV(ji,jj)609 pv_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * pv_ice(ji,jj) & ! previous velocity 610 & + zTauE + zTauO * pv_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 611 & ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 612 & + ( 1._wp - rswitch ) * pv_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 613 & ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) ) & ! v_ice = v_oce if mass < zmmin 614 & ) * zmaskV(ji,jj) 615 615 ! Bouillon 2013 616 !! v_ice(ji,jj) = ( zmV_t(ji,jj) * ( zbeta * v_ice(ji,jj) +v_ice_b(ji,jj) ) &616 !!pv_ice(ji,jj) = ( zmV_t(ji,jj) * ( zbeta * pv_ice(ji,jj) + pv_ice_b(ji,jj) ) & 617 617 !! & + zfV(ji,jj) + zCory(ji,jj) + zTauV_ia(ji,jj) + zTauO * v_oce(ji,jj) + zspgV(ji,jj) & 618 618 !! & ) / MAX( zmV_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchV(ji,jj) 619 619 END DO 620 620 END DO 621 CALL lbc_lnk( v_ice, 'V', -1. )621 CALL lbc_lnk( pv_ice, 'V', -1. ) 622 622 ! 623 623 #if defined key_agrif … … 631 631 IF(ln_ctl) THEN ! Convergence test 632 632 DO jj = 2 , jpjm1 633 zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ), ABS(v_ice(:,jj) - zv_ice(:,jj) ) )633 zresr(:,jj) = MAX( ABS( pu_ice(:,jj) - zu_ice(:,jj) ), ABS( pv_ice(:,jj) - zv_ice(:,jj) ) ) 634 634 END DO 635 635 zresm = MAXVAL( zresr( 1:jpi, 2:jpjm1 ) ) … … 648 648 649 649 ! shear at F points 650 zds(ji,jj) = ( ( u_ice(ji,jj+1) * r1_e1u(ji,jj+1) -u_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj) &651 & + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) -v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) &650 zds(ji,jj) = ( ( pu_ice(ji,jj+1) * r1_e1u(ji,jj+1) - pu_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj) & 651 & + ( pv_ice(ji+1,jj) * r1_e2v(ji+1,jj) - pv_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) & 652 652 & ) * r1_e1e2f(ji,jj) * zfmask(ji,jj) 653 653 … … 660 660 661 661 ! tension**2 at T points 662 zdt = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) -u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj) &663 & - ( v_ice(ji,jj) * r1_e1v(ji,jj) -v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj) &662 zdt = ( ( pu_ice(ji,jj) * r1_e2u(ji,jj) - pu_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj) & 663 & - ( pv_ice(ji,jj) * r1_e1v(ji,jj) - pv_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj) & 664 664 & ) * r1_e1e2t(ji,jj) 665 665 zdt2 = zdt * zdt … … 671 671 672 672 ! shear at T points 673 shear_i(ji,jj) = SQRT( zdt2 + zds2 )673 pshear_i(ji,jj) = SQRT( zdt2 + zds2 ) 674 674 675 675 ! divergence at T points 676 divu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) *u_ice(ji-1,jj) &677 & + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) *v_ice(ji,jj-1) &678 & ) * r1_e1e2t(ji,jj)676 pdivu_i(ji,jj) = ( e2u(ji,jj) * pu_ice(ji,jj) - e2u(ji-1,jj) * pu_ice(ji-1,jj) & 677 & + e1v(ji,jj) * pv_ice(ji,jj) - e1v(ji,jj-1) * pv_ice(ji,jj-1) & 678 & ) * r1_e1e2t(ji,jj) 679 679 680 680 ! delta at T points 681 zdelta = SQRT( divu_i(ji,jj) *divu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 )681 zdelta = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 ) 682 682 rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zdelta ) ) ! 0 if delta=0 683 delta_i(ji,jj) = zdelta + rn_creepl * rswitch683 pdelta_i(ji,jj) = zdelta + rn_creepl * rswitch 684 684 685 685 END DO 686 686 END DO 687 CALL lbc_lnk_multi( shear_i, 'T', 1., divu_i, 'T', 1.,delta_i, 'T', 1. )687 CALL lbc_lnk_multi( pshear_i, 'T', 1., pdivu_i, 'T', 1., pdelta_i, 'T', 1. ) 688 688 689 689 ! --- Store the stress tensor for the next time step --- ! 690 stress1_i (:,:) = zs1 (:,:)691 stress2_i (:,:) = zs2 (:,:)692 stress12_i(:,:) = zs12(:,:)690 pstress1_i (:,:) = zs1 (:,:) 691 pstress2_i (:,:) = zs2 (:,:) 692 pstress12_i(:,:) = zs12(:,:) 693 693 ! 694 694 … … 703 703 704 704 ! --- divergence, shear and strength --- ! 705 IF( iom_use('idive') ) CALL iom_put( "idive" , divu_i(:,:) * zswi(:,:) ) ! divergence706 IF( iom_use('ishear') ) CALL iom_put( "ishear" , shear_i(:,:) * zswi(:,:) ) ! shear705 IF( iom_use('idive') ) CALL iom_put( "idive" , pdivu_i(:,:) * zswi(:,:) ) ! divergence 706 IF( iom_use('ishear') ) CALL iom_put( "ishear" , pshear_i(:,:) * zswi(:,:) ) ! shear 707 707 IF( iom_use('icestr') ) CALL iom_put( "icestr" , strength(:,:) * zswi(:,:) ) ! Ice strength 708 708 … … 714 714 DO jj = 2, jpjm1 715 715 DO ji = 2, jpim1 716 zdum1 = ( zswi(ji-1,jj) * stress12_i(ji-1,jj) + zswi(ji ,jj-1) *stress12_i(ji ,jj-1) + & ! stress12_i at T-point717 & zswi(ji ,jj) * stress12_i(ji ,jj) + zswi(ji-1,jj-1) *stress12_i(ji-1,jj-1) ) &718 & / MAX( 1._wp, zswi(ji-1,jj) + zswi(ji,jj-1) + zswi(ji,jj) + zswi(ji-1,jj-1) )719 720 zshear = SQRT( stress2_i(ji,jj) *stress2_i(ji,jj) + 4._wp * zdum1 * zdum1 ) ! shear stress716 zdum1 = ( zswi(ji-1,jj) * pstress12_i(ji-1,jj) + zswi(ji ,jj-1) * pstress12_i(ji ,jj-1) + & ! stress12_i at T-point 717 & zswi(ji ,jj) * pstress12_i(ji ,jj) + zswi(ji-1,jj-1) * pstress12_i(ji-1,jj-1) ) & 718 & / MAX( 1._wp, zswi(ji-1,jj) + zswi(ji,jj-1) + zswi(ji,jj) + zswi(ji-1,jj-1) ) 719 720 zshear = SQRT( pstress2_i(ji,jj) * pstress2_i(ji,jj) + 4._wp * zdum1 * zdum1 ) ! shear stress 721 721 722 722 zdum2 = zswi(ji,jj) / MAX( 1._wp, strength(ji,jj) ) 723 723 724 !! zsig1(ji,jj) = 0.5_wp * zdum2 * ( stress1_i(ji,jj) + zshear ) ! principal stress (y-direction, see Hunke & Dukowicz 2002)725 !! zsig2(ji,jj) = 0.5_wp * zdum2 * ( stress1_i(ji,jj) - zshear ) ! principal stress (x-direction, see Hunke & Dukowicz 2002)726 !! zsig3(ji,jj) = zdum2**2 * ( ( stress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 ) ! quadratic relation linking compressive stress to shear stress724 !! zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) + zshear ) ! principal stress (y-direction, see Hunke & Dukowicz 2002) 725 !! zsig2(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) - zshear ) ! principal stress (x-direction, see Hunke & Dukowicz 2002) 726 !! zsig3(ji,jj) = zdum2**2 * ( ( pstress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 ) ! quadratic relation linking compressive stress to shear stress 727 727 !! ! (scheme converges if this value is ~1, see Bouillon et al 2009 (eq. 11)) 728 zsig1(ji,jj) = 0.5_wp * zdum2 * ( stress1_i(ji,jj) ) ! compressive stress, see Bouillon et al. 2015728 zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) ) ! compressive stress, see Bouillon et al. 2015 729 729 zsig2(ji,jj) = 0.5_wp * zdum2 * ( zshear ) ! shear stress 730 zsig3(ji,jj) = zdum2**2 * ( ( stress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 )730 zsig3(ji,jj) = zdum2**2 * ( ( pstress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 ) 731 731 END DO 732 732 END DO … … 773 773 774 774 ! 2D ice mass, snow mass, area transport arrays (X, Y) 775 zfac_x = 0.5 * u_ice(ji,jj) * e2u(ji,jj) * rswitch776 zfac_y = 0.5 * v_ice(ji,jj) * e1v(ji,jj) * rswitch775 zfac_x = 0.5 * pu_ice(ji,jj) * e2u(ji,jj) * rswitch 776 zfac_y = 0.5 * pv_ice(ji,jj) * e1v(ji,jj) * rswitch 777 777 778 778 zdiag_xmtrp_ice(ji,jj) = rhoic * zfac_x * ( vt_i(ji+1,jj) + vt_i(ji,jj) ) ! ice mass transport, X-component -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icestp.F90
r8514 r8515 285 285 ! ! Initial sea-ice state 286 286 IF( .NOT. ln_rstart ) THEN ! start from rest: sea-ice deduced from sst 287 CALL ice_istate_init 287 288 CALL ice_istate 288 289 ELSE ! start from a restart file … … 301 302 tn_ice(:,:,:) = t_su(:,:,:) ! initialisation of surface temp for coupled simu 302 303 ! 303 DO jj = 1, jpj 304 DO ji = 1, jpi 305 IF( gphit(ji,jj) > 0._wp ) THEN ; rn_amax_2d(ji,jj) = rn_amax_n ! NH 306 ELSE ; rn_amax_2d(ji,jj) = rn_amax_s ! SH 307 ENDIF 308 END DO 309 END DO 304 WHERE( gphit(:,:) > 0._wp ) ; rn_amax_2d(:,:) = rn_amax_n ! NH 305 ELSEWHERE ; rn_amax_2d(:,:) = rn_amax_s ! SH 306 END WHERE 310 307 ! 311 308 END SUBROUTINE ice_init … … 381 378 !!---------------------------------------------------------------------- 382 379 ! 383 DO jl = 1, jpl 384 DO jj = 2, jpjm1 385 DO ji = 2, jpim1 386 a_i_b (ji,jj,jl) = a_i (ji,jj,jl) ! ice area 387 v_i_b (ji,jj,jl) = v_i (ji,jj,jl) ! ice volume 388 v_s_b (ji,jj,jl) = v_s (ji,jj,jl) ! snow volume 389 smv_i_b(ji,jj,jl) = smv_i(ji,jj,jl) ! salt content 390 oa_i_b (ji,jj,jl) = oa_i (ji,jj,jl) ! areal age content 391 e_s_b (ji,jj,:,jl) = e_s (ji,jj,:,jl) ! snow thermal energy 392 e_i_b (ji,jj,:,jl) = e_i (ji,jj,:,jl) ! ice thermal energy 393 ! ! ice thickness 394 rswitch = MAX( 0._wp , SIGN( 1._wp, a_i_b(ji,jj,jl) - epsi20 ) ) ! 0 if no ice and 1 if yes 395 ht_i_b(ji,jj,jl) = v_i_b (ji,jj,jl) / MAX( a_i_b(ji,jj,jl) , epsi20 ) * rswitch 396 ht_s_b(ji,jj,jl) = v_s_b (ji,jj,jl) / MAX( a_i_b(ji,jj,jl) , epsi20 ) * rswitch 397 END DO 398 END DO 399 END DO 400 CALL lbc_lnk_multi( a_i_b, 'T', 1., v_i_b , 'T', 1., ht_i_b, 'T', 1., smv_i_b, 'T', 1., & 401 & oa_i_b, 'T', 1., v_s_b , 'T', 1., ht_s_b, 'T', 1. ) 402 CALL lbc_lnk( e_i_b, 'T', 1. ) 403 CALL lbc_lnk( e_s_b, 'T', 1. ) 380 a_i_b (:,:,:) = a_i (:,:,:) ! ice area 381 v_i_b (:,:,:) = v_i (:,:,:) ! ice volume 382 v_s_b (:,:,:) = v_s (:,:,:) ! snow volume 383 smv_i_b(:,:,:) = smv_i(:,:,:) ! salt content 384 oa_i_b (:,:,:) = oa_i (:,:,:) ! areal age content 385 e_s_b (:,:,:,:) = e_s (:,:,:,:) ! snow thermal energy 386 e_i_b (:,:,:,:) = e_i (:,:,:,:) ! ice thermal energy 387 WHERE( a_i_b(:,:,:) >= epsi20 ) 388 ht_i_b(:,:,:) = v_i_b (:,:,:) / a_i_b(:,:,:) ! ice thickness 389 ht_s_b(:,:,:) = v_s_b (:,:,:) / a_i_b(:,:,:) ! snw thickness 390 ELSEWHERE 391 ht_i_b(:,:,:) = 0._wp 392 ht_s_b(:,:,:) = 0._wp 393 END WHERE 404 394 405 !!gm Question: here , a_i_b, u_ice and v_ice are defined over the whole domain,406 !!gm so why not just a copy over the whole domain and no lbc_lnk ?407 !!gm that is some thing like:408 ! at_i_b(:,:) = SUM( a_i_b(:,:,:), dim=3 )409 ! u_ice_b(:,:) = u_ice(:,:)410 ! v_ice_b(:,:) = v_ice(:,:)411 ! idem for the loop above412 !!gm413 395 ! ice velocities & total concentration 414 DO jj = 2, jpjm1 415 DO ji = 2, jpim1 416 at_i_b(ji,jj) = SUM( a_i_b(ji,jj,:) ) 417 u_ice_b(ji,jj) = u_ice(ji,jj) 418 v_ice_b(ji,jj) = v_ice(ji,jj) 419 END DO 420 END DO 421 CALL lbc_lnk_multi( at_i_b, 'T', 1., u_ice_b , 'U', -1., v_ice_b , 'V', -1. ) 396 at_i_b(:,:) = SUM( a_i_b(:,:,:), dim=3 ) 397 u_ice_b(:,:) = u_ice(:,:) 398 v_ice_b(:,:) = v_ice(:,:) 422 399 ! 423 400 END SUBROUTINE store_fields … … 433 410 INTEGER :: ji, jj ! dummy loop index 434 411 !!---------------------------------------------------------------------- 435 DO jj = 1, jpj 436 DO ji = 1, jpi 437 sfx (ji,jj) = 0._wp ; 438 sfx_bri(ji,jj) = 0._wp ; sfx_lam(ji,jj) = 0._wp 439 sfx_sni(ji,jj) = 0._wp ; sfx_opw(ji,jj) = 0._wp 440 sfx_bog(ji,jj) = 0._wp ; sfx_dyn(ji,jj) = 0._wp 441 sfx_bom(ji,jj) = 0._wp ; sfx_sum(ji,jj) = 0._wp 442 sfx_res(ji,jj) = 0._wp ; sfx_sub(ji,jj) = 0._wp 443 ! 444 wfx_snw(ji,jj) = 0._wp ; wfx_ice(ji,jj) = 0._wp 445 wfx_sni(ji,jj) = 0._wp ; wfx_opw(ji,jj) = 0._wp 446 wfx_bog(ji,jj) = 0._wp ; wfx_dyn(ji,jj) = 0._wp 447 wfx_bom(ji,jj) = 0._wp ; wfx_sum(ji,jj) = 0._wp 448 wfx_res(ji,jj) = 0._wp ; wfx_sub(ji,jj) = 0._wp 449 wfx_spr(ji,jj) = 0._wp ; wfx_lam(ji,jj) = 0._wp 450 wfx_snw_dyn(ji,jj) = 0._wp ; wfx_snw_sum(ji,jj) = 0._wp 451 wfx_snw_sub(ji,jj) = 0._wp ; wfx_ice_sub(ji,jj) = 0._wp 452 wfx_snw_sni(ji,jj) = 0._wp 453 ! MV MP 2016 454 wfx_pnd(ji,jj) = 0._wp 455 ! END MV MP 2016 456 457 hfx_thd(ji,jj) = 0._wp ; 458 hfx_snw(ji,jj) = 0._wp ; hfx_opw(ji,jj) = 0._wp 459 hfx_bog(ji,jj) = 0._wp ; hfx_dyn(ji,jj) = 0._wp 460 hfx_bom(ji,jj) = 0._wp ; hfx_sum(ji,jj) = 0._wp 461 hfx_res(ji,jj) = 0._wp ; hfx_sub(ji,jj) = 0._wp 462 hfx_spr(ji,jj) = 0._wp ; hfx_dif(ji,jj) = 0._wp 463 hfx_err(ji,jj) = 0._wp ; hfx_err_rem(ji,jj) = 0._wp 464 hfx_err_dif(ji,jj) = 0._wp 465 wfx_err_sub(ji,jj) = 0._wp 466 ! 467 afx_tot(ji,jj) = 0._wp ; 468 ! 469 diag_heat(ji,jj) = 0._wp ; diag_smvi(ji,jj) = 0._wp 470 diag_vice(ji,jj) = 0._wp ; diag_vsnw(ji,jj) = 0._wp 471 472 ! SIMIP diagnostics 473 diag_fc_bo(ji,jj) = 0._wp ; diag_fc_su(ji,jj) = 0._wp 474 475 tau_icebfr(ji,jj) = 0._wp; ! landfast ice param only (clem: important to keep the init here) 476 END DO 477 END DO 412 sfx (:,:) = 0._wp ; 413 sfx_bri(:,:) = 0._wp ; sfx_lam(:,:) = 0._wp 414 sfx_sni(:,:) = 0._wp ; sfx_opw(:,:) = 0._wp 415 sfx_bog(:,:) = 0._wp ; sfx_dyn(:,:) = 0._wp 416 sfx_bom(:,:) = 0._wp ; sfx_sum(:,:) = 0._wp 417 sfx_res(:,:) = 0._wp ; sfx_sub(:,:) = 0._wp 418 ! 419 wfx_snw(:,:) = 0._wp ; wfx_ice(:,:) = 0._wp 420 wfx_sni(:,:) = 0._wp ; wfx_opw(:,:) = 0._wp 421 wfx_bog(:,:) = 0._wp ; wfx_dyn(:,:) = 0._wp 422 wfx_bom(:,:) = 0._wp ; wfx_sum(:,:) = 0._wp 423 wfx_res(:,:) = 0._wp ; wfx_sub(:,:) = 0._wp 424 wfx_spr(:,:) = 0._wp ; wfx_lam(:,:) = 0._wp 425 wfx_snw_dyn(:,:) = 0._wp ; wfx_snw_sum(:,:) = 0._wp 426 wfx_snw_sub(:,:) = 0._wp ; wfx_ice_sub(:,:) = 0._wp 427 wfx_snw_sni(:,:) = 0._wp 428 ! MV MP 2016 429 wfx_pnd(:,:) = 0._wp 430 ! END MV MP 2016 431 432 hfx_thd(:,:) = 0._wp ; 433 hfx_snw(:,:) = 0._wp ; hfx_opw(:,:) = 0._wp 434 hfx_bog(:,:) = 0._wp ; hfx_dyn(:,:) = 0._wp 435 hfx_bom(:,:) = 0._wp ; hfx_sum(:,:) = 0._wp 436 hfx_res(:,:) = 0._wp ; hfx_sub(:,:) = 0._wp 437 hfx_spr(:,:) = 0._wp ; hfx_dif(:,:) = 0._wp 438 hfx_err(:,:) = 0._wp ; hfx_err_rem(:,:) = 0._wp 439 hfx_err_dif(:,:) = 0._wp 440 wfx_err_sub(:,:) = 0._wp 441 ! 442 afx_tot(:,:) = 0._wp ; 443 ! 444 diag_heat(:,:) = 0._wp ; diag_smvi(:,:) = 0._wp 445 diag_vice(:,:) = 0._wp ; diag_vsnw(:,:) = 0._wp 446 447 ! SIMIP diagnostics 448 diag_fc_bo(:,:) = 0._wp ; diag_fc_su(:,:) = 0._wp 449 450 tau_icebfr(:,:) = 0._wp; ! landfast ice param only (clem: important to keep the init here) 478 451 479 452 END SUBROUTINE ice_diag0 -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd.F90
r8514 r8515 539 539 INTEGER :: ios ! Local integer output status for namelist read 540 540 !! 541 NAMELIST/namice_thd/ ln_icethd, rn_kappa_i, ln_cndi_U64, ln_cndi_P07, ln_dqns_i, rn_cnd_s, 542 & ln_icedH, rn_blow_s,&543 & ln_icedA, rn_beta, rn_dmin,&544 & ln_icedO, rn_hinew, ln_frazil, rn_maxfraz, rn_vfraz, rn_Cfraz, rn_himin,&545 & nn_iceflx541 NAMELIST/namice_thd/ ln_icethd, rn_kappa_i, ln_cndi_U64, ln_cndi_P07, ln_dqns_i, rn_cnd_s, & 542 & ln_icedH, rn_blow_s, & 543 & ln_icedA, rn_beta, rn_dmin, & 544 & ln_icedO, rn_hinew, ln_frazil, rn_maxfraz, rn_vfraz, rn_Cfraz, & 545 & nn_iceflx 546 546 !!------------------------------------------------------------------- 547 547 ! … … 581 581 WRITE(numout,*) ' Threshold relative drift speed for collection of frazil rn_vfraz = ', rn_vfraz 582 582 WRITE(numout,*) ' Squeezing coefficient for collection of frazil rn_Cfraz = ', rn_Cfraz 583 WRITE(numout,*) ' -- iceitd --'584 WRITE(numout,*) ' minimum ice thickness rn_himin = ', rn_himin585 583 WRITE(numout,*) ' -- icestp --' 586 584 WRITE(numout,*) ' Multicategory heat flux formulation nn_iceflx = ', nn_iceflx
Note: See TracChangeset
for help on using the changeset viewer.