Changeset 2715 for trunk/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90
- Timestamp:
- 2011-03-30T17:58:35+02:00 (13 years ago)
- File:
-
- 1 edited
-
trunk/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90 (modified) (64 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90
r2528 r2715 2 2 !!====================================================================== 3 3 !! *** MODULE limitd_me *** 4 !! Mechanical impact on ice thickness distribution 5 !! computation of changes in g(h) 4 !! LIM-3 : Mechanical impact on ice thickness distribution 6 5 !!====================================================================== 7 6 !! History : LIM ! 2006-02 (M. Vancoppenolle) Original code 8 7 !! 3.2 ! 2009-07 (M. Vancoppenolle, Y. Aksenov, G. Madec) bug correction in smsw & fsalt_rpo 8 !! 4.0 ! 2011-02 (G. Madec) dynamical allocation 9 9 !!---------------------------------------------------------------------- 10 10 #if defined key_lim3 … … 12 12 !! 'key_lim3' : LIM3 sea-ice model 13 13 !!---------------------------------------------------------------------- 14 USE dom_ice15 14 USE par_oce ! ocean parameters 16 USE dom_oce 17 USE lbclnk 15 USE dom_oce ! ocean domain 18 16 USE phycst ! physical constants (ocean directory) 19 USE sbc_oce ! Surface boundary condition: ocean fields 20 USE thd_ice 21 USE in_out_manager 22 USE ice 23 USE par_ice 24 USE limthd_lac 25 USE limvar 26 USE limcons 17 USE sbc_oce ! surface boundary condition: ocean fields 18 USE thd_ice ! LIM thermodynamics 19 USE ice ! LIM variables 20 USE par_ice ! LIM parameters 21 USE dom_ice ! LIM domain 22 USE limthd_lac ! LIM 23 USE limvar ! LIM 24 USE limcons ! LIM 25 USE in_out_manager ! I/O manager 26 USE lbclnk ! lateral boundary condition - MPP exchanges 27 USE lib_mpp ! MPP library 27 28 USE prtctl ! Print control 28 USE lib_mpp29 USE wrk_nemo ! workspace manager 29 30 30 31 IMPLICIT NONE 31 32 PRIVATE 32 33 33 !! * Routine accessibility 34 PUBLIC lim_itd_me ! called by ice_stp 35 PUBLIC lim_itd_me_icestrength 36 PUBLIC lim_itd_me_ridgeprep 37 PUBLIC lim_itd_me_ridgeshift 38 PUBLIC lim_itd_me_asumr 39 PUBLIC lim_itd_me_init 40 PUBLIC lim_itd_me_zapsmall 41 42 !! * Module variables 43 REAL(wp) :: & ! constant values 44 epsi20 = 1e-20 , & 45 epsi13 = 1e-13 , & 46 epsi11 = 1e-11 , & 47 zzero = 0.e0 , & 48 zone = 1.e0 34 PUBLIC lim_itd_me ! called by ice_stp 35 PUBLIC lim_itd_me_icestrength 36 PUBLIC lim_itd_me_init 37 PUBLIC lim_itd_me_zapsmall 38 PUBLIC lim_itd_me_alloc ! called by nemogcm.F90 39 40 REAL(wp) :: epsi11 = 1.e-11_wp ! constant values 41 REAL(wp) :: epsi10 = 1.e-10_wp ! constant values 42 REAL(wp) :: epsi06 = 1.e-06_wp ! constant values 49 43 50 44 !----------------------------------------------------------------------- 51 45 ! Variables shared among ridging subroutines 52 46 !----------------------------------------------------------------------- 53 REAL(wp), DIMENSION (jpi,jpj) :: & 54 asum , & ! sum of total ice and open water area 55 aksum ! ratio of area removed to area ridged 56 57 REAL(wp), DIMENSION(jpi,jpj,0:jpl) :: & 58 athorn ! participation function; fraction of ridging/ 59 ! closing associated w/ category n 60 61 REAL(wp), DIMENSION(jpi,jpj,jpl) :: & 62 hrmin , & ! minimum ridge thickness 63 hrmax , & ! maximum ridge thickness 64 hraft , & ! thickness of rafted ice 65 krdg , & ! mean ridge thickness/thickness of ridging ice 66 aridge , & ! participating ice ridging 67 araft ! participating ice rafting 68 69 REAL(wp), PARAMETER :: & 70 krdgmin = 1.1, & ! min ridge thickness multiplier 71 kraft = 2.0 ! rafting multipliyer 72 73 REAL(wp) :: & 74 Cp 47 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: asum ! sum of total ice and open water area 48 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: aksum ! ratio of area removed to area ridged 49 50 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: athorn ! participation function; fraction of ridging/ 51 ! ! closing associated w/ category n 52 53 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hrmin ! minimum ridge thickness 54 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hrmax ! maximum ridge thickness 55 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hraft ! thickness of rafted ice 56 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: krdg ! mean ridge thickness/thickness of ridging ice 57 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: aridge ! participating ice ridging 58 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: araft ! participating ice rafting 59 60 REAL(wp), PARAMETER :: krdgmin = 1.1_wp ! min ridge thickness multiplier 61 REAL(wp), PARAMETER :: kraft = 2.0_wp ! rafting multipliyer 62 63 REAL(wp) :: Cp ! 75 64 ! 76 65 !----------------------------------------------------------------------- 77 66 ! Ridging diagnostic arrays for history files 78 67 !----------------------------------------------------------------------- 79 ! 80 REAL (wp), DIMENSION(jpi,jpj) :: & 81 dardg1dt , & ! rate of fractional area loss by ridging ice (1/s) 82 dardg2dt , & ! rate of fractional area gain by new ridges (1/s) 83 dvirdgdt , & ! rate of ice volume ridged (m/s) 84 opening ! rate of opening due to divergence/shear (1/s) 85 68 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: dardg1dt ! rate of fractional area loss by ridging ice (1/s) 69 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: dardg2dt ! rate of fractional area gain by new ridges (1/s) 70 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: dvirdgdt ! rate of ice volume ridged (m/s) 71 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: opening ! rate of opening due to divergence/shear (1/s) 86 72 87 73 !!---------------------------------------------------------------------- 88 74 !! NEMO/LIM3 3.3 , UCL - NEMO Consortium (2010) 89 75 !! $Id$ 90 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)76 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 91 77 !!---------------------------------------------------------------------- 92 93 94 78 CONTAINS 95 79 96 !!-----------------------------------------------------------------------------! 97 !!-----------------------------------------------------------------------------! 98 99 SUBROUTINE lim_itd_me ! (subroutine 1/6) 80 INTEGER FUNCTION lim_itd_me_alloc() 81 !!---------------------------------------------------------------------! 82 !! *** ROUTINE lim_itd_me_alloc *** 83 !!---------------------------------------------------------------------! 84 ALLOCATE( & 85 !* Variables shared among ridging subroutines 86 & asum (jpi,jpj) , athorn(jpi,jpj,0:jpl) , & 87 & aksum(jpi,jpj) , & 88 ! 89 & hrmin(jpi,jpj,jpl) , hraft(jpi,jpj,jpl) , aridge(jpi,jpj,jpl) , & 90 & hrmax(jpi,jpj,jpl) , krdg (jpi,jpj,jpl) , araft (jpi,jpj,jpl) , & 91 ! 92 !* Ridging diagnostic arrays for history files 93 & dardg1dt(jpi,jpj) , dardg2dt(jpi,jpj) , & 94 & dvirdgdt(jpi,jpj) , opening(jpi,jpj) , STAT=lim_itd_me_alloc ) 95 ! 96 IF( lim_itd_me_alloc /= 0 ) CALL ctl_warn( 'lim_itd_me_alloc: failed to allocate arrays' ) 97 ! 98 END FUNCTION lim_itd_me_alloc 99 100 101 SUBROUTINE lim_itd_me 100 102 !!---------------------------------------------------------------------! 101 103 !! *** ROUTINE lim_itd_me *** 102 !! ** Purpose : 103 !! This routine computes the mechanical redistribution 104 !! of ice thickness 105 !! 106 !! ** Method : a very simple method :-) 107 !! 108 !! ** Arguments : 109 !! kideb , kiut : Starting and ending points on which the 110 !! the computation is applied 111 !! 112 !! ** Inputs / Ouputs : (global commons) 113 !! 114 !! ** External : 115 !! 116 !! ** Steps : 117 !! 1) Thickness categories boundaries, ice / o.w. concentrations 118 !! Ridge preparation 119 !! 2) Dynamical inputs (closing rate, divu_adv, opning) 120 !! 3) Ridging iteration 121 !! 4) Ridging diagnostics 122 !! 5) Heat, salt and freshwater fluxes 123 !! 6) Compute increments of tate variables and come back to old values 124 !! 125 !! ** References : There are a lot of references and can be difficult / 126 !! boring to read 127 !! 128 !! Flato, G. M., and W. D. Hibler III, 1995: Ridging and strength 129 !! in modeling the thickness distribution of Arctic sea ice, 130 !! J. Geophys. Res., 100, 18,611-18,626. 131 !! 132 !! Hibler, W. D. III, 1980: Modeling a variable thickness sea ice 133 !! cover, Mon. Wea. Rev., 108, 1943-1973, 1980. 134 !! 135 !! Rothrock, D. A., 1975: The energetics of the plastic deformation of 136 !! pack ice by ridging, J. Geophys. Res., 80, 4514-4519. 137 !! 138 !! Thorndike, A. S., D. A. Rothrock, G. A. Maykut, and R. Colony, 139 !! 1975: The thickness distribution of sea ice, J. Geophys. Res., 140 !! 80, 4501-4513. 141 !! 142 !! Bitz et al., JGR 2001 143 !! 144 !! Amundrud and Melling, JGR 2005 145 !! 146 !! Babko et al., JGR 2002 104 !! 105 !! ** Purpose : computes the mechanical redistribution of ice thickness 106 !! 107 !! ** Method : Steps : 108 !! 1) Thickness categories boundaries, ice / o.w. concentrations 109 !! Ridge preparation 110 !! 2) Dynamical inputs (closing rate, divu_adv, opning) 111 !! 3) Ridging iteration 112 !! 4) Ridging diagnostics 113 !! 5) Heat, salt and freshwater fluxes 114 !! 6) Compute increments of tate variables and come back to old values 115 !! 116 !! References : Flato, G. M., and W. D. Hibler III, 1995, JGR, 100, 18,611-18,626. 117 !! Hibler, W. D. III, 1980, MWR, 108, 1943-1973, 1980. 118 !! Rothrock, D. A., 1975: JGR, 80, 4514-4519. 119 !! Thorndike et al., 1975, JGR, 80, 4501-4513. 120 !! Bitz et al., JGR, 2001 121 !! Amundrud and Melling, JGR 2005 122 !! Babko et al., JGR 2002 147 123 !! 148 124 !! This routine is based on CICE code and authors William H. Lipscomb, 149 125 !! and Elizabeth C. Hunke, LANL are gratefully acknowledged 150 126 !!--------------------------------------------------------------------! 151 !! * Arguments 152 153 !! * Local variables 154 INTEGER :: ji, & ! spatial dummy loop index 155 jj, & ! spatial dummy loop index 156 jk, & ! vertical layering dummy loop index 157 jl, & ! ice category dummy loop index 158 niter, & ! iteration counter 159 nitermax = 20 ! max number of ridging iterations 160 161 REAL(wp) :: & ! constant values 162 zeps = 1.0e-10, & 163 epsi10 = 1.0e-10, & 164 epsi06 = 1.0e-6 165 166 REAL(wp), DIMENSION(jpi,jpj) :: & 167 closing_net, & ! net rate at which area is removed (1/s) 168 ! (ridging ice area - area of new ridges) / dt 169 divu_adv , & ! divu as implied by transport scheme (1/s) 170 opning , & ! rate of opening due to divergence/shear 171 closing_gross, & ! rate at which area removed, not counting 172 ! area of new ridges 173 msnow_mlt , & ! mass of snow added to ocean (kg m-2) 174 esnow_mlt ! energy needed to melt snow in ocean (J m-2) 175 176 REAL(wp) :: & 177 w1, & ! temporary variable 178 tmpfac, & ! factor by which opening/closing rates are cut 179 dti ! 1 / dt 180 181 LOGICAL :: & 182 asum_error ! flag for asum .ne. 1 183 184 INTEGER :: iterate_ridging ! if true, repeat the ridging 185 186 REAL(wp) :: & 187 big = 1.0e8 188 189 REAL (wp), DIMENSION(jpi,jpj) :: & ! 190 vt_i_init, vt_i_final ! ice volume summed over categories 191 192 CHARACTER (len = 15) :: fieldid 193 194 !!-- End of declarations 195 !-----------------------------------------------------------------------------! 196 197 IF( numit == nstart ) CALL lim_itd_me_init ! Initialization (first time-step only) 127 USE wrk_nemo, ONLY: closing_net => wrk_2d_1 ! net rate at which area is removed (1/s) 128 ! ! (ridging ice area - area of new ridges) / dt 129 USE wrk_nemo, ONLY: divu_adv => wrk_2d_2 ! divu as implied by transport scheme (1/s) 130 USE wrk_nemo, ONLY: opning => wrk_2d_3 ! rate of opening due to divergence/shear 131 USE wrk_nemo, ONLY: closing_gross => wrk_2d_4 ! rate at which area removed, not counting area of new ridges 132 USE wrk_nemo, ONLY: msnow_mlt => wrk_2d_5 ! mass of snow added to ocean (kg m-2) 133 USE wrk_nemo, ONLY: esnow_mlt => wrk_2d_6 ! energy needed to melt snow in ocean (J m-2) 134 USE wrk_nemo, ONLY: vt_i_init => wrk_2d_7 ! ice volume summed over 135 USE wrk_nemo, ONLY: vt_i_final => wrk_2d_8 ! categories 136 ! 137 INTEGER :: ji, jj, jk, jl ! dummy loop index 138 INTEGER :: niter, nitermax = 20 ! local integer 139 LOGICAL :: asum_error ! flag for asum .ne. 1 140 INTEGER :: iterate_ridging ! if true, repeat the ridging 141 REAL(wp) :: w1, tmpfac, dti ! local scalar 142 CHARACTER (len = 15) :: fieldid 143 !!----------------------------------------------------------------------------- 144 145 IF( wrk_in_use(2, 1,2,3,4,5,6,7,8) ) THEN 146 CALL ctl_stop('lim_itd_me: requested workspace arrays unavailable') ; RETURN 147 ENDIF 148 149 IF( numit == nstart ) CALL lim_itd_me_init ! Initialization (first time-step only) 198 150 199 151 IF(ln_ctl) THEN … … 210 162 hi_max(jpl) = 999.99 211 163 212 Cp = 0.5 * grav * (rau0-rhoic)*rhoic/rau0! proport const for PE164 Cp = 0.5 * grav * (rau0-rhoic) * rhoic / rau0 ! proport const for PE 213 165 CALL lim_itd_me_ridgeprep ! prepare ridging 214 166 215 ! conservation check 216 IF ( con_i) CALL lim_column_sum (jpl, v_i, vt_i_init) 217 218 ! Initialize arrays. 219 DO jj = 1, jpj 167 IF( con_i) CALL lim_column_sum( jpl, v_i, vt_i_init ) ! conservation check 168 169 DO jj = 1, jpj ! Initialize arrays. 220 170 DO ji = 1, jpi 221 222 msnow_mlt(ji,jj) = 0.0 223 esnow_mlt(ji,jj) = 0.0 224 dardg1dt(ji,jj) = 0.0 225 dardg2dt(ji,jj) = 0.0 226 dvirdgdt(ji,jj) = 0.0 227 opening (ji,jj) = 0.0 171 msnow_mlt(ji,jj) = 0._wp 172 esnow_mlt(ji,jj) = 0._wp 173 dardg1dt (ji,jj) = 0._wp 174 dardg2dt (ji,jj) = 0._wp 175 dvirdgdt (ji,jj) = 0._wp 176 opening (ji,jj) = 0._wp 228 177 229 178 !-----------------------------------------------------------------------------! … … 246 195 ! (thick, newly ridged ice). 247 196 248 closing_net(ji,jj) = & 249 Cs*0.5*(Delta_i(ji,jj)-ABS(divu_i(ji,jj))) - MIN(divu_i(ji,jj),0.0) 197 closing_net(ji,jj) = Cs * 0.5 * ( Delta_i(ji,jj) - ABS( divu_i(ji,jj) ) ) - MIN( divu_i(ji,jj), 0._wp ) 250 198 251 199 ! 2.2 divu_adv … … 258 206 ! to give asum = 1.0 after ridging. 259 207 260 divu_adv(ji,jj) = (1.0-asum(ji,jj)) / rdt_ice ! asum found in ridgeprep 261 262 IF (divu_adv(ji,jj) .LT. 0.0) & 263 closing_net(ji,jj) = max(closing_net(ji,jj), -divu_adv(ji,jj)) 208 divu_adv(ji,jj) = ( 1._wp - asum(ji,jj) ) / rdt_ice ! asum found in ridgeprep 209 210 IF( divu_adv(ji,jj) < 0._wp ) closing_net(ji,jj) = MAX( closing_net(ji,jj), -divu_adv(ji,jj) ) 264 211 265 212 ! 2.3 opning … … 268 215 ! asum = 1.0 after ridging. 269 216 opning(ji,jj) = closing_net(ji,jj) + divu_adv(ji,jj) 270 271 217 END DO 272 218 END DO … … 275 221 ! 3) Ridging iteration 276 222 !-----------------------------------------------------------------------------! 277 niter = 1 ! iteration counter223 niter = 1 ! iteration counter 278 224 iterate_ridging = 1 279 280 225 281 226 DO WHILE ( iterate_ridging > 0 .AND. niter < nitermax ) … … 315 260 DO jj = 1, jpj 316 261 DO ji = 1, jpi 317 IF ( a_i(ji,jj,jl) .GT. epsi11 .AND. athorn(ji,jj,jl) .GT. 0.0 )THEN262 IF ( a_i(ji,jj,jl) > epsi11 .AND. athorn(ji,jj,jl) > 0._wp )THEN 318 263 w1 = athorn(ji,jj,jl) * closing_gross(ji,jj) * rdt_ice 319 IF ( w1 .GT.a_i(ji,jj,jl) ) THEN264 IF ( w1 > a_i(ji,jj,jl) ) THEN 320 265 tmpfac = a_i(ji,jj,jl) / w1 321 266 closing_gross(ji,jj) = closing_gross(ji,jj) * tmpfac 322 opning (ji,jj) = opning(ji,jj) * tmpfac267 opning (ji,jj) = opning (ji,jj) * tmpfac 323 268 ENDIF 324 269 ENDIF … … 330 275 !-----------------------------------------------------------------------------! 331 276 332 CALL lim_itd_me_ridgeshift (opning, closing_gross, & 333 msnow_mlt, esnow_mlt) 277 CALL lim_itd_me_ridgeshift( opning, closing_gross, msnow_mlt, esnow_mlt ) 334 278 335 279 ! 3.4 Compute total area of ice plus open water after ridging. … … 348 292 DO ji = 1, jpi 349 293 IF (ABS(asum(ji,jj) - 1.0) .LT. epsi11) THEN 350 closing_net(ji,jj) = 0. 0351 opning (ji,jj) = 0.0294 closing_net(ji,jj) = 0._wp 295 opning (ji,jj) = 0._wp 352 296 ELSE 353 297 iterate_ridging = 1 354 divu_adv (ji,jj) = (1.0- asum(ji,jj)) / rdt_ice355 closing_net(ji,jj) = MAX( 0.0, -divu_adv(ji,jj))356 opning (ji,jj) = MAX(0.0, divu_adv(ji,jj))298 divu_adv (ji,jj) = (1._wp - asum(ji,jj)) / rdt_ice 299 closing_net(ji,jj) = MAX( 0._wp, -divu_adv(ji,jj) ) 300 opning (ji,jj) = MAX( 0._wp, divu_adv(ji,jj) ) 357 301 ENDIF 358 302 END DO 359 303 END DO 360 304 361 IF( lk_mpp ) CALL mpp_max(iterate_ridging)305 IF( lk_mpp ) CALL mpp_max( iterate_ridging ) 362 306 363 307 ! Repeat if necessary. … … 368 312 niter = niter + 1 369 313 370 IF (iterate_ridging == 1) THEN371 IF (niter .GT. nitermax) THEN314 IF( iterate_ridging == 1 ) THEN 315 IF( niter .GT. nitermax ) THEN 372 316 WRITE(numout,*) ' ALERTE : non-converging ridging scheme ' 373 317 WRITE(numout,*) ' niter, iterate_ridging ', niter, iterate_ridging … … 384 328 ! Update fresh water and heat fluxes due to snow melt. 385 329 386 dti = 1. 0/rdt_ice330 dti = 1._wp / rdt_ice 387 331 388 332 asum_error = .false. … … 401 345 ! 5) Heat, salt and freshwater fluxes 402 346 !-----------------------------------------------------------------------------! 403 ! fresh water source for ocean 404 fmmec(ji,jj) = fmmec(ji,jj) + msnow_mlt(ji,jj)*dti 405 406 ! heat sink for ocean 407 fhmec(ji,jj) = fhmec(ji,jj) + esnow_mlt(ji,jj)*dti 347 fmmec(ji,jj) = fmmec(ji,jj) + msnow_mlt(ji,jj) * dti ! fresh water source for ocean 348 fhmec(ji,jj) = fhmec(ji,jj) + esnow_mlt(ji,jj) * dti ! heat sink for ocean 408 349 409 350 END DO … … 446 387 !----------------- 447 388 448 d_u_ice_dyn(:,:) = u_ice(:,:) - old_u_ice(:,:) 449 d_v_ice_dyn(:,:) = v_ice(:,:) - old_v_ice(:,:) 450 d_a_i_trp(:,:,:) = a_i(:,:,:) - old_a_i(:,:,:) 451 d_v_s_trp(:,:,:) = v_s(:,:,:) - old_v_s(:,:,:) 452 d_v_i_trp(:,:,:) = v_i(:,:,:) - old_v_i(:,:,:) 453 d_e_s_trp(:,:,:,:) = e_s(:,:,:,:) - old_e_s(:,:,:,:) 454 d_e_i_trp(:,:,:,:) = e_i(:,:,:,:) - old_e_i(:,:,:,:) 455 d_oa_i_trp(:,:,:) = oa_i(:,:,:) - old_oa_i(:,:,:) 456 d_smv_i_trp(:,:,:) = 0.0 457 IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) & 458 d_smv_i_trp(:,:,:) = smv_i(:,:,:) - old_smv_i(:,:,:) 389 d_u_ice_dyn(:,:) = u_ice(:,:) - old_u_ice(:,:) 390 d_v_ice_dyn(:,:) = v_ice(:,:) - old_v_ice(:,:) 391 d_a_i_trp (:,:,:) = a_i (:,:,:) - old_a_i (:,:,:) 392 d_v_s_trp (:,:,:) = v_s (:,:,:) - old_v_s (:,:,:) 393 d_v_i_trp (:,:,:) = v_i (:,:,:) - old_v_i (:,:,:) 394 d_e_s_trp (:,:,:,:) = e_s (:,:,:,:) - old_e_s (:,:,:,:) 395 d_e_i_trp (:,:,:,:) = e_i (:,:,:,:) - old_e_i (:,:,:,:) 396 d_oa_i_trp (:,:,:) = oa_i (:,:,:) - old_oa_i (:,:,:) 397 d_smv_i_trp(:,:,:) = 0._wp 398 IF( num_sal == 2 .OR. num_sal == 4 ) d_smv_i_trp(:,:,:) = smv_i(:,:,:) - old_smv_i(:,:,:) 459 399 460 400 IF(ln_ctl) THEN ! Control print … … 503 443 e_i(:,:,:,:) = old_e_i(:,:,:,:) 504 444 oa_i(:,:,:) = old_oa_i(:,:,:) 505 IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) & 506 smv_i(:,:,:) = old_smv_i(:,:,:) 445 IF( num_sal == 2 .OR. num_sal == 4 ) smv_i(:,:,:) = old_smv_i(:,:,:) 507 446 508 447 !----------------------------------------------------! … … 518 457 DO jj = 1, jpj 519 458 DO ji = 1, jpi 520 IF ( ( old_v_i(ji,jj,jl) .LT.epsi06 ) .AND. &521 ( d_v_i_trp(ji,jj,jl) .GT.epsi06 ) ) THEN459 IF ( ( old_v_i(ji,jj,jl) < epsi06 ) .AND. & 460 ( d_v_i_trp(ji,jj,jl) > epsi06 ) ) THEN 522 461 old_e_i(ji,jj,jk,jl) = d_e_i_trp(ji,jj,jk,jl) 523 d_e_i_trp(ji,jj,jk,jl) = 0. 0462 d_e_i_trp(ji,jj,jk,jl) = 0._wp 524 463 ENDIF 525 464 END DO … … 531 470 DO jj = 1, jpj 532 471 DO ji = 1, jpi 533 IF ( ( old_v_i(ji,jj,jl) .LT.epsi06 ) .AND. &534 ( d_v_i_trp(ji,jj,jl) .GT.epsi06 ) ) THEN472 IF ( ( old_v_i(ji,jj,jl) < epsi06 ) .AND. & 473 ( d_v_i_trp(ji,jj,jl) > epsi06 ) ) THEN 535 474 old_v_i(ji,jj,jl) = d_v_i_trp(ji,jj,jl) 536 d_v_i_trp(ji,jj,jl) = 0. 0475 d_v_i_trp(ji,jj,jl) = 0._wp 537 476 old_a_i(ji,jj,jl) = d_a_i_trp(ji,jj,jl) 538 d_a_i_trp(ji,jj,jl) = 0. 0477 d_a_i_trp(ji,jj,jl) = 0._wp 539 478 old_v_s(ji,jj,jl) = d_v_s_trp(ji,jj,jl) 540 d_v_s_trp(ji,jj,jl) = 0. 0479 d_v_s_trp(ji,jj,jl) = 0._wp 541 480 old_e_s(ji,jj,1,jl) = d_e_s_trp(ji,jj,1,jl) 542 d_e_s_trp(ji,jj,1,jl) = 0. 0481 d_e_s_trp(ji,jj,1,jl) = 0._wp 543 482 old_oa_i(ji,jj,jl) = d_oa_i_trp(ji,jj,jl) 544 d_oa_i_trp(ji,jj,jl) = 0.0 545 IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) & 546 old_smv_i(ji,jj,jl) = d_smv_i_trp(ji,jj,jl) 547 d_smv_i_trp(ji,jj,jl) = 0.0 483 d_oa_i_trp(ji,jj,jl) = 0._wp 484 IF( num_sal == 2 .OR. num_sal == 4 ) old_smv_i(ji,jj,jl) = d_smv_i_trp(ji,jj,jl) 485 d_smv_i_trp(ji,jj,jl) = 0._wp 548 486 ENDIF 549 487 END DO … … 551 489 END DO 552 490 491 IF( wrk_not_released(2, 1,2,3,4,5,6,7,8) ) CALL ctl_stop('lim_itd_me: failed to release workspace arrays') 492 ! 553 493 END SUBROUTINE lim_itd_me 554 494 555 !=============================================================================== 556 557 SUBROUTINE lim_itd_me_icestrength (kstrngth) ! (subroutine 2/6) 558 495 496 SUBROUTINE lim_itd_me_icestrength( kstrngth ) 559 497 !!---------------------------------------------------------------------- 560 498 !! *** ROUTINE lim_itd_me_icestrength *** 561 !! ** Purpose : 562 !! This routine computes ice strength used in dynamics routines 563 !! of ice thickness 564 !! 565 !! ** Method : 566 !! Compute the strength of the ice pack, defined as the energy (J m-2) 567 !! dissipated per unit area removed from the ice pack under compression, 568 !! and assumed proportional to the change in potential energy caused 569 !! by ridging. Note that only Hibler's formulation is stable and that 570 !! ice strength has to be smoothed 499 !! 500 !! ** Purpose : computes ice strength used in dynamics routines of ice thickness 501 !! 502 !! ** Method : Compute the strength of the ice pack, defined as the energy (J m-2) 503 !! dissipated per unit area removed from the ice pack under compression, 504 !! and assumed proportional to the change in potential energy caused 505 !! by ridging. Note that only Hibler's formulation is stable and that 506 !! ice strength has to be smoothed 571 507 !! 572 508 !! ** Inputs / Ouputs : kstrngth (what kind of ice strength we are using) 573 !!574 !! ** External :575 !!576 !! ** References :577 !!578 509 !!---------------------------------------------------------------------- 579 !! * Arguments 580 581 INTEGER, INTENT(in) :: & 582 kstrngth ! = 1 for Rothrock formulation, 0 for Hibler (1979) 583 584 INTEGER :: & 585 ji,jj, & !: horizontal indices 586 jl, & !: thickness category index 587 ksmooth, & !: smoothing the resistance to deformation 588 numts_rm !: number of time steps for the P smoothing 589 590 REAL(wp) :: & 591 hi, & !: ice thickness (m) 592 zw1, & !: temporary variable 593 zp, & !: temporary ice strength 594 zdummy 595 596 REAL(wp), DIMENSION(jpi,jpj) :: & 597 zworka !: temporary array used here 510 USE wrk_nemo, ONLY: zworka => wrk_2d_1 ! 2D workspace 511 ! 512 INTEGER, INTENT(in) :: kstrngth ! = 1 for Rothrock formulation, 0 for Hibler (1979) 513 514 INTEGER :: ji,jj, jl ! dummy loop indices 515 INTEGER :: ksmooth ! smoothing the resistance to deformation 516 INTEGER :: numts_rm ! number of time steps for the P smoothing 517 518 REAL(wp) :: hi, zw1, zp, zdummy, zzc, z1_3 ! local scalars 519 !!---------------------------------------------------------------------- 520 521 IF( wrk_in_use(2, 1) ) THEN 522 CALL ctl_stop('lim_itd_me_icestrength : requested workspace array unavailable') ; RETURN 523 ENDIF 598 524 599 525 !------------------------------------------------------------------------------! 600 526 ! 1) Initialize 601 527 !------------------------------------------------------------------------------! 602 strength(:,:) = 0. 0528 strength(:,:) = 0._wp 603 529 604 530 !------------------------------------------------------------------------------! … … 610 536 ! 3) Rothrock(1975)'s method 611 537 !------------------------------------------------------------------------------! 612 IF (kstrngth == 1) then613 538 IF( kstrngth == 1 ) THEN 539 z1_3 = 1._wp / 3._wp 614 540 DO jl = 1, jpl 615 541 DO jj= 1, jpj 616 542 DO ji = 1, jpi 617 618 IF( ( a_i(ji,jj,jl) .GT. epsi11 )&619 .AND. ( athorn(ji,jj,jl) .GT. 0.0 )) THEN543 ! 544 IF( a_i(ji,jj,jl) > epsi11 .AND. & 545 athorn(ji,jj,jl) > 0._wp ) THEN 620 546 hi = v_i(ji,jj,jl) / a_i(ji,jj,jl) 621 547 !---------------------------- 622 548 ! PE loss from deforming ice 623 549 !---------------------------- 624 strength(ji,jj) = strength(ji,jj) - athorn(ji,jj,jl) * & 625 hi * hi 550 strength(ji,jj) = strength(ji,jj) - athorn(ji,jj,jl) * hi * hi 626 551 627 552 !-------------------------- 628 553 ! PE gain from rafting ice 629 554 !-------------------------- 630 strength(ji,jj) = strength(ji,jj) + 2.0 * araft(ji,jj,jl) & 631 * hi * hi 555 strength(ji,jj) = strength(ji,jj) + 2._wp * araft(ji,jj,jl) * hi * hi 632 556 633 557 !---------------------------- 634 558 ! PE gain from ridging ice 635 559 !---------------------------- 636 strength(ji,jj) = strength(ji,jj) & 637 + aridge(ji,jj,jl)/krdg(ji,jj,jl) & 638 * 1.0/3.0 * (hrmax(ji,jj,jl)**3 - hrmin(ji,jj,jl)**3) & 639 / (hrmax(ji,jj,jl)-hrmin(ji,jj,jl)) 560 strength(ji,jj) = strength(ji,jj) + aridge(ji,jj,jl)/krdg(ji,jj,jl) & 561 * z1_3 * (hrmax(ji,jj,jl)**3 - hrmin(ji,jj,jl)**3) / ( hrmax(ji,jj,jl)-hrmin(ji,jj,jl) ) 562 !!gm Optimization: (a**3-b**3)/(a-b) = a*a+ab+b*b ==> less costly operations even if a**3 is replaced by a*a*a... 640 563 ENDIF ! aicen > epsi11 641 564 ! 642 565 END DO ! ji 643 566 END DO !jj 644 567 END DO !jl 645 568 646 DO jj = 1, jpj 647 DO ji = 1, jpi 648 strength(ji,jj) = Cf * Cp * strength(ji,jj) / aksum(ji,jj) 649 ! Cp = (g/2)*(rhow-rhoi)*(rhoi/rhow) 650 ! Cf accounts for frictional dissipation 651 652 END DO ! j 653 END DO ! i 569 zzc = Cf * Cp ! where Cp = (g/2)*(rhow-rhoi)*(rhoi/rhow) and Cf accounts for frictional dissipation 570 strength(:,:) = zzc * strength(:,:) / aksum(:,:) 654 571 655 572 ksmooth = 1 … … 659 576 !------------------------------------------------------------------------------! 660 577 ELSE ! kstrngth ne 1: Hibler (1979) form 661 662 DO jj = 1, jpj 663 DO ji = 1, jpi 664 strength(ji,jj) = Pstar*vt_i(ji,jj)*exp(-C_rhg*(1.0-at_i(ji,jj))) 665 END DO ! j 666 END DO ! i 667 578 ! 579 strength(:,:) = Pstar * vt_i(:,:) * EXP( - C_rhg * ( 1._wp - at_i(:,:) ) ) 580 ! 668 581 ksmooth = 1 669 582 ! 670 583 ENDIF ! kstrngth 671 584 … … 676 589 ! CAN BE REMOVED 677 590 ! 678 IF ( brinstren_swi .EQ.1 ) THEN591 IF ( brinstren_swi == 1 ) THEN 679 592 680 593 DO jj = 1, jpj … … 699 612 ! Spatial smoothing 700 613 !------------------- 701 IF ( ksmooth .EQ.1 ) THEN614 IF ( ksmooth == 1 ) THEN 702 615 703 616 CALL lbc_lnk( strength, 'T', 1. ) … … 713 626 + strength(ji,jj+1) * tms(ji,jj+1) 714 627 715 zw1 = 4.0 + tms(ji-1,jj) + tms(ji+1,jj) & 716 + tms(ji,jj-1) + tms(ji,jj+1) 628 zw1 = 4.0 + tms(ji-1,jj) + tms(ji+1,jj) + tms(ji,jj-1) + tms(ji,jj+1) 717 629 zworka(ji,jj) = zworka(ji,jj) / zw1 718 630 ELSE … … 734 646 ! Temporal smoothing 735 647 !-------------------- 736 IF ( numit .EQ.nit000 + nn_fsbc - 1 ) THEN648 IF ( numit == nit000 + nn_fsbc - 1 ) THEN 737 649 strp1(:,:) = 0.0 738 650 strp2(:,:) = 0.0 739 651 ENDIF 740 652 741 IF ( ksmooth .EQ.2 ) THEN653 IF ( ksmooth == 2 ) THEN 742 654 743 655 … … 746 658 DO jj = 1, jpj - 1 747 659 DO ji = 1, jpi - 1 748 IF ( ( asum(ji,jj) - ato_i(ji,jj) ) .GT. epsi11) THEN ! ice is 749 ! present 660 IF ( ( asum(ji,jj) - ato_i(ji,jj) ) .GT. epsi11) THEN ! ice is present 750 661 numts_rm = 1 ! number of time steps for the running mean 751 662 IF ( strp1(ji,jj) .GT. 0.0 ) numts_rm = numts_rm + 1 752 663 IF ( strp2(ji,jj) .GT. 0.0 ) numts_rm = numts_rm + 1 753 zp = ( strength(ji,jj) + strp1(ji,jj) + strp2(ji,jj) ) / & 754 numts_rm 664 zp = ( strength(ji,jj) + strp1(ji,jj) + strp2(ji,jj) ) / numts_rm 755 665 strp2(ji,jj) = strp1(ji,jj) 756 666 strp1(ji,jj) = strength(ji,jj) … … 763 673 ENDIF ! ksmooth 764 674 765 ! Boundary conditions 766 CALL lbc_lnk( strength, 'T', 1. ) 767 675 CALL lbc_lnk( strength, 'T', 1. ) ! Boundary conditions 676 677 IF( wrk_not_released(2, 1) ) CALL ctl_stop('lim_itd_me_icestrength: failed to release workspace array') 678 ! 768 679 END SUBROUTINE lim_itd_me_icestrength 769 680 770 !=============================================================================== 771 772 SUBROUTINE lim_itd_me_ridgeprep !(subroutine 3/6) 773 681 682 SUBROUTINE lim_itd_me_ridgeprep 774 683 !!---------------------------------------------------------------------! 775 684 !! *** ROUTINE lim_itd_me_ridgeprep *** 776 !! ** Purpose : 777 !! preparation for ridging and strength calculations 778 !! 779 !! ** Method : 780 !! Compute the thickness distribution of the ice and open water 781 !! participating in ridging and of the resulting ridges. 782 !! 783 !! ** Arguments : 784 !! 785 !! ** External : 786 !! 685 !! 686 !! ** Purpose : preparation for ridging and strength calculations 687 !! 688 !! ** Method : Compute the thickness distribution of the ice and open water 689 !! participating in ridging and of the resulting ridges. 787 690 !!---------------------------------------------------------------------! 788 !! * Arguments 789 790 INTEGER :: & 791 ji,jj, & ! horizontal indices 792 jl, & ! thickness category index 793 krdg_index ! which participation function using 794 795 REAL(wp) :: & 796 Gstari, & ! = 1.0/Gstar 797 astari ! = 1.0/astar 798 799 REAL(wp), DIMENSION(jpi,jpj,-1:jpl) :: & 800 Gsum ! Gsum(n) = sum of areas in categories 0 to n 801 802 REAL(wp) :: & 803 hi, & ! ice thickness for each cat (m) 804 hrmean ! mean ridge thickness (m) 805 806 REAL(wp), DIMENSION(jpi,jpj) :: & 807 zworka ! temporary array used here 808 809 REAL(wp) :: & 810 zdummy, & 811 epsi06 = 1.0e-6 812 691 INTEGER :: ji,jj, jl ! dummy loop indices 692 INTEGER :: krdg_index ! 693 694 REAL(wp) :: Gstari, astari, hi, hrmean, zdummy ! local scalar 695 696 REAL(wp), DIMENSION(jpi,jpj,-1:jpl) :: Gsum ! Gsum(n) = sum of areas in categories 0 to n 697 698 REAL(wp), DIMENSION(jpi,jpj) :: zworka ! temporary array used here 813 699 !------------------------------------------------------------------------------! 814 700 … … 841 727 ! initial value (in h = 0) equals open water area 842 728 843 Gsum(:,:,-1) = 0. 0729 Gsum(:,:,-1) = 0._wp 844 730 845 731 DO jj = 1, jpj 846 732 DO ji = 1, jpi 847 IF (ato_i(ji,jj) .GT. epsi11) THEN 848 Gsum(ji,jj,0) = ato_i(ji,jj) 849 ELSE 850 Gsum(ji,jj,0) = 0.0 733 IF( ato_i(ji,jj) > epsi11 ) THEN ; Gsum(ji,jj,0) = ato_i(ji,jj) 734 ELSE ; Gsum(ji,jj,0) = 0._wp 851 735 ENDIF 852 736 END DO … … 857 741 DO jj = 1, jpj 858 742 DO ji = 1, jpi 859 IF ( a_i(ji,jj,jl) .GT. epsi11 ) THEN 860 Gsum(ji,jj,jl) = Gsum(ji,jj,jl-1) + a_i(ji,jj,jl) 861 ELSE 862 Gsum(ji,jj,jl) = Gsum(ji,jj,jl-1) 743 IF( a_i(ji,jj,jl) .GT. epsi11 ) THEN ; Gsum(ji,jj,jl) = Gsum(ji,jj,jl-1) + a_i(ji,jj,jl) 744 ELSE ; Gsum(ji,jj,jl) = Gsum(ji,jj,jl-1) 863 745 ENDIF 864 746 END DO … … 867 749 868 750 ! Normalize the cumulative distribution to 1 869 DO jj = 1, jpj 870 DO ji = 1, jpi 871 zworka(ji,jj) = 1.0 / Gsum(ji,jj,jpl) 872 END DO 873 END DO 874 751 zworka(:,:) = 1._wp / Gsum(:,:,jpl) 875 752 DO jl = 0, jpl 876 DO jj = 1, jpj 877 DO ji = 1, jpi 878 Gsum(ji,jj,jl) = Gsum(ji,jj,jl) * zworka(ji,jj) 879 END DO 880 END DO 753 Gsum(:,:,jl) = Gsum(:,:,jl) * zworka(:,:) 881 754 END DO 882 755 … … 895 768 krdg_index = 1 896 769 897 IF ( krdg_index .EQ. 0 ) THEN 898 899 !--- Linear formulation (Thorndike et al., 1975) 900 DO jl = 0, ice_cat_bounds(1,2) ! only undeformed ice participates 770 IF( krdg_index == 0 ) THEN !--- Linear formulation (Thorndike et al., 1975) 771 DO jl = 0, ice_cat_bounds(1,2) ! only undeformed ice participates 901 772 DO jj = 1, jpj 902 773 DO ji = 1, jpi 903 IF (Gsum(ji,jj,jl) < Gstar) THEN774 IF( Gsum(ji,jj,jl) < Gstar) THEN 904 775 athorn(ji,jj,jl) = Gstari * (Gsum(ji,jj,jl)-Gsum(ji,jj,jl-1)) * & 905 776 (2.0 - (Gsum(ji,jj,jl-1)+Gsum(ji,jj,jl))*Gstari) … … 914 785 END DO ! jl 915 786 916 ELSE ! krdg_index = 1 917 918 !--- Exponential, more stable formulation (Lipscomb et al, 2007) 919 ! precompute exponential terms using Gsum as a work array 920 zdummy = 1.0 / (1.0-EXP(-astari)) 787 ELSE !--- Exponential, more stable formulation (Lipscomb et al, 2007) 788 ! 789 zdummy = 1._wp / ( 1._wp - EXP(-astari) ) ! precompute exponential terms using Gsum as a work array 921 790 922 791 DO jl = -1, jpl 923 DO jj = 1, jpj 924 DO ji = 1, jpi 925 Gsum(ji,jj,jl) = EXP(-Gsum(ji,jj,jl)*astari)*zdummy 926 END DO !ji 927 END DO !jj 792 Gsum(:,:,jl) = EXP( -Gsum(:,:,jl) * astari ) * zdummy 928 793 END DO !jl 929 930 ! compute athorn931 794 DO jl = 0, ice_cat_bounds(1,2) 932 DO jj = 1, jpj 933 DO ji = 1, jpi 934 athorn(ji,jj,jl) = Gsum(ji,jj,jl-1) - Gsum(ji,jj,jl) 935 END DO !ji 936 END DO ! jj 937 END DO !jl 938 795 athorn(:,:,jl) = Gsum(:,:,jl-1) - Gsum(:,:,jl) 796 END DO 797 ! 939 798 ENDIF ! krdg_index 940 799 941 ! Ridging and rafting ice participation functions 942 IF ( raftswi .EQ. 1 ) THEN 943 800 IF( raftswi == 1 ) THEN ! Ridging and rafting ice participation functions 801 ! 944 802 DO jl = 1, jpl 945 803 DO jj = 1, jpj 946 804 DO ji = 1, jpi 947 IF ( athorn(ji,jj,jl) .GT. 0.0 ) THEN 948 aridge(ji,jj,jl) = ( TANH ( Craft * ( ht_i(ji,jj,jl) - & 949 hparmeter ) ) + 1.0 ) / 2.0 * & 950 athorn(ji,jj,jl) 951 araft (ji,jj,jl) = ( TANH ( - Craft * ( ht_i(ji,jj,jl) - & 952 hparmeter ) ) + 1.0 ) / 2.0 * & 953 athorn(ji,jj,jl) 954 IF ( araft(ji,jj,jl) .LT. epsi06 ) araft(ji,jj,jl) = 0.0 955 aridge(ji,jj,jl) = MAX( athorn(ji,jj,jl) - araft(ji,jj,jl), 0.0) 805 IF ( athorn(ji,jj,jl) .GT. 0._wp ) THEN 806 !!gm TANH( -X ) = - TANH( X ) so can be computed only 1 time.... 807 aridge(ji,jj,jl) = ( TANH ( Craft * ( ht_i(ji,jj,jl) - hparmeter ) ) + 1.0 ) * 0.5 * athorn(ji,jj,jl) 808 araft (ji,jj,jl) = ( TANH ( -Craft * ( ht_i(ji,jj,jl) - hparmeter ) ) + 1.0 ) * 0.5 * athorn(ji,jj,jl) 809 IF ( araft(ji,jj,jl) < epsi06 ) araft(ji,jj,jl) = 0._wp 810 aridge(ji,jj,jl) = MAX( athorn(ji,jj,jl) - araft(ji,jj,jl), 0.0 ) 956 811 ENDIF ! athorn 957 812 END DO ! ji … … 960 815 961 816 ELSE ! raftswi = 0 962 817 ! 963 818 DO jl = 1, jpl 964 DO jj = 1, jpj 965 DO ji = 1, jpi 966 aridge(ji,jj,jl) = 1.0*athorn(ji,jj,jl) 967 END DO 968 END DO 969 END DO 970 819 aridge(:,:,jl) = athorn(:,:,jl) 820 END DO 821 ! 971 822 ENDIF 972 823 973 IF ( raftswi .EQ.1 ) THEN824 IF ( raftswi == 1 ) THEN 974 825 975 826 IF( MAXVAL(aridge + araft - athorn(:,:,1:jpl)) .GT. epsi11 ) THEN … … 1043 894 1044 895 ! Normalization factor : aksum, ensures mass conservation 1045 DO jj = 1, jpj1046 DO ji = 1, jpi1047 aksum(ji,jj) = athorn(ji,jj,0)1048 END DO896 aksum(:,:) = athorn(ji,jj,0) 897 DO jl = 1, jpl 898 aksum(:,:) = aksum(:,:) + aridge(:,:,jl) * ( 1._wp - 1._wp / krdg(:,:,jl) ) & 899 & + araft (:,:,jl) * ( 1._wp - 1._wp / kraft ) 1049 900 END DO 1050 1051 DO jl = 1, jpl 1052 DO jj = 1, jpj 1053 DO ji = 1, jpi 1054 aksum(ji,jj) = aksum(ji,jj) & 1055 + aridge(ji,jj,jl) * (1.0 - 1.0/krdg(ji,jj,jl)) & 1056 + araft (ji,jj,jl) * (1.0 - 1.0/kraft) 1057 END DO 1058 END DO 1059 END DO 1060 901 ! 1061 902 END SUBROUTINE lim_itd_me_ridgeprep 1062 903 1063 !=============================================================================== 1064 1065 SUBROUTINE lim_itd_me_ridgeshift(opning, closing_gross, & 1066 msnow_mlt, esnow_mlt) ! (subroutine 4/6) 1067 1068 !!----------------------------------------------------------------------------- 904 905 SUBROUTINE lim_itd_me_ridgeshift( opning, closing_gross, msnow_mlt, esnow_mlt ) 906 !!---------------------------------------------------------------------- 1069 907 !! *** ROUTINE lim_itd_me_icestrength *** 1070 !! ** Purpose : 1071 !! This routine shift ridging ice among thickness categories 1072 !! of ice thickness 1073 !! 1074 !! ** Method : 1075 !! Remove area, volume, and energy from each ridging category 1076 !! and add to thicker ice categories. 1077 !! 1078 !! ** Arguments : 1079 !! 1080 !! ** Inputs / Ouputs : 1081 !! 1082 !! ** External : 1083 !! 1084 1085 REAL (wp), DIMENSION(jpi,jpj), INTENT(IN) :: & 1086 opning, & ! rate of opening due to divergence/shear 1087 closing_gross ! rate at which area removed, not counting 1088 ! area of new ridges 1089 1090 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: & 1091 msnow_mlt, & ! mass of snow added to ocean (kg m-2) 1092 esnow_mlt ! energy needed to melt snow in ocean (J m-2) 1093 1094 INTEGER :: & 1095 ji, jj, & ! horizontal indices 1096 jl, jl1, jl2, & ! thickness category indices 1097 jk, & ! ice layer index 1098 ij, & ! horizontal index, combines i and j loops 1099 icells ! number of cells with aicen > puny 1100 1101 INTEGER, DIMENSION(1:(jpi+1)*(jpj+1)) :: & 1102 indxi, indxj ! compressed indices 1103 1104 REAL(wp), DIMENSION(jpi,jpj) :: & 1105 vice_init, vice_final, & ! ice volume summed over categories 1106 eice_init, eice_final ! ice energy summed over layers 1107 1108 REAL(wp), DIMENSION(jpi,jpj,jpl) :: & 1109 aicen_init, & ! ice area before ridging 1110 vicen_init, & ! ice volume before ridging 1111 vsnon_init, & ! snow volume before ridging 1112 esnon_init, & ! snow energy before ridging 1113 smv_i_init, & ! ice salinity before ridging 1114 oa_i_init ! ice age before ridging 1115 1116 REAL(wp), DIMENSION(jpi,jpj,jkmax,jpl) :: & 1117 eicen_init ! ice energy before ridging 1118 1119 REAL(wp), DIMENSION(jpi,jpj) :: & 1120 afrac , & ! fraction of category area ridged 1121 ardg1 , & ! area of ice ridged 1122 ardg2 , & ! area of new ridges 1123 vsrdg , & ! snow volume of ridging ice 1124 esrdg , & ! snow energy of ridging ice 1125 oirdg1 , & ! areal age content of ridged ice 1126 oirdg2 , & ! areal age content of ridging ice 1127 dhr , & ! hrmax - hrmin 1128 dhr2 , & ! hrmax^2 - hrmin^2 1129 fvol ! fraction of new ridge volume going to n2 1130 1131 REAL(wp), DIMENSION(jpi,jpj) :: & 1132 vrdg1 , & ! volume of ice ridged 1133 vrdg2 , & ! volume of new ridges 1134 vsw , & ! volume of seawater trapped into ridges 1135 srdg1 , & ! sal*volume of ice ridged 1136 srdg2 , & ! sal*volume of new ridges 1137 smsw ! sal*volume of water trapped into ridges 1138 1139 REAL(wp), DIMENSION(jpi,jpj) :: & 1140 afrft , & ! fraction of category area rafted 1141 arft1 , & ! area of ice rafted 1142 arft2 , & ! area of new rafted zone 1143 virft , & ! ice volume of rafting ice 1144 vsrft , & ! snow volume of rafting ice 1145 esrft , & ! snow energy of rafting ice 1146 smrft , & ! salinity of rafting ice 1147 oirft1 , & ! areal age content of rafted ice 1148 oirft2 ! areal age content of rafting ice 1149 1150 REAL(wp), DIMENSION(jpi,jpj,jkmax) :: & 1151 eirft , & ! ice energy of rafting ice 1152 erdg1 , & ! enth*volume of ice ridged 1153 erdg2 , & ! enth*volume of new ridges 1154 ersw ! enth of water trapped into ridges 1155 1156 REAL(wp) :: & 1157 hL, hR , & ! left and right limits of integration 1158 farea , & ! fraction of new ridge area going to n2 1159 zdummy , & ! dummy argument 1160 zdummy0 , & ! dummy argument 1161 ztmelts ! ice melting point 1162 1163 REAL(wp) :: zsrdg2 1164 1165 CHARACTER (len=80) :: & 1166 fieldid ! field identifier 1167 1168 LOGICAL, PARAMETER :: & 1169 l_conservation_check = .true. ! if true, check conservation 1170 ! (useful for debugging) 1171 LOGICAL :: & 1172 neg_ato_i , & ! flag for ato_i(i,j) < -puny 1173 large_afrac , & ! flag for afrac > 1 1174 large_afrft ! flag for afrac > 1 1175 1176 REAL(wp) :: & 1177 zeps , & 1178 epsi10 , & 1179 zindb ! switch for the presence of ridge poros or not 1180 1181 !---------------------------------------------------------------------------- 908 !! 909 !! ** Purpose : shift ridging ice among thickness categories of ice thickness 910 !! 911 !! ** Method : Remove area, volume, and energy from each ridging category 912 !! and add to thicker ice categories. 913 !!---------------------------------------------------------------------- 914 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: opning ! rate of opening due to divergence/shear 915 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: closing_gross ! rate at which area removed, excluding area of new ridges 916 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: msnow_mlt ! mass of snow added to ocean (kg m-2) 917 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: esnow_mlt ! energy needed to melt snow in ocean (J m-2) 918 ! 919 CHARACTER (len=80) :: fieldid ! field identifier 920 LOGICAL, PARAMETER :: l_conservation_check = .true. ! if true, check conservation (useful for debugging) 921 ! 922 LOGICAL :: neg_ato_i ! flag for ato_i(i,j) < -puny 923 LOGICAL :: large_afrac ! flag for afrac > 1 924 LOGICAL :: large_afrft ! flag for afrac > 1 925 INTEGER :: ji, jj, jl, jl1, jl2, jk ! dummy loop indices 926 INTEGER :: ij ! horizontal index, combines i and j loops 927 INTEGER :: icells ! number of cells with aicen > puny 928 REAL(wp) :: zeps, zindb, zsrdg2 ! local scalar 929 REAL(wp) :: hL, hR, farea, zdummy, zdummy0, ztmelts ! left and right limits of integration 930 931 INTEGER, DIMENSION(1:(jpi+1)*(jpj+1)) :: indxi, indxj ! compressed indices 932 933 REAL(wp), DIMENSION(jpi,jpj) :: vice_init, vice_final ! ice volume summed over categories 934 REAL(wp), DIMENSION(jpi,jpj) :: eice_init, eice_final ! ice energy summed over layers 935 936 REAL(wp), DIMENSION(jpi,jpj,jpl) :: aicen_init, vicen_init ! ice area & volume before ridging 937 REAL(wp), DIMENSION(jpi,jpj,jpl) :: vsnon_init, esnon_init ! snow volume & energy before ridging 938 REAL(wp), DIMENSION(jpi,jpj,jpl) :: smv_i_init, oa_i_init ! ice salinity & age before ridging 939 940 REAL(wp), DIMENSION(jpi,jpj,jkmax,jpl) :: eicen_init ! ice energy before ridging 941 942 REAL(wp), DIMENSION(jpi,jpj) :: afrac , fvol ! fraction of category area ridged & new ridge volume going to n2 943 REAL(wp), DIMENSION(jpi,jpj) :: ardg1 , ardg2 ! area of ice ridged & new ridges 944 REAL(wp), DIMENSION(jpi,jpj) :: vsrdg , esrdg ! snow volume & energy of ridging ice 945 REAL(wp), DIMENSION(jpi,jpj) :: oirdg1, oirdg2 ! areal age content of ridged & rifging ice 946 REAL(wp), DIMENSION(jpi,jpj) :: dhr , dhr2 ! hrmax - hrmin & hrmax^2 - hrmin^2 947 948 REAL(wp), DIMENSION(jpi,jpj) :: vrdg1 ! volume of ice ridged 949 REAL(wp), DIMENSION(jpi,jpj) :: vrdg2 ! volume of new ridges 950 REAL(wp), DIMENSION(jpi,jpj) :: vsw ! volume of seawater trapped into ridges 951 REAL(wp), DIMENSION(jpi,jpj) :: srdg1 ! sal*volume of ice ridged 952 REAL(wp), DIMENSION(jpi,jpj) :: srdg2 ! sal*volume of new ridges 953 REAL(wp), DIMENSION(jpi,jpj) :: smsw ! sal*volume of water trapped into ridges 954 955 REAL(wp), DIMENSION(jpi,jpj) :: afrft ! fraction of category area rafted 956 REAL(wp), DIMENSION(jpi,jpj) :: arft1 , arft2 ! area of ice rafted and new rafted zone 957 REAL(wp), DIMENSION(jpi,jpj) :: virft , vsrft ! ice & snow volume of rafting ice 958 REAL(wp), DIMENSION(jpi,jpj) :: esrft , smrft ! snow energy & salinity of rafting ice 959 REAL(wp), DIMENSION(jpi,jpj) :: oirft1, oirft2 ! areal age content of rafted ice & rafting ice 960 961 REAL(wp), DIMENSION(jpi,jpj,jkmax) :: eirft ! ice energy of rafting ice 962 REAL(wp), DIMENSION(jpi,jpj,jkmax) :: erdg1 ! enth*volume of ice ridged 963 REAL(wp), DIMENSION(jpi,jpj,jkmax) :: erdg2 ! enth*volume of new ridges 964 REAL(wp), DIMENSION(jpi,jpj,jkmax) :: ersw ! enth of water trapped into ridges 965 !!---------------------------------------------------------------------- 1182 966 1183 967 ! Conservation check 1184 eice_init(:,:) = 0. 01185 1186 IF 968 eice_init(:,:) = 0._wp 969 970 IF( con_i ) THEN 1187 971 CALL lim_column_sum (jpl, v_i, vice_init ) 1188 972 WRITE(numout,*) ' vice_init : ', vice_init(jiindx,jjindx) … … 1191 975 ENDIF 1192 976 1193 zeps = 1.0d-20 1194 epsi10 = 1.0d-10 977 zeps = 1.e-20_wp 1195 978 1196 979 !------------------------------------------------------------------------------- … … 1202 985 DO jj = 1, jpj 1203 986 DO ji = 1, jpi 1204 ato_i(ji,jj) = ato_i(ji,jj) & 1205 - athorn(ji,jj,0)*closing_gross(ji,jj)*rdt_ice & 1206 + opning(ji,jj)*rdt_ice 1207 IF (ato_i(ji,jj) .LT. -epsi11) THEN 1208 neg_ato_i = .true. 1209 ELSEIF (ato_i(ji,jj) .LT. 0.0) THEN ! roundoff error 1210 ato_i(ji,jj) = 0.0 987 ato_i(ji,jj) = ato_i(ji,jj) - athorn(ji,jj,0) * closing_gross(ji,jj) * rdt_ice & 988 & + opning(ji,jj) * rdt_ice 989 IF( ato_i(ji,jj) < -epsi11 ) THEN 990 neg_ato_i = .TRUE. 991 ELSEIF( ato_i(ji,jj) < 0._wp ) THEN ! roundoff error 992 ato_i(ji,jj) = 0._wp 1211 993 ENDIF 1212 994 END DO !jj … … 1214 996 1215 997 ! if negative open water area alert it 1216 IF (neg_ato_i) THEN ! there is a bug998 IF( neg_ato_i ) THEN ! there is a bug 1217 999 DO jj = 1, jpj 1218 1000 DO ji = 1, jpi 1219 IF (ato_i(ji,jj) .LT. -epsi11) THEN1001 IF( ato_i(ji,jj) < -epsi11 ) THEN 1220 1002 WRITE(numout,*) '' 1221 1003 WRITE(numout,*) 'Ridging error: ato_i < 0' 1222 1004 WRITE(numout,*) 'ato_i : ', ato_i(ji,jj) 1223 1005 ENDIF ! ato_i < -epsi11 1224 END DO ! ji1225 END DO ! jj1226 ENDIF ! neg_ato_i1006 END DO 1007 END DO 1008 ENDIF 1227 1009 1228 1010 !----------------------------------------------------------------- … … 1231 1013 1232 1014 DO jl = 1, jpl 1233 DO jj = 1, jpj 1234 DO ji = 1, jpi 1235 aicen_init(ji,jj,jl) = a_i(ji,jj,jl) 1236 vicen_init(ji,jj,jl) = v_i(ji,jj,jl) 1237 vsnon_init(ji,jj,jl) = v_s(ji,jj,jl) 1238 1239 smv_i_init(ji,jj,jl) = smv_i(ji,jj,jl) 1240 oa_i_init (ji,jj,jl) = oa_i(ji,jj,jl) 1241 END DO !ji 1242 END DO ! jj 1015 aicen_init(:,:,jl) = a_i(:,:,jl) 1016 vicen_init(:,:,jl) = v_i(:,:,jl) 1017 vsnon_init(:,:,jl) = v_s(:,:,jl) 1018 ! 1019 smv_i_init(:,:,jl) = smv_i(:,:,jl) 1020 oa_i_init (:,:,jl) = oa_i (:,:,jl) 1243 1021 END DO !jl 1244 1022 … … 1247 1025 DO jl = 1, jpl 1248 1026 DO jk = 1, nlay_i 1249 DO jj = 1, jpj 1250 DO ji = 1, jpi 1251 eicen_init(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) 1252 END DO !ji 1253 END DO !jj 1254 END DO !jk 1255 END DO !jl 1027 eicen_init(:,:,jk,jl) = e_i(:,:,jk,jl) 1028 END DO 1029 END DO 1256 1030 1257 1031 ! … … 1324 1098 ! / rafting category n1. 1325 1099 !-------------------------------------------------------------------------- 1326 vrdg1(ji,jj) = vicen_init(ji,jj,jl1) * afrac(ji,jj) / & 1327 ( 1.0 + ridge_por ) 1100 vrdg1(ji,jj) = vicen_init(ji,jj,jl1) * afrac(ji,jj) / ( 1._wp + ridge_por ) 1328 1101 vrdg2(ji,jj) = vrdg1(ji,jj) * ( 1. + ridge_por ) 1329 1102 vsw (ji,jj) = vrdg1(ji,jj) * ridge_por … … 1331 1104 vsrdg(ji,jj) = vsnon_init(ji,jj,jl1) * afrac(ji,jj) 1332 1105 esrdg(ji,jj) = esnon_init(ji,jj,jl1) * afrac(ji,jj) 1333 srdg1(ji,jj) = smv_i_init(ji,jj,jl1) * afrac(ji,jj) / & 1334 ( 1. + ridge_por ) 1106 srdg1(ji,jj) = smv_i_init(ji,jj,jl1) * afrac(ji,jj) / ( 1._wp + ridge_por ) 1335 1107 srdg2(ji,jj) = smv_i_init(ji,jj,jl1) * afrac(ji,jj) 1336 1108 … … 1371 1143 ! ij looping 1-icells 1372 1144 1373 dardg1dt (ji,jj) = dardg1dt(ji,jj) + ardg1(ji,jj) + arft1(ji,jj)1374 dardg2dt (ji,jj) = dardg2dt(ji,jj) + ardg2(ji,jj) + arft2(ji,jj)1145 dardg1dt (ji,jj) = dardg1dt(ji,jj) + ardg1(ji,jj) + arft1(ji,jj) 1146 dardg2dt (ji,jj) = dardg2dt(ji,jj) + ardg2(ji,jj) + arft2(ji,jj) 1375 1147 diag_dyn_gr(ji,jj) = diag_dyn_gr(ji,jj) + ( vrdg2(ji,jj) + virft(ji,jj) ) / rdt_ice 1376 opening (ji,jj)= opening (ji,jj) + opning(ji,jj)*rdt_ice1377 1378 IF (con_i)vice_init(ji,jj) = vice_init(ji,jj) + vrdg2(ji,jj) - vrdg1(ji,jj)1148 opening (ji,jj) = opening (ji,jj) + opning(ji,jj)*rdt_ice 1149 1150 IF( con_i ) vice_init(ji,jj) = vice_init(ji,jj) + vrdg2(ji,jj) - vrdg1(ji,jj) 1379 1151 1380 1152 !------------------------------------------ … … 1390 1162 ! ij looping 1-icells 1391 1163 1392 msnow_mlt(ji,jj) = msnow_mlt(ji,jj) & 1393 + rhosn*vsrdg(ji,jj)*(1.0-fsnowrdg) & 1394 !rafting included 1395 + rhosn*vsrft(ji,jj)*(1.0-fsnowrft) 1396 1397 esnow_mlt(ji,jj) = esnow_mlt(ji,jj) & 1398 + esrdg(ji,jj)*(1.0-fsnowrdg) & 1399 !rafting included 1400 + esrft(ji,jj)*(1.0-fsnowrft) 1164 msnow_mlt(ji,jj) = msnow_mlt(ji,jj) + rhosn*vsrdg(ji,jj)*(1.0-fsnowrdg) & ! rafting included 1165 & + rhosn*vsrft(ji,jj)*(1.0-fsnowrft) 1166 1167 esnow_mlt(ji,jj) = esnow_mlt(ji,jj) + esrdg(ji,jj)*(1.0-fsnowrdg) & !rafting included 1168 & + esrft(ji,jj)*(1.0-fsnowrft) 1401 1169 1402 1170 !----------------------------------------------------------------- … … 1409 1177 1410 1178 dhr(ji,jj) = hrmax(ji,jj,jl1) - hrmin(ji,jj,jl1) 1411 dhr2(ji,jj) = hrmax(ji,jj,jl1) * hrmax(ji,jj,jl1) & 1412 - hrmin(ji,jj,jl1) * hrmin(ji,jj,jl1) 1179 dhr2(ji,jj) = hrmax(ji,jj,jl1) * hrmax(ji,jj,jl1) - hrmin(ji,jj,jl1) * hrmin(ji,jj,jl1) 1413 1180 1414 1181 … … 1425 1192 jj = indxj(ij) 1426 1193 ! heat content of ridged ice 1427 erdg1(ji,jj,jk) = eicen_init(ji,jj,jk,jl1) * afrac(ji,jj) / & 1428 ( 1.0 + ridge_por ) 1194 erdg1(ji,jj,jk) = eicen_init(ji,jj,jk,jl1) * afrac(ji,jj) / ( 1._wp + ridge_por ) 1429 1195 eirft(ji,jj,jk) = eicen_init(ji,jj,jk,jl1) * afrft(ji,jj) 1430 e_i(ji,jj,jk,jl1) = e_i(ji,jj,jk,jl1) & 1431 - erdg1(ji,jj,jk) & 1432 - eirft(ji,jj,jk) 1196 e_i (ji,jj,jk,jl1) = e_i(ji,jj,jk,jl1) - erdg1(ji,jj,jk) - eirft(ji,jj,jk) 1433 1197 ! sea water heat content 1434 1198 ztmelts = - tmut * sss_m(ji,jj) + rtt … … 1437 1201 1438 1202 ! corrected sea water salinity 1439 zindb = MAX( 0.0, SIGN( 1.0, vsw(ji,jj) - zeps ) ) 1440 zdummy = zindb * ( srdg1(ji,jj) - srdg2(ji,jj) ) / & 1441 MAX( ridge_por * vsw(ji,jj), zeps ) 1203 zindb = MAX( 0._wp , SIGN( 1._wp , vsw(ji,jj) - zeps ) ) 1204 zdummy = zindb * ( srdg1(ji,jj) - srdg2(ji,jj) ) / MAX( ridge_por * vsw(ji,jj), zeps ) 1442 1205 1443 1206 ztmelts = - tmut * zdummy + rtt … … 1445 1208 1446 1209 ! heat flux 1447 fheat_rpo(ji,jj) = fheat_rpo(ji,jj) + ( ersw(ji,jj,jk) - zdummy0 ) / & 1448 rdt_ice 1210 fheat_rpo(ji,jj) = fheat_rpo(ji,jj) + ( ersw(ji,jj,jk) - zdummy0 ) / rdt_ice 1449 1211 1450 1212 ! Correct dimensions to avoid big values 1451 ersw(ji,jj,jk) = ersw(ji,jj,jk) / 1.0d+091213 ersw(ji,jj,jk) = ersw(ji,jj,jk) * 1.e-09 1452 1214 1453 1215 ! Mutliply by ice volume, and divide by number of layers to get heat content in 10^9 J 1454 ersw(ji,jj,jk) = ersw(ji,jj,jk) * & 1455 area(ji,jj) * vsw(ji,jj) / & 1456 nlay_i 1216 ersw (ji,jj,jk) = ersw(ji,jj,jk) * area(ji,jj) * vsw(ji,jj) / nlay_i 1457 1217 1458 1218 erdg2(ji,jj,jk) = erdg1(ji,jj,jk) + ersw(ji,jj,jk) … … 1461 1221 1462 1222 1463 IF 1223 IF( con_i ) THEN 1464 1224 DO jk = 1, nlay_i 1465 1225 !CDIR NODEP … … 1467 1227 ji = indxi(ij) 1468 1228 jj = indxj(ij) 1469 eice_init(ji,jj) = eice_init(ji,jj) + erdg2(ji,jj,jk) - & 1470 erdg1(ji,jj,jk) 1229 eice_init(ji,jj) = eice_init(ji,jj) + erdg2(ji,jj,jk) - erdg1(ji,jj,jk) 1471 1230 END DO ! ij 1472 1231 END DO !jk 1473 1232 ENDIF 1474 1233 1475 IF (large_afrac) THEN! there is a bug1234 IF( large_afrac ) THEN ! there is a bug 1476 1235 !CDIR NODEP 1477 1236 DO ij = 1, icells 1478 1237 ji = indxi(ij) 1479 1238 jj = indxj(ij) 1480 IF 1239 IF( afrac(ji,jj) > 1.0 + epsi11 ) THEN 1481 1240 WRITE(numout,*) '' 1482 1241 WRITE(numout,*) ' ardg > a_i' 1483 WRITE(numout,*) ' ardg, aicen_init : ', & 1484 ardg1(ji,jj), aicen_init(ji,jj,jl1) 1485 ENDIF ! afrac > 1 + puny 1486 ENDDO ! if 1487 ENDIF ! large_afrac 1488 IF (large_afrft) THEN ! there is a bug 1242 WRITE(numout,*) ' ardg, aicen_init : ', ardg1(ji,jj), aicen_init(ji,jj,jl1) 1243 ENDIF 1244 END DO 1245 ENDIF 1246 IF( large_afrft ) THEN ! there is a bug 1489 1247 !CDIR NODEP 1490 1248 DO ij = 1, icells 1491 1249 ji = indxi(ij) 1492 1250 jj = indxj(ij) 1493 IF 1251 IF( afrft(ji,jj) > 1.0 + epsi11 ) THEN 1494 1252 WRITE(numout,*) '' 1495 1253 WRITE(numout,*) ' arft > a_i' 1496 WRITE(numout,*) ' arft, aicen_init : ', & 1497 arft1(ji,jj), aicen_init(ji,jj,jl1) 1498 ENDIF ! afrft > 1 + puny 1499 ENDDO ! if 1500 ENDIF ! large_afrft 1254 WRITE(numout,*) ' arft, aicen_init : ', arft1(ji,jj), aicen_init(ji,jj,jl1) 1255 ENDIF 1256 END DO 1257 ENDIF 1501 1258 1502 1259 !------------------------------------------------------------------------------- … … 1528 1285 fvol(ji,jj) = (hR*hR - hL*hL) / dhr2(ji,jj) 1529 1286 1530 a_i(ji,jj,jl2) = a_i(ji,jj,jl2) + farea * ardg2(ji,jj) 1531 v_i(ji,jj,jl2) = v_i(ji,jj,jl2) + fvol(ji,jj) * vrdg2(ji,jj) 1532 v_s(ji,jj,jl2) = v_s(ji,jj,jl2) & 1533 + fvol(ji,jj) * vsrdg(ji,jj) * fsnowrdg 1534 e_s(ji,jj,1,jl2) = e_s(ji,jj,1,jl2) & 1535 + fvol(ji,jj) * esrdg(ji,jj) * fsnowrdg 1536 smv_i(ji,jj,jl2) = smv_i(ji,jj,jl2) + fvol(ji,jj) * srdg2(ji,jj) 1537 oa_i(ji,jj,jl2) = oa_i(ji,jj,jl2) + farea * oirdg2(ji,jj) 1287 a_i (ji,jj,jl2) = a_i (ji,jj,jl2) + ardg2 (ji,jj) * farea 1288 v_i (ji,jj,jl2) = v_i (ji,jj,jl2) + vrdg2 (ji,jj) * fvol(ji,jj) 1289 v_s (ji,jj,jl2) = v_s (ji,jj,jl2) + vsrdg (ji,jj) * fvol(ji,jj) * fsnowrdg 1290 e_s (ji,jj,1,jl2) = e_s (ji,jj,1,jl2) + esrdg (ji,jj) * fvol(ji,jj) * fsnowrdg 1291 smv_i(ji,jj,jl2) = smv_i(ji,jj,jl2) + srdg2 (ji,jj) * fvol(ji,jj) 1292 oa_i (ji,jj,jl2) = oa_i (ji,jj,jl2) + oirdg2(ji,jj) * farea 1538 1293 1539 1294 END DO ! ij … … 1545 1300 ji = indxi(ij) 1546 1301 jj = indxj(ij) 1547 e_i(ji,jj,jk,jl2) = e_i(ji,jj,jk,jl2) & 1548 + fvol(ji,jj)*erdg2(ji,jj,jk) 1549 END DO ! ij 1550 END DO !jk 1551 1552 1302 e_i(ji,jj,jk,jl2) = e_i(ji,jj,jk,jl2) + fvol(ji,jj)*erdg2(ji,jj,jk) 1303 END DO 1304 END DO 1305 ! 1553 1306 END DO ! jl2 (new ridges) 1554 1307 1555 DO jl2 1308 DO jl2 = ice_cat_bounds(1,1), ice_cat_bounds(1,2) 1556 1309 1557 1310 !CDIR NODEP … … 1566 1319 a_i(ji,jj,jl2) = a_i(ji,jj,jl2) + arft2(ji,jj) 1567 1320 v_i(ji,jj,jl2) = v_i(ji,jj,jl2) + virft(ji,jj) 1568 v_s(ji,jj,jl2) = v_s(ji,jj,jl2) & 1569 + vsrft(ji,jj)*fsnowrft 1570 e_s(ji,jj,1,jl2) = e_s(ji,jj,1,jl2) & 1571 + esrft(ji,jj)*fsnowrft 1572 smv_i(ji,jj,jl2) = smv_i(ji,jj,jl2) & 1573 + smrft(ji,jj) 1574 oa_i(ji,jj,jl2) = oa_i(ji,jj,jl2) & 1575 + oirft2(ji,jj) 1321 v_s(ji,jj,jl2) = v_s(ji,jj,jl2) + vsrft(ji,jj)*fsnowrft 1322 e_s(ji,jj,1,jl2) = e_s(ji,jj,1,jl2) + esrft(ji,jj)*fsnowrft 1323 smv_i(ji,jj,jl2) = smv_i(ji,jj,jl2) + smrft(ji,jj) 1324 oa_i(ji,jj,jl2) = oa_i(ji,jj,jl2) + oirft2(ji,jj) 1576 1325 ENDIF ! hraft 1577 1326 … … 1586 1335 IF (hraft(ji,jj,jl1) .LE. hi_max(jl2) .AND. & 1587 1336 hraft(ji,jj,jl1) .GT. hi_max(jl2-1)) THEN 1588 e_i(ji,jj,jk,jl2) = e_i(ji,jj,jk,jl2) & 1589 + eirft(ji,jj,jk) 1337 e_i(ji,jj,jk,jl2) = e_i(ji,jj,jk,jl2) + eirft(ji,jj,jk) 1590 1338 ENDIF 1591 1339 END DO ! ij … … 1610 1358 WRITE(numout,*) ' eice_final : ', eice_final(jiindx,jjindx) 1611 1359 ENDIF 1612 1360 ! 1613 1361 END SUBROUTINE lim_itd_me_ridgeshift 1614 1362 1615 !============================================================================== 1616 1617 SUBROUTINE lim_itd_me_asumr !(subroutine 5/6) 1618 1363 1364 SUBROUTINE lim_itd_me_asumr 1619 1365 !!----------------------------------------------------------------------------- 1620 1366 !! *** ROUTINE lim_itd_me_asumr *** 1621 !! ** Purpose : 1622 !! This routine finds total fractional area 1623 !! 1624 !! ** Method : 1625 !! Find the total area of ice plus open water in each grid cell. 1626 !! 1627 !! This is similar to the aggregate_area subroutine except that the 1628 !! total area can be greater than 1, so the open water area is 1629 !! included in the sum instead of being computed as a residual. 1630 !! 1631 !! ** Arguments : 1632 1633 INTEGER :: ji, jj, jl 1634 1635 !----------------------------------------------------------------- 1636 ! open water 1637 !----------------------------------------------------------------- 1638 1639 DO jj = 1, jpj 1640 DO ji = 1, jpi 1641 asum(ji,jj) = ato_i(ji,jj) 1642 END DO 1367 !! 1368 !! ** Purpose : finds total fractional area 1369 !! 1370 !! ** Method : Find the total area of ice plus open water in each grid cell. 1371 !! This is similar to the aggregate_area subroutine except that the 1372 !! total area can be greater than 1, so the open water area is 1373 !! included in the sum instead of being computed as a residual. 1374 !!----------------------------------------------------------------------------- 1375 INTEGER :: jl ! dummy loop index 1376 !!----------------------------------------------------------------------------- 1377 ! 1378 asum(:,:) = ato_i(:,:) ! open water 1379 DO jl = 1, jpl ! ice categories 1380 asum(:,:) = asum(:,:) + a_i(:,:,jl) 1643 1381 END DO 1644 1645 !----------------------------------------------------------------- 1646 ! ice categories 1647 !----------------------------------------------------------------- 1648 1649 DO jl = 1, jpl 1650 DO jj= 1, jpj 1651 DO ji = 1, jpi 1652 asum(ji,jj) = asum(ji,jj) + a_i(ji,jj,jl) 1653 END DO !ji 1654 END DO !jj 1655 END DO ! jl 1656 1382 ! 1657 1383 END SUBROUTINE lim_itd_me_asumr 1658 1384 1659 !============================================================================== 1660 1661 SUBROUTINE lim_itd_me_init ! (subroutine 6/6) 1385 1386 SUBROUTINE lim_itd_me_init 1662 1387 !!------------------------------------------------------------------- 1663 1388 !! *** ROUTINE lim_itd_me_init *** … … 1671 1396 !! 1672 1397 !! ** input : Namelist namiceitdme 1673 !!1674 !! history :1675 !! 9.0, LIM3.0 - 02-2006 (M. Vancoppenolle) original code1676 1398 !!------------------------------------------------------------------- 1677 1399 NAMELIST/namiceitdme/ ridge_scheme_swi, Cs, Cf, fsnowrdg, fsnowrft,& … … 1681 1403 brinstren_swi 1682 1404 !!------------------------------------------------------------------- 1683 1684 ! Define the initial parameters 1685 ! ------------------------- 1686 REWIND( numnam_ice ) 1405 ! 1406 REWIND( numnam_ice ) ! read namiceitdme namelist 1687 1407 READ ( numnam_ice , namiceitdme) 1688 IF (lwp) THEN 1408 ! 1409 IF (lwp) THEN ! control print 1689 1410 WRITE(numout,*) 1690 1411 WRITE(numout,*)' lim_itd_me_init : ice parameters for mechanical ice redistribution ' … … 1707 1428 WRITE(numout,*)' Switch for including brine volume in ice strength comp. brinstren_swi ', brinstren_swi 1708 1429 ENDIF 1709 1430 ! 1710 1431 END SUBROUTINE lim_itd_me_init 1711 1432 1712 !==============================================================================1713 1433 1714 1434 SUBROUTINE lim_itd_me_zapsmall … … 1717 1437 !! 1718 1438 !! ** Purpose : Remove too small sea ice areas and correct salt fluxes 1719 !!1720 1439 !! 1721 1440 !! history : … … 1726 1445 !! 9.0, LIM3.0 - 02-2006 (M. Vancoppenolle) original code 1727 1446 !!------------------------------------------------------------------- 1728 1729 INTEGER :: & 1730 ji,jj, & ! horizontal indices 1731 jl, & ! ice category index 1732 jk, & ! ice layer index 1733 ! ij, & ! combined i/j horizontal index 1734 icells ! number of cells with ice to zap 1735 1736 ! INTEGER, DIMENSION(1:(jpi+1)*(jpj+1)) :: & 1737 ! indxi, & ! compressed indices for i/j directions 1738 ! indxj 1739 1740 INTEGER, DIMENSION(jpi,jpj) :: zmask 1741 1742 1743 REAL(wp) :: & 1744 xtmp ! temporary variable 1447 INTEGER :: ji, jj, jl, jk ! dummy loop indices 1448 INTEGER :: icells ! number of cells with ice to zap 1449 1450 REAL(wp), DIMENSION(jpi,jpj) :: zmask ! 2D workspace 1451 1452 !!gm REAL(wp) :: xtmp ! temporary variable 1453 !!------------------------------------------------------------------- 1745 1454 1746 1455 DO jl = 1, jpl … … 1763 1472 1764 1473 icells = 0 1765 zmask = 0.e01474 zmask = 0._wp 1766 1475 DO jj = 1, jpj 1767 1476 DO ji = 1, jpi 1768 IF ( ( a_i(ji,jj,jl) .GE. -epsi11 .AND. a_i(ji,jj,jl) .LT. 0.0) & 1769 .OR. & 1770 ( a_i(ji,jj,jl) .GT. 0.0 .AND. a_i(ji,jj,jl) .LE. 1.0e-11 ) & 1771 .OR. & 1772 !new line 1773 ( v_i(ji,jj,jl) .EQ. 0.0 .AND. a_i(ji,jj,jl) .GT. 0.0 ) & 1774 .OR. & 1775 ( v_i(ji,jj,jl) .GT. 0.0 .AND. v_i(ji,jj,jl) .LT. 1.e-12 ) ) THEN 1776 zmask(ji,jj) = 1 1777 ENDIF 1778 END DO 1779 END DO 1780 IF( ln_nicep ) WRITE(numout,*) SUM(zmask), ' cells of ice zapped in the ocean ' 1477 IF( ( a_i(ji,jj,jl) .GE. -epsi11 .AND. a_i(ji,jj,jl) .LT. 0._wp ) .OR. & 1478 & ( a_i(ji,jj,jl) .GT. 0._wp .AND. a_i(ji,jj,jl) .LE. 1.0e-11 ) .OR. & 1479 & ( v_i(ji,jj,jl) == 0._wp .AND. a_i(ji,jj,jl) .GT. 0._wp ) .OR. & 1480 & ( v_i(ji,jj,jl) .GT. 0._wp .AND. v_i(ji,jj,jl) .LT. 1.e-12 ) ) zmask(ji,jj) = 1._wp 1481 END DO 1482 END DO 1483 IF( ln_nicep ) WRITE(numout,*) SUM(zmask), ' cells of ice zapped in the ocean ' 1781 1484 1782 1485 !----------------------------------------------------------------- … … 1787 1490 DO jj = 1 , jpj 1788 1491 DO ji = 1 , jpi 1789 1790 xtmp = e_i(ji,jj,jk,jl) / area(ji,jj) / rdt_ice 1791 xtmp = xtmp * unit_fac 1792 ! fheat_res(ji,jj) = fheat_res(ji,jj) - xtmp 1492 !!gm xtmp = e_i(ji,jj,jk,jl) / area(ji,jj) / rdt_ice 1493 !!gm xtmp = xtmp * unit_fac 1494 ! fheat_res(ji,jj) = fheat_res(ji,jj) - xtmp 1793 1495 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * ( 1 - zmask(ji,jj) ) 1794 END DO ! ji1795 END DO ! jj1796 END DO ! jk1496 END DO 1497 END DO 1498 END DO 1797 1499 1798 1500 DO jj = 1 , jpj … … 1802 1504 ! Zap snow energy and use ocean heat to melt snow 1803 1505 !----------------------------------------------------------------- 1804 1805 1506 ! xtmp = esnon(i,j,n) / dt ! < 0 1806 1507 ! fhnet(i,j) = fhnet(i,j) + xtmp … … 1809 1510 ! fluxes are positive to the ocean 1810 1511 ! here the flux has to be negative for the ocean 1811 xtmp = ( rhosn*cpic*( rtt-t_s(ji,jj,1,jl) ) + rhosn*lfus ) / rdt_ice1512 !!gm xtmp = ( rhosn*cpic*( rtt-t_s(ji,jj,1,jl) ) + rhosn*lfus ) / rdt_ice 1812 1513 ! fheat_res(ji,jj) = fheat_res(ji,jj) - xtmp 1813 1514 1814 xtmp = ( rhosn*cpic*( rtt-t_s(ji,jj,1,jl) ) + rhosn*lfus ) / rdt_ice !RB ???????1515 !!gm xtmp = ( rhosn*cpic*( rtt-t_s(ji,jj,1,jl) ) + rhosn*lfus ) / rdt_ice !RB ??????? 1815 1516 1816 1517 t_s(ji,jj,1,jl) = rtt * zmask(ji,jj) + t_s(ji,jj,1,jl) * ( 1 - zmask(ji,jj) ) … … 1833 1534 ! fsalt_hist(i,j) = fsalt_hist(i,j) + xtmp 1834 1535 1835 ato_i(ji,jj) = a_i(ji,jj,jl) * zmask(ji,jj)+ ato_i(ji,jj)1836 a_i (ji,jj,jl) = a_i(ji,jj,jl) * ( 1 - zmask(ji,jj) )1837 v_i (ji,jj,jl) = v_i(ji,jj,jl) * ( 1 - zmask(ji,jj) )1838 v_s (ji,jj,jl) = v_s(ji,jj,jl) * ( 1 - zmask(ji,jj) )1839 t_su (ji,jj,jl) = t_su(ji,jj,jl) * (1 -zmask(ji,jj) ) + t_bo(ji,jj) * zmask(ji,jj)1840 oa_i (ji,jj,jl) = oa_i(ji,jj,jl) * ( 1 - zmask(ji,jj) )1536 ato_i(ji,jj) = a_i (ji,jj,jl) * zmask(ji,jj) + ato_i(ji,jj) 1537 a_i (ji,jj,jl) = a_i (ji,jj,jl) * ( 1 - zmask(ji,jj) ) 1538 v_i (ji,jj,jl) = v_i (ji,jj,jl) * ( 1 - zmask(ji,jj) ) 1539 v_s (ji,jj,jl) = v_s (ji,jj,jl) * ( 1 - zmask(ji,jj) ) 1540 t_su (ji,jj,jl) = t_su (ji,jj,jl) * ( 1 - zmask(ji,jj) ) + t_bo(ji,jj) * zmask(ji,jj) 1541 oa_i (ji,jj,jl) = oa_i (ji,jj,jl) * ( 1 - zmask(ji,jj) ) 1841 1542 smv_i(ji,jj,jl) = smv_i(ji,jj,jl) * ( 1 - zmask(ji,jj) ) 1842 1843 END DO ! ji1844 END DO ! jj1845 1543 ! 1544 END DO 1545 END DO 1546 ! 1846 1547 END DO ! jl 1847 1548 ! 1848 1549 END SUBROUTINE lim_itd_me_zapsmall 1849 1550 1850 1551 #else 1851 !!====================================================================== 1852 !! *** MODULE limitd_me *** 1853 !! no sea ice model 1854 !!====================================================================== 1855 1552 !!---------------------------------------------------------------------- 1553 !! Default option Empty module NO LIM-3 sea-ice model 1554 !!---------------------------------------------------------------------- 1856 1555 CONTAINS 1857 1858 1556 SUBROUTINE lim_itd_me ! Empty routines 1859 1557 END SUBROUTINE lim_itd_me 1860 1558 SUBROUTINE lim_itd_me_icestrength 1861 1559 END SUBROUTINE lim_itd_me_icestrength 1862 SUBROUTINE lim_itd_me_ridgeprep1863 END SUBROUTINE lim_itd_me_ridgeprep1864 SUBROUTINE lim_itd_me_ridgeshift1865 END SUBROUTINE lim_itd_me_ridgeshift1866 SUBROUTINE lim_itd_me_asumr1867 END SUBROUTINE lim_itd_me_asumr1868 1560 SUBROUTINE lim_itd_me_sort 1869 1561 END SUBROUTINE lim_itd_me_sort … … 1873 1565 END SUBROUTINE lim_itd_me_zapsmall 1874 1566 #endif 1567 !!====================================================================== 1875 1568 END MODULE limitd_me
Note: See TracChangeset
for help on using the changeset viewer.
