Changeset 4332 for branches/2013/dev_r4028_CNRS_LIM3
- Timestamp:
- 2013-12-11T15:38:42+01:00 (10 years ago)
- Location:
- branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM
- Files:
-
- 1 deleted
- 30 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/CONFIG/cfg.txt
r4099 r4332 7 7 ORCA2_LIM_PISCES OPA_SRC LIM_SRC_2 NST_SRC TOP_SRC 8 8 AMM12 OPA_SRC 9 ORCA2_LIM3 OPA_SRC LIM_SRC_3 NST_SRC 9 10 ORCA2_LIM OPA_SRC LIM_SRC_2 NST_SRC 10 ORCA2_LIM3 OPA_SRC LIM_SRC_3 NST_SRC -
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 -
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/BDY/bdy_oce.F90
r3651 r4332 8 8 !! 3.3 ! 2010-09 (D. Storkey) add ice boundary conditions 9 9 !! 3.4 ! 2011 (D. Storkey) rewrite in preparation for OBC-BDY merge 10 !! - ! 2012-01 (C. Rousset) add ice boundary conditions for lim3 10 11 !!---------------------------------------------------------------------- 11 12 #if defined key_bdy … … 45 46 REAL, POINTER, DIMENSION(:) :: hicif 46 47 REAL, POINTER, DIMENSION(:) :: hsnif 48 #elif defined key_lim3 49 REAL, POINTER, DIMENSION(:,:) :: a_i !: now ice leads fraction climatology 50 REAL, POINTER, DIMENSION(:,:) :: ht_i !: Now ice thickness climatology 51 REAL, POINTER, DIMENSION(:,:) :: ht_s !: now snow thickness 47 52 #endif 48 53 END TYPE OBC_DATA … … 78 83 REAL, DIMENSION(jp_bdy) :: rn_time_dmp !: Damping time scale in days 79 84 80 #if defined key_lim281 INTEGER, DIMENSION(jp_bdy) :: nn_ice_lim 2! Choice of boundary condition for sea ice variables82 INTEGER, DIMENSION(jp_bdy) :: nn_ice_lim 2_dta!: = 0 use the initial state as bdy dta ;85 #if ( defined key_lim2 || defined key_lim3 ) 86 INTEGER, DIMENSION(jp_bdy) :: nn_ice_lim ! Choice of boundary condition for sea ice variables 87 INTEGER, DIMENSION(jp_bdy) :: nn_ice_lim_dta !: = 0 use the initial state as bdy dta ; 83 88 !: = 1 read it in a NetCDF file 84 89 #endif -
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90
r3909 r4332 11 11 !! 3.3 ! 2010-09 (D.Storkey) add ice boundary conditions 12 12 !! 3.4 ! 2011 (D. Storkey) rewrite in preparation for OBC-BDY merge 13 !! - ! 2012-01 (C. Rousset) add ice boundary conditions for lim3 13 14 !!---------------------------------------------------------------------- 14 15 #if defined key_bdy … … 31 32 #if defined key_lim2 32 33 USE ice_2 34 #elif defined key_lim3 35 USE par_ice 36 USE ice 37 USE limcat_1D ! redistribute ice input into categories 33 38 #endif 34 39 USE sbcapr … … 49 54 50 55 TYPE(MAP_POINTER), ALLOCATABLE, DIMENSION(:) :: nbmap_ptr ! array of pointers to nbmap 56 57 #if defined key_lim3 58 LOGICAL :: ll_bdylim3 ! determine whether ice input is lim2 (F) or lim3 (T) type 59 INTEGER :: jfld_hti, jfld_hts, jfld_ai ! indices of ice thickness, snow thickness and concentration in bf structure 60 #endif 51 61 52 62 # include "domzgr_substitute.h90" … … 77 87 ! etc. 78 88 !! 79 INTEGER :: ib_bdy, jfld, jstart, jend, ib, ii, ij, ik, igrd ! local indices89 INTEGER :: ib_bdy, jfld, jstart, jend, ib, ii, ij, ik, igrd, jl ! local indices 80 90 INTEGER, DIMENSION(jpbgrd) :: ilen1 81 91 INTEGER, POINTER, DIMENSION(:) :: nblen, nblenrim ! short cuts … … 164 174 165 175 #if defined key_lim2 166 IF( nn_ice_lim 2(ib_bdy) .gt. 0 .and. nn_ice_lim2_dta(ib_bdy) .eq. 0 ) THEN176 IF( nn_ice_lim(ib_bdy) .gt. 0 .and. nn_ice_lim_dta(ib_bdy) .eq. 0 ) THEN 167 177 ilen1(:) = nblen(:) 168 178 igrd = 1 ! Everything is at T-points here … … 170 180 ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 171 181 ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 172 dta_bdy(ib_bdy)%frld (ib) =frld(ii,ij) * tmask(ii,ij,1)182 dta_bdy(ib_bdy)%frld (ib) = frld(ii,ij) * tmask(ii,ij,1) 173 183 dta_bdy(ib_bdy)%hicif(ib) = hicif(ii,ij) * tmask(ii,ij,1) 174 184 dta_bdy(ib_bdy)%hsnif(ib) = hsnif(ii,ij) * tmask(ii,ij,1) 175 185 END DO 186 ENDIF 187 #elif defined key_lim3 188 IF( nn_ice_lim(ib_bdy) .gt. 0 .and. nn_ice_lim_dta(ib_bdy) .eq. 0 ) THEN 189 ilen1(:) = nblen(:) 190 igrd = 1 ! Everything is at T-points here 191 DO jl = 1, jpl 192 DO ib = 1, ilen1(igrd) 193 ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 194 ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 195 dta_bdy(ib_bdy)%a_i (ib,jl) = a_i(ii,ij,jl) * tmask(ii,ij,1) 196 dta_bdy(ib_bdy)%ht_i(ib,jl) = ht_i(ii,ij,jl) * tmask(ii,ij,1) 197 dta_bdy(ib_bdy)%ht_s(ib,jl) = ht_s(ii,ij,jl) * tmask(ii,ij,1) 198 END DO 199 END DO 176 200 ENDIF 177 201 #endif … … 319 343 ENDIF 320 344 ENDIF 345 #if defined key_lim3 346 IF( .NOT. ll_bdylim3 .AND. nn_ice_lim(ib_bdy) > 0 .AND. nn_ice_lim_dta(ib_bdy) == 1 ) THEN ! bdy ice input (case input is lim2 type) 347 CALL lim_cat_1D ( bf(jfld_hti)%fnow(:,1,1), bf(jfld_hts)%fnow(:,1,1), bf(jfld_ai)%fnow(:,1,1), & 348 & dta_bdy(ib_bdy)%ht_i, dta_bdy(ib_bdy)%ht_s, dta_bdy(ib_bdy)%a_i ) 349 ENDIF 350 #endif 351 321 352 ENDIF 322 353 jstart = jend+1 … … 366 397 INTEGER, ALLOCATABLE, DIMENSION(:) :: igrid ! index for grid type (1,2,3 = T,U,V) 367 398 INTEGER, POINTER, DIMENSION(:) :: nblen, nblenrim ! short cuts 399 #if defined key_lim3 400 INTEGER, DIMENSION(3) :: zdimsz ! number of elements in each of the 4 dimensions (i.e. i,j,t,ice-cat) for an array 401 INTEGER :: zndims ! number of dimensions in an array (i.e. 3 = wo ice cat; 4 = w ice cat) 402 INTEGER :: inum,id1 ! local integer 403 #endif 368 404 TYPE(FLD_N), ALLOCATABLE, DIMENSION(:) :: blf_i ! array of namelist information structures 369 405 TYPE(FLD_N) :: bn_tem, bn_sal, bn_u3d, bn_v3d ! 370 406 TYPE(FLD_N) :: bn_ssh, bn_u2d, bn_v2d ! informations about the fields to be read 371 407 #if defined key_lim2 372 TYPE(FLD_N) :: bn_frld, bn_hicif, bn_hsnif ! 408 TYPE(FLD_N) :: bn_frld, bn_hicif, bn_hsnif 409 #elif defined key_lim3 410 TYPE(FLD_N) :: bn_a_i, bn_ht_i, bn_ht_s 373 411 #endif 374 412 NAMELIST/nambdy_dta/ cn_dir, bn_tem, bn_sal, bn_u3d, bn_v3d, bn_ssh, bn_u2d, bn_v2d 375 413 #if defined key_lim2 376 414 NAMELIST/nambdy_dta/ bn_frld, bn_hicif, bn_hsnif 415 #elif defined key_lim3 416 NAMELIST/nambdy_dta/ bn_a_i, bn_ht_i, bn_ht_s 377 417 #endif 378 418 NAMELIST/nambdy_dta/ ln_full_vel … … 388 428 ! Set nn_dta 389 429 DO ib_bdy = 1, nb_bdy 390 nn_dta(ib_bdy) = MAX( 391 392 393 #if defined key_lim2394 ,nn_ice_lim2_dta(ib_bdy) &430 nn_dta(ib_bdy) = MAX( nn_dyn2d_dta(ib_bdy) & 431 ,nn_dyn3d_dta(ib_bdy) & 432 ,nn_tra_dta(ib_bdy) & 433 #if ( defined key_lim2 || defined key_lim3 ) 434 ,nn_ice_lim_dta(ib_bdy) & 395 435 #endif 396 436 ) … … 412 452 nb_bdy_fld(ib_bdy) = nb_bdy_fld(ib_bdy) + 2 413 453 ENDIF 414 #if defined key_lim2415 IF( nn_ice_lim 2(ib_bdy) .gt. 0 .and. nn_ice_lim2_dta(ib_bdy) .eq. 1 ) THEN454 #if ( defined key_lim2 || defined key_lim3 ) 455 IF( nn_ice_lim(ib_bdy) .gt. 0 .and. nn_ice_lim_dta(ib_bdy) .eq. 1 ) THEN 416 456 nb_bdy_fld(ib_bdy) = nb_bdy_fld(ib_bdy) + 3 417 457 ENDIF … … 448 488 ln_full_vel = .false. 449 489 ! ... default values (NB: frequency positive => hours, negative => months) 450 ! ! file ! frequency ! variable! time intep ! clim ! 'yearly' or ! weights ! rotation !451 ! ! name ! hours ! name ! (T/F)! (T/F) ! 'monthly' ! filename ! pairs !452 bn_ssh = FLD_N( 'bdy_ssh' , 24 , 'sossheig' , .false., .false. , 'yearly' , '' , '' )453 bn_u2d = FLD_N( 'bdy_vel2d_u' , 24 , 'vobtcrtx' , .false., .false. , 'yearly' , '' , '' )454 bn_v2d = FLD_N( 'bdy_vel2d_v' , 24 , 'vobtcrty' , .false., .false. , 'yearly' , '' , '' )455 bn_u3d = FLD_N( 'bdy_vel3d_u' , 24 , 'vozocrtx' , .false., .false. , 'yearly' , '' , '' )456 bn_v3d = FLD_N( 'bdy_vel3d_v' , 24 , 'vomecrty' , .false., .false. , 'yearly' , '' , '' )457 bn_tem = FLD_N( 'bdy_tem' , 24 , 'votemper' , .false., .false. , 'yearly' , '' , '' )458 bn_sal = FLD_N( 'bdy_sal' , 24 , 'vosaline' , .false., .false. , 'yearly' , '' , '' )490 ! ! file ! frequency ! variable ! time intep ! clim ! 'yearly' or ! weights ! rotation ! 491 ! ! name ! hours ! name ! (T/F) ! (T/F) ! 'monthly' ! filename ! pairs ! 492 bn_ssh = FLD_N( 'bdy_ssh' , 24 , 'sossheig' , .false. , .false. , 'yearly' , '' , '' ) 493 bn_u2d = FLD_N( 'bdy_vel2d_u' , 24 , 'vobtcrtx' , .false. , .false. , 'yearly' , '' , '' ) 494 bn_v2d = FLD_N( 'bdy_vel2d_v' , 24 , 'vobtcrty' , .false. , .false. , 'yearly' , '' , '' ) 495 bn_u3d = FLD_N( 'bdy_vel3d_u' , 24 , 'vozocrtx' , .false. , .false. , 'yearly' , '' , '' ) 496 bn_v3d = FLD_N( 'bdy_vel3d_v' , 24 , 'vomecrty' , .false. , .false. , 'yearly' , '' , '' ) 497 bn_tem = FLD_N( 'bdy_tem' , 24 , 'votemper' , .false. , .false. , 'yearly' , '' , '' ) 498 bn_sal = FLD_N( 'bdy_sal' , 24 , 'vosaline' , .false. , .false. , 'yearly' , '' , '' ) 459 499 #if defined key_lim2 460 bn_frld = FLD_N( 'bdy_frld' , 24 , 'ildsconc' , .false. , .false. , 'yearly' , '' , '' ) 461 bn_hicif = FLD_N( 'bdy_hicif' , 24 , 'iicethic' , .false. , .false. , 'yearly' , '' , '' ) 462 bn_hsnif = FLD_N( 'bdy_hsnif' , 24 , 'isnothic' , .false. , .false. , 'yearly' , '' , '' ) 500 bn_frld = FLD_N( 'bdy_frld' , 24 , 'ildsconc' , .false. , .false. , 'yearly' , '' , '' ) 501 bn_hicif = FLD_N( 'bdy_hicif' , 24 , 'iicethic' , .false. , .false. , 'yearly' , '' , '' ) 502 bn_hsnif = FLD_N( 'bdy_hsnif' , 24 , 'isnothic' , .false. , .false. , 'yearly' , '' , '' ) 503 #elif defined key_lim3 504 bn_a_i = FLD_N( 'bdy_a_i' , 24 , 'ileadfra' , .false. , .false. , 'yearly' , '' , '' ) 505 bn_ht_i = FLD_N( 'bdy_ht_i' , 24 , 'iicethic' , .false. , .false. , 'yearly' , '' , '' ) 506 bn_ht_s = FLD_N( 'bdy_ht_s' , 24 , 'isnowthi' , .false. , .false. , 'yearly' , '' , '' ) 463 507 #endif 464 508 … … 545 589 #if defined key_lim2 546 590 ! sea ice 547 IF( nn_ice_lim 2(ib_bdy) .gt. 0 .and. nn_ice_lim2_dta(ib_bdy) .eq. 1 ) THEN591 IF( nn_ice_lim(ib_bdy) .gt. 0 .and. nn_ice_lim_dta(ib_bdy) .eq. 1 ) THEN 548 592 549 593 jfld = jfld + 1 … … 567 611 ilen1(jfld) = nblen(igrid(jfld)) 568 612 ilen3(jfld) = 1 613 614 ENDIF 615 #elif defined key_lim3 616 ! sea ice 617 IF( nn_ice_lim(ib_bdy) .gt. 0 .and. nn_ice_lim_dta(ib_bdy) .eq. 1 ) THEN 618 619 ! Test for types of ice input (lim2 or lim3) 620 CALL iom_open ( bn_a_i%clname, inum ) 621 id1 = iom_varid ( inum, bn_a_i%clvar, kdimsz=zdimsz, kndims=zndims, ldstop = .FALSE. ) 622 CALL iom_close ( inum ) 623 !CALL fld_clopn ( bn_a_i, nyear, nmonth, nday, ldstop=.TRUE. ) 624 !CALL iom_open ( bn_a_i%clname, inum ) 625 !id1 = iom_varid ( bn_a_i%num, bn_a_i%clvar, kdimsz=zdimsz, kndims=zndims, ldstop = .FALSE. ) 626 IF ( zndims == 4 ) THEN 627 ll_bdylim3 = .TRUE. ! lim3 input 628 ELSE 629 ll_bdylim3 = .FALSE. ! lim2 input 630 ENDIF 631 ! End test 632 633 jfld = jfld + 1 634 blf_i(jfld) = bn_a_i 635 ibdy(jfld) = ib_bdy 636 igrid(jfld) = 1 637 ilen1(jfld) = nblen(igrid(jfld)) 638 IF ( ll_bdylim3 ) THEN ; ilen3(jfld)=jpl ; ELSE ; ilen3(jfld)=1 ; ENDIF 639 640 jfld = jfld + 1 641 blf_i(jfld) = bn_ht_i 642 ibdy(jfld) = ib_bdy 643 igrid(jfld) = 1 644 ilen1(jfld) = nblen(igrid(jfld)) 645 IF ( ll_bdylim3 ) THEN ; ilen3(jfld)=jpl ; ELSE ; ilen3(jfld)=1 ; ENDIF 646 647 jfld = jfld + 1 648 blf_i(jfld) = bn_ht_s 649 ibdy(jfld) = ib_bdy 650 igrid(jfld) = 1 651 ilen1(jfld) = nblen(igrid(jfld)) 652 IF ( ll_bdylim3 ) THEN ; ilen3(jfld)=jpl ; ELSE ; ilen3(jfld)=1 ; ENDIF 569 653 570 654 ENDIF … … 662 746 663 747 #if defined key_lim2 664 IF (nn_ice_lim 2(ib_bdy) .gt. 0) THEN665 IF( nn_ice_lim 2_dta(ib_bdy) .eq. 0 ) THEN748 IF (nn_ice_lim(ib_bdy) .gt. 0) THEN 749 IF( nn_ice_lim_dta(ib_bdy) .eq. 0 ) THEN 666 750 ilen0(1:3) = nblen(1:3) 667 751 ALLOCATE( dta_bdy(ib_bdy)%frld(ilen0(1)) ) … … 675 759 jfld = jfld + 1 676 760 dta_bdy(ib_bdy)%hsnif => bf(jfld)%fnow(:,1,1) 761 ENDIF 762 ENDIF 763 #elif defined key_lim3 764 IF (nn_ice_lim(ib_bdy) .gt. 0) THEN 765 ilen0(1:3) = nblen(1:3) 766 IF( nn_ice_lim_dta(ib_bdy) .eq. 0 ) THEN 767 ALLOCATE( dta_bdy(ib_bdy)%a_i (ilen0(1),jpl) ) 768 ALLOCATE( dta_bdy(ib_bdy)%ht_i(ilen0(1),jpl) ) 769 ALLOCATE( dta_bdy(ib_bdy)%ht_s(ilen0(1),jpl) ) 770 ELSE 771 IF ( ll_bdylim3 ) THEN ! case input is lim3 type 772 jfld = jfld + 1 773 dta_bdy(ib_bdy)%a_i => bf(jfld)%fnow(:,1,:) 774 jfld = jfld + 1 775 dta_bdy(ib_bdy)%ht_i => bf(jfld)%fnow(:,1,:) 776 jfld = jfld + 1 777 dta_bdy(ib_bdy)%ht_s => bf(jfld)%fnow(:,1,:) 778 ELSE ! case input is lim2 type 779 jfld_ai = jfld + 1 780 jfld_hti = jfld + 2 781 jfld_hts = jfld + 3 782 jfld = jfld + 3 783 ALLOCATE( dta_bdy(ib_bdy)%a_i (ilen0(1),jpl) ) 784 ALLOCATE( dta_bdy(ib_bdy)%ht_i(ilen0(1),jpl) ) 785 ALLOCATE( dta_bdy(ib_bdy)%ht_s(ilen0(1),jpl) ) 786 dta_bdy(ib_bdy)%a_i (:,:) = 0.0 787 dta_bdy(ib_bdy)%ht_i(:,:) = 0.0 788 dta_bdy(ib_bdy)%ht_s(:,:) = 0.0 789 ENDIF 790 677 791 ENDIF 678 792 ENDIF -
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90
r3909 r4332 96 96 & nn_dyn3d, nn_dyn3d_dta, nn_tra, nn_tra_dta, & 97 97 & ln_tra_dmp, ln_dyn3d_dmp, rn_time_dmp, & 98 #if defined key_lim299 & nn_ice_lim 2, nn_ice_lim2_dta, &98 #if ( defined key_lim2 || defined key_lim3 ) 99 & nn_ice_lim, nn_ice_lim_dta, & 100 100 #endif 101 101 & ln_vol, nn_volctl, nn_rimwidth … … 137 137 ln_dyn3d_dmp(:) = .false. 138 138 rn_time_dmp(:) = 1. 139 #if defined key_lim2140 nn_ice_lim 2(:) = 0141 nn_ice_lim 2_dta(:)= -1 ! uninitialised flag139 #if ( defined key_lim2 || defined key_lim3 ) 140 nn_ice_lim(:) = 0 141 nn_ice_lim_dta(:)= -1 ! uninitialised flag 142 142 #endif 143 143 ln_vol = .false. … … 161 161 162 162 DO ib_bdy = 1,nb_bdy 163 icount = 0 ! flag to set max rimwidth to 1 if no relaxation 163 164 IF(lwp) WRITE(numout,*) ' ' 164 165 IF(lwp) WRITE(numout,*) '------ Open boundary data set ',ib_bdy,'------' … … 175 176 CASE(jp_none) ; IF(lwp) WRITE(numout,*) ' no open boundary condition' 176 177 CASE(jp_frs) ; IF(lwp) WRITE(numout,*) ' Flow Relaxation Scheme' 178 icount = icount + 1 177 179 CASE(jp_flather) ; IF(lwp) WRITE(numout,*) ' Flather radiation condition' 178 180 CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_dyn2d' ) … … 196 198 CASE(jp_none) ; IF(lwp) WRITE(numout,*) ' no open boundary condition' 197 199 CASE(jp_frs) ; IF(lwp) WRITE(numout,*) ' Flow Relaxation Scheme' 200 icount = icount + 1 198 201 CASE( 2 ) ; IF(lwp) WRITE(numout,*) ' Specified value' 199 202 CASE( 3 ) ; IF(lwp) WRITE(numout,*) ' Zero baroclinic velocities (runoff case)' … … 215 218 CALL ctl_stop( 'Use FRS OR relaxation' ) 216 219 ELSE 220 icount = icount + 1 217 221 IF(lwp) WRITE(numout,*) ' + baroclinic velocities relaxation zone' 218 222 IF(lwp) WRITE(numout,*) ' Damping time scale: ',rn_time_dmp(ib_bdy),' days' … … 228 232 CASE(jp_none) ; IF(lwp) WRITE(numout,*) ' no open boundary condition' 229 233 CASE(jp_frs) ; IF(lwp) WRITE(numout,*) ' Flow Relaxation Scheme' 234 icount = icount + 1 230 235 CASE( 2 ) ; IF(lwp) WRITE(numout,*) ' Specified value' 231 236 CASE( 3 ) ; IF(lwp) WRITE(numout,*) ' Neumann conditions' … … 248 253 CALL ctl_stop( 'Use FRS OR relaxation' ) 249 254 ELSE 255 icount = icount + 1 250 256 IF(lwp) WRITE(numout,*) ' + T/S relaxation zone' 251 257 IF(lwp) WRITE(numout,*) ' Damping time scale: ',rn_time_dmp(ib_bdy),' days' … … 257 263 IF(lwp) WRITE(numout,*) 258 264 259 #if defined key_lim2265 #if ( defined key_lim2 || defined key_lim3 ) 260 266 IF(lwp) WRITE(numout,*) 'Boundary conditions for sea ice: ' 261 SELECT CASE( nn_ice_lim 2(ib_bdy) )267 SELECT CASE( nn_ice_lim(ib_bdy) ) 262 268 CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' no open boundary condition' 263 269 CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' Flow Relaxation Scheme' 270 icount = icount + 1 264 271 CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_tra' ) 265 272 END SELECT 266 IF( nn_ice_lim 2(ib_bdy) .gt. 0 ) THEN267 SELECT CASE( nn_ice_lim 2_dta(ib_bdy) ) !273 IF( nn_ice_lim(ib_bdy) .gt. 0 ) THEN 274 SELECT CASE( nn_ice_lim_dta(ib_bdy) ) ! 268 275 CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' initial state used for bdy data' 269 276 CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' boundary data taken from file' 270 CASE DEFAULT ; CALL ctl_stop( 'nn_ice_lim 2_dta must be 0 or 1' )277 CASE DEFAULT ; CALL ctl_stop( 'nn_ice_lim_dta must be 0 or 1' ) 271 278 END SELECT 272 279 ENDIF 273 280 IF(lwp) WRITE(numout,*) 274 281 #endif 275 276 IF(lwp) WRITE(numout,*) ' Width of relaxation zone = ', nn_rimwidth(ib_bdy) 277 IF(lwp) WRITE(numout,*) 282 IF ( icount>0 ) THEN 283 IF(lwp) WRITE(numout,*) ' Width of relaxation zone = ', nn_rimwidth(ib_bdy) 284 IF(lwp) WRITE(numout,*) 285 ELSE 286 nn_rimwidth(ib_bdy) = 1 ! no relaxation 287 ENDIF 278 288 279 289 ENDDO … … 387 397 DO igrd = 1, jpbgrd 388 398 id_dummy = iom_varid( inum, 'nbi'//cgrid(igrd), kdimsz=kdimsz ) 389 nblendta(igrd,ib_bdy) = kdimsz(1) 390 jpbdtau = MAX(jpbdtau, kdimsz(1)) 399 !clem nblendta(igrd,ib_bdy) = kdimsz(1) 400 !clem jpbdtau = MAX(jpbdtau, kdimsz(1)) 401 nblendta(igrd,ib_bdy) = MAXVAL(kdimsz) 402 jpbdtau = MAX(jpbdtau, MAXVAL(kdimsz)) 391 403 ENDDO 392 404 CALL iom_close( inum ) 393 394 405 ENDIF 395 406 … … 398 409 IF (nb_bdy>0) THEN 399 410 jpbdta = MAXVAL(nblendta(1:jpbgrd,1:nb_bdy)) 400 401 411 ! Allocate arrays 402 412 !--------------- … … 446 456 ENDIF 447 457 448 ENDDO 458 ENDDO 449 459 450 460 ! 2. Now fill indices corresponding to straight open boundary arrays: … … 752 762 ! check if point is in local domain 753 763 IF( nbidta(ib,igrd,ib_bdy) >= iw .AND. nbidta(ib,igrd,ib_bdy) <= ie .AND. & 754 & nbjdta(ib,igrd,ib_bdy) >= is .AND. nbjdta(ib,igrd,ib_bdy) <= in ) THEN 764 & nbjdta(ib,igrd,ib_bdy) >= is .AND. nbjdta(ib,igrd,ib_bdy) <= in .AND. & 765 & nbrdta(ib,igrd,ib_bdy) <= nn_rimwidth(ib_bdy) ) THEN 755 766 ! 756 767 icount = icount + 1 … … 765 776 ! Allocate index arrays for this boundary set 766 777 !-------------------------------------------- 767 ilen1 = MAXVAL(idx_bdy(ib_bdy)%nblen(:)) 778 779 ilen1 = MAXVAL(idx_bdy(ib_bdy)%nblen(1:jpbgrd)) 780 ilen1 = MAX(1,ilen1) 768 781 ALLOCATE( idx_bdy(ib_bdy)%nbi(ilen1,jpbgrd) ) 769 782 ALLOCATE( idx_bdy(ib_bdy)%nbj(ilen1,jpbgrd) ) … … 773 786 ALLOCATE( idx_bdy(ib_bdy)%nbw(ilen1,jpbgrd) ) 774 787 ALLOCATE( idx_bdy(ib_bdy)%flagu(ilen1) ) 775 ALLOCATE( idx_bdy(ib_bdy)%flagv(ilen1) ) 788 ALLOCATE( idx_bdy(ib_bdy)%flagv(ilen1) ) 776 789 777 790 ! Dispatch mapping indices and discrete distances on each processor … … 1008 1021 ! bdytmask = 1 on the computational domain AND on open boundaries 1009 1022 ! = 0 elsewhere 1010 1023 bdytmask(:,:) = 1.e0 1024 bdyumask(:,:) = 1.e0 1025 bdyvmask(:,:) = 1.e0 1026 1011 1027 IF( ln_mask_file ) THEN 1012 1028 CALL iom_open( cn_mask_file, inum ) … … 1053 1069 1054 1070 bdytmask(:,:) = tmask(:,:,1) 1055 IF( .not. ln_mask_file ) THEN1056 ! If .not. ln_mask_file then we need to derive mask on U and V grid1057 ! from mask on T grid here.1058 bdyumask(:,:) = 0.e01059 bdyvmask(:,:) = 0.e01060 DO ij=1, jpjm11061 DO ii=1, jpim11062 bdyumask(ii,ij)=bdytmask(ii,ij)*bdytmask(ii+1, ij )1063 bdyvmask(ii,ij)=bdytmask(ii,ij)*bdytmask(ii ,ij+1)1064 END DO1065 END DO1066 CALL lbc_lnk( bdyumask(:,:), 'U', 1. ) ; CALL lbc_lnk( bdyvmask(:,:), 'V', 1. ) ! Lateral boundary cond.1067 ENDIF1068 1071 1069 1072 ! bdy masks and bmask are now set to zero on boundary points: -
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90
r4220 r4332 237 237 ENDIF 238 238 ! 239 IF( lk_lim3 ) THEN 239 IF( lk_lim3 ) THEN 240 240 CALL iom_get( numror, jpdom_autoglo, 'iatte' , iatte ) ! clem modif 241 241 CALL iom_get( numror, jpdom_autoglo, 'oatte' , oatte ) ! clem modif -
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90
r4040 r4332 724 724 Cd_n10(:,:) = cdn_wave 725 725 ELSE 726 Cd_n10 = 1 E-3 * ( 2.7/dU10 + 0.142 + dU10/13.09 ) ! L & Y eq. (6a)726 Cd_n10 = 1.e-3 * ( 2.7/dU10 + 0.142 + dU10/13.09 ) ! L & Y eq. (6a) 727 727 ENDIF 728 728 sqrt_Cd_n10 = sqrt(Cd_n10) 729 Ce_n10 = 1 E-3 * ( 34.6 * sqrt_Cd_n10 ) ! L & Y eq. (6b)730 Ch_n10 = 1 E-3*sqrt_Cd_n10*(18*stab + 32.7*(1-stab)) ! L & Y eq. (6c), (6d)729 Ce_n10 = 1.e-3 * ( 34.6 * sqrt_Cd_n10 ) ! L & Y eq. (6b) 730 Ch_n10 = 1.e-3*sqrt_Cd_n10*(18.*stab + 32.7*(1.-stab)) ! L & Y eq. (6c), (6d) 731 731 !! 732 732 !! Initializing transfert coefficients with their first guess neutral equivalents : … … 755 755 756 756 !! Updating the neutral 10m transfer coefficients : 757 Cd_n10 = 1 E-3 * (2.7/U_n10 + 0.142 + U_n10/13.09) ! L & Y eq. (6a)757 Cd_n10 = 1.e-3 * (2.7/U_n10 + 0.142 + U_n10/13.09) ! L & Y eq. (6a) 758 758 sqrt_Cd_n10 = sqrt(Cd_n10) 759 Ce_n10 = 1 E-3 * (34.6 * sqrt_Cd_n10) ! L & Y eq. (6b)759 Ce_n10 = 1.e-3 * (34.6 * sqrt_Cd_n10) ! L & Y eq. (6b) 760 760 stab = 0.5 + sign(0.5,zeta) 761 Ch_n10 = 1 E-3*sqrt_Cd_n10*(18.*stab + 32.7*(1-stab)) ! L & Y eq. (6c), (6d)761 Ch_n10 = 1.e-3*sqrt_Cd_n10*(18.*stab + 32.7*(1.-stab)) ! L & Y eq. (6c), (6d) 762 762 763 763 !! Shifting the neutral 10m transfer coefficients to ( zu , zeta ) : … … 858 858 Cd_n10(:,:) = cdn_wave 859 859 ELSE 860 Cd_n10 = 1 E-3*( 2.7/dU10 + 0.142 + dU10/13.09 )860 Cd_n10 = 1.e-3*( 2.7/dU10 + 0.142 + dU10/13.09 ) 861 861 ENDIF 862 862 sqrt_Cd_n10 = sqrt(Cd_n10) 863 Ce_n10 = 1 E-3*( 34.6 * sqrt_Cd_n10 )864 Ch_n10 = 1 E-3*sqrt_Cd_n10*(18*stab + 32.7*(1- stab))863 Ce_n10 = 1.e-3*( 34.6 * sqrt_Cd_n10 ) 864 Ch_n10 = 1.e-3*sqrt_Cd_n10*(18.*stab + 32.7*(1. - stab)) 865 865 866 866 !! Initializing transf. coeff. with their first guess neutral equivalents : … … 904 904 ELSE 905 905 !! Updating the neutral 10m transfer coefficients : 906 Cd_n10 = 1 E-3 * (2.7/U_n10 + 0.142 + U_n10/13.09) ! L & Y eq. (6a)906 Cd_n10 = 1.e-3 * (2.7/U_n10 + 0.142 + U_n10/13.09) ! L & Y eq. (6a) 907 907 sqrt_Cd_n10 = sqrt(Cd_n10) 908 Ce_n10 = 1 E-3 * (34.6 * sqrt_Cd_n10) ! L & Y eq. (6b)908 Ce_n10 = 1.e-3 * (34.6 * sqrt_Cd_n10) ! L & Y eq. (6b) 909 909 stab = 0.5 + sign(0.5,zeta_u) 910 Ch_n10 = 1 E-3*sqrt_Cd_n10*(18.*stab + 32.7*(1-stab)) ! L & Y eq. (6c-6d)910 Ch_n10 = 1.e-3*sqrt_Cd_n10*(18.*stab + 32.7*(1.-stab)) ! L & Y eq. (6c-6d) 911 911 !! 912 912 !! -
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90
r4220 r4332 27 27 USE iceini ! LIM-3: ice initialisation 28 28 USE dom_ice ! LIM-3: ice domain 29 USE domvvl ! domain: variable volume level(fse3t_m) 29 30 30 31 USE sbc_oce ! Surface boundary condition: ocean fields … … 58 59 USE in_out_manager ! I/O manager 59 60 USE prtctl ! Print control 61 USE lib_fortran ! 60 62 61 63 #if defined key_bdy … … 167 169 ! 168 170 IF( ln_nicep ) THEN ! control print at a given point 169 jiindx = 6 ; jjindx = 47170 WRITE(numout,*) ' The debugging point is : jiindx : ',jiindx, ' jjindx : ',jjindx171 jiindx = 177 ; jjindx = 112 172 IF(lwp) WRITE(numout,*) ' The debugging point is : jiindx : ',jiindx, ' jjindx : ',jjindx 171 173 ENDIF 172 174 ENDIF … … 303 305 fhbri (:,:) = 0._wp ; fheat_mec(:,:) = 0._wp ; fheat_res(:,:) = 0._wp 304 306 fhmec (:,:) = 0._wp ; 305 fmmec (:,:) = 0._wp 307 fmmec (:,:) = 0._wp 308 fmmflx (:,:) = 0._wp 306 309 focea2D(:,:) = 0._wp 307 310 fsup2D (:,:) = 0._wp … … 328 331 CALL lim_rst_opn( kt ) ! Open Ice restart file 329 332 ! 330 IF( ln_nicep ) CALL lim_prt_state( jiindx, jjindx, 1, ' - Beginning the time step - ' ) ! control print 331 ! 332 IF( .NOT. lk_c1d ) THEN ! Ice dynamics & transport (except in 1D case) 333 IF( ln_nicep ) CALL lim_prt_state( kt, jiindx, jjindx, 1, ' - Beginning the time step - ' ) ! control print 334 ! ---------------------------------------------- 335 ! ice dynamics and transport (except in 1D case) 336 ! ---------------------------------------------- 337 IF( .NOT. lk_c1d ) THEN 333 338 CALL lim_dyn( kt ) ! Ice dynamics ( rheology/dynamics ) 334 339 CALL lim_trp( kt ) ! Ice transport ( Advection/diffusion ) 335 340 CALL lim_var_glo2eqv ! equivalent variables, requested for rafting 336 IF( ln_nicep ) CALL lim_prt_state( jiindx, jjindx,-1, ' - ice dyn & trp - ' ) ! control print341 IF( ln_nicep ) CALL lim_prt_state( kt, jiindx, jjindx,-1, ' - ice dyn & trp - ' ) ! control print 337 342 CALL lim_itd_me ! Mechanical redistribution ! (ridging/rafting) 338 343 CALL lim_var_agg( 1 ) 344 #if defined key_bdy 345 ! bdy ice thermo 346 CALL lim_var_glo2eqv ! equivalent variables 347 CALL bdy_ice_lim( kt ) 348 CALL lim_itd_me_zapsmall 349 CALL lim_var_agg(1) 350 IF( ln_nicep ) CALL lim_prt_state( kt, jiindx, jjindx, 1, ' - ice thermo bdy - ' ) ! control print 351 #endif 339 352 CALL lim_update1 340 353 ENDIF … … 349 362 old_oa_i(:,:,:) = oa_i(:,:,:) 350 363 old_smv_i(:,:,:) = smv_i (:,:,:) 351 ! ! Ice thermodynamics 364 365 ! ---------------------------------------------- 366 ! ice thermodynamic 367 ! ---------------------------------------------- 352 368 CALL lim_var_glo2eqv ! equivalent variables 353 369 CALL lim_var_agg(1) ! aggregate ice categories 370 ! previous lead fraction and ice volume for flux calculations 371 pfrld(:,:) = 1._wp - at_i(:,:) 372 phicif(:,:) = vt_i(:,:) 373 ! 354 374 CALL lim_var_bv ! bulk brine volume (diag) 355 375 CALL lim_thd( kt ) ! Ice thermodynamics … … 357 377 oa_i(:,:,:) = oa_i(:,:,:) + a_i(:,:,:) * zcoef 358 378 CALL lim_var_glo2eqv ! this CALL is maybe not necessary (Martin) 359 IF( ln_nicep ) CALL lim_prt_state( jiindx, jjindx, 1, ' - ice thermodyn. - ' ) ! control print379 IF( ln_nicep ) CALL lim_prt_state( kt, jiindx, jjindx, 1, ' - ice thermodyn. - ' ) ! control print 360 380 CALL lim_itd_th( kt ) ! Remap ice categories, lateral accretion ! 361 !362 ! ! Global variables update363 381 CALL lim_var_agg( 1 ) ! requested by limupdate 364 CALL lim_update2 ! Global variables update 365 #if defined key_bdy 366 CALL lim_var_glo2eqv ! 367 CALL bdy_ice_lim( kt ) ! clem modif: bdy ice 368 #endif 382 CALL lim_update2 ! Global variables update 383 369 384 CALL lim_var_glo2eqv ! equivalent variables (outputs) 370 385 CALL lim_var_agg(2) ! aggregate ice thickness categories 371 IF( ln_nicep ) CALL lim_prt_state( jiindx, jjindx, 2, ' - Final state - ' ) ! control print386 IF( ln_nicep ) CALL lim_prt_state( kt, jiindx, jjindx, 2, ' - Final state - ' ) ! control print 372 387 ! 373 388 CALL lim_sbc_flx( kt ) ! Update surface ocean mass, heat and salt fluxes 374 389 ! 375 IF( ln_nicep ) CALL lim_prt_state( jiindx, jjindx, 3, ' - Final state lim_sbc - ' ) ! control print390 IF( ln_nicep ) CALL lim_prt_state( kt, jiindx, jjindx, 3, ' - Final state lim_sbc - ' ) ! control print 376 391 ! 377 392 ! ! Diagnostics and outputs … … 386 401 CALL lim_var_glo2eqv ! ??? 387 402 ! 388 IF( ln_nicep ) CALL lim_ctl 403 IF( ln_nicep ) CALL lim_ctl( kt ) ! alerts in case of model crash 389 404 ! 390 405 ENDIF ! End sea-ice time step only … … 413 428 414 429 415 SUBROUTINE lim_ctl 430 SUBROUTINE lim_ctl( kt ) 416 431 !!----------------------------------------------------------------------- 417 432 !! *** ROUTINE lim_ctl *** … … 419 434 !! ** Purpose : Alerts in case of model crash 420 435 !!------------------------------------------------------------------- 436 INTEGER, INTENT(in) :: kt ! ocean time step 421 437 INTEGER :: ji, jj, jk, jl ! dummy loop indices 422 438 INTEGER :: inb_altests ! number of alert tests (max 20) … … 438 454 DO ji = 1, jpi 439 455 IF( v_i(ji,jj,jl) /= 0._wp .AND. a_i(ji,jj,jl) == 0._wp ) THEN 440 WRITE(numout,*) ' ALERTE 2 : Incompatible volume and concentration '441 WRITE(numout,*) ' at_i ', at_i(ji,jj)442 WRITE(numout,*) ' Point - category', ji, jj, jl443 WRITE(numout,*) ' a_i *** a_i_old ', a_i (ji,jj,jl), old_a_i (ji,jj,jl)444 WRITE(numout,*) ' v_i *** v_i_old ', v_i (ji,jj,jl), old_v_i (ji,jj,jl)445 WRITE(numout,*) ' d_a_i_thd/trp ', d_a_i_thd(ji,jj,jl), d_a_i_trp(ji,jj,jl)446 WRITE(numout,*) ' d_v_i_thd/trp ', d_v_i_thd(ji,jj,jl), d_v_i_trp(ji,jj,jl)456 !WRITE(numout,*) ' ALERTE 2 : Incompatible volume and concentration ' 457 !WRITE(numout,*) ' at_i ', at_i(ji,jj) 458 !WRITE(numout,*) ' Point - category', ji, jj, jl 459 !WRITE(numout,*) ' a_i *** a_i_old ', a_i (ji,jj,jl), old_a_i (ji,jj,jl) 460 !WRITE(numout,*) ' v_i *** v_i_old ', v_i (ji,jj,jl), old_v_i (ji,jj,jl) 461 !WRITE(numout,*) ' d_a_i_thd/trp ', d_a_i_thd(ji,jj,jl), d_a_i_trp(ji,jj,jl) 462 !WRITE(numout,*) ' d_v_i_thd/trp ', d_v_i_thd(ji,jj,jl), d_v_i_trp(ji,jj,jl) 447 463 inb_alp(ialert_id) = inb_alp(ialert_id) + 1 448 464 ENDIF … … 458 474 DO ji = 1, jpi 459 475 IF( ht_i(ji,jj,jl) > 50._wp ) THEN 460 CALL lim_prt_state(ji, jj, 2, ' ALERTE 3 : Very thick ice ' )476 !CALL lim_prt_state( kt, ji, jj, 2, ' ALERTE 3 : Very thick ice ' ) 461 477 inb_alp(ialert_id) = inb_alp(ialert_id) + 1 462 478 ENDIF … … 469 485 DO jj = 1, jpj 470 486 DO ji = 1, jpi 471 IF( MAX( ABS( u_ice(ji,jj) ), ABS( v_ice(ji,jj) ) ) > 0.5 .AND. &487 IF( MAX( ABS( u_ice(ji,jj) ), ABS( v_ice(ji,jj) ) ) > 1.5 .AND. & 472 488 & at_i(ji,jj) > 0._wp ) THEN 473 CALL lim_prt_state(ji, jj, 1, ' ALERTE 4 : Very fast ice ' )474 WRITE(numout,*) ' ice strength : ', strength(ji,jj)475 WRITE(numout,*) ' oceanic stress utau : ', utau(ji,jj)476 WRITE(numout,*) ' oceanic stress vtau : ', vtau(ji,jj)477 WRITE(numout,*) ' sea-ice stress utau_ice : ', utau_ice(ji,jj)478 WRITE(numout,*) ' sea-ice stress vtau_ice : ', vtau_ice(ji,jj)479 WRITE(numout,*) ' oceanic speed u : ', u_oce(ji,jj)480 WRITE(numout,*) ' oceanic speed v : ', v_oce(ji,jj)481 WRITE(numout,*) ' sst : ', sst_m(ji,jj)482 WRITE(numout,*) ' sss : ', sss_m(ji,jj)483 WRITE(numout,*)489 !CALL lim_prt_state( kt, ji, jj, 1, ' ALERTE 4 : Very fast ice ' ) 490 !WRITE(numout,*) ' ice strength : ', strength(ji,jj) 491 !WRITE(numout,*) ' oceanic stress utau : ', utau(ji,jj) 492 !WRITE(numout,*) ' oceanic stress vtau : ', vtau(ji,jj) 493 !WRITE(numout,*) ' sea-ice stress utau_ice : ', utau_ice(ji,jj) 494 !WRITE(numout,*) ' sea-ice stress vtau_ice : ', vtau_ice(ji,jj) 495 !WRITE(numout,*) ' oceanic speed u : ', u_oce(ji,jj) 496 !WRITE(numout,*) ' oceanic speed v : ', v_oce(ji,jj) 497 !WRITE(numout,*) ' sst : ', sst_m(ji,jj) 498 !WRITE(numout,*) ' sss : ', sss_m(ji,jj) 499 !WRITE(numout,*) 484 500 inb_alp(ialert_id) = inb_alp(ialert_id) + 1 485 501 ENDIF … … 493 509 DO ji = 1, jpi 494 510 IF( tms(ji,jj) <= 0._wp .AND. at_i(ji,jj) > 0._wp ) THEN 495 CALL lim_prt_state(ji, jj, 1, ' ALERTE 6 : Ice on continents ' )496 WRITE(numout,*) ' masks s, u, v : ', tms(ji,jj), tmu(ji,jj), tmv(ji,jj)497 WRITE(numout,*) ' sst : ', sst_m(ji,jj)498 WRITE(numout,*) ' sss : ', sss_m(ji,jj)499 WRITE(numout,*) ' at_i(ji,jj) : ', at_i(ji,jj)500 WRITE(numout,*) ' v_ice(ji,jj) : ', v_ice(ji,jj)501 WRITE(numout,*) ' v_ice(ji,jj-1) : ', v_ice(ji,jj-1)502 WRITE(numout,*) ' u_ice(ji-1,jj) : ', u_ice(ji-1,jj)503 WRITE(numout,*) ' u_ice(ji,jj) : ', v_ice(ji,jj)511 !CALL lim_prt_state( kt, ji, jj, 1, ' ALERTE 6 : Ice on continents ' ) 512 !WRITE(numout,*) ' masks s, u, v : ', tms(ji,jj), tmu(ji,jj), tmv(ji,jj) 513 !WRITE(numout,*) ' sst : ', sst_m(ji,jj) 514 !WRITE(numout,*) ' sss : ', sss_m(ji,jj) 515 !WRITE(numout,*) ' at_i(ji,jj) : ', at_i(ji,jj) 516 !WRITE(numout,*) ' v_ice(ji,jj) : ', v_ice(ji,jj) 517 !WRITE(numout,*) ' v_ice(ji,jj-1) : ', v_ice(ji,jj-1) 518 !WRITE(numout,*) ' u_ice(ji-1,jj) : ', u_ice(ji-1,jj) 519 !WRITE(numout,*) ' u_ice(ji,jj) : ', v_ice(ji,jj) 504 520 ! 505 521 inb_alp(ialert_id) = inb_alp(ialert_id) + 1 … … 515 531 DO jj = 1, jpj 516 532 DO ji = 1, jpi 517 !!gm test twice sm_i ... ???? bug? 518 IF( ( ( ABS( sm_i(ji,jj,jl) ) < 0.5 ) .OR. & 519 ( ABS( sm_i(ji,jj,jl) ) < 0.5 ) ) .AND. & 520 ( a_i(ji,jj,jl) > 0._wp ) ) THEN 521 ! CALL lim_prt_state(ji,jj,1, ' ALERTE 7 : Very fresh ice ' ) 533 IF( sm_i(ji,jj,jl) < 0.1 .AND. a_i(ji,jj,jl) > 0._wp ) THEN 534 ! CALL lim_prt_state(kt,ji,jj,1, ' ALERTE 7 : Very fresh ice ' ) 522 535 ! WRITE(numout,*) ' sst : ', sst_m(ji,jj) 523 536 ! WRITE(numout,*) ' sss : ', sss_m(ji,jj) … … 540 553 ( ABS( o_i(ji,jj,jl) ) < 0._wp) ) .AND. & 541 554 ( a_i(ji,jj,jl) > 0._wp ) ) THEN 542 CALL lim_prt_state(ji, jj, 1, ' ALERTE 9 : Wrong ice age ')555 !CALL lim_prt_state( kt, ji, jj, 1, ' ALERTE 9 : Wrong ice age ') 543 556 inb_alp(ialert_id) = inb_alp(ialert_id) + 1 544 557 ENDIF … … 552 565 DO jj = 1, jpj 553 566 DO ji = 1, jpi 554 IF( ABS( sfx (ji,jj) ) .GT. 1.0e-2 ) THEN 555 CALL lim_prt_state(ji, jj, 3, ' ALERTE 5 : High salt flux ' )556 DO jl = 1, jpl557 WRITE(numout,*) ' Category no: ', jl558 WRITE(numout,*) ' a_i : ', a_i (ji,jj,jl) , ' old_a_i : ', old_a_i (ji,jj,jl)559 WRITE(numout,*) ' d_a_i_trp : ', d_a_i_trp(ji,jj,jl) , ' d_a_i_thd : ', d_a_i_thd(ji,jj,jl)560 WRITE(numout,*) ' v_i : ', v_i (ji,jj,jl) , ' old_v_i : ', old_v_i (ji,jj,jl)561 WRITE(numout,*) ' d_v_i_trp : ', d_v_i_trp(ji,jj,jl) , ' d_v_i_thd : ', d_v_i_thd(ji,jj,jl)562 WRITE(numout,*) ' '563 END DO567 IF( ABS( sfx (ji,jj) ) .GT. 1.0e-2 ) THEN ! = 1 psu/day for 1m ocean depth 568 !CALL lim_prt_state( kt, ji, jj, 3, ' ALERTE 5 : High salt flux ' ) 569 !DO jl = 1, jpl 570 !WRITE(numout,*) ' Category no: ', jl 571 !WRITE(numout,*) ' a_i : ', a_i (ji,jj,jl) , ' old_a_i : ', old_a_i (ji,jj,jl) 572 !WRITE(numout,*) ' d_a_i_trp : ', d_a_i_trp(ji,jj,jl) , ' d_a_i_thd : ', d_a_i_thd(ji,jj,jl) 573 !WRITE(numout,*) ' v_i : ', v_i (ji,jj,jl) , ' old_v_i : ', old_v_i (ji,jj,jl) 574 !WRITE(numout,*) ' d_v_i_trp : ', d_v_i_trp(ji,jj,jl) , ' d_v_i_thd : ', d_v_i_thd(ji,jj,jl) 575 !WRITE(numout,*) ' ' 576 !END DO 564 577 inb_alp(ialert_id) = inb_alp(ialert_id) + 1 565 578 ENDIF … … 572 585 DO jj = 1, jpj 573 586 DO ji = 1, jpi 574 IF( ABS( qns(ji,jj) ) > 1500._wp .AND. at_i(ji,jj) > 0._wp )THEN587 IF( ABS( qns(ji,jj) ) > 1500._wp .AND. at_i(ji,jj) > 0._wp ) THEN 575 588 ! 576 WRITE(numout,*) ' ALERTE 8 : Very high non-solar heat flux'577 WRITE(numout,*) ' ji, jj : ', ji, jj578 WRITE(numout,*) ' qns : ', qns(ji,jj)579 WRITE(numout,*) ' sst : ', sst_m(ji,jj)580 WRITE(numout,*) ' sss : ', sss_m(ji,jj)581 WRITE(numout,*) ' qcmif : ', qcmif(ji,jj)582 WRITE(numout,*) ' qldif : ', qldif(ji,jj)583 WRITE(numout,*) ' qcmif : ', qcmif(ji,jj) / rdt_ice584 WRITE(numout,*) ' qldif : ', qldif(ji,jj) / rdt_ice585 WRITE(numout,*) ' qfvbq : ', qfvbq(ji,jj)586 WRITE(numout,*) ' qdtcn : ', qdtcn(ji,jj)587 WRITE(numout,*) ' qfvbq / dt: ', qfvbq(ji,jj) / rdt_ice588 WRITE(numout,*) ' qdtcn / dt: ', qdtcn(ji,jj) / rdt_ice589 WRITE(numout,*) ' fdtcn : ', fdtcn(ji,jj)590 WRITE(numout,*) ' fhmec : ', fhmec(ji,jj)591 WRITE(numout,*) ' fheat_mec : ', fheat_mec(ji,jj)592 WRITE(numout,*) ' fheat_res : ', fheat_res(ji,jj)593 WRITE(numout,*) ' fhbri : ', fhbri(ji,jj)589 !WRITE(numout,*) ' ALERTE 8 : Very high non-solar heat flux' 590 !WRITE(numout,*) ' ji, jj : ', ji, jj 591 !WRITE(numout,*) ' qns : ', qns(ji,jj) 592 !WRITE(numout,*) ' sst : ', sst_m(ji,jj) 593 !WRITE(numout,*) ' sss : ', sss_m(ji,jj) 594 !WRITE(numout,*) ' qcmif : ', qcmif(ji,jj) 595 !WRITE(numout,*) ' qldif : ', qldif(ji,jj) 596 !WRITE(numout,*) ' qcmif : ', qcmif(ji,jj) / rdt_ice 597 !WRITE(numout,*) ' qldif : ', qldif(ji,jj) / rdt_ice 598 !WRITE(numout,*) ' qfvbq : ', qfvbq(ji,jj) 599 !WRITE(numout,*) ' qdtcn : ', qdtcn(ji,jj) 600 !WRITE(numout,*) ' qfvbq / dt: ', qfvbq(ji,jj) / rdt_ice 601 !WRITE(numout,*) ' qdtcn / dt: ', qdtcn(ji,jj) / rdt_ice 602 !WRITE(numout,*) ' fdtcn : ', fdtcn(ji,jj) 603 !WRITE(numout,*) ' fhmec : ', fhmec(ji,jj) 604 !WRITE(numout,*) ' fheat_mec : ', fheat_mec(ji,jj) 605 !WRITE(numout,*) ' fheat_res : ', fheat_res(ji,jj) 606 !WRITE(numout,*) ' fhbri : ', fhbri(ji,jj) 594 607 ! 595 CALL lim_prt_state(ji, jj, 2, ' ')608 !CALL lim_prt_state( kt, ji, jj, 2, ' ') 596 609 inb_alp(ialert_id) = inb_alp(ialert_id) + 1 597 610 ! … … 610 623 DO ji = 1, jpi 611 624 ztmelts = -tmut * s_i(ji,jj,jk,jl) + rtt 612 IF( t_i(ji,jj,jk,jl) >= ztmelts .AND. v_i(ji,jj,jl) > 1.e- 6&625 IF( t_i(ji,jj,jk,jl) >= ztmelts .AND. v_i(ji,jj,jl) > 1.e-10 & 613 626 & .AND. a_i(ji,jj,jl) > 0._wp ) THEN 614 WRITE(numout,*) ' ALERTE 10 : Very warm ice'615 WRITE(numout,*) ' ji, jj, jk, jl : ', ji, jj, jk, jl616 WRITE(numout,*) ' t_i : ', t_i(ji,jj,jk,jl)617 WRITE(numout,*) ' e_i : ', e_i(ji,jj,jk,jl)618 WRITE(numout,*) ' s_i : ', s_i(ji,jj,jk,jl)619 WRITE(numout,*) ' ztmelts : ', ztmelts627 !WRITE(numout,*) ' ALERTE 10 : Very warm ice' 628 !WRITE(numout,*) ' ji, jj, jk, jl : ', ji, jj, jk, jl 629 !WRITE(numout,*) ' t_i : ', t_i(ji,jj,jk,jl) 630 !WRITE(numout,*) ' e_i : ', e_i(ji,jj,jk,jl) 631 !WRITE(numout,*) ' s_i : ', s_i(ji,jj,jk,jl) 632 !WRITE(numout,*) ' ztmelts : ', ztmelts 620 633 inb_alp(ialert_id) = inb_alp(ialert_id) + 1 621 634 ENDIF … … 625 638 END DO 626 639 627 ialert_id = 1 ! reference number of this alert 628 cl_alname(ialert_id) = ' NO alerte 1 ' ! name of the alert 629 WRITE(numout,*) 630 WRITE(numout,*) ' All alerts at the end of ice model ' 631 DO ialert_id = 1, inb_altests 632 WRITE(numout,*) ialert_id, cl_alname(ialert_id)//' : ', inb_alp(ialert_id), ' times ! ' 633 END DO 634 ! 640 ! sum of the alerts on all processors 641 IF( lk_mpp ) THEN 642 DO ialert_id = 1, inb_altests 643 CALL mpp_sum(inb_alp(ialert_id)) 644 END DO 645 ENDIF 646 647 ! print alerts 648 IF( lwp ) THEN 649 ialert_id = 1 ! reference number of this alert 650 cl_alname(ialert_id) = ' NO alerte 1 ' ! name of the alert 651 WRITE(numout,*) ' time step ',kt 652 WRITE(numout,*) ' All alerts at the end of ice model ' 653 DO ialert_id = 1, inb_altests 654 WRITE(numout,*) ialert_id, cl_alname(ialert_id)//' : ', inb_alp(ialert_id), ' times ! ' 655 END DO 656 ENDIF 657 ! 635 658 END SUBROUTINE lim_ctl 636 659 637 660 638 SUBROUTINE lim_prt_state( k i, kj, kn, cd1 )661 SUBROUTINE lim_prt_state( kt, ki, kj, kn, cd1 ) 639 662 !!----------------------------------------------------------------------- 640 663 !! *** ROUTINE lim_prt_state *** … … 650 673 !! n : number of the option 651 674 !!------------------------------------------------------------------- 675 INTEGER , INTENT(in) :: kt ! ocean time step 652 676 INTEGER , INTENT(in) :: ki, kj, kn ! ocean gridpoint indices 653 677 CHARACTER(len=*), INTENT(in) :: cd1 ! 654 678 !! 655 INTEGER :: jl 679 INTEGER :: jl, ji, jj 656 680 !!------------------------------------------------------------------- 657 681 658 WRITE(numout,*) cd1 ! print title 659 660 !---------------- 661 ! Simple state 662 !---------------- 663 664 IF ( kn == 1 .OR. kn == -1 ) THEN 665 WRITE(numout,*) ' lim_prt_state - Point : ',ki,kj 666 WRITE(numout,*) ' ~~~~~~~~~~~~~~ ' 667 WRITE(numout,*) ' Simple state ' 668 WRITE(numout,*) ' masks s,u,v : ', tms(ki,kj), tmu(ki,kj), tmv(ki,kj) 669 WRITE(numout,*) ' lat - long : ', gphit(ki,kj), glamt(ki,kj) 670 WRITE(numout,*) ' Time step : ', numit 671 WRITE(numout,*) ' - Ice drift ' 672 WRITE(numout,*) ' ~~~~~~~~~~~ ' 673 WRITE(numout,*) ' u_ice(i-1,j) : ', u_ice(ki-1,kj) 674 WRITE(numout,*) ' u_ice(i ,j) : ', u_ice(ki,kj) 675 WRITE(numout,*) ' v_ice(i ,j-1): ', v_ice(ki,kj-1) 676 WRITE(numout,*) ' v_ice(i ,j) : ', v_ice(ki,kj) 677 WRITE(numout,*) ' strength : ', strength(ki,kj) 678 WRITE(numout,*) 679 WRITE(numout,*) ' - Cell values ' 680 WRITE(numout,*) ' ~~~~~~~~~~~ ' 681 WRITE(numout,*) ' cell area : ', area(ki,kj) 682 WRITE(numout,*) ' at_i : ', at_i(ki,kj) 683 WRITE(numout,*) ' vt_i : ', vt_i(ki,kj) 684 WRITE(numout,*) ' vt_s : ', vt_s(ki,kj) 685 DO jl = 1, jpl 686 WRITE(numout,*) ' - Category (', jl,')' 687 WRITE(numout,*) ' a_i : ', a_i(ki,kj,jl) 688 WRITE(numout,*) ' ht_i : ', ht_i(ki,kj,jl) 689 WRITE(numout,*) ' ht_s : ', ht_s(ki,kj,jl) 690 WRITE(numout,*) ' v_i : ', v_i(ki,kj,jl) 691 WRITE(numout,*) ' v_s : ', v_s(ki,kj,jl) 692 WRITE(numout,*) ' e_s : ', e_s(ki,kj,1,jl)/1.0e9 693 WRITE(numout,*) ' e_i : ', e_i(ki,kj,1:nlay_i,jl)/1.0e9 694 WRITE(numout,*) ' t_su : ', t_su(ki,kj,jl) 695 WRITE(numout,*) ' t_snow : ', t_s(ki,kj,1,jl) 696 WRITE(numout,*) ' t_i : ', t_i(ki,kj,1:nlay_i,jl) 697 WRITE(numout,*) ' sm_i : ', sm_i(ki,kj,jl) 698 WRITE(numout,*) ' smv_i : ', smv_i(ki,kj,jl) 699 WRITE(numout,*) 700 WRITE(numout,*) ' Pathological case : ', patho_case(ki,kj,jl) 701 END DO 702 ENDIF 703 IF( kn == -1 ) THEN 704 WRITE(numout,*) ' Mechanical Check ************** ' 705 WRITE(numout,*) ' Check what means ice divergence ' 706 WRITE(numout,*) ' Total ice concentration ', at_i (ki,kj) 707 WRITE(numout,*) ' Total lead fraction ', ato_i(ki,kj) 708 WRITE(numout,*) ' Sum of both ', ato_i(ki,kj) + at_i(ki,kj) 709 WRITE(numout,*) ' Sum of both minus 1 ', ato_i(ki,kj) + at_i(ki,kj) - 1.00 710 ENDIF 711 712 713 !-------------------- 714 ! Exhaustive state 715 !-------------------- 716 717 IF ( kn .EQ. 2 ) THEN 718 WRITE(numout,*) ' lim_prt_state - Point : ',ki,kj 719 WRITE(numout,*) ' ~~~~~~~~~~~~~~ ' 720 WRITE(numout,*) ' Exhaustive state ' 721 WRITE(numout,*) ' lat - long ', gphit(ki,kj), glamt(ki,kj) 722 WRITE(numout,*) ' Time step ', numit 723 WRITE(numout,*) 724 WRITE(numout,*) ' - Cell values ' 725 WRITE(numout,*) ' ~~~~~~~~~~~ ' 726 WRITE(numout,*) ' cell area : ', area(ki,kj) 727 WRITE(numout,*) ' at_i : ', at_i(ki,kj) 728 WRITE(numout,*) ' vt_i : ', vt_i(ki,kj) 729 WRITE(numout,*) ' vt_s : ', vt_s(ki,kj) 730 WRITE(numout,*) ' u_ice(i-1,j) : ', u_ice(ki-1,kj) 731 WRITE(numout,*) ' u_ice(i ,j) : ', u_ice(ki,kj) 732 WRITE(numout,*) ' v_ice(i ,j-1): ', v_ice(ki,kj-1) 733 WRITE(numout,*) ' v_ice(i ,j) : ', v_ice(ki,kj) 734 WRITE(numout,*) ' strength : ', strength(ki,kj) 735 WRITE(numout,*) ' d_u_ice_dyn : ', d_u_ice_dyn(ki,kj), ' d_v_ice_dyn : ', d_v_ice_dyn(ki,kj) 736 WRITE(numout,*) ' old_u_ice : ', old_u_ice(ki,kj) , ' old_v_ice : ', old_v_ice(ki,kj) 737 WRITE(numout,*) 738 739 DO jl = 1, jpl 740 WRITE(numout,*) ' - Category (',jl,')' 741 WRITE(numout,*) ' ~~~~~~~~ ' 742 WRITE(numout,*) ' ht_i : ', ht_i(ki,kj,jl) , ' ht_s : ', ht_s(ki,kj,jl) 743 WRITE(numout,*) ' t_i : ', t_i(ki,kj,1:nlay_i,jl) 744 WRITE(numout,*) ' t_su : ', t_su(ki,kj,jl) , ' t_s : ', t_s(ki,kj,1,jl) 745 WRITE(numout,*) ' sm_i : ', sm_i(ki,kj,jl) , ' o_i : ', o_i(ki,kj,jl) 746 WRITE(numout,*) ' a_i : ', a_i(ki,kj,jl) , ' old_a_i : ', old_a_i(ki,kj,jl) 747 WRITE(numout,*) ' d_a_i_trp : ', d_a_i_trp(ki,kj,jl) , ' d_a_i_thd : ', d_a_i_thd(ki,kj,jl) 748 WRITE(numout,*) ' v_i : ', v_i(ki,kj,jl) , ' old_v_i : ', old_v_i(ki,kj,jl) 749 WRITE(numout,*) ' d_v_i_trp : ', d_v_i_trp(ki,kj,jl) , ' d_v_i_thd : ', d_v_i_thd(ki,kj,jl) 750 WRITE(numout,*) ' v_s : ', v_s(ki,kj,jl) , ' old_v_s : ', old_v_s(ki,kj,jl) 751 WRITE(numout,*) ' d_v_s_trp : ', d_v_s_trp(ki,kj,jl) , ' d_v_s_thd : ', d_v_s_thd(ki,kj,jl) 752 WRITE(numout,*) ' e_i1 : ', e_i(ki,kj,1,jl)/1.0e9 , ' old_ei1 : ', old_e_i(ki,kj,1,jl)/1.0e9 753 WRITE(numout,*) ' de_i1_trp : ', d_e_i_trp(ki,kj,1,jl)/1.0e9, ' de_i1_thd : ', d_e_i_thd(ki,kj,1,jl)/1.0e9 754 WRITE(numout,*) ' e_i2 : ', e_i(ki,kj,2,jl)/1.0e9 , ' old_ei2 : ', old_e_i(ki,kj,2,jl)/1.0e9 755 WRITE(numout,*) ' de_i2_trp : ', d_e_i_trp(ki,kj,2,jl)/1.0e9, ' de_i2_thd : ', d_e_i_thd(ki,kj,2,jl)/1.0e9 756 WRITE(numout,*) ' e_snow : ', e_s(ki,kj,1,jl) , ' old_e_snow : ', old_e_s(ki,kj,1,jl) 757 WRITE(numout,*) ' d_e_s_trp : ', d_e_s_trp(ki,kj,1,jl) , ' d_e_s_thd : ', d_e_s_thd(ki,kj,1,jl) 758 WRITE(numout,*) ' smv_i : ', smv_i(ki,kj,jl) , ' old_smv_i : ', old_smv_i(ki,kj,jl) 759 WRITE(numout,*) ' d_smv_i_trp: ', d_smv_i_trp(ki,kj,jl) , ' d_smv_i_thd: ', d_smv_i_thd(ki,kj,jl) 760 WRITE(numout,*) ' oa_i : ', oa_i(ki,kj,jl) , ' old_oa_i : ', old_oa_i(ki,kj,jl) 761 WRITE(numout,*) ' d_oa_i_trp : ', d_oa_i_trp(ki,kj,jl) , ' d_oa_i_thd : ', d_oa_i_thd(ki,kj,jl) 762 WRITE(numout,*) ' Path. case : ', patho_case(ki,kj,jl) 763 END DO !jl 764 765 WRITE(numout,*) 766 WRITE(numout,*) ' - Heat / FW fluxes ' 767 WRITE(numout,*) ' ~~~~~~~~~~~~~~~~ ' 768 ! WRITE(numout,*) ' sfx_bri : ', sfx_bri (ki,kj) 769 ! WRITE(numout,*) ' sfx : ', sfx (ki,kj) 770 ! WRITE(numout,*) ' sfx_res : ', sfx_res(ki,kj) 771 WRITE(numout,*) ' fmmec : ', fmmec (ki,kj) 772 WRITE(numout,*) ' fhmec : ', fhmec (ki,kj) 773 WRITE(numout,*) ' fhbri : ', fhbri (ki,kj) 774 WRITE(numout,*) ' fheat_mec : ', fheat_mec(ki,kj) 775 WRITE(numout,*) 776 WRITE(numout,*) ' sst : ', sst_m(ki,kj) 777 WRITE(numout,*) ' sss : ', sss_m(ki,kj) 778 WRITE(numout,*) 779 WRITE(numout,*) ' - Stresses ' 780 WRITE(numout,*) ' ~~~~~~~~ ' 781 WRITE(numout,*) ' utau_ice : ', utau_ice(ki,kj) 782 WRITE(numout,*) ' vtau_ice : ', vtau_ice(ki,kj) 783 WRITE(numout,*) ' utau : ', utau (ki,kj) 784 WRITE(numout,*) ' vtau : ', vtau (ki,kj) 785 WRITE(numout,*) ' oc. vel. u : ', u_oce (ki,kj) 786 WRITE(numout,*) ' oc. vel. v : ', v_oce (ki,kj) 787 ENDIF 788 789 !--------------------- 790 ! Salt / heat fluxes 791 !--------------------- 792 793 IF ( kn .EQ. 3 ) THEN 794 WRITE(numout,*) ' lim_prt_state - Point : ',ki,kj 795 WRITE(numout,*) ' ~~~~~~~~~~~~~~ ' 796 WRITE(numout,*) ' - Salt / Heat Fluxes ' 797 WRITE(numout,*) ' ~~~~~~~~~~~~~~~~ ' 798 WRITE(numout,*) ' lat - long ', gphit(ki,kj), glamt(ki,kj) 799 WRITE(numout,*) ' Time step ', numit 800 WRITE(numout,*) 801 WRITE(numout,*) ' - Heat fluxes at bottom interface ***' 802 WRITE(numout,*) ' qsr : ', qsr(ki,kj) 803 WRITE(numout,*) ' qns : ', qns(ki,kj) 804 WRITE(numout,*) 805 WRITE(numout,*) ' - Salt fluxes at bottom interface ***' 806 WRITE(numout,*) ' emp : ', emp (ki,kj) 807 WRITE(numout,*) ' sfx_bri : ', sfx_bri(ki,kj) 808 WRITE(numout,*) ' sfx : ', sfx (ki,kj) 809 WRITE(numout,*) ' sfx_res : ', sfx_res(ki,kj) 810 WRITE(numout,*) ' sfx_mec : ', sfx_mec(ki,kj) 811 WRITE(numout,*) ' - Heat fluxes at bottom interface ***' 812 WRITE(numout,*) ' fheat_res : ', fheat_res(ki,kj) 813 WRITE(numout,*) 814 WRITE(numout,*) ' - Momentum fluxes ' 815 WRITE(numout,*) ' utau : ', utau(ki,kj) 816 WRITE(numout,*) ' vtau : ', vtau(ki,kj) 817 ENDIF 818 WRITE(numout,*) ' ' 819 ! 682 DO ji = mi0(ki), mi1(ki) 683 DO jj = mj0(kj), mj1(kj) 684 685 WRITE(numout,*) ' time step ',kt,' ',cd1 ! print title 686 687 !---------------- 688 ! Simple state 689 !---------------- 690 691 IF ( kn == 1 .OR. kn == -1 ) THEN 692 WRITE(numout,*) ' lim_prt_state - Point : ',ji,jj 693 WRITE(numout,*) ' ~~~~~~~~~~~~~~ ' 694 WRITE(numout,*) ' Simple state ' 695 WRITE(numout,*) ' masks s,u,v : ', tms(ji,jj), tmu(ji,jj), tmv(ji,jj) 696 WRITE(numout,*) ' lat - long : ', gphit(ji,jj), glamt(ji,jj) 697 WRITE(numout,*) ' Time step : ', numit 698 WRITE(numout,*) ' - Ice drift ' 699 WRITE(numout,*) ' ~~~~~~~~~~~ ' 700 WRITE(numout,*) ' u_ice(i-1,j) : ', u_ice(ji-1,jj) 701 WRITE(numout,*) ' u_ice(i ,j) : ', u_ice(ji,jj) 702 WRITE(numout,*) ' v_ice(i ,j-1): ', v_ice(ji,jj-1) 703 WRITE(numout,*) ' v_ice(i ,j) : ', v_ice(ji,jj) 704 WRITE(numout,*) ' strength : ', strength(ji,jj) 705 WRITE(numout,*) 706 WRITE(numout,*) ' - Cell values ' 707 WRITE(numout,*) ' ~~~~~~~~~~~ ' 708 WRITE(numout,*) ' cell area : ', area(ji,jj) 709 WRITE(numout,*) ' at_i : ', at_i(ji,jj) 710 WRITE(numout,*) ' vt_i : ', vt_i(ji,jj) 711 WRITE(numout,*) ' vt_s : ', vt_s(ji,jj) 712 DO jl = 1, jpl 713 WRITE(numout,*) ' - Category (', jl,')' 714 WRITE(numout,*) ' a_i : ', a_i(ji,jj,jl) 715 WRITE(numout,*) ' ht_i : ', ht_i(ji,jj,jl) 716 WRITE(numout,*) ' ht_s : ', ht_s(ji,jj,jl) 717 WRITE(numout,*) ' v_i : ', v_i(ji,jj,jl) 718 WRITE(numout,*) ' v_s : ', v_s(ji,jj,jl) 719 WRITE(numout,*) ' e_s : ', e_s(ji,jj,1,jl)/1.0e9 720 WRITE(numout,*) ' e_i : ', e_i(ji,jj,1:nlay_i,jl)/1.0e9 721 WRITE(numout,*) ' t_su : ', t_su(ji,jj,jl) 722 WRITE(numout,*) ' t_snow : ', t_s(ji,jj,1,jl) 723 WRITE(numout,*) ' t_i : ', t_i(ji,jj,1:nlay_i,jl) 724 WRITE(numout,*) ' sm_i : ', sm_i(ji,jj,jl) 725 WRITE(numout,*) ' smv_i : ', smv_i(ji,jj,jl) 726 WRITE(numout,*) 727 END DO 728 ENDIF 729 IF( kn == -1 ) THEN 730 WRITE(numout,*) ' Mechanical Check ************** ' 731 WRITE(numout,*) ' Check what means ice divergence ' 732 WRITE(numout,*) ' Total ice concentration ', at_i (ji,jj) 733 WRITE(numout,*) ' Total lead fraction ', ato_i(ji,jj) 734 WRITE(numout,*) ' Sum of both ', ato_i(ji,jj) + at_i(ji,jj) 735 WRITE(numout,*) ' Sum of both minus 1 ', ato_i(ji,jj) + at_i(ji,jj) - 1.00 736 ENDIF 737 738 739 !-------------------- 740 ! Exhaustive state 741 !-------------------- 742 743 IF ( kn .EQ. 2 ) THEN 744 WRITE(numout,*) ' lim_prt_state - Point : ',ji,jj 745 WRITE(numout,*) ' ~~~~~~~~~~~~~~ ' 746 WRITE(numout,*) ' Exhaustive state ' 747 WRITE(numout,*) ' lat - long ', gphit(ji,jj), glamt(ji,jj) 748 WRITE(numout,*) ' Time step ', numit 749 WRITE(numout,*) 750 WRITE(numout,*) ' - Cell values ' 751 WRITE(numout,*) ' ~~~~~~~~~~~ ' 752 WRITE(numout,*) ' cell area : ', area(ji,jj) 753 WRITE(numout,*) ' at_i : ', at_i(ji,jj) 754 WRITE(numout,*) ' vt_i : ', vt_i(ji,jj) 755 WRITE(numout,*) ' vt_s : ', vt_s(ji,jj) 756 WRITE(numout,*) ' u_ice(i-1,j) : ', u_ice(ji-1,jj) 757 WRITE(numout,*) ' u_ice(i ,j) : ', u_ice(ji,jj) 758 WRITE(numout,*) ' v_ice(i ,j-1): ', v_ice(ji,jj-1) 759 WRITE(numout,*) ' v_ice(i ,j) : ', v_ice(ji,jj) 760 WRITE(numout,*) ' strength : ', strength(ji,jj) 761 WRITE(numout,*) ' d_u_ice_dyn : ', d_u_ice_dyn(ji,jj), ' d_v_ice_dyn : ', d_v_ice_dyn(ji,jj) 762 WRITE(numout,*) ' old_u_ice : ', old_u_ice(ji,jj) , ' old_v_ice : ', old_v_ice(ji,jj) 763 WRITE(numout,*) 764 765 DO jl = 1, jpl 766 WRITE(numout,*) ' - Category (',jl,')' 767 WRITE(numout,*) ' ~~~~~~~~ ' 768 WRITE(numout,*) ' ht_i : ', ht_i(ji,jj,jl) , ' ht_s : ', ht_s(ji,jj,jl) 769 WRITE(numout,*) ' t_i : ', t_i(ji,jj,1:nlay_i,jl) 770 WRITE(numout,*) ' t_su : ', t_su(ji,jj,jl) , ' t_s : ', t_s(ji,jj,1,jl) 771 WRITE(numout,*) ' sm_i : ', sm_i(ji,jj,jl) , ' o_i : ', o_i(ji,jj,jl) 772 WRITE(numout,*) ' a_i : ', a_i(ji,jj,jl) , ' old_a_i : ', old_a_i(ji,jj,jl) 773 WRITE(numout,*) ' d_a_i_trp : ', d_a_i_trp(ji,jj,jl) , ' d_a_i_thd : ', d_a_i_thd(ji,jj,jl) 774 WRITE(numout,*) ' v_i : ', v_i(ji,jj,jl) , ' old_v_i : ', old_v_i(ji,jj,jl) 775 WRITE(numout,*) ' d_v_i_trp : ', d_v_i_trp(ji,jj,jl) , ' d_v_i_thd : ', d_v_i_thd(ji,jj,jl) 776 WRITE(numout,*) ' v_s : ', v_s(ji,jj,jl) , ' old_v_s : ', old_v_s(ji,jj,jl) 777 WRITE(numout,*) ' d_v_s_trp : ', d_v_s_trp(ji,jj,jl) , ' d_v_s_thd : ', d_v_s_thd(ji,jj,jl) 778 WRITE(numout,*) ' e_i1 : ', e_i(ji,jj,1,jl)/1.0e9 , ' old_ei1 : ', old_e_i(ji,jj,1,jl)/1.0e9 779 WRITE(numout,*) ' de_i1_trp : ', d_e_i_trp(ji,jj,1,jl)/1.0e9, ' de_i1_thd : ', d_e_i_thd(ji,jj,1,jl)/1.0e9 780 WRITE(numout,*) ' e_i2 : ', e_i(ji,jj,2,jl)/1.0e9 , ' old_ei2 : ', old_e_i(ji,jj,2,jl)/1.0e9 781 WRITE(numout,*) ' de_i2_trp : ', d_e_i_trp(ji,jj,2,jl)/1.0e9, ' de_i2_thd : ', d_e_i_thd(ji,jj,2,jl)/1.0e9 782 WRITE(numout,*) ' e_snow : ', e_s(ji,jj,1,jl) , ' old_e_snow : ', old_e_s(ji,jj,1,jl) 783 WRITE(numout,*) ' d_e_s_trp : ', d_e_s_trp(ji,jj,1,jl) , ' d_e_s_thd : ', d_e_s_thd(ji,jj,1,jl) 784 WRITE(numout,*) ' smv_i : ', smv_i(ji,jj,jl) , ' old_smv_i : ', old_smv_i(ji,jj,jl) 785 WRITE(numout,*) ' d_smv_i_trp: ', d_smv_i_trp(ji,jj,jl) , ' d_smv_i_thd: ', d_smv_i_thd(ji,jj,jl) 786 WRITE(numout,*) ' oa_i : ', oa_i(ji,jj,jl) , ' old_oa_i : ', old_oa_i(ji,jj,jl) 787 WRITE(numout,*) ' d_oa_i_trp : ', d_oa_i_trp(ji,jj,jl) , ' d_oa_i_thd : ', d_oa_i_thd(ji,jj,jl) 788 END DO !jl 789 790 WRITE(numout,*) 791 WRITE(numout,*) ' - Heat / FW fluxes ' 792 WRITE(numout,*) ' ~~~~~~~~~~~~~~~~ ' 793 WRITE(numout,*) ' emp : ', emp (ji,jj) 794 WRITE(numout,*) ' sfx : ', sfx (ji,jj) 795 WRITE(numout,*) ' sfx_thd : ', sfx_thd(ji,jj) 796 WRITE(numout,*) ' sfx_bri : ', sfx_bri (ji,jj) 797 WRITE(numout,*) ' sfx_mec : ', sfx_mec (ji,jj) 798 WRITE(numout,*) ' sfx_res : ', sfx_res(ji,jj) 799 WRITE(numout,*) ' fmmec : ', fmmec (ji,jj) 800 WRITE(numout,*) ' fhmec : ', fhmec (ji,jj) 801 WRITE(numout,*) ' fhbri : ', fhbri (ji,jj) 802 WRITE(numout,*) ' fheat_mec : ', fheat_mec(ji,jj) 803 WRITE(numout,*) 804 WRITE(numout,*) ' sst : ', sst_m(ji,jj) 805 WRITE(numout,*) ' sss : ', sss_m(ji,jj) 806 WRITE(numout,*) 807 WRITE(numout,*) ' - Stresses ' 808 WRITE(numout,*) ' ~~~~~~~~ ' 809 WRITE(numout,*) ' utau_ice : ', utau_ice(ji,jj) 810 WRITE(numout,*) ' vtau_ice : ', vtau_ice(ji,jj) 811 WRITE(numout,*) ' utau : ', utau (ji,jj) 812 WRITE(numout,*) ' vtau : ', vtau (ji,jj) 813 WRITE(numout,*) ' oc. vel. u : ', u_oce (ji,jj) 814 WRITE(numout,*) ' oc. vel. v : ', v_oce (ji,jj) 815 ENDIF 816 817 !--------------------- 818 ! Salt / heat fluxes 819 !--------------------- 820 821 IF ( kn .EQ. 3 ) THEN 822 WRITE(numout,*) ' lim_prt_state - Point : ',ji,jj 823 WRITE(numout,*) ' ~~~~~~~~~~~~~~ ' 824 WRITE(numout,*) ' - Salt / Heat Fluxes ' 825 WRITE(numout,*) ' ~~~~~~~~~~~~~~~~ ' 826 WRITE(numout,*) ' lat - long ', gphit(ji,jj), glamt(ji,jj) 827 WRITE(numout,*) ' Time step ', numit 828 WRITE(numout,*) 829 WRITE(numout,*) ' - Heat fluxes at bottom interface ***' 830 WRITE(numout,*) ' qsr : ', qsr(ji,jj) 831 WRITE(numout,*) ' qns : ', qns(ji,jj) 832 WRITE(numout,*) ' fdtcn : ', fdtcn(ji,jj) 833 WRITE(numout,*) ' qcmif : ', qcmif(ji,jj) * r1_rdtice 834 WRITE(numout,*) ' qldif : ', qldif(ji,jj) * r1_rdtice 835 WRITE(numout,*) 836 WRITE(numout,*) ' - Salt fluxes at bottom interface ***' 837 WRITE(numout,*) ' emp : ', emp (ji,jj) 838 WRITE(numout,*) ' sfx_bri : ', sfx_bri(ji,jj) 839 WRITE(numout,*) ' sfx : ', sfx (ji,jj) 840 WRITE(numout,*) ' sfx_res : ', sfx_res(ji,jj) 841 WRITE(numout,*) ' sfx_mec : ', sfx_mec(ji,jj) 842 WRITE(numout,*) ' - Heat fluxes at bottom interface ***' 843 WRITE(numout,*) ' fheat_res : ', fheat_res(ji,jj) 844 WRITE(numout,*) 845 WRITE(numout,*) ' - Momentum fluxes ' 846 WRITE(numout,*) ' utau : ', utau(ji,jj) 847 WRITE(numout,*) ' vtau : ', vtau(ji,jj) 848 ENDIF 849 WRITE(numout,*) ' ' 850 ! 851 END DO 852 END DO 853 820 854 END SUBROUTINE lim_prt_state 821 855 -
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90
r4220 r4332 335 335 ! ! (update freshwater fluxes) 336 336 !RBbug do not understand why see ticket 667 337 CALL lbc_lnk( emp, 'T', 1. )337 !clem-bugsal CALL lbc_lnk( emp, 'T', 1. ) 338 338 ! 339 339 IF( kt == nit000 ) THEN ! set the forcing field at nit000 - 1 ! -
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90
r4220 r4332 33 33 USE timing ! Timing 34 34 USE sbc_ice, ONLY : lk_lim3 35 36 35 37 36 IMPLICIT NONE -
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/step.F90
r4220 r4332 89 89 ! Update data, open boundaries, surface boundary condition (including sea-ice) 90 90 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 91 CALL sbc ( kstp ) ! Sea Boundary Condition (including sea-ice)92 93 91 IF( lk_tide.AND.(kstp /= nit000 )) CALL tide_init ( kstp ) 94 92 IF( lk_tide ) CALL sbc_tide( kstp ) … … 97 95 IF( lk_bdy ) CALL bdy_dta( kstp, time_offset=+1 ) ! update dynamic and tracer data at open boundaries 98 96 97 CALL sbc ( kstp ) ! Sea Boundary Condition (including sea-ice) 98 ! clem: moved here for bdy ice purpose 99 99 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 100 100 ! Ocean dynamics : ssh, wn, hdiv, rot ! -
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/SETTE/sette.sh
r3708 r4332 131 131 #- 132 132 # Compiler among those in NEMOGCM/ARCH 133 COMPILER= PW6_VARGAS133 COMPILER=x3750_ADA 134 134 export BATCH_COMMAND_PAR="llsubmit" 135 135 export BATCH_COMMAND_SEQ=$BATCH_COMMAND_PAR
Note: See TracChangeset
for help on using the changeset viewer.