Changeset 4332 for branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3
- Timestamp:
- 2013-12-11T15:38:42+01:00 (10 years ago)
- Location:
- branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3
- Files:
-
- 1 deleted
- 19 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/ice.F90
r4220 r4332 292 292 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: dh_i_surf2D, dh_i_bott2D, fstbif, fsup2D, focea2D, q_s 293 293 294 INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: patho_case ! number of the pathological case (if any, of course)295 296 294 !!-------------------------------------------------------------------------- 297 295 !! * Ice global state variables … … 396 394 LOGICAL , PUBLIC :: ln_limdyn = .TRUE. !: flag for ice dynamics (T) or not (F) 397 395 LOGICAL , PUBLIC :: ln_nicep = .FALSE. !: flag for sea-ice points output (T) or not (F) 398 REAL(wp) , PUBLIC :: cai = 1.4 0e-3 !: atmospheric drag over sea ice399 REAL(wp) , PUBLIC :: cao = 1.0 0e-3 !: atmospheric drag over ocean396 REAL(wp) , PUBLIC :: cai = 1.4e-3 !: atmospheric drag over sea ice 397 REAL(wp) , PUBLIC :: cao = 1.0e-3 !: atmospheric drag over ocean 400 398 REAL(wp) , PUBLIC :: amax = 0.99 !: maximum ice concentration 401 399 ! … … 465 463 & fsup2D (jpi,jpj) , focea2D (jpi,jpj) , q_s (jpi,jpj) , STAT=ierr(ii) ) 466 464 467 ii = ii + 1468 ALLOCATE( patho_case(jpi, jpj, jpl) , STAT=ierr(ii) )469 470 465 ! * Ice global state variables 471 466 ii = ii + 1 -
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/iceini.F90
r4220 r4332 88 88 ! 89 89 ! ! Initial sea-ice state 90 IF( .NOT. ln_rstart ) THEN ! start from rest90 IF( .NOT. ln_rstart ) THEN ! start from rest 91 91 numit = 0 92 92 numit = nit000 - 1 … … 101 101 ENDIF 102 102 ! 103 CALL lim_sbc_init ! ice surface boundary condition 104 ! 105 fr_i(:,:) = at_i(:,:) ! initialisation of sea-ice fraction 106 tn_ice(:,:,:) = t_su(:,:,:) 103 hi_max(jpl) = 99._wp ! Set hi_max(jpl) to a big value to ensure that all ice is thinner than hi_max(jpl) 104 ! 105 CALL lim_sbc_init ! ice surface boundary condition 106 ! 107 fr_i(:,:) = at_i(:,:) ! initialisation of sea-ice fraction 108 tn_ice(:,:,:) = t_su(:,:,:) ! initialisation of surface temp for coupled simu 107 109 ! 108 110 nstart = numit + nn_fsbc … … 132 134 READ ( numnam_ice , namicerun ) 133 135 ! 134 IF( lk_mpp .AND. ln_nicep ) THEN135 ln_nicep = .FALSE.136 CALL ctl_warn( 'ice_run : specific control print for LIM3 desactivated with MPI' )137 ENDIF136 !IF( lk_mpp .AND. ln_nicep ) THEN 137 ! ln_nicep = .FALSE. 138 ! CALL ctl_warn( 'ice_run : specific control print for LIM3 desactivated with MPI' ) 139 !ENDIF 138 140 ! 139 141 IF(lwp) THEN ! control print … … 159 161 !! ** Purpose : Initializes the ice thickness distribution 160 162 !! ** Method : ... 163 !! Note : hi_max(jpl) is here set up to a value close to 7 m for 164 !! limistate (only) and is changed to 99 m in ice_init 161 165 !!------------------------------------------------------------------ 162 166 INTEGER :: jl, jm ! dummy loop index -
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limcat_1D.F90
r4099 r4332 44 44 REAL(wp), DIMENSION(:), INTENT(in) :: zhti, zhts, zai ! input ice/snow variables 45 45 REAL(wp), DIMENSION(:,:), INTENT(inout) :: zht_i, zht_s, za_i ! output ice/snow variables 46 REAL(wp) :: epsi06 = 1.0e-6 47 46 REAL(wp) :: epsi06 = 1.0e-6 47 REAL(wp) :: zc1, zc2, zc3, zx1, zdh ! local scalars 48 REAL(wp), DIMENSION(0:jpl) :: zhi_max !:Boundary of ice thickness categories in thickness space 49 48 50 IF( nn_timing == 1 ) CALL timing_start('limcat_1D') 49 51 !-------------------------------------------------------------------- … … 51 53 !-------------------------------------------------------------------- 52 54 ijpij = SIZE(zhti,1) 53 zht_i(1:ijpij,1:jpl) = 0. d054 zht_s(1:ijpij,1:jpl) = 0. d055 za_i (1:ijpij,1:jpl) = 0. d055 zht_i(1:ijpij,1:jpl) = 0._wp 56 zht_s(1:ijpij,1:jpl) = 0._wp 57 za_i (1:ijpij,1:jpl) = 0._wp 56 58 57 59 !------------------------------------------------------------------------------------ … … 66 68 ! fulfills ice volume concervation between input and output (ztests=4) 67 69 !-------------------------------------------------------------------------------------- 70 71 !- Thickness categories boundaries 72 ! hi_max is calculated in iceini.F90 but since limcat_1D.F90 routine 73 ! is called before (in bdydta.F90), one must recalculate it. 74 ! Note clem: there may be a way of doing things cleaner 75 !---------------------------------- 76 zhi_max(:) = 0._wp 77 zc1 = 3._wp / REAL( jpl , wp ) ; zc2 = 10._wp * zc1 ; zc3 = 3._wp 78 DO jl = 1, jpl 79 zx1 = REAL( jl-1 , wp ) / REAL( jpl , wp ) 80 zhi_max(jl) = zhi_max(jl-1) + zc1 + zc2 * ( 1._wp + TANH( zc3 * ( zx1 - 1._wp ) ) ) 81 END DO 68 82 69 83 ! ---------------------------------------- … … 71 85 ! ---------------------------------------- 72 86 DO ji = 1, ijpij 73 ! snow thickness in each category74 zht_s(ji,1:jpl) = zhts(ji)75 87 88 IF( zhti(ji) > 0._wp ) THEN 89 76 90 ! initialisation of tests 77 91 ztest_1 = 0 … … 87 101 88 102 ! initialisation of ice variables for each try 89 zht_i(ji,1:jpl) = 0. d090 za_i (ji,1:jpl) = 0. d0103 zht_i(ji,1:jpl) = 0._wp 104 za_i (ji,1:jpl) = 0._wp 91 105 92 106 ! *** case very thin ice: fill only category 1 … … 94 108 zht_i(ji,1) = zhti(ji) 95 109 za_i (ji,1) = zai (ji) 96 97 110 ! *** case ice is thicker: fill categories >1 98 111 ELSE … … 100 113 ! Fill ice thicknesses except the last one (i_fill) by (hmax-hmin)/2 101 114 DO jl = 1, i_fill - 1 102 zht_i(ji,jl) = ( hi_max(jl) + hi_max(jl-1) ) / 2.115 zht_i(ji,jl) = ( zhi_max(jl) + zhi_max(jl-1) ) * 0.5_wp 103 116 END DO 104 117 … … 106 119 jl0 = i_fill 107 120 DO jl = 1, i_fill 108 IF ( ( zhti(ji) >= hi_max(jl-1) ) .AND. ( zhti(ji) <hi_max(jl) ) ) THEN121 IF ( ( zhti(ji) >= zhi_max(jl-1) ) .AND. ( zhti(ji) < zhi_max(jl) ) ) THEN 109 122 jl0 = jl 110 123 CYCLE … … 116 129 DO jl = 1, i_fill - 1 117 130 IF ( jl == jl0 ) CYCLE 118 zarg = ( zht_i(ji,jl) - zhti(ji) ) / ( zhti(ji) / 2.)131 zarg = ( zht_i(ji,jl) - zhti(ji) ) / ( zhti(ji) * 0.5_wp ) 119 132 za_i(ji,jl) = za_i (ji,jl0) * EXP(-zarg**2) 120 133 END DO … … 141 154 142 155 ! Test 3: thickness of the last category is in-bounds ? 143 IF ( zht_i(ji,i_fill) >= hi_max(i_fill-1) ) ztest_3 = 1156 IF ( zht_i(ji,i_fill) >= zhi_max(i_fill-1) ) ztest_3 = 1 144 157 145 158 ! Test 4: positivity of ice concentrations 146 159 ztest_4 = 1 147 160 DO jl = 1, i_fill 148 IF ( za_i(ji,jl) < 0. 0d0) ztest_4 = 0161 IF ( za_i(ji,jl) < 0._wp ) ztest_4 = 0 149 162 END DO 150 163 … … 164 177 ! WRITE(numout,*) 'za_i(ji,jpl)=',za_i(ji,:) 165 178 !ENDIF 166 179 180 ENDIF ! if zhti > 0 181 167 182 END DO ! i loop 183 184 ! ------------------------------------------------ 185 ! Adding Snow in each category where za_i is not 0 186 ! ------------------------------------------------ 187 DO jl = 1, jpl 188 DO ji = 1, ijpij 189 IF( za_i(ji,jl) > 0._wp ) THEN 190 zht_s(ji,jl) = zht_i(ji,jl) * ( zhts(ji) / zhti(ji) ) 191 ! In case snow load is in excess that would lead to transformation from snow to ice 192 ! Then, transfer the snow excess into the ice (different from limthd_dh) 193 zdh = MAX( 0._wp, ( rhosn * zht_s(ji,jl) + ( rhoic - rau0 ) * zht_i(ji,jl) ) * r1_rau0 ) 194 ! recompute ht_i, ht_s avoiding out of bounds values 195 zht_i(ji,jl) = MIN( zhi_max(jl), zht_i(ji,jl) + zdh ) 196 zht_s(ji,jl) = MAX( 0._wp, zht_s(ji,jl) - zdh * rhoic / rhosn ) 197 ENDIF 198 ENDDO 199 ENDDO 168 200 169 201 IF( nn_timing == 1 ) CALL timing_stop('limcat_1D') -
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limhdf.F90
r4045 r4332 29 29 30 30 LOGICAL :: linit = .TRUE. ! initialization flag (set to flase after the 1st call) 31 REAL(wp) :: epsi04 = 1 e-04 ! constant31 REAL(wp) :: epsi04 = 1.e-04 ! constant 32 32 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: efact ! metric coefficient 33 33 … … 54 54 REAL(wp), DIMENSION(jpi,jpj), INTENT( inout ) :: ptab ! Field on which the diffusion is applied 55 55 ! 56 INTEGER :: ji, jj ! dummy loop indices57 INTEGER :: its, iter, ierr ! local integers58 REAL(wp) :: zalfa, zrlxint, zconv , zeps! local scalars56 INTEGER :: ji, jj ! dummy loop indices 57 INTEGER :: its, iter, ierr ! local integers 58 REAL(wp) :: zalfa, zrlxint, zconv ! local scalars 59 59 REAL(wp), POINTER, DIMENSION(:,:) :: zrlx, zflu, zflv, zdiv0, zdiv, ztab0 60 60 CHARACTER(lc) :: charout ! local character … … 79 79 zalfa = 0.5_wp ! =1.0/0.5/0.0 = implicit/Cranck-Nicholson/explicit 80 80 its = 100 ! Maximum number of iteration 81 zeps = 2._wp * epsi0482 81 ! 83 82 ztab0(:, : ) = ptab(:,:) ! Arrays initialization … … 94 93 iter = 0 95 94 ! 96 DO WHILE( zconv > zeps.AND. iter <= its ) ! Sub-time step loop95 DO WHILE( zconv > ( 2._wp * epsi04 ) .AND. iter <= its ) ! Sub-time step loop 97 96 ! 98 97 iter = iter + 1 ! incrementation of the sub-time step number -
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90
r4099 r4332 76 76 !! 6) Lateral boundary conditions 77 77 !! 78 !! ** Notes : o_i, t_su, t_s, t_i, s_i must be filled everywhere, even 79 !! where there is no ice (clem: I do not know why but it is mandatory) 80 !! 78 81 !! History : 79 82 !! 2.0 ! 01-04 (C. Ethe, G. Madec) Original code … … 84 87 !! * Local variables 85 88 INTEGER :: ji, jj, jk, jl ! dummy loop indices 86 REAL(wp) :: epsi 06, epsi20, ztmelts89 REAL(wp) :: epsi20, ztmelts, zdh 87 90 INTEGER :: i_hemis, i_fill, jl0 88 91 REAL(wp) :: ztest_1, ztest_2, ztest_3, ztest_4, ztests, zsigma, zarg, zA, zV, zA_cons, zV_cons, zconv … … 98 101 CALL wrk_alloc( 2, zhm_i_ini, zat_i_ini, zvt_i_ini, zhm_s_ini, zsm_i_ini ) 99 102 100 epsi06 = 1.0e-6101 103 epsi20 = 1.0e-20 102 104 IF(lwp) WRITE(numout,*) … … 308 310 !--------------------------------------------------------------------- 309 311 310 ! Ice concentration, thickness and volume, snow depth, ice 311 ! salinity, ice age, surface temperature 312 ! Ice concentration, thickness and volume, ice salinity, ice age, surface temperature 312 313 DO jl = 1, jpl ! loop over categories 313 314 DO jj = 1, jpj … … 315 316 a_i(ji,jj,jl) = zidto(ji,jj) * za_i_ini (jl,zhemis(ji,jj)) ! concentration 316 317 ht_i(ji,jj,jl) = zidto(ji,jj) * zht_i_ini(jl,zhemis(ji,jj)) ! ice thickness 317 ht_s(ji,jj,jl) = zidto(ji,jj) * zhm_s_ini(zhemis(ji,jj))! snow depth318 sm_i(ji,jj,jl) = zidto(ji,jj) * zsm_i_ini(zhemis(ji,jj)) 318 ht_s(ji,jj,jl) = ht_i(ji,jj,jl) * ( zhm_s_ini( zhemis(ji,jj) ) / zhm_i_ini( zhemis(ji,jj) ) ) ! snow depth 319 sm_i(ji,jj,jl) = zidto(ji,jj) * zsm_i_ini(zhemis(ji,jj)) + ( 1._wp - zidto(ji,jj) ) * s_i_min ! salinity 319 320 o_i(ji,jj,jl) = zidto(ji,jj) * 1._wp + ( 1._wp - zidto(ji,jj) ) ! age 320 t_su(ji,jj,jl) = zidto(ji,jj) * 270.0 + ( 1._wp - zidto(ji,jj) ) * t_bo(ji,jj) ! surf temp 321 322 ! ice volume, snow volume, salt content, age content 321 t_su(ji,jj,jl) = zidto(ji,jj) * 270.0 + ( 1._wp - zidto(ji,jj) ) * 270.0 ! surf temp 322 323 ! This case below should not be used if (ht_s/ht_i) is ok in namelist 324 ! In case snow load is in excess that would lead to transformation from snow to ice 325 ! Then, transfer the snow excess into the ice (different from limthd_dh) 326 zdh = MAX( 0._wp, ( rhosn * ht_s(ji,jj,jl) + ( rhoic - rau0 ) * ht_i(ji,jj,jl) ) * r1_rau0 ) 327 ! recompute ht_i, ht_s avoiding out of bounds values 328 ht_i(ji,jj,jl) = MIN( hi_max(jl), ht_i(ji,jj,jl) + zdh ) 329 ht_s(ji,jj,jl) = MAX( 0._wp, ht_s(ji,jj,jl) - zdh * rhoic / rhosn ) 330 331 ! ice volume, salt content, age content 323 332 v_i(ji,jj,jl) = ht_i(ji,jj,jl) * a_i(ji,jj,jl) ! ice volume 324 333 v_s(ji,jj,jl) = ht_s(ji,jj,jl) * a_i(ji,jj,jl) ! snow volume -
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90
r4072 r4332 43 43 PUBLIC lim_itd_me_alloc ! called by iceini.F90 44 44 45 REAL(wp) :: epsi 11 = 1.e-11_wp ! constant values45 REAL(wp) :: epsi20 = 1.e-20_wp ! constant values 46 46 REAL(wp) :: epsi10 = 1.e-10_wp ! constant values 47 47 REAL(wp) :: epsi06 = 1.e-06_wp ! constant values … … 50 50 ! Variables shared among ridging subroutines 51 51 !----------------------------------------------------------------------- 52 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: asum ! sum of total ice and open water area 53 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: aksum ! ratio of area removed to area ridged 54 ! 55 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: athorn ! participation function; fraction of ridging/ 56 ! ! closing associated w/ category n 57 ! 58 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hrmin ! minimum ridge thickness 59 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hrmax ! maximum ridge thickness 60 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hraft ! thickness of rafted ice 61 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: krdg ! mean ridge thickness/thickness of ridging ice 62 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: aridge ! participating ice ridging 63 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: araft ! participating ice rafting 52 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: asum ! sum of total ice and open water area 53 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: aksum ! ratio of area removed to area ridged 54 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: athorn ! participation function; fraction of ridging/ 55 ! ! closing associated w/ category n 56 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: hrmin ! minimum ridge thickness 57 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: hrmax ! maximum ridge thickness 58 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: hraft ! thickness of rafted ice 59 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: krdg ! mean ridge thickness/thickness of ridging ice 60 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: aridge ! participating ice ridging 61 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: araft ! participating ice rafting 64 62 65 63 REAL(wp), PARAMETER :: krdgmin = 1.1_wp ! min ridge thickness multiplier 66 64 REAL(wp), PARAMETER :: kraft = 2.0_wp ! rafting multipliyer 67 REAL(wp), PARAMETER :: kamax = 1.0 65 REAL(wp), PARAMETER :: kamax = 1.0_wp ! max of ice area authorized (clem: scheme is not stable if kamax <= 0.99) 68 66 69 67 REAL(wp) :: Cp ! … … 183 181 ! 1) Thickness categories boundaries, ice / o.w. concentrations, init_ons 184 182 !-----------------------------------------------------------------------------! 185 ! Set hi_max(ncat) to a big value to ensure that all ridged ice is thinner than hi_max(ncat).186 187 hi_max(jpl) = 999.99188 189 183 Cp = 0.5 * grav * (rau0-rhoic) * rhoic / rau0 ! proport const for PE 190 184 ! … … 265 259 ! Reduce the closing rate if more than 100% of the open water 266 260 ! would be removed. Reduce the opening rate proportionately. 267 IF ( ato_i(ji,jj) .GT. epsi1 1.AND. athorn(ji,jj,0) .GT. 0.0 ) THEN261 IF ( ato_i(ji,jj) .GT. epsi10 .AND. athorn(ji,jj,0) .GT. 0.0 ) THEN 268 262 w1 = athorn(ji,jj,0) * closing_gross(ji,jj) * rdt_ice 269 263 IF ( w1 .GT. ato_i(ji,jj)) THEN … … 285 279 DO jj = 1, jpj 286 280 DO ji = 1, jpi 287 IF ( a_i(ji,jj,jl) > epsi1 1.AND. athorn(ji,jj,jl) > 0._wp )THEN281 IF ( a_i(ji,jj,jl) > epsi10 .AND. athorn(ji,jj,jl) > 0._wp )THEN 288 282 w1 = athorn(ji,jj,jl) * closing_gross(ji,jj) * rdt_ice 289 283 IF ( w1 > a_i(ji,jj,jl) ) THEN … … 316 310 DO jj = 1, jpj 317 311 DO ji = 1, jpi 318 IF (ABS(asum(ji,jj) - kamax ) .LT. epsi1 1) THEN312 IF (ABS(asum(ji,jj) - kamax ) .LT. epsi10) THEN 319 313 closing_net(ji,jj) = 0._wp 320 314 opning (ji,jj) = 0._wp … … 358 352 DO ji = 1, jpi 359 353 360 IF(ABS(asum(ji,jj) - kamax) > epsi1 1) asum_error = .true.354 IF(ABS(asum(ji,jj) - kamax) > epsi10 ) asum_error = .true. 361 355 362 356 dardg1dt(ji,jj) = dardg1dt(ji,jj) * r1_rdtice … … 377 371 DO jj = 1, jpj 378 372 DO ji = 1, jpi 379 IF( ABS( asum(ji,jj) - kamax) > epsi1 1) THEN ! there is a bug373 IF( ABS( asum(ji,jj) - kamax) > epsi10 ) THEN ! there is a bug 380 374 WRITE(numout,*) ' ' 381 375 WRITE(numout,*) ' ALERTE : Ridging error: total area = ', asum(ji,jj) … … 400 394 401 395 !-----------------------------------------------------------------------------! 402 ! 6) Updating state variables and trend terms 396 ! 6) Updating state variables and trend terms (done in limupdate) 403 397 !-----------------------------------------------------------------------------! 404 405 398 CALL lim_var_glo2eqv 406 399 CALL lim_itd_me_zapsmall … … 419 412 END DO 420 413 END DO 421 422 !-----------------423 ! Trend terms424 !-----------------425 426 d_u_ice_dyn(:,:) = u_ice(:,:) - old_u_ice(:,:)427 d_v_ice_dyn(:,:) = v_ice(:,:) - old_v_ice(:,:)428 d_a_i_trp (:,:,:) = a_i (:,:,:) - old_a_i (:,:,:)429 d_v_s_trp (:,:,:) = v_s (:,:,:) - old_v_s (:,:,:)430 d_v_i_trp (:,:,:) = v_i (:,:,:) - old_v_i (:,:,:)431 d_e_s_trp (:,:,:,:) = e_s (:,:,:,:) - old_e_s (:,:,:,:)432 d_e_i_trp (:,:,:,:) = e_i (:,:,:,:) - old_e_i (:,:,:,:)433 d_oa_i_trp (:,:,:) = oa_i (:,:,:) - old_oa_i (:,:,:)434 d_smv_i_trp(:,:,:) = 0._wp435 IF( num_sal == 2 ) d_smv_i_trp(:,:,:) = smv_i(:,:,:) - old_smv_i(:,:,:)436 414 437 415 IF(ln_ctl) THEN ! Control print … … 491 469 ! ------------------------------- 492 470 493 !-------------------------!494 ! Back to initial values495 !-------------------------!496 497 ! update of fields will be made later in lim update498 u_ice(:,:) = old_u_ice(:,:)499 v_ice(:,:) = old_v_ice(:,:)500 a_i(:,:,:) = old_a_i(:,:,:)501 v_s(:,:,:) = old_v_s(:,:,:)502 v_i(:,:,:) = old_v_i(:,:,:)503 e_s(:,:,:,:) = old_e_s(:,:,:,:)504 e_i(:,:,:,:) = old_e_i(:,:,:,:)505 oa_i(:,:,:) = old_oa_i(:,:,:)506 IF( num_sal == 2 ) smv_i(:,:,:) = old_smv_i(:,:,:)507 508 !----------------------------------------------------!509 ! Advection of ice in a free cell, newly ridged ice510 !----------------------------------------------------!511 512 ! to allow for thermodynamics to melt new ice513 ! we immediately advect ice in free cells514 515 ! heat content has to be corrected before ice volume516 !clem@order517 ! DO jl = 1, jpl518 ! DO jk = 1, nlay_i519 ! DO jj = 1, jpj520 ! DO ji = 1, jpi521 ! IF ( ( old_v_i(ji,jj,jl) < epsi06 ) .AND. &522 ! ( d_v_i_trp(ji,jj,jl) > epsi06 ) ) THEN523 ! old_e_i(ji,jj,jk,jl) = d_e_i_trp(ji,jj,jk,jl)524 ! d_e_i_trp(ji,jj,jk,jl) = 0._wp525 ! ENDIF526 ! END DO527 ! END DO528 ! END DO529 ! END DO530 !531 ! DO jl = 1, jpl532 ! DO jj = 1, jpj533 ! DO ji = 1, jpi534 ! IF( old_v_i (ji,jj,jl) < epsi06 .AND. &535 ! d_v_i_trp(ji,jj,jl) > epsi06 ) THEN536 ! old_v_i (ji,jj,jl) = d_v_i_trp(ji,jj,jl)537 ! d_v_i_trp (ji,jj,jl) = 0._wp538 ! old_a_i (ji,jj,jl) = d_a_i_trp(ji,jj,jl)539 ! d_a_i_trp (ji,jj,jl) = 0._wp540 ! old_v_s (ji,jj,jl) = d_v_s_trp(ji,jj,jl)541 ! d_v_s_trp (ji,jj,jl) = 0._wp542 ! old_e_s (ji,jj,1,jl) = d_e_s_trp(ji,jj,1,jl)543 ! d_e_s_trp (ji,jj,1,jl) = 0._wp544 ! old_oa_i (ji,jj,jl) = d_oa_i_trp(ji,jj,jl)545 ! d_oa_i_trp(ji,jj,jl) = 0._wp546 ! IF( num_sal == 2 ) old_smv_i(ji,jj,jl) = d_smv_i_trp(ji,jj,jl)547 ! d_smv_i_trp(ji,jj,jl) = 0._wp548 ! ENDIF549 ! END DO550 ! END DO551 ! END DO552 !clem@order553 471 ENDIF ! ln_limdyn=.true. 554 472 ! … … 605 523 DO ji = 1, jpi 606 524 ! 607 IF( a_i(ji,jj,jl) > epsi11 .AND. & 608 athorn(ji,jj,jl) > 0._wp ) THEN 525 IF( a_i(ji,jj,jl) > epsi10 .AND. athorn(ji,jj,jl) > 0._wp ) THEN 609 526 hi = v_i(ji,jj,jl) / a_i(ji,jj,jl) 610 527 !---------------------------- … … 624 541 * z1_3 * (hrmax(ji,jj,jl)**3 - hrmin(ji,jj,jl)**3) / ( hrmax(ji,jj,jl)-hrmin(ji,jj,jl) ) 625 542 !!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... 626 ENDIF ! aicen > epsi1 1543 ENDIF ! aicen > epsi10 627 544 ! 628 545 END DO ! ji … … 681 598 DO jj = 2, jpj - 1 682 599 DO ji = 2, jpi - 1 683 IF ( ( asum(ji,jj) - ato_i(ji,jj) ) .GT. epsi1 1) THEN ! ice is600 IF ( ( asum(ji,jj) - ato_i(ji,jj) ) .GT. epsi10) THEN ! ice is 684 601 ! present 685 602 zworka(ji,jj) = 4.0 * strength(ji,jj) & … … 721 638 DO jj = 1, jpj - 1 722 639 DO ji = 1, jpi - 1 723 IF ( ( asum(ji,jj) - ato_i(ji,jj) ) .GT. epsi1 1) THEN ! ice is present640 IF ( ( asum(ji,jj) - ato_i(ji,jj) ) .GT. epsi10) THEN ! ice is present 724 641 numts_rm = 1 ! number of time steps for the running mean 725 642 IF ( strp1(ji,jj) .GT. 0.0 ) numts_rm = numts_rm + 1 … … 794 711 DO jj = 1, jpj 795 712 DO ji = 1, jpi 796 IF( ato_i(ji,jj) > epsi1 1) THEN ; Gsum(ji,jj,0) = ato_i(ji,jj)713 IF( ato_i(ji,jj) > epsi10 ) THEN ; Gsum(ji,jj,0) = ato_i(ji,jj) 797 714 ELSE ; Gsum(ji,jj,0) = 0._wp 798 715 ENDIF … … 804 721 DO jj = 1, jpj 805 722 DO ji = 1, jpi 806 IF( a_i(ji,jj,jl) .GT. epsi1 1) THEN ; Gsum(ji,jj,jl) = Gsum(ji,jj,jl-1) + a_i(ji,jj,jl)723 IF( a_i(ji,jj,jl) .GT. epsi10 ) THEN ; Gsum(ji,jj,jl) = Gsum(ji,jj,jl-1) + a_i(ji,jj,jl) 807 724 ELSE ; Gsum(ji,jj,jl) = Gsum(ji,jj,jl-1) 808 725 ENDIF … … 887 804 IF ( raftswi == 1 ) THEN 888 805 889 IF( MAXVAL(aridge + araft - athorn(:,:,1:jpl)) .GT. epsi1 1) THEN806 IF( MAXVAL(aridge + araft - athorn(:,:,1:jpl)) .GT. epsi10 ) THEN 890 807 DO jl = 1, jpl 891 808 DO jj = 1, jpj 892 809 DO ji = 1, jpi 893 IF ( aridge(ji,jj,jl) + araft(ji,jj,jl) - athorn(ji,jj,jl) .GT. & 894 epsi11 ) THEN 810 IF ( aridge(ji,jj,jl) + araft(ji,jj,jl) - athorn(ji,jj,jl) .GT. epsi10 ) THEN 895 811 WRITE(numout,*) ' ALERTE 96 : wrong participation function ... ' 896 812 WRITE(numout,*) ' ji, jj, jl : ', ji, jj, jl … … 938 854 DO ji = 1, jpi 939 855 940 IF (a_i(ji,jj,jl) .GT. epsi1 1.AND. athorn(ji,jj,jl) .GT. 0.0 ) THEN856 IF (a_i(ji,jj,jl) .GT. epsi10 .AND. athorn(ji,jj,jl) .GT. 0.0 ) THEN 941 857 hi = v_i(ji,jj,jl) / a_i(ji,jj,jl) 942 858 hrmean = MAX(SQRT(Hstar*hi), hi*krdgmin) … … 992 908 INTEGER :: ij ! horizontal index, combines i and j loops 993 909 INTEGER :: icells ! number of cells with aicen > puny 994 REAL(wp) :: z eps, zindb, zsrdg2 ! local scalar910 REAL(wp) :: zindb, zsrdg2 ! local scalar 995 911 REAL(wp) :: hL, hR, farea, zdummy, zdummy0, ztmelts ! left and right limits of integration 996 912 … … 1044 960 1045 961 IF( con_i ) THEN 1046 CALL lim_column_sum (jpl, v_i, vice_init ) 1047 WRITE(numout,*) ' vice_init : ', vice_init(jiindx,jjindx) 962 CALL lim_column_sum (jpl, v_i, vice_init ) 1048 963 CALL lim_column_sum_energy (jpl, nlay_i, e_i, eice_init ) 1049 WRITE(numout,*) ' eice_init : ', eice_init(jiindx,jjindx) 964 DO ji = mi0(jiindx), mi1(jiindx) 965 DO jj = mj0(jjindx), mj1(jjindx) 966 WRITE(numout,*) ' vice_init : ', vice_init(ji,jj) 967 WRITE(numout,*) ' eice_init : ', eice_init(ji,jj) 968 END DO 969 END DO 1050 970 ENDIF 1051 1052 zeps = 1.e-20_wp1053 971 1054 972 !------------------------------------------------------------------------------- … … 1062 980 ato_i(ji,jj) = ato_i(ji,jj) - athorn(ji,jj,0) * closing_gross(ji,jj) * rdt_ice & 1063 981 & + opning(ji,jj) * rdt_ice 1064 IF( ato_i(ji,jj) < -epsi1 1) THEN982 IF( ato_i(ji,jj) < -epsi10 ) THEN 1065 983 neg_ato_i = .TRUE. 1066 984 ELSEIF( ato_i(ji,jj) < 0._wp ) THEN ! roundoff error … … 1074 992 DO jj = 1, jpj 1075 993 DO ji = 1, jpi 1076 IF( ato_i(ji,jj) < -epsi1 1) THEN994 IF( ato_i(ji,jj) < -epsi10 ) THEN 1077 995 WRITE(numout,*) '' 1078 996 WRITE(numout,*) 'Ridging error: ato_i < 0' 1079 997 WRITE(numout,*) 'ato_i : ', ato_i(ji,jj) 1080 ENDIF ! ato_i < -epsi1 1998 ENDIF ! ato_i < -epsi10 1081 999 END DO 1082 1000 END DO … … 1120 1038 DO jj = 1, jpj 1121 1039 DO ji = 1, jpi 1122 IF( aicen_init(ji,jj,jl1) > epsi11 .AND. athorn(ji,jj,jl1) > 0._wp&1123 .AND. closing_gross(ji,jj) > 0._wp ) THEN1040 IF( aicen_init(ji,jj,jl1) > epsi10 .AND. athorn(ji,jj,jl1) > 0._wp & 1041 & .AND. closing_gross(ji,jj) > 0._wp ) THEN 1124 1042 icells = icells + 1 1125 1043 indxi(icells) = ji … … 1158 1076 afrft(ji,jj) = arft1(ji,jj) / aicen_init(ji,jj,jl1) !rafting 1159 1077 1160 IF (afrac(ji,jj) > kamax + epsi1 1) THEN !riging1078 IF (afrac(ji,jj) > kamax + epsi10) THEN !riging 1161 1079 large_afrac = .true. 1162 1080 ELSEIF (afrac(ji,jj) > kamax) THEN ! roundoff error 1163 1081 afrac(ji,jj) = kamax 1164 1082 ENDIF 1165 IF (afrft(ji,jj) > kamax + epsi1 1) THEN !rafting1083 IF (afrft(ji,jj) > kamax + epsi10) THEN !rafting 1166 1084 large_afrft = .true. 1167 1085 ELSEIF (afrft(ji,jj) > kamax) THEN ! roundoff error … … 1222 1140 dardg1dt (ji,jj) = dardg1dt(ji,jj) + ardg1(ji,jj) + arft1(ji,jj) 1223 1141 dardg2dt (ji,jj) = dardg2dt(ji,jj) + ardg2(ji,jj) + arft2(ji,jj) 1224 !clem diag_dyn_gr(ji,jj) = diag_dyn_gr(ji,jj) + ( vrdg2(ji,jj) + virft(ji,jj) ) * r1_rdtice1225 1142 opening (ji,jj) = opening (ji,jj) + opning(ji,jj) * rdt_ice 1226 1143 … … 1276 1193 1277 1194 ! corrected sea water salinity 1278 zindb = MAX( 0._wp , SIGN( 1._wp , vsw(ji,jj) - zeps) )1279 zdummy = zindb * ( srdg1(ji,jj) - srdg2(ji,jj) ) / MAX( ridge_por * vsw(ji,jj), zeps)1195 zindb = MAX( 0._wp , SIGN( 1._wp , vsw(ji,jj) - epsi20 ) ) 1196 zdummy = zindb * ( srdg1(ji,jj) - srdg2(ji,jj) ) / MAX( ridge_por * vsw(ji,jj), epsi20 ) 1280 1197 1281 1198 ztmelts = - tmut * zdummy + rtt … … 1312 1229 ji = indxi(ij) 1313 1230 jj = indxj(ij) 1314 IF( afrac(ji,jj) > kamax + epsi1 1) THEN1231 IF( afrac(ji,jj) > kamax + epsi10 ) THEN 1315 1232 WRITE(numout,*) '' 1316 1233 WRITE(numout,*) ' ardg > a_i' … … 1324 1241 ji = indxi(ij) 1325 1242 jj = indxj(ij) 1326 IF( afrft(ji,jj) > kamax + epsi1 1) THEN1243 IF( afrft(ji,jj) > kamax + epsi10 ) THEN 1327 1244 WRITE(numout,*) '' 1328 1245 WRITE(numout,*) ' arft > a_i' … … 1424 1341 fieldid = ' v_i : limitd_me ' 1425 1342 CALL lim_cons_check (vice_init, vice_final, 1.0e-6, fieldid) 1426 WRITE(numout,*) ' vice_init : ', vice_init(jiindx,jjindx)1427 WRITE(numout,*) ' vice_final : ', vice_final(jiindx,jjindx)1428 1343 1429 1344 CALL lim_column_sum_energy (jpl, nlay_i, e_i, eice_final ) 1430 1345 fieldid = ' e_i : limitd_me ' 1431 1346 CALL lim_cons_check (eice_init, eice_final, 1.0e-2, fieldid) 1432 WRITE(numout,*) ' eice_init : ', eice_init(jiindx,jjindx) 1433 WRITE(numout,*) ' eice_final : ', eice_final(jiindx,jjindx) 1347 1348 DO ji = mi0(jiindx), mi1(jiindx) 1349 DO jj = mj0(jjindx), mj1(jjindx) 1350 WRITE(numout,*) ' vice_init : ', vice_init (ji,jj) 1351 WRITE(numout,*) ' vice_final : ', vice_final(ji,jj) 1352 WRITE(numout,*) ' eice_init : ', eice_init (ji,jj) 1353 WRITE(numout,*) ' eice_final : ', eice_final(ji,jj) 1354 END DO 1355 END DO 1434 1356 ENDIF 1435 1357 ! … … 1533 1455 1534 1456 REAL(wp), POINTER, DIMENSION(:,:) :: zmask ! 2D workspace 1535 1457 REAL(wp) :: zmask_glo 1536 1458 !!gm REAL(wp) :: xtmp ! temporary variable 1537 1459 !!------------------------------------------------------------------- … … 1545 1467 ! Abort model in case of negative area. 1546 1468 !----------------------------------------------------------------- 1547 IF( MINVAL(a_i(:,:,jl)) .LT. -epsi11 .AND. ln_nicep ) THEN1548 DO jj = 1, jpj1549 DO ji = 1, jpi1550 IF ( a_i(ji,jj,jl) .LT. -epsi11 ) THEN1551 WRITE (numout,*) ' ALERTE 98 '1552 WRITE (numout,*) ' Negative ice area: ji, jj, jl: ', ji, jj,jl1553 WRITE (numout,*) ' a_i ', a_i(ji,jj,jl)1554 ENDIF1555 END DO1556 END DO1557 ENDIF1558 1559 1469 icells = 0 1560 zmask = 0._wp1470 zmask(:,:) = 0._wp 1561 1471 DO jj = 1, jpj 1562 1472 DO ji = 1, jpi 1563 IF( ( a_i(ji,jj,jl) .GE. -epsi11 .AND. a_i(ji,jj,jl) .LT. 0._wp ) .OR. & 1564 & ( a_i(ji,jj,jl) .GT. 0._wp .AND. a_i(ji,jj,jl) .LE. 1.0e-11 ) .OR. & 1565 & ( v_i(ji,jj,jl) == 0._wp .AND. a_i(ji,jj,jl) .GT. 0._wp ) .OR. & 1566 & ( v_i(ji,jj,jl) .GT. 0._wp .AND. v_i(ji,jj,jl) .LT. 1.e-12 ) ) zmask(ji,jj) = 1._wp 1567 END DO 1568 END DO 1569 IF( ln_nicep ) WRITE(numout,*) SUM(zmask), ' cells of ice zapped in the ocean ' 1473 IF( ( a_i(ji,jj,jl) >= -epsi10 .AND. a_i(ji,jj,jl) < 0._wp ) .OR. & 1474 & ( a_i(ji,jj,jl) > 0._wp .AND. a_i(ji,jj,jl) <= epsi10 ) .OR. & 1475 & ( v_i(ji,jj,jl) == 0._wp .AND. a_i(ji,jj,jl) > 0._wp ) .OR. & 1476 & ( v_i(ji,jj,jl) > 0._wp .AND. v_i(ji,jj,jl) <= epsi10 ) ) zmask(ji,jj) = 1._wp 1477 END DO 1478 END DO 1479 zmask_glo = glob_sum(zmask) 1480 !IF( ln_nicep .AND. lwp ) WRITE(numout,*) zmask_glo, ' cells of ice zapped in the ocean ' 1570 1481 1571 1482 !----------------------------------------------------------------- … … 1579 1490 !!gm xtmp = xtmp * unit_fac 1580 1491 ! fheat_res(ji,jj) = fheat_res(ji,jj) - xtmp 1581 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * ( 1 - zmask(ji,jj) )1492 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * ( 1._wp - zmask(ji,jj) ) 1582 1493 END DO 1583 1494 END DO … … 1601 1512 !!gm xtmp = ( rhosn*cpic*( rtt-t_s(ji,jj,1,jl) ) + rhosn*lfus ) * r1_rdtice !RB ??????? 1602 1513 1603 t_s(ji,jj,1,jl) = rtt * zmask(ji,jj) + t_s(ji,jj,1,jl) * ( 1 - zmask(ji,jj) )1514 t_s(ji,jj,1,jl) = rtt * zmask(ji,jj) + t_s(ji,jj,1,jl) * ( 1._wp - zmask(ji,jj) ) 1604 1515 1605 1516 !----------------------------------------------------------------- … … 1615 1526 1616 1527 ato_i(ji,jj) = a_i (ji,jj,jl) * zmask(ji,jj) + ato_i(ji,jj) 1617 a_i (ji,jj,jl) = a_i (ji,jj,jl) * ( 1 - zmask(ji,jj) )1618 v_i (ji,jj,jl) = v_i (ji,jj,jl) * ( 1 - zmask(ji,jj) )1619 v_s (ji,jj,jl) = v_s (ji,jj,jl) * ( 1 - zmask(ji,jj) )1620 t_su (ji,jj,jl) = t_su (ji,jj,jl) * ( 1 - zmask(ji,jj) ) + t_bo(ji,jj) * zmask(ji,jj)1621 oa_i (ji,jj,jl) = oa_i (ji,jj,jl) * ( 1 - zmask(ji,jj) )1622 smv_i(ji,jj,jl) = smv_i(ji,jj,jl) * ( 1 - zmask(ji,jj) )1528 a_i (ji,jj,jl) = a_i (ji,jj,jl) * ( 1._wp - zmask(ji,jj) ) 1529 v_i (ji,jj,jl) = v_i (ji,jj,jl) * ( 1._wp - zmask(ji,jj) ) 1530 v_s (ji,jj,jl) = v_s (ji,jj,jl) * ( 1._wp - zmask(ji,jj) ) 1531 t_su (ji,jj,jl) = t_su (ji,jj,jl) * ( 1._wp - zmask(ji,jj) ) + t_bo(ji,jj) * zmask(ji,jj) 1532 oa_i (ji,jj,jl) = oa_i (ji,jj,jl) * ( 1._wp - zmask(ji,jj) ) 1533 smv_i(ji,jj,jl) = smv_i(ji,jj,jl) * ( 1._wp - zmask(ji,jj) ) 1623 1534 ! 1624 1535 END DO -
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90
r4155 r4332 45 45 PUBLIC lim_itd_shiftice 46 46 47 REAL(wp) :: epsi20 = 1.e-20_wp ! constant values48 REAL(wp) :: epsi13 = 1.e-13_wp !49 47 REAL(wp) :: epsi10 = 1.e-10_wp ! 48 REAL(wp) :: epsi06 = 1.e-6_wp ! 50 49 51 50 !!---------------------------------------------------------------------- … … 110 109 CALL lim_var_glo2eqv ! only for info 111 110 112 !---------------------------------------------------------------------------------------- 113 ! 4) Computation of trend terms and get back to old values 114 !---------------------------------------------------------------------------------------- 115 116 !- Trend terms 117 d_a_i_thd(:,:,:) = a_i(:,:,:) - old_a_i(:,:,:) 118 d_v_s_thd(:,:,:) = v_s(:,:,:) - old_v_s(:,:,:) 119 d_v_i_thd(:,:,:) = v_i(:,:,:) - old_v_i(:,:,:) 120 d_e_s_thd(:,:,:,:) = e_s(:,:,:,:) - old_e_s(:,:,:,:) 121 d_e_i_thd(:,:,:,:) = e_i(:,:,:,:) - old_e_i(:,:,:,:) 122 !?? d_oa_i_thd(:,:,:) = oa_i (:,:,:) - old_oa_i (:,:,:) 123 d_smv_i_thd(:,:,:) = 0._wp 124 IF( num_sal == 2 ) d_smv_i_thd(:,:,:) = smv_i(:,:,:) - old_smv_i(:,:,:) 125 126 ! diag only (clem) 127 dv_dt_thd(:,:,:) = d_v_i_thd(:,:,:) * r1_rdtice * rday 128 129 IF(ln_ctl) THEN ! Control print 111 IF(ln_ctl) THEN ! Control print 130 112 CALL prt_ctl_info(' ') 131 113 CALL prt_ctl_info(' - Cell values : ') … … 183 165 ! ------------------------------- 184 166 ! 185 !- Recover Old values 186 a_i(:,:,:) = old_a_i (:,:,:) 187 v_s(:,:,:) = old_v_s (:,:,:) 188 v_i(:,:,:) = old_v_i (:,:,:) 189 e_s(:,:,:,:) = old_e_s (:,:,:,:) 190 e_i(:,:,:,:) = old_e_i (:,:,:,:) 191 !?? oa_i(:,:,:) = old_oa_i(:,:,:) 192 IF( num_sal == 2 ) smv_i(:,:,:) = old_smv_i(:,:,:) 193 ! 194 IF( nn_timing == 1 ) CALL timing_stop('limitd_th') 167 IF( nn_timing == 1 ) CALL timing_stop('limitd_th') 195 168 END SUBROUTINE lim_itd_th 196 169 ! … … 281 254 DO jj = 1, jpj 282 255 DO ji = 1, jpi 283 zindb = 1.0 -MAX(0.0,SIGN(1.0,-a_i(ji,jj,jl)+epsi10)) !0 if no ice and 1 if yes284 ht_i(ji,jj,jl) = v_i(ji,jj,jl) / MAX( a_i(ji,jj,jl),epsi10) * zindb285 zindb = 1.0 -MAX(0.0,SIGN(1.0,-old_a_i(ji,jj,jl)+epsi10)) !0 if no ice and 1 if yes286 zht_i_o(ji,jj,jl) = old_v_i(ji,jj,jl) / MAX( old_a_i(ji,jj,jl),epsi10) * zindb287 IF( a_i(ji,jj,jl) > 1e-6 ) zdhice(ji,jj,jl) = ht_i(ji,jj,jl) - zht_i_o(ji,jj,jl)256 zindb = 1.0 - MAX( 0.0, SIGN( 1.0, - a_i(ji,jj,jl) + epsi10 ) ) !0 if no ice and 1 if yes 257 ht_i(ji,jj,jl) = v_i(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi10 ) * zindb 258 zindb = 1.0 - MAX( 0.0, SIGN( 1.0, - old_a_i(ji,jj,jl) + epsi10) ) !0 if no ice and 1 if yes 259 zht_i_o(ji,jj,jl) = old_v_i(ji,jj,jl) / MAX( old_a_i(ji,jj,jl), epsi10 ) * zindb 260 IF( a_i(ji,jj,jl) > epsi06 ) zdhice(ji,jj,jl) = ht_i(ji,jj,jl) - zht_i_o(ji,jj,jl) 288 261 END DO 289 262 END DO … … 321 294 ! Tricky trick see limitd_me.F90 322 295 ! will be soon removed, CT 323 ! hi_max(kubnd) = 99 9.99296 ! hi_max(kubnd) = 99. 324 297 zhbnew(:,:,:) = 0._wp 325 298 … … 329 302 ij = nind_j(ji) 330 303 ! 331 IF ( ( zht_i_o(ii,ij,jl) .GT.epsi10 ) .AND. &332 ( zht_i_o(ii,ij,jl+1) .GT.epsi10 ) ) THEN304 IF ( ( zht_i_o(ii,ij,jl) .GT. epsi10 ) .AND. & 305 ( zht_i_o(ii,ij,jl+1) .GT. epsi10 ) ) THEN 333 306 !interpolate between adjacent category growth rates 334 307 zslope = ( zdhice(ii,ij,jl+1) - zdhice(ii,ij,jl) ) / & … … 402 375 zhbnew(ji,jj,kubnd) = 3._wp * ht_i(ji,jj,kubnd) - 2._wp * zhbnew(ji,jj,kubnd-1) 403 376 ELSE 404 zhbnew(ji,jj,kubnd) = hi_max(kubnd) 377 zhbnew(ji,jj,kubnd) = hi_max(kubnd) 378 !!? clem bug: since hi_max(jpl)=99, this limit is very high 379 !!? but I think it is erased in fitline subroutine 405 380 ENDIF 406 381 … … 447 422 * a_i(ii,ij,klbnd) / ( a_i(ii,ij,klbnd) - zda0 ) 448 423 a_i(ii,ij,klbnd) = a_i(ii,ij,klbnd) - zda0 449 v_i(ii,ij,klbnd) = a_i(ii,ij,klbnd)*ht_i(ii,ij,klbnd) ! clem @useless ?424 v_i(ii,ij,klbnd) = a_i(ii,ij,klbnd)*ht_i(ii,ij,klbnd) ! clem-useless ? 450 425 ENDIF ! zetamax > 0 451 426 ! ji, a_i > epsi10 … … 539 514 a_i(ii,ij,1) = a_i(ii,ij,1) * ht_i(ii,ij,1) / hiclim 540 515 ht_i(ii,ij,1) = hiclim 541 v_i(ii,ij,1) = a_i(ii,ij,1) * ht_i(ii,ij,1) !clem @useless516 v_i(ii,ij,1) = a_i(ii,ij,1) * ht_i(ii,ij,1) !clem-useless 542 517 ENDIF 543 518 END DO !ji … … 602 577 REAL(wp) :: zdhr ! 1 / (hR - hL) 603 578 REAL(wp) :: zwk1, zwk2 ! temporary variables 604 REAL(wp) :: zacrith ! critical minimum concentration in an ice category 605 !!------------------------------------------------------------------ 606 ! 607 zacrith = 1.0e-6 579 !!------------------------------------------------------------------ 580 ! 608 581 ! 609 582 DO jj = 1, jpj 610 583 DO ji = 1, jpi 611 584 ! 612 IF( zremap_flag(ji,jj) == 1 .AND. a_i(ji,jj,num_cat) > zacrith&585 IF( zremap_flag(ji,jj) == 1 .AND. a_i(ji,jj,num_cat) > epsi10 & 613 586 & .AND. hice(ji,jj) > 0._wp ) THEN 614 587 … … 1009 982 !zdaice(ji,jj,jl) = a_i(ji,jj,jl) 1010 983 !zdvice(ji,jj,jl) = v_i(ji,jj,jl) 1011 zdaice(ji,jj,jl) = a_i(ji,jj,jl)/21012 zdvice(ji,jj,jl) = v_i(ji,jj,jl)-zdaice(ji,jj,jl)*(hi_max(jl)+hi_max(jl-1))/2984 !zdaice(ji,jj,jl) = a_i(ji,jj,jl) * 0.5_wp 985 !zdvice(ji,jj,jl) = v_i(ji,jj,jl)-zdaice(ji,jj,jl)*(hi_max(jl)+hi_max(jl-1)) * 0.5_wp 1013 986 ! end TECLIM change 987 ! clem: how much of a_i you send in cat sup is somewhat arbitrary 988 zdaice(ji,jj,jl) = a_i(ji,jj,jl) * ( ht_i(ji,jj,jl) - hi_max(jl) + epsi10 ) / ht_i(ji,jj,jl) 989 zdvice(ji,jj,jl) = v_i(ji,jj,jl) - ( a_i(ji,jj,jl) - zdaice(ji,jj,jl) ) * ( hi_max(jl) - epsi10 ) 1014 990 ENDIF 1015 991 END DO ! ji … … 1038 1014 zshiftflag = 0 1039 1015 1016 !clem-change 1017 ! DO jj = 1, jpj 1018 ! DO ji = 1, jpi 1019 ! IF( a_i(ji,jj,jl+1) > epsi10 .AND. & 1020 ! ht_i(ji,jj,jl+1) <= hi_max(jl) ) THEN 1021 ! ! 1022 ! zshiftflag = 1 1023 ! zdonor(ji,jj,jl) = jl + 1 1024 ! zdaice(ji,jj,jl) = a_i(ji,jj,jl+1) 1025 ! zdvice(ji,jj,jl) = v_i(ji,jj,jl+1) 1026 ! ENDIF 1027 ! END DO ! ji 1028 ! END DO ! jj 1029 ! 1030 ! IF(lk_mpp) CALL mpp_max( zshiftflag ) 1031 ! 1032 ! IF( zshiftflag == 1 ) THEN ! Shift ice between categories 1033 ! CALL lim_itd_shiftice( klbnd, kubnd, zdonor, zdaice, zdvice ) 1034 ! ! Reset shift parameters 1035 ! zdonor(:,:,jl) = 0 1036 ! zdaice(:,:,jl) = 0._wp 1037 ! zdvice(:,:,jl) = 0._wp 1038 ! ENDIF 1039 !clem-change 1040 1041 ! clem-change begin: why not doing that? 1040 1042 DO jj = 1, jpj 1041 1043 DO ji = 1, jpi 1042 1044 IF( a_i(ji,jj,jl+1) > epsi10 .AND. & 1043 1045 ht_i(ji,jj,jl+1) <= hi_max(jl) ) THEN 1044 ! 1045 zshiftflag = 1 1046 zdonor(ji,jj,jl) = jl + 1 1047 zdaice(ji,jj,jl) = a_i(ji,jj,jl+1) 1048 zdvice(ji,jj,jl) = v_i(ji,jj,jl+1) 1046 ht_i(ji,jj,jl+1) = hi_max(jl) + epsi10 1047 a_i (ji,jj,jl+1) = v_i(ji,jj,jl+1) / ht_i(ji,jj,jl+1) 1049 1048 ENDIF 1050 1049 END DO ! ji 1051 1050 END DO ! jj 1052 1053 IF(lk_mpp) CALL mpp_max( zshiftflag ) 1054 1055 IF( zshiftflag == 1 ) THEN ! Shift ice between categories 1056 CALL lim_itd_shiftice( klbnd, kubnd, zdonor, zdaice, zdvice ) 1057 ! Reset shift parameters 1058 zdonor(:,:,jl) = 0 1059 zdaice(:,:,jl) = 0._wp 1060 zdvice(:,:,jl) = 0._wp 1061 ENDIF 1051 ! clem-change end 1062 1052 1063 1053 END DO ! jl -
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90
r4220 r4332 50 50 PUBLIC lim_rhg ! routine called by lim_dyn (or lim_dyn_2) 51 51 52 REAL(wp) :: epsi10 = 1.e-10_wp ! 52 53 REAL(wp) :: rzero = 0._wp ! constant values 53 54 REAL(wp) :: rone = 1._wp ! constant values … … 486 487 CALL lbc_lnk( zs12(:,:), 'F', 1. ) 487 488 488 !#if defined key_bdy489 ! ! clem: change zs1, zs2, zs12 at the boundary for each iteration490 ! CALL bdy_ice_lim_dyn( 2, zs1, zs2, zs12 )491 ! CALL lbc_lnk( zs1 (:,:), 'T', 1. )492 ! CALL lbc_lnk( zs2 (:,:), 'T', 1. )493 ! CALL lbc_lnk( zs12(:,:), 'F', 1. )494 !#endif495 496 489 ! Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) 497 490 DO jj = k_j1+1, k_jpj-1 … … 544 537 CALL agrif_rhg_lim2( jter, nevp, 'U' ) 545 538 #endif 539 #if defined key_bdy 540 ! clem: change u_ice and v_ice at the boundary for each iteration 541 CALL bdy_ice_lim_dyn( 'U' ) 542 #endif 546 543 547 544 !CDIR NOVERRCHK … … 572 569 CALL agrif_rhg_lim2( jter, nevp, 'V' ) 573 570 #endif 571 #if defined key_bdy 572 ! clem: change u_ice and v_ice at the boundary for each iteration 573 CALL bdy_ice_lim_dyn( 'V' ) 574 #endif 574 575 575 576 ELSE … … 599 600 CALL lbc_lnk( v_ice(:,:), 'V', -1. ) 600 601 #if defined key_agrif && defined key_lim2 601 CALL agrif_rhg_lim2( jter, nevp , 'V' ) 602 #endif 602 CALL agrif_rhg_lim2( jter, nevp, 'V' ) 603 #endif 604 #if defined key_bdy 605 ! clem: change u_ice and v_ice at the boundary for each iteration 606 CALL bdy_ice_lim_dyn( 'V' ) 607 #endif 603 608 604 609 !CDIR NOVERRCHK … … 632 637 CALL agrif_rhg_lim2( jter, nevp, 'U' ) 633 638 #endif 639 #if defined key_bdy 640 ! clem: change u_ice and v_ice at the boundary for each iteration 641 CALL bdy_ice_lim_dyn( 'U' ) 642 #endif 634 643 635 644 ENDIF 636 645 637 !#if defined key_bdy638 ! ! clem: change u_ice and v_ice at the boundary for each iteration639 ! CALL bdy_ice_lim_dyn( 1 )640 !#endif641 642 646 IF(ln_ctl) THEN 643 647 !--- Convergence test. … … 666 670 !CDIR NOVERRCHK 667 671 DO ji = fs_2, fs_jpim1 668 zindb = MAX( 0.0, SIGN( 1.0, at_i(ji,jj) - 1.0e-6) )669 !zdummy = zindb * vt_i(ji,jj) / MAX(at_i(ji,jj) , 1.0e-06)672 zindb = MAX( 0.0, SIGN( 1.0, at_i(ji,jj) - epsi10 ) ) 673 !zdummy = zindb * vt_i(ji,jj) / MAX(at_i(ji,jj) , epsi10 ) 670 674 zdummy = vt_i(ji,jj) 671 675 IF ( zdummy .LE. hminrhg ) THEN … … 684 688 #if defined key_bdy 685 689 ! clem: change u_ice and v_ice at the boundary 686 CALL bdy_ice_lim_dyn( 1 ) 690 CALL bdy_ice_lim_dyn( 'U' ) 691 CALL bdy_ice_lim_dyn( 'V' ) 687 692 #endif 688 693 689 694 DO jj = k_j1+1, k_jpj-1 690 695 DO ji = fs_2, fs_jpim1 691 zindb = MAX( 0.0, SIGN( 1.0, at_i(ji,jj) - 1.0e-6) )692 !zdummy = zindb * vt_i(ji,jj) / MAX(at_i(ji,jj) , 1.0e-06)696 zindb = MAX( 0.0, SIGN( 1.0, at_i(ji,jj) - epsi10 ) ) 697 !zdummy = zindb * vt_i(ji,jj) / MAX(at_i(ji,jj) , epsi10 ) 693 698 zdummy = vt_i(ji,jj) 694 699 IF ( zdummy .LE. hminrhg ) THEN … … 714 719 !- zdd(:,:), zdt(:,:): divergence and tension at centre 715 720 !- zds(:,:): shear on northeast corner of grid cells 716 zindb = MAX( 0.0, SIGN( 1.0, at_i(ji,jj) - 1.0e-6) )717 !zdummy = zindb * vt_i(ji,jj) / MAX(at_i(ji,jj) , 1.0e-06)721 zindb = MAX( 0.0, SIGN( 1.0, at_i(ji,jj) - epsi10 ) ) 722 !zdummy = zindb * vt_i(ji,jj) / MAX(at_i(ji,jj) , epsi10 ) 718 723 zdummy = vt_i(ji,jj) 719 724 IF ( zdummy .LE. hminrhg ) THEN -
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90
r4220 r4332 49 49 PUBLIC lim_sbc_tau ! called by sbc_ice_lim 50 50 51 REAL(wp) :: epsi16 = 1.e-16_wp ! constant values52 51 REAL(wp) :: rzero = 0._wp 53 52 REAL(wp) :: rone = 1._wp … … 145 144 146 145 ! computation the solar flux at ocean surface 147 IF (lk_cpl) THEN ! be carfeful: not be ingtested yet146 IF (lk_cpl) THEN ! be carfeful: not been tested yet 148 147 ! original line 149 148 !zfcm1 = qsr_tot(ji,jj) + fstric(ji,jj) * at_i(ji,jj) … … 188 187 qns(ji,jj) = zfcm2 - fdtcn(ji,jj) ! non solar heat flux 189 188 ! ! fdtcn : turbulent oceanic heat flux 190 191 !!gm this IF prevents the vertorisation of the whole loop192 ! IF ( ( ji == jiindx ) .AND. ( jj == jjindx) ) THEN193 ! WRITE(numout,*) ' lim_sbc : heat fluxes '194 ! WRITE(numout,*) ' qsr : ', qsr(jiindx,jjindx)195 ! WRITE(numout,*) ' pfrld : ', pfrld(jiindx,jjindx)196 ! WRITE(numout,*) ' fstric : ', fstric (jiindx,jjindx)197 ! WRITE(numout,*)198 ! WRITE(numout,*) ' qns : ', qns(jiindx,jjindx)199 ! WRITE(numout,*) ' fdtcn : ', fdtcn(jiindx,jjindx)200 ! WRITE(numout,*) ' ifral : ', ifral201 ! WRITE(numout,*) ' ial : ', ial202 ! WRITE(numout,*) ' qcmif : ', qcmif(jiindx,jjindx)203 ! WRITE(numout,*) ' qldif : ', qldif(jiindx,jjindx)204 ! !WRITE(numout,*) ' qcmif / dt: ', qcmif(jiindx,jjindx) * r1_rdtice205 ! !WRITE(numout,*) ' qldif / dt: ', qldif(jiindx,jjindx) * r1_rdtice206 ! WRITE(numout,*) ' ifrdv : ', ifrdv207 ! WRITE(numout,*) ' qfvbq : ', qfvbq(jiindx,jjindx)208 ! WRITE(numout,*) ' qdtcn : ', qdtcn(jiindx,jjindx)209 ! !WRITE(numout,*) ' qfvbq / dt: ', qfvbq(jiindx,jjindx) * r1_rdtice210 ! !WRITE(numout,*) ' qdtcn / dt: ', qdtcn(jiindx,jjindx) * r1_rdtice211 ! WRITE(numout,*) ' '212 ! WRITE(numout,*) ' fdtcn : ', fdtcn(jiindx,jjindx)213 ! WRITE(numout,*) ' fhmec : ', fhmec(jiindx,jjindx)214 ! WRITE(numout,*) ' fheat_mec : ', fheat_mec(jiindx,jjindx)215 ! WRITE(numout,*) ' fhbri : ', fhbri(jiindx,jjindx)216 ! WRITE(numout,*) ' fheat_res : ', fheat_res(jiindx,jjindx)217 ! ENDIF218 !!gm end219 189 END DO 220 190 END DO -
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90
r4220 r4332 50 50 PUBLIC lim_thd_init ! called by iceini module 51 51 52 REAL(wp) :: epsi20 = 1e-20_wp ! constant values 53 REAL(wp) :: epsi16 = 1e-16_wp ! 54 REAL(wp) :: epsi10 = 1e-10_wp ! 55 REAL(wp) :: epsi06 = 1e-06_wp ! 56 REAL(wp) :: epsi04 = 1e-04_wp ! 52 REAL(wp) :: epsi10 = 1.e-10_wp ! 57 53 REAL(wp) :: zzero = 0._wp ! 58 54 REAL(wp) :: zone = 1._wp ! … … 126 122 DO ji = 1, jpi 127 123 !Energy of melting q(S,T) [J.m-3] 128 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_i(ji,jj,jl) , epsi 06) ) * REAL( nlay_i )124 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_i(ji,jj,jl) , epsi10 ) ) * REAL( nlay_i ) 129 125 !0 if no ice and 1 if yes 130 zindb = 1.0 - MAX( 0.0 , SIGN( 1.0 , - ht_i(ji,jj,jl) ) )126 zindb = 1.0 - MAX( 0.0 , SIGN( 1.0 , - v_i(ji,jj,jl) + epsi10 ) ) 131 127 !convert units ! very important that this line is here 132 128 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * unit_fac * zindb … … 138 134 DO ji = 1, jpi 139 135 !Energy of melting q(S,T) [J.m-3] 140 e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_s(ji,jj,jl) , epsi 06) ) * REAL( nlay_s )136 e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_s(ji,jj,jl) , epsi10 ) ) * REAL( nlay_s ) 141 137 !0 if no ice and 1 if yes 142 zindb = 1.0 - MAX( 0.0 , SIGN( 1.0 , - ht_s(ji,jj,jl) ) )138 zindb = 1.0 - MAX( 0.0 , SIGN( 1.0 , - v_s(ji,jj,jl) + epsi10 ) ) 143 139 !convert units ! very important that this line is here 144 140 e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * unit_fac * zindb … … 147 143 END DO 148 144 END DO 149 150 !-----------------------------151 ! 1.3) Set some dummies to 0152 !-----------------------------153 !clem rdvosif(:,:) = 0.e0 ! variation of ice volume at surface154 !clem rdvobif(:,:) = 0.e0 ! variation of ice volume at bottom155 !clem fdvolif(:,:) = 0.e0 ! total variation of ice volume156 !clem rdvonif(:,:) = 0.e0 ! lateral variation of ice volume157 !clem fstric (:,:) = 0.e0 ! part of solar radiation transmitted through the ice158 !clem ffltbif(:,:) = 0.e0 ! linked with fstric159 !clem qfvbq (:,:) = 0.e0 ! linked with fstric160 !clem rdm_snw(:,:) = 0.e0 ! variation of snow mass per unit area161 !clem rdm_ice(:,:) = 0.e0 ! variation of ice mass per unit area162 !clem hicifp (:,:) = 0.e0 ! daily thermodynamic ice production.163 !clem sfx_bri(:,:) = 0.e0 ! brine flux contribution to salt flux to the ocean164 !clem fhbri (:,:) = 0.e0 ! brine flux contribution to heat flux to the ocean165 !clem sfx_thd(:,:) = 0.e0 ! equivalent salt flux to the ocean due to ice/growth decay166 145 167 146 !----------------------------------- … … 182 161 !CDIR NOVERRCHK 183 162 DO ji = 1, jpi 184 phicif(ji,jj) = vt_i(ji,jj) 185 pfrld(ji,jj) = 1.0 - at_i(ji,jj) 186 zinda = tms(ji,jj) * ( 1.0 - MAX( zzero , SIGN( zone , - at_i(ji,jj) ) ) ) 163 zinda = tms(ji,jj) * ( 1.0 - MAX( zzero , SIGN( zone , - at_i(ji,jj) + epsi10 ) ) ) 187 164 ! 188 165 ! ! solar irradiance transmission at the mixed layer bottom and used in the lead heat budget … … 195 172 196 173 ! here the drag will depend on ice thickness and type (0.006) 197 fdtcn(ji,jj) = zinda * rau0 * rcp * 0.006 * zfric_u * ( ( sst_m(ji,jj) + rt0) - t_bo(ji,jj) )174 fdtcn(ji,jj) = zinda * rau0 * rcp * 0.006 * zfric_u * ( ( sst_m(ji,jj) + rt0 ) - t_bo(ji,jj) ) 198 175 ! also category dependent 199 176 ! !-- Energy from the turbulent oceanic heat flux heat flux coming in the lead 200 qdtcn(ji,jj) = zinda * fdtcn(ji,jj) * ( 1.0 - at_i(ji,jj)) * rdt_ice177 qdtcn(ji,jj) = zinda * fdtcn(ji,jj) * ( 1.0 - at_i(ji,jj) ) * rdt_ice 201 178 ! 202 179 ! !-- Lead heat budget, qldif (part 1, next one is in limthd_dh) … … 214 191 zpareff = 1.0 - zinda * zfntlat 215 192 != 0 if ice and positive heat budget and 1 if one of those two is false 216 zqlbsbq(ji,jj) = qldif(ji,jj) * ( 1.0 - zpareff ) / MAX( at_i(ji,jj) * rdt_ice , epsi16)193 zqlbsbq(ji,jj) = qldif(ji,jj) * ( 1.0 - zpareff ) / ( rdt_ice * MAX( at_i(ji,jj), epsi10 ) ) 217 194 ! 218 195 ! Heat budget of the lead, energy transferred from ice to ocean … … 221 198 ! 222 199 ! Energy needed to bring ocean surface layer until its freezing (qcmif, limflx) 223 qcmif (ji,jj) = rau0 * rcp * fse3t_m(ji,jj,1) * ( t_bo(ji,jj) - ( sst_m(ji,jj) + rt0) )200 qcmif (ji,jj) = rau0 * rcp * fse3t_m(ji,jj,1) * ( t_bo(ji,jj) - ( sst_m(ji,jj) + rt0 ) ) 224 201 ! 225 202 ! oceanic heat flux (limthd_dh) 226 fbif (ji,jj) = zinda * ( fsbbq(ji,jj) / MAX( at_i(ji,jj) , epsi 20 ) + fdtcn(ji,jj) )203 fbif (ji,jj) = zinda * ( fsbbq(ji,jj) / MAX( at_i(ji,jj) , epsi10 ) + fdtcn(ji,jj) ) 227 204 ! 228 205 END DO … … 240 217 ENDIF 241 218 242 zareamin = 1.e-10219 zareamin = epsi10 243 220 nbpb = 0 244 221 DO jj = 1, jpj … … 248 225 npb(nbpb) = (jj - 1) * jpi + ji 249 226 ENDIF 250 ! debug point to follow 251 IF ( (ji.eq.jiindx).AND.(jj.eq.jjindx) ) THEN 252 jiindex_1d = nbpb 253 ENDIF 254 END DO 255 END DO 227 END DO 228 END DO 229 230 ! debug point to follow 231 jiindex_1d = 0 232 IF( ln_nicep ) THEN 233 DO ji = mi0(jiindx), mi1(jiindx) 234 DO jj = mj0(jjindx), mj1(jjindx) 235 jiindex_1d = (jj - 1) * jpi + ji 236 END DO 237 END DO 238 ENDIF 256 239 257 240 !------------------------------------------------------------------------------! … … 315 298 !-------------------------------- 316 299 317 IF( con_i ) CALL lim_thd_enmelt( 1, nbpb ) ! computes sea ice energy of melting318 IF( con_i ) CALL lim_thd_glohec( qt_i_in, qt_s_in, q_i_layer_in, 1, nbpb, jl )300 IF( con_i .AND. jiindex_1d > 0 ) CALL lim_thd_enmelt( 1, nbpb ) ! computes sea ice energy of melting 301 IF( con_i .AND. jiindex_1d > 0 ) CALL lim_thd_glohec( qt_i_in, qt_s_in, q_i_layer_in, 1, nbpb, jl ) 319 302 320 303 ! !---------------------------------! … … 324 307 CALL lim_thd_enmelt( 1, nbpb ) ! computes sea ice energy of melting compulsory for limthd_dh 325 308 326 IF( con_i ) CALL lim_thd_glohec ( qt_i_fin, qt_s_fin, q_i_layer_fin, 1, nbpb, jl )327 IF( con_i ) CALL lim_thd_con_dif( 1 , nbpb , jl )309 IF( con_i .AND. jiindex_1d > 0 ) CALL lim_thd_glohec ( qt_i_fin, qt_s_fin, q_i_layer_fin, 1, nbpb, jl ) 310 IF( con_i .AND. jiindex_1d > 0 ) CALL lim_thd_con_dif( 1 , nbpb , jl ) 328 311 329 312 ! !---------------------------------! … … 340 323 341 324 ! CALL lim_thd_enmelt(1,nbpb) ! computes sea ice energy of melting 342 IF( con_i ) CALL lim_thd_glohec( qt_i_fin, qt_s_fin, q_i_layer_fin, 1, nbpb, jl )343 IF( con_i ) CALL lim_thd_con_dh ( 1 , nbpb , jl )325 IF( con_i .AND. jiindex_1d > 0 ) CALL lim_thd_glohec( qt_i_fin, qt_s_fin, q_i_layer_fin, 1, nbpb, jl ) 326 IF( con_i .AND. jiindex_1d > 0 ) CALL lim_thd_con_dh ( 1 , nbpb , jl ) 344 327 345 328 !-------------------------------- … … 431 414 !clem@mv-to-itd dv_dt_thd(:,:,:) = d_v_i_thd(:,:,:) * r1_rdtice * rday 432 415 433 IF( con_i ) fbif(:,:) = fbif(:,:) + zqlbsbq(:,:)416 IF( con_i .AND. jiindex_1d > 0 ) fbif(:,:) = fbif(:,:) + zqlbsbq(:,:) 434 417 435 418 IF(ln_ctl) THEN ! Control print … … 522 505 END DO 523 506 ! 524 IF(lwp)WRITE(numout,*) ' lim_thd_glohec '525 IF(lwp)WRITE(numout,*) ' qt_i_in : ', eti(jiindex_1d,jl) * r1_rdtice526 IF(lwp)WRITE(numout,*) ' qt_s_in : ', ets(jiindex_1d,jl) * r1_rdtice527 IF(lwp)WRITE(numout,*) ' qt_in : ', ( eti(jiindex_1d,jl) + ets(jiindex_1d,jl) ) * r1_rdtice507 WRITE(numout,*) ' lim_thd_glohec ' 508 WRITE(numout,*) ' qt_i_in : ', eti(jiindex_1d,jl) * r1_rdtice 509 WRITE(numout,*) ' qt_s_in : ', ets(jiindex_1d,jl) * r1_rdtice 510 WRITE(numout,*) ' qt_in : ', ( eti(jiindex_1d,jl) + ets(jiindex_1d,jl) ) * r1_rdtice 528 511 ! 529 512 END SUBROUTINE lim_thd_glohec … … 611 594 WRITE(numout,*) ' Number of points where there is a surf err gt than surf_err : ', numce, numit 612 595 613 IF (jiindex_1D.GT.0)WRITE(numout,*) ' fc_su : ', fc_su(jiindex_1d)614 IF (jiindex_1D.GT.0)WRITE(numout,*) ' fatm : ', fatm(jiindex_1d,jl)615 IF (jiindex_1D.GT.0)WRITE(numout,*) ' t_su : ', t_su_b(jiindex_1d)596 WRITE(numout,*) ' fc_su : ', fc_su(jiindex_1d) 597 WRITE(numout,*) ' fatm : ', fatm(jiindex_1d,jl) 598 WRITE(numout,*) ' t_su : ', t_su_b(jiindex_1d) 616 599 617 600 !--------------------------------------- -
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90
r4099 r4332 32 32 PUBLIC lim_thd_dh ! called by lim_thd 33 33 34 REAL(wp) :: epsi20 = 1 e-20 ! constant values35 REAL(wp) :: epsi1 3 = 1e-13!36 REAL(wp) :: epsi1 6 = 1e-16!37 REAL(wp) :: zzero = 0. e0!38 REAL(wp) :: zone = 1. e0!34 REAL(wp) :: epsi20 = 1.e-20 ! constant values 35 REAL(wp) :: epsi10 = 1.e-10 ! 36 REAL(wp) :: epsi13 = 1.e-13 ! 37 REAL(wp) :: zzero = 0._wp ! 38 REAL(wp) :: zone = 1._wp ! 39 39 40 40 !!---------------------------------------------------------------------- … … 125 125 INTEGER :: num_iter_max, numce_dh 126 126 REAL(wp) :: meance_dh 127 REAL(wp) :: zinda 127 128 REAL(wp), POINTER, DIMENSION(:) :: zinnermelt 128 129 REAL(wp), POINTER, DIMENSION(:) :: zfbase, zdq_i … … 288 289 END DO 289 290 290 ! !-------------------291 IF( con_i ) THEN ! Conservation test292 ! !-------------------291 ! !------------------- 292 IF( con_i .AND. jiindex_1d > 0 ) THEN ! Conservation test 293 ! !------------------- 293 294 numce_dh = 0 294 295 meance_dh = 0._wp … … 494 495 END DO ! jk 495 496 496 ! !-------------------497 IF( con_i ) THEN ! Conservation test498 ! !-------------------497 ! !------------------- 498 IF( con_i .AND. jiindex_1d > 0 ) THEN ! Conservation test 499 ! !------------------- 499 500 DO ji = kideb, kiut 500 501 IF( ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) >= 0.e0 ) THEN … … 544 545 sfx_thd_1d(ji) = sfx_thd_1d(ji) - sm_i_b(ji) * a_i_b(ji) * zdvres * rhoic * r1_rdtice 545 546 ! ! excessive energy is sent to lateral ablation 546 fsup (ji) = rhoic * lfus * at_i_b(ji) / MAX( 1.0 - at_i_b(ji) , epsi13 ) * zdvres * r1_rdtice 547 zinda = MAX( 0._wp, SIGN( 1._wp , 1.0 - at_i_b(ji) - epsi10 ) ) 548 fsup(ji) = zinda * rhoic * lfus * at_i_b(ji) / MAX( 1.0 - at_i_b(ji) , epsi10 ) * zdvres * r1_rdtice 547 549 END DO 548 550 … … 577 579 578 580 ! energy of melting of remaining snow 579 zqt_s(ji) = ( 1. - zihg ) * zqt_s(ji) / MAX( zhni, epsi13 ) 581 zinda = MAX( 0._wp, SIGN( 1._wp , zhni - epsi10 ) ) 582 zqt_s(ji) = ( 1. - zihg ) * zqt_s(ji) / MAX( zhni, epsi10 ) * zinda 580 583 zdhnm = - ( 1. - zihg ) * ( 1. - zihgnew ) * zfdt_final(ji) / MAX( zqt_s(ji) , epsi13 ) 581 584 zhnfi = zhni + zdhnm … … 613 616 ENDIF 614 617 ! 615 ! diagnostic ( bottom ice growth )618 ! diagnostic 616 619 ii = MOD( npb(ji) - 1, jpi ) + 1 617 620 ij = ( npb(ji) - 1 ) / jpi + 1 … … 689 692 ! clem: new salinity difference stored (to be used in limthd_ent.F90) 690 693 IF ( num_sal == 2 ) THEN 691 i_ice_switch = 1.0 - MAX( 0.e0 , SIGN( 1.0 , - zhgnew(ji) + epsi13) )694 i_ice_switch = MAX( 0._wp , SIGN( 1._wp , zhgnew(ji) - epsi10 ) ) 692 695 ! salinity dif due to snow-ice formation 693 dsm_i_si_1d(ji) = ( zsm_snowice - sm_i_b(ji) ) * dh_snowice(ji) / MAX( zhgnew(ji), epsi1 3) * i_ice_switch696 dsm_i_si_1d(ji) = ( zsm_snowice - sm_i_b(ji) ) * dh_snowice(ji) / MAX( zhgnew(ji), epsi10 ) * i_ice_switch 694 697 ! salinity dif due to bottom growth 695 698 IF ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) < 0._wp ) THEN 696 dsm_i_se_1d(ji) = ( s_i_new(ji) - sm_i_b(ji) ) * dh_i_bott(ji) / MAX( zhgnew(ji), epsi1 3) * i_ice_switch699 dsm_i_se_1d(ji) = ( s_i_new(ji) - sm_i_b(ji) ) * dh_i_bott(ji) / MAX( zhgnew(ji), epsi10 ) * i_ice_switch 697 700 ENDIF 698 701 ENDIF -
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90
r4220 r4332 31 31 PUBLIC lim_thd_dif ! called by lim_thd 32 32 33 REAL(wp) :: epsi20 = 1e-20 ! constant values 34 REAL(wp) :: epsi13 = 1e-13 ! constant values 35 33 REAL(wp) :: epsi10 = 1.e-10_wp ! 36 34 !!---------------------------------------------------------------------- 37 35 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) … … 107 105 INTEGER, DIMENSION(kiut) :: numeqmax ! reference number of bottom equation 108 106 INTEGER, DIMENSION(kiut) :: isnow ! switch for presence (1) or absence (0) of snow 109 REAL(wp) :: zeps = 1.e-10_wp !110 107 REAL(wp) :: zg1s = 2._wp ! for the tridiagonal system 111 108 REAL(wp) :: zg1 = 2._wp ! … … 114 111 REAL(wp) :: zraext_s = 1.e+8_wp ! extinction coefficient of radiation in the snow 115 112 REAL(wp) :: zkimin = 0.10_wp ! minimum ice thermal conductivity 116 REAL(wp) :: zht_smin = 1.e-4_wp ! minimum snow depth117 113 REAL(wp) :: ztmelt_i ! ice melting temperature 118 114 REAL(wp) :: zerritmax ! current maximal error on temperature … … 312 308 IF( thcon_i_swi == 0 ) THEN ! Untersteiner (1964) formula 313 309 DO ji = kideb , kiut 314 ztcond_i(ji,0) = rcdic + zbeta*s_i_b(ji,1) / MIN(- zeps,t_i_b(ji,1)-rtt)310 ztcond_i(ji,0) = rcdic + zbeta*s_i_b(ji,1) / MIN(-epsi10,t_i_b(ji,1)-rtt) 315 311 ztcond_i(ji,0) = MAX(ztcond_i(ji,0),zkimin) 316 312 END DO … … 318 314 DO ji = kideb , kiut 319 315 ztcond_i(ji,layer) = rcdic + zbeta*( s_i_b(ji,layer) + s_i_b(ji,layer+1) ) / & 320 MIN(-2.0_wp * zeps, t_i_b(ji,layer)+t_i_b(ji,layer+1) - 2.0_wp * rtt)316 MIN(-2.0_wp * epsi10, t_i_b(ji,layer)+t_i_b(ji,layer+1) - 2.0_wp * rtt) 321 317 ztcond_i(ji,layer) = MAX(ztcond_i(ji,layer),zkimin) 322 318 END DO … … 326 322 IF( thcon_i_swi == 1 ) THEN ! Pringle et al formula included: 2.11 + 0.09 S/T - 0.011.T 327 323 DO ji = kideb , kiut 328 ztcond_i(ji,0) = rcdic + 0.090_wp * s_i_b(ji,1) / MIN( - zeps, t_i_b(ji,1)-rtt ) &324 ztcond_i(ji,0) = rcdic + 0.090_wp * s_i_b(ji,1) / MIN( -epsi10, t_i_b(ji,1)-rtt ) & 329 325 & - 0.011_wp * ( t_i_b(ji,1) - rtt ) 330 326 ztcond_i(ji,0) = MAX( ztcond_i(ji,0), zkimin ) … … 333 329 DO ji = kideb , kiut 334 330 ztcond_i(ji,layer) = rcdic + 0.090_wp * ( s_i_b(ji,layer) + s_i_b(ji,layer+1) ) & 335 & / MIN(-2.0_wp * zeps, t_i_b(ji,layer)+t_i_b(ji,layer+1) - 2.0_wp * rtt) &331 & / MIN(-2.0_wp * epsi10, t_i_b(ji,layer)+t_i_b(ji,layer+1) - 2.0_wp * rtt) & 336 332 & - 0.0055_wp* ( t_i_b(ji,layer) + t_i_b(ji,layer+1) - 2.0*rtt ) 337 333 ztcond_i(ji,layer) = MAX( ztcond_i(ji,layer), zkimin ) … … 339 335 END DO 340 336 DO ji = kideb , kiut 341 ztcond_i(ji,nlay_i) = rcdic + 0.090_wp * s_i_b(ji,nlay_i) / MIN(- zeps,t_bo_b(ji)-rtt) &337 ztcond_i(ji,nlay_i) = rcdic + 0.090_wp * s_i_b(ji,nlay_i) / MIN(-epsi10,t_bo_b(ji)-rtt) & 342 338 & - 0.011_wp * ( t_bo_b(ji) - rtt ) 343 339 ztcond_i(ji,nlay_i) = MAX( ztcond_i(ji,nlay_i), zkimin ) … … 352 348 353 349 !-- Snow kappa factors 354 zkappa_s(ji,0) = rcdsn / MAX( zeps,zh_s(ji))355 zkappa_s(ji,nlay_s) = rcdsn / MAX( zeps,zh_s(ji))350 zkappa_s(ji,0) = rcdsn / MAX(epsi10,zh_s(ji)) 351 zkappa_s(ji,nlay_s) = rcdsn / MAX(epsi10,zh_s(ji)) 356 352 END DO 357 353 … … 359 355 DO ji = kideb , kiut 360 356 zkappa_s(ji,layer) = 2.0 * rcdsn / & 361 MAX( zeps,2.0*zh_s(ji))357 MAX(epsi10,2.0*zh_s(ji)) 362 358 END DO 363 359 END DO … … 367 363 !-- Ice kappa factors 368 364 zkappa_i(ji,layer) = 2.0*ztcond_i(ji,layer)/ & 369 MAX( zeps,2.0*zh_i(ji))370 END DO 371 END DO 372 373 DO ji = kideb , kiut 374 zkappa_i(ji,0) = ztcond_i(ji,0)/MAX( zeps,zh_i(ji))375 zkappa_i(ji,nlay_i) = ztcond_i(ji,nlay_i) / MAX( zeps,zh_i(ji))365 MAX(epsi10,2.0*zh_i(ji)) 366 END DO 367 END DO 368 369 DO ji = kideb , kiut 370 zkappa_i(ji,0) = ztcond_i(ji,0)/MAX(epsi10,zh_i(ji)) 371 zkappa_i(ji,nlay_i) = ztcond_i(ji,nlay_i) / MAX(epsi10,zh_i(ji)) 376 372 !-- Interface 377 zkappa_s(ji,nlay_s) = 2.0*rcdsn*ztcond_i(ji,0)/MAX( zeps, &373 zkappa_s(ji,nlay_s) = 2.0*rcdsn*ztcond_i(ji,0)/MAX(epsi10, & 378 374 (ztcond_i(ji,0)*zh_s(ji) + rcdsn*zh_i(ji))) 379 375 zkappa_i(ji,0) = zkappa_s(ji,nlay_s)*REAL( isnow(ji) ) & … … 389 385 ztitemp(ji,layer) = t_i_b(ji,layer) 390 386 zspeche_i(ji,layer) = cpic + zgamma*s_i_b(ji,layer)/ & 391 MAX((t_i_b(ji,layer)-rtt)*(ztiold(ji,layer)-rtt), zeps)387 MAX((t_i_b(ji,layer)-rtt)*(ztiold(ji,layer)-rtt),epsi10) 392 388 zeta_i(ji,layer) = rdt_ice / MAX(rhoic*zspeche_i(ji,layer)*zh_i(ji), & 393 zeps)389 epsi10) 394 390 END DO 395 391 END DO … … 398 394 DO ji = kideb , kiut 399 395 ztstemp(ji,layer) = t_s_b(ji,layer) 400 zeta_s(ji,layer) = rdt_ice / MAX(rhosn*cpic*zh_s(ji), zeps)396 zeta_s(ji,layer) = rdt_ice / MAX(rhosn*cpic*zh_s(ji),epsi10) 401 397 END DO 402 398 END DO … … 707 703 END DO ! End of the do while iterative procedure 708 704 709 IF( ln_nicep ) THEN705 IF( ln_nicep .AND. lwp ) THEN 710 706 WRITE(numout,*) ' zerritmax : ', zerritmax 711 707 WRITE(numout,*) ' nconv : ', nconv … … 732 728 ! Heat conservation ! 733 729 !-------------------------! 734 IF( con_i ) THEN730 IF( con_i .AND. jiindex_1d > 0 ) THEN 735 731 DO ji = kideb, kiut 736 732 ! Upper snow value -
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_ent.F90
r4072 r4332 36 36 PUBLIC lim_thd_ent ! called by lim_thd 37 37 38 REAL(wp) :: epsi20 = 1e-20_wp ! constant values 39 REAL(wp) :: epsi13 = 1e-13_wp ! 40 REAL(wp) :: epsi10 = 1e-10_wp ! 41 REAL(wp) :: epsi06 = 1e-06_wp ! 38 REAL(wp) :: epsi20 = 1.e-20_wp ! constant values 39 REAL(wp) :: epsi10 = 1.e-10_wp ! 42 40 REAL(wp) :: zzero = 0._wp ! 43 41 REAL(wp) :: zone = 1._wp ! … … 122 120 REAL(wp), POINTER, DIMENSION(:,:) :: zhl0 ! old and new layer thicknesses 123 121 REAL(wp), POINTER, DIMENSION(:,:) :: zrl01 122 123 REAL(wp) :: zinda 124 124 !!------------------------------------------------------------------- 125 125 … … 391 391 DO layer1 = ntop1, nbot1 392 392 DO ji = kideb, kiut 393 zrl01(layer1,layer0) = MAX(0.0,( MIN(zm0(ji,layer0),z_s(ji,layer1)) & 393 zinda = MAX( 0._wp, SIGN( 1._wp , zhl0(ji,layer0) - epsi10 ) ) 394 zrl01(layer1,layer0) = zinda * MAX(0.0,( MIN(zm0(ji,layer0),z_s(ji,layer1)) & 394 395 & - MAX(zm0(ji,layer0-1), z_s(ji,layer1-1))) / MAX(zhl0(ji,layer0),epsi10)) 395 396 q_s_b(ji,layer1) = q_s_b(ji,layer1) + zrl01(layer1,layer0)*qm0(ji,layer0) & … … 407 408 END DO 408 409 409 IF ( con_i ) THEN410 IF ( con_i .AND. jiindex_1d > 0 ) THEN 410 411 DO ji = kideb, kiut 411 412 IF ( ABS ( zqts_in(ji) - zqts_fin(ji) ) * r1_rdtice > 1.0e-6 ) THEN … … 429 430 DO jk = 1, nlay_s 430 431 DO ji = kideb, kiut 431 q_s_b(ji,jk) = q_s_b(ji,jk) / MAX( zh_s(ji) , epsi20 ) 432 zinda = MAX( 0._wp, SIGN( 1._wp , zh_s(ji) - epsi10 ) ) 433 q_s_b(ji,jk) = zinda * q_s_b(ji,jk) / MAX( zh_s(ji) , epsi10 ) 432 434 END DO !ji 433 435 END DO !jk … … 555 557 DO ji = kideb, kiut 556 558 ! Heat conservation 557 zqti_in(ji) = zqti_in(ji) + qm0(ji,jk) * MAX( 0.0 , SIGN(1.0,ht_i_b(ji)-epsi 06) ) &559 zqti_in(ji) = zqti_in(ji) + qm0(ji,jk) * MAX( 0.0 , SIGN(1.0,ht_i_b(ji)-epsi10) ) & 558 560 & * MAX( 0.0 , SIGN( 1. , REAL(nbot0(ji) - jk) ) ) 559 561 END DO … … 601 603 DO layer1 = ntop1, nbot1 602 604 DO ji = kideb, kiut 603 zrl01(layer1,layer0) = MAX(0.0,( MIN(zm0(ji,layer0),z_i(ji,layer1)) & 605 zinda = MAX( 0._wp, SIGN( 1._wp , zhl0(ji,layer0) - epsi10 ) ) 606 zrl01(layer1,layer0) = zinda * MAX(0.0,( MIN(zm0(ji,layer0),z_i(ji,layer1)) & 604 607 - MAX(zm0(ji,layer0-1), z_i(ji,layer1-1)))/MAX(zhl0(ji,layer0),epsi10)) 605 608 q_i_b(ji,layer1) = q_i_b(ji,layer1) & 606 609 + zrl01(layer1,layer0)*qm0(ji,layer0) & 607 * MAX(0.0,SIGN(1.0,ht_i_b(ji)-epsi 06)) &610 * MAX(0.0,SIGN(1.0,ht_i_b(ji)-epsi10)) & 608 611 * MAX(0.0,SIGN(1.0,REAL(nbot0(ji)-layer0))) 609 612 END DO … … 621 624 END DO 622 625 ! 623 IF ( con_i ) THEN626 IF ( con_i .AND. jiindex_1d > 0 ) THEN 624 627 DO ji = kideb, kiut 625 628 IF ( ABS ( zqti_in(ji) - zqti_fin(ji) ) * r1_rdtice > 1.0e-6 ) THEN … … 646 649 DO jk = 1, nlay_i 647 650 DO ji = kideb, kiut 648 q_i_b(ji,jk) = q_i_b(ji,jk) / MAX( zh_i(ji) , epsi20 ) 651 zinda = MAX( 0._wp, SIGN( 1._wp , zh_i(ji) - epsi10 ) ) 652 q_i_b(ji,jk) = zinda * q_i_b(ji,jk) / MAX( zh_i(ji) , epsi10 ) 649 653 END DO !ji 650 654 END DO !jk -
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90
r4045 r4332 36 36 PUBLIC lim_thd_lac ! called by lim_thd 37 37 38 REAL(wp) :: epsi20 = 1e-20_wp ! constant values 39 REAL(wp) :: epsi13 = 1e-13_wp ! 40 REAL(wp) :: epsi11 = 1e-11_wp ! 41 REAL(wp) :: epsi10 = 1e-10_wp ! 42 REAL(wp) :: epsi06 = 1e-06_wp ! 43 REAL(wp) :: epsi03 = 1e-03_wp ! 38 REAL(wp) :: epsi10 = 1.e-10_wp ! 44 39 REAL(wp) :: zzero = 0._wp ! 45 40 REAL(wp) :: zone = 1._wp ! … … 160 155 !Energy of melting q(S,T) [J.m-3] 161 156 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / MAX( area(ji,jj) * v_i(ji,jj,jl) , epsi10 ) * REAL( nlay_i ) 162 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , -v_i(ji,jj,jl) ) ) !0 if no ice and 1 if yes157 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , -v_i(ji,jj,jl) + epsi10 ) ) !0 if no ice and 1 if yes 163 158 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * unit_fac * zindb 164 159 END DO … … 286 281 nbpac = nbpac + 1 287 282 npac( nbpac ) = (jj - 1) * jpi + ji 288 IF( ji == jiindx .AND. jj == jjindx ) jiindex_1d = nbpac289 283 ENDIF 290 284 END DO 291 285 END DO 292 286 293 IF( ln_nicep ) WRITE(numout,*) 'lim_thd_lac : nbpac = ', nbpac 287 ! debug point to follow 288 jiindex_1d = 0 289 IF( ln_nicep ) THEN 290 DO ji = mi0(jiindx), mi1(jiindx) 291 DO jj = mj0(jjindx), mj1(jjindx) 292 IF ( tms(ji,jj) * ( qcmif(ji,jj) - qldif(ji,jj) ) > 0._wp ) THEN 293 jiindex_1d = (jj - 1) * jpi + ji 294 ENDIF 295 END DO 296 END DO 297 ENDIF 298 299 IF( ln_nicep ) WRITE(numout,*) 'lim_thd_lac : nbpac = ', nbpac 294 300 295 301 !------------------------------ … … 499 505 END DO 500 506 501 IF( ln_nicep ) WRITE(numout,*) ' zv_i_ac : ', zv_i_ac(jiindx, 1:jpl)507 IF( ln_nicep .AND. jiindex_1d > 0 ) WRITE(numout,*) ' zv_i_ac : ', zv_i_ac(jiindex_1d, 1:jpl) 502 508 DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2) 503 509 DO ji = 1, nbpac 504 510 zindb = MAX( 0._wp, SIGN( 1._wp , zdv_res(ji) ) ) 505 zinda = MAX( 0._wp, SIGN( 1._wp , zat_i_lev(ji) - epsi 06) ) ! clem506 zv_i_ac(ji,jl) = zv_i_ac(ji,jl) + zindb * zinda * zdv_res(ji) * za_i_ac(ji,jl) / MAX( zat_i_lev(ji) , epsi 06)507 END DO 508 END DO 509 IF( ln_nicep ) WRITE(numout,*) ' zv_i_ac : ', zv_i_ac(jiindx, 1:jpl)511 zinda = MAX( 0._wp, SIGN( 1._wp , zat_i_lev(ji) - epsi10 ) ) ! clem 512 zv_i_ac(ji,jl) = zv_i_ac(ji,jl) + zindb * zinda * zdv_res(ji) * za_i_ac(ji,jl) / MAX( zat_i_lev(ji) , epsi10 ) 513 END DO 514 END DO 515 IF( ln_nicep .AND. jiindex_1d > 0 ) WRITE(numout,*) ' zv_i_ac : ', zv_i_ac(jiindex_1d, 1:jpl) 510 516 511 517 !--------------------------------- … … 649 655 ! CALL lim_cons_check (et_s_init, et_s_final, 1.0e-3, fieldid) 650 656 IF( ln_nicep ) THEN 651 WRITE(numout,*) ' vt_i_init : ', vt_i_init(jiindx,jjindx) 652 WRITE(numout,*) ' vt_i_final: ', vt_i_final(jiindx,jjindx) 653 WRITE(numout,*) ' et_i_init : ', et_i_init(jiindx,jjindx) 654 WRITE(numout,*) ' et_i_final: ', et_i_final(jiindx,jjindx) 657 DO ji = mi0(jiindx), mi1(jiindx) 658 DO jj = mj0(jjindx), mj1(jjindx) 659 WRITE(numout,*) ' vt_i_init : ', vt_i_init (ji,jj) 660 WRITE(numout,*) ' vt_i_final: ', vt_i_final(ji,jj) 661 WRITE(numout,*) ' et_i_init : ', et_i_init (ji,jj) 662 WRITE(numout,*) ' et_i_final: ', et_i_final(ji,jj) 663 END DO 664 END DO 655 665 ENDIF 656 666 ! -
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90
r4220 r4332 36 36 PUBLIC lim_trp ! called by ice_step 37 37 38 REAL(wp) :: epsi06 = 1.e-06_wp ! constant values39 REAL(wp) :: epsi03 = 1.e-03_wp40 38 REAL(wp) :: epsi10 = 1.e-10_wp 41 REAL(wp) :: epsi16 = 1.e-16_wp42 REAL(wp) :: epsi20 = 1.e-20_wp43 39 REAL(wp) :: rzero = 0._wp 44 40 REAL(wp) :: rone = 1._wp 45 46 REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:,:,:) :: zs0e47 41 48 42 !! * Substitution … … 448 442 449 443 ! Switches and dummy variables 450 zusvosn = 1.0/MAX( v_s(ji,jj,jl) , epsi1 6)451 zusvoic = 1.0/MAX( v_i(ji,jj,jl) , epsi1 6)444 zusvosn = 1.0/MAX( v_s(ji,jj,jl) , epsi10 ) 445 zusvoic = 1.0/MAX( v_i(ji,jj,jl) , epsi10 ) 452 446 zindsn = MAX( rzero, SIGN( rone, v_s(ji,jj,jl) - epsi10 ) ) 453 447 zindic = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi10 ) ) … … 460 454 ENDIF 461 455 462 zage = MAX( MIN( zbigval, zs0oi(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi1 6) ), 0._wp ) * a_i(ji,jj,jl)456 zage = MAX( MIN( zbigval, zs0oi(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi10 ) ), 0._wp ) * a_i(ji,jj,jl) 463 457 oa_i (ji,jj,jl) = zindic * zage 464 458 … … 500 494 at_i(ji,jj) = at_i(ji,jj) + a_i(ji,jj,jl) ! ice concentration 501 495 ! 502 zinda = MAX( rzero , SIGN( rone , at_i(ji,jj) - epsi1 6) )503 icethi(ji,jj) = vt_i(ji,jj) / MAX( at_i(ji,jj) , epsi1 6) * zinda ! ice thickness496 zinda = MAX( rzero , SIGN( rone , at_i(ji,jj) - epsi10 ) ) 497 icethi(ji,jj) = vt_i(ji,jj) / MAX( at_i(ji,jj) , epsi10 ) * zinda ! ice thickness 504 498 END DO 505 499 END DO -
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90
r4099 r4332 48 48 PUBLIC lim_update1 ! routine called by ice_step 49 49 50 REAL(wp) :: epsi06 = 1.e-06_wp ! module constants51 REAL(wp) :: epsi04 = 1.e-04_wp ! - -52 REAL(wp) :: epsi03 = 1.e-03_wp ! - -53 50 REAL(wp) :: epsi10 = 1.e-10_wp ! - - 54 REAL(wp) :: epsi16 = 1.e-16_wp ! - -55 REAL(wp) :: epsi20 = 1.e-20_wp ! - -56 51 REAL(wp) :: rzero = 0._wp ! - - 57 52 REAL(wp) :: rone = 1._wp ! - - … … 108 103 !------------------------------------------------------------------------------ 109 104 110 !--------------------- 111 ! Ice dynamics 112 !--------------------- 113 u_ice(:,:) = u_ice(:,:) + d_u_ice_dyn(:,:) 114 v_ice(:,:) = v_ice(:,:) + d_v_ice_dyn(:,:) 115 116 !----------------------------- 117 ! Update ice and snow volumes 118 !----------------------------- 119 DO jl = 1, jpl 120 v_i(:,:,jl) = v_i(:,:,jl) + d_v_i_trp(:,:,jl) 121 v_s(:,:,jl) = v_s(:,:,jl) + d_v_s_trp(:,:,jl) 122 END DO 123 124 125 !--------------------------------------------- 126 ! Ice concentration and ice heat content 127 !--------------------------------------------- 128 a_i (:,:,:) = a_i (:,:,:) + d_a_i_trp(:,:,:) 129 e_i(:,:,:,:) = e_i(:,:,:,:) + d_e_i_trp(:,:,:,:) 130 131 !------------------------------ 132 ! Snow temperature and ice age 133 !------------------------------ 134 e_s (:,:,:,:) = e_s (:,:,:,:) + d_e_s_trp (:,:,:,:) 135 oa_i(:,:,:) = oa_i(:,:,:) + d_oa_i_trp(:,:,:) 136 137 !-------------- 138 ! Ice salinity 139 !-------------- 140 IF( num_sal == 2 ) THEN ! general case 141 smv_i(:,:,:) = smv_i(:,:,:) + d_smv_i_trp(:,:,:) 142 ENDIF 105 !----------------- 106 ! Trend terms 107 !----------------- 108 d_u_ice_dyn(:,:) = u_ice(:,:) - old_u_ice(:,:) 109 d_v_ice_dyn(:,:) = v_ice(:,:) - old_v_ice(:,:) 110 d_a_i_trp (:,:,:) = a_i (:,:,:) - old_a_i (:,:,:) 111 d_v_s_trp (:,:,:) = v_s (:,:,:) - old_v_s (:,:,:) 112 d_v_i_trp (:,:,:) = v_i (:,:,:) - old_v_i (:,:,:) 113 d_e_s_trp (:,:,:,:) = e_s (:,:,:,:) - old_e_s (:,:,:,:) 114 d_e_i_trp (:,:,:,:) = e_i (:,:,:,:) - old_e_i (:,:,:,:) 115 d_oa_i_trp (:,:,:) = oa_i (:,:,:) - old_oa_i (:,:,:) 116 d_smv_i_trp(:,:,:) = 0._wp 117 IF( num_sal == 2 ) d_smv_i_trp(:,:,:) = smv_i(:,:,:) - old_smv_i(:,:,:) 143 118 144 119 ! mass and salt flux init (clem) … … 160 135 CALL lim_var_glo2eqv 161 136 162 !---------------------------------163 ! Classify the pathological cases164 !---------------------------------165 ! (1) v_i (new) > 0; d_v_i_thd + v_i(old) > 0 (easy case)166 ! (3) v_i (new) < 0; d_v_i_thd + v_i(old) > 0 (combined total ablation)167 ! (5) v_i (old) = 0; d_v_i_trp > 0 (advection of ice in a free-cell)168 !169 DO jl = 1, jpl170 DO jj = 1, jpj171 DO ji = 1, jpi172 patho_case(ji,jj,jl) = 1173 IF( v_i(ji,jj,jl) .LT. 0.0 ) THEN174 patho_case(ji,jj,jl) = 3175 ENDIF176 IF( ( old_v_i(ji,jj,jl) .LE. epsi10 ) .AND. ( d_v_i_trp(ji,jj,jl) .GT. epsi06 ) ) THEN177 patho_case(ji,jj,jl) = 5 ! advection of ice in an ice-free cell178 ENDIF179 END DO180 END DO181 END DO182 183 137 !-------------------------------------- 184 138 ! 2. Review of all pathological cases 185 139 !-------------------------------------- 186 140 141 ! clem: useless now 187 142 !------------------------------------------- 188 143 ! 2.1) Advection of ice in an ice-free cell 189 144 !------------------------------------------- 190 145 ! should be removed since it is treated after dynamics now 191 192 ! !IF( ln_nicep ) THEN 193 ! WRITE(numout,*) ' limupdate1 - before h correction ' 194 ! WRITE(numout,*) ' a_i : ', a_i(12, 44, 1:jpl) 195 ! WRITE(numout,*) ' at_i : ', at_i(12,44) 196 ! WRITE(numout,*) ' v_i : ', v_i(12, 44, 1:jpl) 197 ! WRITE(numout,*) ' ht_i : ', ht_i(12, 44, 1:jpl) 198 ! !ENDIF 146 ! zhimax = 5._wp 147 ! ! first category 148 ! DO jj = 1, jpj 149 ! DO ji = 1, jpi 150 ! !--- the thickness of such an ice is often out of bounds 151 ! !--- thus we recompute a new area while conserving ice volume 152 ! zat_i_old = SUM( old_a_i(ji,jj,:) ) 153 ! zindb = MAX( 0._wp, SIGN( 1._wp, ABS( d_a_i_trp(ji,jj,1) ) - epsi10 ) ) 154 ! IF( ( ABS( d_v_i_trp(ji,jj,1) ) / MAX( ABS( d_a_i_trp(ji,jj,1) ), epsi10 ) * zindb .GT. zhimax ) & 155 ! & .AND.( ( v_i(ji,jj,1) / MAX( a_i(ji,jj,1), epsi10 ) * zindb ) .GT. zhimax ) & 156 ! & .AND.( zat_i_old .LT. 1.e-6 ) ) THEN ! new line 157 ! ht_i(ji,jj,1) = hi_max(1) * 0.5_wp 158 ! a_i (ji,jj,1) = v_i(ji,jj,1) / ht_i(ji,jj,1) 159 ! ENDIF 160 ! END DO 161 ! END DO 199 162 ! 200 zhimax = .3_wp 201 ! first category 202 DO jj = 1, jpj 203 DO ji = 1, jpi 204 !--- the thickness of such an ice is often out of bounds 205 !--- thus we recompute a new area while conserving ice volume 206 zat_i_old = SUM( old_a_i(ji,jj,:) ) 207 zindb = MAX( 0._wp, SIGN( 1._wp, ABS( d_a_i_trp(ji,jj,1)) - epsi10 ) ) 208 IF ( ( ABS(d_v_i_trp(ji,jj,1))/MAX(ABS(d_a_i_trp(ji,jj,1)),epsi10)*zindb .GT. zhimax) & 209 .AND.( ( v_i(ji,jj,1)/MAX(a_i(ji,jj,1),epsi10)*zindb) .GT. zhimax ) & 210 .AND.( zat_i_old .LT. 1.e-6 ) ) THEN ! new line 211 ht_i(ji,jj,1) = hi_max(1) / 2.0 212 a_i (ji,jj,1) = v_i(ji,jj,1) / ht_i(ji,jj,1) 213 ENDIF 214 END DO 215 END DO 216 217 218 ! !IF( ln_nicep ) THEN 219 ! at_i(:,:) = 0._wp 220 ! DO jl = 1, jpl 221 ! at_i(:,:) = a_i(:,:,jl) + at_i(:,:) 222 ! END DO 223 ! WRITE(numout,*) ' limupdate1 - after h correction 1 ' 224 ! WRITE(numout,*) ' a_i : ', a_i(12, 44, 1:jpl) 225 ! WRITE(numout,*) ' at_i : ', at_i(12,44) 226 ! WRITE(numout,*) ' v_i : ', v_i(12, 44, 1:jpl) 227 ! WRITE(numout,*) ' ht_i : ', ht_i(12, 44, 1:jpl) 228 ! !ENDIF 229 230 zhimax = 1._wp 231 ! other categories 232 DO jl = 2, jpl 233 jm = ice_types(jl) 234 DO jj = 1, jpj 235 DO ji = 1, jpi 236 zindb = MAX( rzero, SIGN( rone, ABS(d_a_i_trp(ji,jj,jl)) - epsi10 ) ) 237 ! this correction is very tricky... sometimes, advection gets wrong i don't know why 238 ! it makes problems when the advected volume and concentration do not seem to be 239 ! related with each other 240 ! the new thickness is sometimes very big! 241 ! and sometimes d_a_i_trp and d_v_i_trp have different sign 242 ! which of course is plausible 243 ! but fuck! it fucks everything up :) 244 IF ( (ABS(d_v_i_trp(ji,jj,jl))/MAX(ABS(d_a_i_trp(ji,jj,jl)),epsi10)*zindb .GT. zhimax) & 245 .AND.(v_i(ji,jj,jl)/MAX(a_i(ji,jj,jl),epsi10)*zindb) .GT. zhimax ) THEN 246 ht_i(ji,jj,jl) = ( hi_max_typ(jl-ice_cat_bounds(jm,1),jm) + hi_max_typ(jl-ice_cat_bounds(jm,1)+1,jm) ) / 2.0 247 a_i (ji,jj,jl) = v_i(ji,jj,jl) / ht_i(ji,jj,jl) 248 ENDIF 249 END DO ! ji 250 END DO !jj 251 END DO !jl 163 ! zhimax = 20._wp 164 ! ! other categories 165 ! DO jl = 2, jpl 166 ! jm = ice_types(jl) 167 ! DO jj = 1, jpj 168 ! DO ji = 1, jpi 169 ! zindb = MAX( rzero, SIGN( rone, ABS( d_a_i_trp(ji,jj,jl) ) - epsi10 ) ) 170 ! ! this correction is very tricky... sometimes, advection gets wrong i don't know why 171 ! ! it makes problems when the advected volume and concentration do not seem to be 172 ! ! related with each other 173 ! ! the new thickness is sometimes very big! 174 ! ! and sometimes d_a_i_trp and d_v_i_trp have different sign 175 ! ! which of course is plausible 176 ! ! but fuck! it fucks everything up :) 177 ! IF ( ( ABS( d_v_i_trp(ji,jj,jl) ) / MAX( ABS( d_a_i_trp(ji,jj,jl) ), epsi10 ) * zindb .GT. zhimax ) & 178 ! & .AND. ( v_i(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi10 ) * zindb ) .GT. zhimax ) THEN 179 ! ht_i(ji,jj,jl) = ( hi_max_typ(jl-ice_cat_bounds(jm,1),jm) + hi_max_typ(jl-ice_cat_bounds(jm,1)+1,jm) ) * 0.5_wp 180 ! a_i (ji,jj,jl) = v_i(ji,jj,jl) / ht_i(ji,jj,jl) 181 ! ENDIF 182 ! END DO ! ji 183 ! END DO !jj 184 ! END DO !jl 252 185 253 186 at_i(:,:) = 0._wp … … 255 188 at_i(:,:) = a_i(:,:,jl) + at_i(:,:) 256 189 END DO 257 258 ! !IF( ln_nicep ) THEN259 ! WRITE(numout,*) ' limupdate1 - after h correction 2 '260 ! WRITE(numout,*) ' a_i : ', a_i(12, 44, 1:jpl)261 ! WRITE(numout,*) ' at_i : ', at_i(12,44)262 ! WRITE(numout,*) ' v_i : ', v_i(12, 44, 1:jpl)263 ! WRITE(numout,*) ' ht_i : ', ht_i(12, 44, 1:jpl)264 ! !ENDIF265 190 266 191 !---------------------------------------------------- … … 285 210 286 211 !switches 287 zindb = MAX( rzero, SIGN( rone, a_i(ji,jj,jl) - epsi 06) )212 zindb = MAX( rzero, SIGN( rone, a_i(ji,jj,jl) - epsi10 ) ) 288 213 !switch = 1 if a_i > 1e-06 and 0 if not 289 zindsn = MAX( rzero, SIGN( rone, v_s(ji,jj,jl) - epsi 06 ) ) !=1 if hs > 1e-6and 0 if not290 zindic = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi 04 ) ) !=1 if hi > 1e-3and 0 if not214 zindsn = MAX( rzero, SIGN( rone, v_s(ji,jj,jl) - epsi10 ) ) !=1 if hs > 1e-10 and 0 if not 215 zindic = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi10 ) ) !=1 if hi > 1e-10 and 0 if not 291 216 ! bug fix 25 avril 2007 292 217 zindb = zindb*zindic … … 332 257 !-------------------------------------------- 333 258 ! if greater than 0, ice concentration cannot be smaller than 1e-10 334 !clem a_i(ji,jj,jl) = zindb * MAX(zindsn, zindic) * MAX( a_i(ji,jj,jl), epsi06 )335 259 a_i(ji,jj,jl) = zindb * a_i(ji,jj,jl) 336 260 … … 352 276 DO jj = 1, jpj 353 277 DO ji = 1, jpi 354 zindic = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi 04) )278 zindic = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi10 ) ) 355 279 e_i(ji,jj,jk,jl)= zindic * ( MIN ( MAX ( 0.0, e_i(ji,jj,jk,jl) ), zbigvalue ) ) 356 280 END DO ! ji … … 370 294 DO ji = 1, jpi 371 295 z_da_ex = MAX( at_i(ji,jj) - amax , 0.0 ) 372 zindb = MAX( rzero, SIGN( rone, at_i(ji,jj) - epsi 06) )296 zindb = MAX( rzero, SIGN( rone, at_i(ji,jj) - epsi10 ) ) 373 297 DO jl = 1, jpl 374 z_da_i = a_i(ji,jj,jl) * z_da_ex / MAX( at_i(ji,jj), epsi 06) * zindb298 z_da_i = a_i(ji,jj,jl) * z_da_ex / MAX( at_i(ji,jj), epsi10 ) * zindb 375 299 a_i(ji,jj,jl) = MAX( 0._wp, a_i(ji,jj,jl) - z_da_i ) 376 300 ! 377 zinda = MAX( rzero, SIGN( rone, a_i(ji,jj,jl) - epsi 06) )378 ht_i(ji,jj,jl) = v_i(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi 06) * zinda301 zinda = MAX( rzero, SIGN( rone, a_i(ji,jj,jl) - epsi10 ) ) 302 ht_i(ji,jj,jl) = v_i(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi10 ) * zinda 379 303 !v_i(ji,jj,jl) = ht_i(ji,jj,jl) * a_i(ji,jj,jl) ! makes ice shrinken but should not be used 380 304 END DO … … 412 336 !clem smv_i(ji,jj,jl) = MAX(MIN((rhoic-rhosn)/rhoic*sss_m(ji,jj),smv_i(ji,jj,jl)),0.1 * v_i(ji,jj,jl) ) 413 337 smv_i(ji,jj,jl) = MAX( MIN( s_i_max * v_i(ji,jj,jl), smv_i(ji,jj,jl) ), s_i_min * v_i(ji,jj,jl) ) 414 i_ice_switch = 1._wp - MAX( 0._wp, SIGN( 1._wp, -v_i(ji,jj,jl) + epsi20) )415 smv_i(ji,jj,jl) = i_ice_switch * smv_i(ji,jj,jl) + s_i_min * ( 1._wp - i_ice_switch ) * v_i(ji,jj,jl)338 i_ice_switch = 1._wp - MAX( 0._wp, SIGN( 1._wp, -v_i(ji,jj,jl) ) ) 339 smv_i(ji,jj,jl) = i_ice_switch * smv_i(ji,jj,jl) !+ s_i_min * ( 1._wp - i_ice_switch ) * v_i(ji,jj,jl) 416 340 END DO ! ji 417 341 END DO ! jj … … 512 436 CALL prt_ctl(tab2d_1=old_oa_i (:,:,jl) , clinfo1= ' lim_update1 : old_oa_i : ') 513 437 CALL prt_ctl(tab2d_1=d_oa_i_trp (:,:,jl) , clinfo1= ' lim_update1 : d_oa_i_trp : ') 514 CALL prt_ctl(tab2d_1=REAL(patho_case(:,:,jl)) , clinfo1= ' lim_update1 : Path. case : ')515 438 516 439 DO jk = 1, nlay_i -
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90
r4099 r4332 45 45 PUBLIC lim_update2 ! routine called by ice_step 46 46 47 REAL(wp) :: epsi06 = 1.e-06_wp ! module constants48 REAL(wp) :: epsi04 = 1.e-04_wp ! - -49 47 REAL(wp) :: epsi10 = 1.e-10_wp ! - - 50 REAL(wp) :: epsi16 = 1.e-16_wp ! - -51 REAL(wp) :: epsi20 = 1.e-20_wp ! - -52 48 REAL(wp) :: rzero = 0._wp ! - - 53 49 REAL(wp) :: rone = 1._wp ! - - … … 102 98 CALL wrk_alloc( jpi,jpj,jpl,zviold, zvsold, zsmvold ) ! clem 103 99 104 !------------------------------------------------------------------------------ 105 ! 1. Update of Global variables | 106 !------------------------------------------------------------------------------ 107 108 !----------------------------- 109 ! Update ice and snow volumes 110 !----------------------------- 111 DO jl = 1, jpl 112 v_i(:,:,jl) = v_i(:,:,jl) + d_v_i_thd(:,:,jl) 113 v_s(:,:,jl) = v_s(:,:,jl) + d_v_s_thd(:,:,jl) 114 END DO 115 116 !--------------------------------------------- 117 ! Ice concentration and ice heat content 118 !--------------------------------------------- 119 a_i (:,:,:) = a_i (:,:,:) + d_a_i_thd(:,:,:) 120 e_i(:,:,:,:) = e_i(:,:,:,:) + d_e_i_thd(:,:,:,:) 121 122 !------------------------------ 123 ! Snow temperature and ice age 124 !------------------------------ 125 e_s (:,:,:,:) = e_s (:,:,:,:) + d_e_s_thd (:,:,:,:) 126 oa_i(:,:,:) = oa_i(:,:,:) + d_oa_i_thd(:,:,:) 127 128 !-------------- 129 ! Ice salinity 130 !-------------- 131 IF( num_sal == 2 ) THEN ! general case 132 smv_i(:,:,:) = smv_i(:,:,:) + d_smv_i_thd(:,:,:) 133 ENDIF 100 !---------------------------------------------------------------------------------------- 101 ! 1. Computation of trend terms 102 !---------------------------------------------------------------------------------------- 103 !- Trend terms 104 d_a_i_thd(:,:,:) = a_i(:,:,:) - old_a_i(:,:,:) 105 d_v_s_thd(:,:,:) = v_s(:,:,:) - old_v_s(:,:,:) 106 d_v_i_thd(:,:,:) = v_i(:,:,:) - old_v_i(:,:,:) 107 d_e_s_thd(:,:,:,:) = e_s(:,:,:,:) - old_e_s(:,:,:,:) 108 d_e_i_thd(:,:,:,:) = e_i(:,:,:,:) - old_e_i(:,:,:,:) 109 !?? d_oa_i_thd(:,:,:) = oa_i (:,:,:) - old_oa_i (:,:,:) 110 d_smv_i_thd(:,:,:) = 0._wp 111 IF( num_sal == 2 ) d_smv_i_thd(:,:,:) = smv_i(:,:,:) - old_smv_i(:,:,:) 112 ! diag only (clem) 113 dv_dt_thd(:,:,:) = d_v_i_thd(:,:,:) * r1_rdtice * rday 134 114 135 115 ! mass and salt flux init (clem) … … 151 131 CALL lim_var_glo2eqv 152 132 153 !---------------------------------154 ! Classify the pathological cases155 !---------------------------------156 ! (1) v_i (new) > 0; d_v_i_thd + v_i(old) > 0 (easy case)157 ! (2) v_i (new) > 0; d_v_i_thd + v_i(old) = 0 (total thermodynamic ablation)158 ! (3) v_i (new) < 0; d_v_i_thd + v_i(old) > 0 (combined total ablation)159 ! (4) v_i (new) < 0; d_v_i_thd + v_i(old) = 0 (total thermodynamic ablation160 ! with negative advection, very pathological )161 !162 DO jl = 1, jpl163 DO jj = 1, jpj164 DO ji = 1, jpi165 patho_case(ji,jj,jl) = 1166 IF( v_i(ji,jj,jl) .GE. 0.0 ) THEN167 IF ( old_v_i(ji,jj,jl) + d_v_i_thd(ji,jj,jl) .LT. epsi10 ) THEN168 patho_case(ji,jj,jl) = 2169 ENDIF170 ELSE171 patho_case(ji,jj,jl) = 3172 IF( old_v_i(ji,jj,jl) + d_v_i_thd(ji,jj,jl) .LT. epsi10 ) THEN173 patho_case(ji,jj,jl) = 4174 ENDIF175 ENDIF176 END DO177 END DO178 END DO179 180 133 !-------------------------------------- 181 134 ! 2. Review of all pathological cases 182 135 !-------------------------------------- 136 137 ! clem: useless now 183 138 !------------------------------------------- 184 139 ! 2.1) Advection of ice in an ice-free cell 185 140 !------------------------------------------- 186 141 ! should be removed since it is treated after dynamics now 187 188 ! !IF( ln_nicep ) THEN 189 ! WRITE(numout,*) ' limupdate2 - before h correction ' 190 ! WRITE(numout,*) ' a_i : ', a_i(12, 44, 1:jpl) 191 ! WRITE(numout,*) ' at_i : ', at_i(12,44) 192 ! WRITE(numout,*) ' v_i : ', v_i(12, 44, 1:jpl) 193 ! WRITE(numout,*) ' ht_i : ', ht_i(12, 44, 1:jpl) 194 ! !ENDIF 195 196 zhimax = .3_wp 197 ! first category 198 DO jj = 1, jpj 199 DO ji = 1, jpi 200 !--- the thickness of such an ice is often out of bounds 201 !--- thus we recompute a new area while conserving ice volume 202 zat_i_old = SUM( old_a_i(ji,jj,:) ) 203 zindb = MAX( 0._wp, SIGN( 1._wp, ABS( d_a_i_thd(ji,jj,1)) - epsi10 ) ) 204 IF ( ( ABS(d_v_i_thd(ji,jj,1))/MAX(ABS(d_a_i_thd(ji,jj,1)),epsi10)*zindb .GT. zhimax) & 205 .AND.( ( v_i(ji,jj,1)/MAX(a_i(ji,jj,1),epsi10)*zindb) .GT. zhimax ) & 206 .AND.( zat_i_old .LT. 1.e-6 ) ) THEN ! new line 207 ht_i(ji,jj,1) = hi_max(1) / 2.0 208 a_i (ji,jj,1) = v_i(ji,jj,1) / ht_i(ji,jj,1) 209 ENDIF 210 END DO 211 END DO 212 213 ! !IF( ln_nicep ) THEN 214 ! at_i(:,:) = 0._wp 215 ! DO jl = 1, jpl 216 ! at_i(:,:) = a_i(:,:,jl) + at_i(:,:) 142 ! zhimax = 5._wp 143 ! ! first category 144 ! DO jj = 1, jpj 145 ! DO ji = 1, jpi 146 ! !--- the thickness of such an ice is often out of bounds 147 ! !--- thus we recompute a new area while conserving ice volume 148 ! zat_i_old = SUM( old_a_i(ji,jj,:) ) 149 ! zindb = MAX( 0._wp, SIGN( 1._wp, ABS( d_a_i_thd(ji,jj,1) ) - epsi10 ) ) 150 ! IF ( ( ABS( d_v_i_thd(ji,jj,1) ) / MAX( ABS( d_a_i_thd(ji,jj,1) ),epsi10 ) * zindb .GT. zhimax ) & 151 ! & .AND. ( ( v_i(ji,jj,1) / MAX( a_i(ji,jj,1), epsi10 ) * zindb ) .GT. zhimax ) & 152 ! & .AND. ( zat_i_old .LT. 1.e-6 ) ) THEN ! new line 153 ! ht_i(ji,jj,1) = hi_max(1) * 0.5_wp 154 ! a_i (ji,jj,1) = v_i(ji,jj,1) / ht_i(ji,jj,1) 155 ! ENDIF 217 156 ! END DO 218 ! WRITE(numout,*) ' limupdate2 - after h correction 1 ' 219 ! WRITE(numout,*) ' a_i : ', a_i(12, 44, 1:jpl) 220 ! WRITE(numout,*) ' at_i : ', at_i(12,44) 221 ! WRITE(numout,*) ' v_i : ', v_i(12, 44, 1:jpl) 222 ! WRITE(numout,*) ' ht_i : ', ht_i(12, 44, 1:jpl) 223 ! !ENDIF 224 ! 225 zhimax = 1._wp 226 ! other categories 227 DO jl = 2, jpl 228 jm = ice_types(jl) 229 DO jj = 1, jpj 230 DO ji = 1, jpi 231 zindb = MAX( rzero, SIGN( rone, ABS(d_a_i_thd(ji,jj,jl)) - epsi10 ) ) 232 ! this correction is very tricky... sometimes, advection gets wrong i don't know why 233 ! it makes problems when the advected volume and concentration do not seem to be 234 ! related with each other 235 ! the new thickness is sometimes very big! 236 ! and sometimes d_a_i_trp and d_v_i_trp have different sign 237 ! which of course is plausible 238 ! but fuck! it fucks everything up :) 239 IF ( (ABS(d_v_i_thd(ji,jj,jl))/MAX(ABS(d_a_i_thd(ji,jj,jl)),epsi10)*zindb .GT. zhimax) & 240 .AND.(v_i(ji,jj,jl)/MAX(a_i(ji,jj,jl),epsi10)*zindb) .GT. zhimax ) THEN 241 ht_i(ji,jj,jl) = ( hi_max_typ(jl-ice_cat_bounds(jm,1),jm) + hi_max_typ(jl-ice_cat_bounds(jm,1)+1,jm) ) / 2.0 242 a_i (ji,jj,jl) = v_i(ji,jj,jl) / ht_i(ji,jj,jl) 243 ENDIF 244 END DO ! ji 245 END DO !jj 246 END DO !jl 157 ! END DO 158 159 ! zhimax = 20._wp 160 ! ! other categories 161 ! DO jl = 2, jpl 162 ! jm = ice_types(jl) 163 ! DO jj = 1, jpj 164 ! DO ji = 1, jpi 165 ! zindb = MAX( rzero, SIGN( rone, ABS( d_a_i_thd(ji,jj,jl)) - epsi10 ) ) 166 ! ! this correction is very tricky... sometimes, advection gets wrong i don't know why 167 ! ! it makes problems when the advected volume and concentration do not seem to be 168 ! ! related with each other 169 ! ! the new thickness is sometimes very big! 170 ! ! and sometimes d_a_i_trp and d_v_i_trp have different sign 171 ! ! which of course is plausible 172 ! ! but fuck! it fucks everything up :) 173 ! IF ( ( ABS( d_v_i_thd(ji,jj,jl) ) / MAX( ABS( d_a_i_thd(ji,jj,jl) ), epsi10 ) * zindb .GT. zhimax ) & 174 ! & .AND. ( v_i(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi10 ) * zindb ) .GT. zhimax ) THEN 175 ! ht_i(ji,jj,jl) = ( hi_max_typ(jl-ice_cat_bounds(jm,1),jm) + hi_max_typ(jl-ice_cat_bounds(jm,1)+1,jm) ) * 0.5_wp 176 ! a_i (ji,jj,jl) = v_i(ji,jj,jl) / ht_i(ji,jj,jl) 177 ! ENDIF 178 ! END DO ! ji 179 ! END DO !jj 180 ! END DO !jl 247 181 248 182 at_i(:,:) = 0._wp … … 250 184 at_i(:,:) = a_i(:,:,jl) + at_i(:,:) 251 185 END DO 252 253 ! !IF( ln_nicep ) THEN254 ! WRITE(numout,*) ' limupdate2 - after h correction 2 '255 ! WRITE(numout,*) ' a_i : ', a_i(12, 44, 1:jpl)256 ! WRITE(numout,*) ' at_i : ', at_i(12,44)257 ! WRITE(numout,*) ' v_i : ', v_i(12, 44, 1:jpl)258 ! WRITE(numout,*) ' ht_i : ', ht_i(12, 44, 1:jpl)259 ! !ENDIF260 186 261 187 !---------------------------------------------------- … … 304 230 zesum = zesum + e_i(ji,jj,jk,jl) 305 231 END DO 306 !clem IF (ind_im .LT. nlay_i ) THEN307 !clem smv_i(ji,jj,jl) = smv_i(ji,jj,jl) / ht_i(ji,jj,jl) * ( ht_i(ji,jj,jl) - REAL(ind_im)*ht_i(ji,jj,jl) / REAL(nlay_i) )308 !clem ENDIF309 232 ht_i(ji,jj,jl) = ht_i(ji,jj,jl) - REAL(ind_im)*ht_i(ji,jj,jl) / REAL(nlay_i) 310 233 v_i(ji,jj,jl) = ht_i(ji,jj,jl) * a_i(ji,jj,jl) … … 368 291 DO ji = 1, jpi 369 292 ! snow energy of melting 370 ze_s = e_s(ji,jj,1,jl) * unit_fac / area(ji,jj) / MAX( v_s(ji,jj,jl), 1.0e-6 ) ! snow energy of melting 293 zinda = MAX( 0._wp, SIGN( 1._wp, v_s(ji,jj,jl) - epsi10 ) ) 294 ze_s = zinda * e_s(ji,jj,1,jl) * unit_fac / area(ji,jj) / MAX( v_s(ji,jj,jl), epsi10 ) ! snow energy of melting 371 295 372 296 ! If snow energy of melting smaller then Lf … … 401 325 402 326 !switches 403 zindb = MAX( rzero, SIGN( rone, a_i(ji,jj,jl) - epsi 06) )327 zindb = MAX( rzero, SIGN( rone, a_i(ji,jj,jl) - epsi10 ) ) 404 328 !switch = 1 if a_i > 1e-06 and 0 if not 405 zindsn = MAX( rzero, SIGN( rone, v_s(ji,jj,jl) - epsi 06 ) ) !=1 if hs > 1e-6and 0 if not406 zindic = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi 04 ) ) !=1 if hi > 1e-3and 0 if not329 zindsn = MAX( rzero, SIGN( rone, v_s(ji,jj,jl) - epsi10 ) ) !=1 if hs > 1e-10 and 0 if not 330 zindic = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi10 ) ) !=1 if hi > 1e-10 and 0 if not 407 331 ! bug fix 25 avril 2007 408 332 zindb = zindb*zindic … … 444 368 !--- 2.7 Correction to ice concentrations 445 369 !-------------------------------------------- 446 !clem a_i(ji,jj,jl) = zindb * MAX(zindsn, zindic) * MAX( a_i(ji,jj,jl), epsi06 )447 370 a_i(ji,jj,jl) = zindb * a_i(ji,jj,jl) 448 371 … … 464 387 DO jj = 1, jpj 465 388 DO ji = 1, jpi 466 zindic = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi 04) )389 zindic = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi10 ) ) 467 390 e_i(ji,jj,jk,jl)= zindic * ( MIN ( MAX ( 0.0, e_i(ji,jj,jk,jl) ), zbigvalue ) ) 468 391 END DO ! ji … … 478 401 !--- 2.12 Constrain the thickness of the smallest category above 5 cm 479 402 !---------------------------------------------------------------------- 480 zindb = MAX( rzero, SIGN( rone, a_i(ji,jj,jl) - epsi 06) )481 ht_i(ji,jj,jl) = zindb*v_i(ji,jj,jl)/MAX(a_i(ji,jj,jl), epsi 06)482 zh = MAX( rone , zindb * hiclim / MAX( ht_i(ji,jj,jl) , epsi 20 ) )403 zindb = MAX( rzero, SIGN( rone, a_i(ji,jj,jl) - epsi10 ) ) 404 ht_i(ji,jj,jl) = zindb*v_i(ji,jj,jl)/MAX(a_i(ji,jj,jl), epsi10) 405 zh = MAX( rone , zindb * hiclim / MAX( ht_i(ji,jj,jl) , epsi10 ) ) 483 406 ht_s(ji,jj,jl) = ht_s(ji,jj,jl)* zh 484 407 ht_i(ji,jj,jl) = ht_i(ji,jj,jl)* zh … … 502 425 DO ji = 1, jpi 503 426 z_da_ex = MAX( at_i(ji,jj) - amax , 0.0 ) 504 zindb = MAX( rzero, SIGN( rone, at_i(ji,jj) - epsi 06) )427 zindb = MAX( rzero, SIGN( rone, at_i(ji,jj) - epsi10 ) ) 505 428 DO jl = 1, jpl 506 z_da_i = a_i(ji,jj,jl) * z_da_ex / MAX( at_i(ji,jj), epsi 06) * zindb429 z_da_i = a_i(ji,jj,jl) * z_da_ex / MAX( at_i(ji,jj), epsi10 ) * zindb 507 430 a_i(ji,jj,jl) = MAX( 0._wp, a_i(ji,jj,jl) - z_da_i ) 508 431 ! 509 zinda = MAX( rzero, SIGN( rone, a_i(ji,jj,jl) - epsi 06) )510 ht_i(ji,jj,jl) = v_i(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi 06) * zinda432 zinda = MAX( rzero, SIGN( rone, a_i(ji,jj,jl) - epsi10 ) ) 433 ht_i(ji,jj,jl) = v_i(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi10 ) * zinda 511 434 !v_i(ji,jj,jl) = ht_i(ji,jj,jl) * a_i(ji,jj,jl) ! makes ice shrinken but should not be used 512 435 END DO … … 542 465 !clem smv_i(ji,jj,jl) = MAX(MIN((rhoic-rhosn)/rhoic*sss_m(ji,jj),smv_i(ji,jj,jl)),0.1 * v_i(ji,jj,jl) ) 543 466 smv_i(ji,jj,jl) = MAX( MIN( s_i_max * v_i(ji,jj,jl), smv_i(ji,jj,jl) ), s_i_min * v_i(ji,jj,jl) ) 544 i_ice_switch = 1._wp - MAX( 0._wp, SIGN( 1._wp, -v_i(ji,jj,jl) + epsi20) )545 smv_i(ji,jj,jl) = i_ice_switch * smv_i(ji,jj,jl) + s_i_min * ( 1._wp - i_ice_switch ) * v_i(ji,jj,jl)467 i_ice_switch = 1._wp - MAX( 0._wp, SIGN( 1._wp, -v_i(ji,jj,jl) ) ) 468 smv_i(ji,jj,jl) = i_ice_switch * smv_i(ji,jj,jl) !+ s_i_min * ( 1._wp - i_ice_switch ) * v_i(ji,jj,jl) 546 469 END DO ! ji 547 470 END DO ! jj … … 663 586 CALL prt_ctl(tab2d_1=old_oa_i (:,:,jl) , clinfo1= ' lim_update2 : old_oa_i : ') 664 587 CALL prt_ctl(tab2d_1=d_oa_i_thd (:,:,jl) , clinfo1= ' lim_update2 : d_oa_i_thd : ') 665 CALL prt_ctl(tab2d_1=REAL(patho_case(:,:,jl)) , clinfo1= ' lim_update2 : Path. case : ')666 588 667 589 DO jk = 1, nlay_i -
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90
r4220 r4332 66 66 PUBLIC lim_var_salprof1d ! 67 67 68 REAL(wp) :: epsi20 = 1.e-20_wp ! module constants69 REAL(wp) :: epsi16 = 1.e-16_wp ! - -70 REAL(wp) :: epsi13 = 1.e-13_wp ! - -71 68 REAL(wp) :: epsi10 = 1.e-10_wp ! - - 72 REAL(wp) :: epsi06 = 1.e-06_wp ! - -73 69 REAL(wp) :: zzero = 0.e0 ! - - 74 70 REAL(wp) :: zone = 1.e0 ! - - … … 117 113 at_i(ji,jj) = at_i(ji,jj) + a_i(ji,jj,jl) ! ice concentration 118 114 ! 119 zinda = MAX( zzero , SIGN( zone , at_i(ji,jj) - epsi1 6) )120 icethi(ji,jj) = vt_i(ji,jj) / MAX( at_i(ji,jj) , epsi1 6) * zinda ! ice thickness115 zinda = MAX( zzero , SIGN( zone , at_i(ji,jj) - epsi10 ) ) 116 icethi(ji,jj) = vt_i(ji,jj) / MAX( at_i(ji,jj) , epsi10 ) * zinda ! ice thickness 121 117 END DO 122 118 END DO … … 138 134 DO jj = 1, jpj 139 135 DO ji = 1, jpi 140 zinda = MAX( zzero , SIGN( zone , vt_i(ji,jj) - epsi1 6) )141 zindb = MAX( zzero , SIGN( zone , at_i(ji,jj) - epsi1 6) )136 zinda = MAX( zzero , SIGN( zone , vt_i(ji,jj) - epsi10 ) ) 137 zindb = MAX( zzero , SIGN( zone , at_i(ji,jj) - epsi10 ) ) 142 138 et_s(ji,jj) = et_s(ji,jj) + e_s(ji,jj,1,jl) ! snow heat content 143 smt_i(ji,jj) = smt_i(ji,jj) + smv_i(ji,jj,jl) / MAX( vt_i(ji,jj) , epsi1 6) * zinda ! ice salinity144 ot_i(ji,jj) = ot_i(ji,jj) + oa_i(ji,jj,jl) / MAX( at_i(ji,jj) , epsi1 6) * zindb ! ice age139 smt_i(ji,jj) = smt_i(ji,jj) + smv_i(ji,jj,jl) / MAX( vt_i(ji,jj) , epsi10 ) * zinda ! ice salinity 140 ot_i(ji,jj) = ot_i(ji,jj) + oa_i(ji,jj,jl) / MAX( at_i(ji,jj) , epsi10 ) * zindb ! ice age 145 141 END DO 146 142 END DO … … 209 205 DO ji = 1, jpi 210 206 ! ! Energy of melting q(S,T) [J.m-3] 211 zq_i = e_i(ji,jj,jk,jl) / area(ji,jj) / MAX( v_i(ji,jj,jl) , epsi 06) * REAL(nlay_i,wp)212 zindb = 1.0 - MAX( 0.0 , SIGN( 1.0 , - v_i(ji,jj,jl) ) ) ! zindb = 0 if no ice and 1 if yes207 zq_i = e_i(ji,jj,jk,jl) / area(ji,jj) / MAX( v_i(ji,jj,jl) , epsi10 ) * REAL(nlay_i,wp) 208 zindb = 1.0 - MAX( 0.0 , SIGN( 1.0 , - v_i(ji,jj,jl) + epsi10 ) ) ! zindb = 0 if no ice and 1 if yes 213 209 zq_i = zq_i * unit_fac * zindb !convert units 214 210 ztmelts = -tmut * s_i(ji,jj,jk,jl) + rtt ! Ice layer melt temperature … … 235 231 DO ji = 1, jpi 236 232 !Energy of melting q(S,T) [J.m-3] 237 zq_s = e_s(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_s(ji,jj,jl) , epsi 06) ) * REAL(nlay_s,wp)238 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - v_s(ji,jj,jl) ) ) ! zindb = 0 if no ice and 1 if yes233 zq_s = e_s(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_s(ji,jj,jl) , epsi10 ) ) * REAL(nlay_s,wp) 234 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - v_s(ji,jj,jl) + epsi10 ) ) ! zindb = 0 if no ice and 1 if yes 239 235 zq_s = zq_s * unit_fac * zindb ! convert units 240 236 ! … … 341 337 zind01 = ( 1._wp - zind0 ) * MAX( 0._wp , SIGN( 1._wp , s_i_1 - sm_i(ji,jj,jl) ) ) 342 338 ! If 2.sm_i GE sss_m then zindbal = 1 339 ! this is to force a constant salinity profile in the Baltic Sea 343 340 zindbal = MAX( 0._wp , SIGN( 1._wp , 2._wp * sm_i(ji,jj,jl) - sss_m(ji,jj) ) ) 344 341 zalpha(ji,jj,jl) = zind0 + zind01 * ( sm_i(ji,jj,jl) * dummy_fac0 + dummy_fac1 ) … … 434 431 DO jj = 1, jpj 435 432 DO ji = 1, jpi 436 zinda = ( 1._wp - MAX( 0._wp , SIGN( 1._wp , (t_i(ji,jj,jk,jl) - rtt) + epsi1 6) ) )437 zindb = ( 1._wp - MAX( 0._wp , SIGN( 1._wp , - vt_i(ji,jj) + epsi1 6) ) )438 zbvi = - zinda * tmut * s_i(ji,jj,jk,jl) / MIN( t_i(ji,jj,jk,jl) - rtt, - epsi1 6) &433 zinda = ( 1._wp - MAX( 0._wp , SIGN( 1._wp , (t_i(ji,jj,jk,jl) - rtt) + epsi10 ) ) ) 434 zindb = ( 1._wp - MAX( 0._wp , SIGN( 1._wp , - vt_i(ji,jj) + epsi10 ) ) ) 435 zbvi = - zinda * tmut * s_i(ji,jj,jk,jl) / MIN( t_i(ji,jj,jk,jl) - rtt, - epsi10 ) & 439 436 & * v_i(ji,jj,jl) / REAL(nlay_i,wp) 440 bv_i(ji,jj) = bv_i(ji,jj) + zindb * zbvi / MAX( vt_i(ji,jj) , epsi1 6)437 bv_i(ji,jj) = bv_i(ji,jj) + zindb * zbvi / MAX( vt_i(ji,jj) , epsi10 ) 441 438 END DO 442 439 END DO … … 498 495 zind01 = ( 1._wp - zind0 ) * MAX( 0._wp , SIGN( 1._wp , s_i_1 - sm_i_b(ji) ) ) 499 496 ! if 2.sm_i GE sss_m then zindbal = 1 497 ! this is to force a constant salinity profile in the Baltic Sea 500 498 zindbal = MAX( 0._wp , SIGN( 1._wp , 2._wp * sm_i_b(ji) - sss_m(ii,ij) ) ) 501 499 ! -
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90
r4045 r4332 52 52 INTEGER , DIMENSION(jpnoumax) :: nc , nca ! switch for saving field ( = 1 ) or not ( = 0 ) 53 53 54 REAL(wp) :: epsi06 = 1 e-6_wp54 REAL(wp) :: epsi06 = 1.e-6_wp 55 55 REAL(wp) :: zzero = 0._wp 56 56 REAL(wp) :: zone = 1._wp … … 192 192 193 193 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 194 IF( ln_nicep ) THEN195 WRITE(numout,*)196 WRITE(numout,*) 'lim_wri : write ice outputs in NetCDF files at time : ', nyear, nmonth, nday, numit197 WRITE(numout,*) '~~~~~~~ '198 WRITE(numout,*) ' kindic = ', kindic199 ENDIF194 !IF( ln_nicep ) THEN 195 ! WRITE(numout,*) 196 ! WRITE(numout,*) 'lim_wri : write ice outputs in NetCDF files at time : ', nyear, nmonth, nday, numit 197 ! WRITE(numout,*) '~~~~~~~ ' 198 ! WRITE(numout,*) ' kindic = ', kindic 199 !ENDIF 200 200 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 201 201 … … 564 564 CALL histdef( kid, "iicevolu", "Ice volume" , "m" , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 565 565 CALL histdef( kid, "iicedive", "Ice divergence" , "10-8s-1", jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 566 CALL histdef( kid, "iicebopr", "Ice bottom production" , "m/s" , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 567 CALL histdef( kid, "iicedypr", "Ice dynamic production" , "m/s" , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 568 CALL histdef( kid, "iicelapr", "Ice open water prod" , "m/s" , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 569 CALL histdef( kid, "iicesipr", "Snow ice production " , "m/s" , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 570 CALL histdef( kid, "iicerepr", "Ice prod from limupdate" , "m/s" , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 571 CALL histdef( kid, "iicebome", "Ice bottom melt" , "m/s" , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 572 CALL histdef( kid, "iicesume", "Ice surface melt" , "m/s" , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 573 CALL histdef( kid, "iisfxthd", "Salt flux from thermo" , "" , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 574 CALL histdef( kid, "iisfxmec", "Salt flux from dynmics" , "" , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 575 CALL histdef( kid, "iisfxres", "Salt flux from limupdate", "" , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 576 566 577 567 578 !CALL histdef( kid, "iice_itd", "Ice concentration by cat", "%" , jpi, jpj, kh_i, jpl, 1, jpl, -99, 32, "inst(x)", rdt, rdt ) … … 586 597 CALL histwrite( kid, "iicedive", kt, divu_i*1.0e8 , jpi*jpj, (/1/) ) 587 598 599 CALL histwrite( kid, "iicebopr", kt, diag_bot_gr , jpi*jpj, (/1/) ) 600 CALL histwrite( kid, "iicedypr", kt, diag_dyn_gr , jpi*jpj, (/1/) ) 601 CALL histwrite( kid, "iicelapr", kt, diag_lat_gr , jpi*jpj, (/1/) ) 602 CALL histwrite( kid, "iicesipr", kt, diag_sni_gr , jpi*jpj, (/1/) ) 603 CALL histwrite( kid, "iicerepr", kt, diag_res_pr , jpi*jpj, (/1/) ) 604 CALL histwrite( kid, "iicebome", kt, diag_bot_me , jpi*jpj, (/1/) ) 605 CALL histwrite( kid, "iicesume", kt, diag_sur_me , jpi*jpj, (/1/) ) 606 CALL histwrite( kid, "iisfxthd", kt, sfx_thd , jpi*jpj, (/1/) ) 607 CALL histwrite( kid, "iisfxmec", kt, sfx_mec , jpi*jpj, (/1/) ) 608 CALL histwrite( kid, "iisfxres", kt, sfx_res , jpi*jpj, (/1/) ) 609 588 610 !CALL histwrite( kid, "iice_itd", kt, a_i , jpi*jpj*jpl, (/1/) ) ! area 589 611 !CALL histwrite( kid, "iice_hid", kt, ht_i , jpi*jpj*jpl, (/1/) ) ! thickness
Note: See TracChangeset
for help on using the changeset viewer.