New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 5053 for branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90 – NEMO

Ignore:
Timestamp:
2015-02-03T18:11:02+01:00 (9 years ago)
Author:
clem
Message:

LIM3 cleaning continues. No change in the physics besides the introduction of the monocategory sea ice

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90

    r5051 r5053  
    3434   USE limthd_sal     ! LIM: thermodynamics, ice salinity 
    3535   USE limthd_ent     ! LIM: thermodynamics, ice enthalpy redistribution 
     36   USE limthd_lac     ! LIM-3 lateral accretion 
     37   USE limitd_th      ! remapping thickness distribution 
    3638   USE limtab         ! LIM: 1D <==> 2D transformation 
    3739   USE limvar         ! LIM: sea-ice variables 
     
    4446   USE timing         ! Timing 
    4547   USE limcons        ! conservation tests 
     48   USE limctl 
    4649 
    4750   IMPLICIT NONE 
     
    8083      !! ** References :  
    8184      !!--------------------------------------------------------------------- 
    82       INTEGER, INTENT(in) ::   kt    ! number of iteration 
     85      INTEGER, INTENT(in) :: kt    ! number of iteration 
    8386      !! 
    8487      INTEGER  :: ji, jj, jk, jl   ! dummy loop indices 
    8588      INTEGER  :: nbpb             ! nb of icy pts for vertical thermo calculations 
    86       INTEGER  :: nbplm            ! nb of icy pts for lateral melting calculations (mono-cat) 
    8789      INTEGER  :: ii, ij           ! temporary dummy loop index 
    88       REAL(wp) :: zfric_umin = 0._wp        ! lower bound for the friction velocity (cice value=5.e-04) 
    89       REAL(wp) :: zch        = 0.0057_wp    ! heat transfer coefficient 
    90       REAL(wp) :: zareamin  
    9190      REAL(wp) :: zfric_u, zqld, zqfr 
    92       ! 
    9391      REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b  
     92      REAL(wp), PARAMETER :: zfric_umin = 0._wp        ! lower bound for the friction velocity (cice value=5.e-04) 
     93      REAL(wp), PARAMETER :: zch        = 0.0057_wp    ! heat transfer coefficient 
    9494      ! 
    9595      REAL(wp), POINTER, DIMENSION(:,:) ::  zqsr, zqns 
     
    106106      !------------------------------------------------------------------------! 
    107107      ftr_ice(:,:,:) = 0._wp  ! part of solar radiation transmitted through the ice 
    108  
    109108 
    110109      !-------------------- 
     
    162161      ENDIF 
    163162 
    164 !CDIR NOVERRCHK 
    165163      DO jj = 1, jpj 
    166 !CDIR NOVERRCHK 
    167164         DO ji = 1, jpi 
    168             rswitch          = tms(ji,jj) * ( 1._wp - MAX( 0._wp , SIGN( 1._wp , - at_i(ji,jj) + epsi10 ) ) ) ! 0 if no ice 
     165            rswitch  = tms(ji,jj) * ( 1._wp - MAX( 0._wp , SIGN( 1._wp , - at_i(ji,jj) + epsi10 ) ) ) ! 0 if no ice 
    169166            ! 
    170167            !           !  solar irradiance transmission at the mixed layer bottom and used in the lead heat budget 
     
    260257         ENDIF 
    261258 
    262          zareamin = epsi10 
    263259         nbpb = 0 
    264260         DO jj = 1, jpj 
    265261            DO ji = 1, jpi 
    266                IF ( a_i(ji,jj,jl) .gt. zareamin ) THEN      
     262               IF ( a_i(ji,jj,jl) > epsi10 ) THEN      
    267263                  nbpb      = nbpb  + 1 
    268264                  npb(nbpb) = (jj - 1) * jpi + ji 
     
    442438              CALL tab_1d_2d( nbpb, hfx_res       , npb, hfx_res_1d(1:nbpb)   , jpi, jpj ) 
    443439              CALL tab_1d_2d( nbpb, hfx_err_rem   , npb, hfx_err_rem_1d(1:nbpb), jpi, jpj ) 
    444  
    445 !clem            IF ( ( ( nn_monocat == 1 ) .OR. ( nn_monocat == 4 ) ) .AND. ( jpl == 1 ) ) THEN 
    446 !clem              CALL tab_1d_2d( nbpb, dh_i_melt(:,:,jl) , npb, dh_i_melt_1d(1:nbpb) , jpi, jpj ) 
    447 !clem            ENDIF 
    448440            ! 
    449441              CALL tab_1d_2d( nbpb, qns_ice(:,:,jl), npb, qns_ice_1d(1:nbpb) , jpi, jpj) 
     
    460452 
    461453      !------------------------ 
    462       ! 5.1) Ice heat content               
     454      ! Ice heat content               
    463455      !------------------------ 
    464456      ! Enthalpies are global variables we have to readjust the units (heat content in Joules) 
     
    470462 
    471463      !------------------------ 
    472       ! 5.2) Snow heat content               
     464      ! Snow heat content               
    473465      !------------------------ 
    474466      ! Enthalpies are global variables we have to readjust the units (heat content in Joules) 
     
    478470         END DO 
    479471      END DO 
     472  
     473      !------------------------ 
     474      ! Ice natural aging               
     475      !------------------------ 
     476      oa_i(:,:,:) = oa_i(:,:,:) + a_i(:,:,:) * rdt_ice /rday 
    480477 
    481478      !---------------------------------- 
    482       ! 5.3) Change thickness to volume 
     479      ! Change thickness to volume 
    483480      !---------------------------------- 
    484481      CALL lim_var_eqv2glo 
     
    487484      ! 5.4) Diagnostic thermodynamic growth rates 
    488485      !-------------------------------------------- 
     486      IF( ln_nicep )   CALL lim_prt( kt, jiindx, jjindx, 1, ' - ice thermodyn. - ' )   ! control print 
     487 
    489488      IF(ln_ctl) THEN            ! Control print 
    490489         CALL prt_ctl_info(' ') 
     
    522521      CALL wrk_dealloc( jpi, jpj, zqsr, zqns ) 
    523522 
     523      !------------------------------------------------------------------------------| 
     524      !  1) Transport of ice between thickness categories.                           | 
     525      !------------------------------------------------------------------------------| 
     526      ! Given thermodynamic growth rates, transport ice between 
     527      ! thickness categories. 
     528      IF( jpl > 1 )   CALL lim_itd_th_rem( 1, jpl, kt ) 
     529      ! 
     530      CALL lim_var_glo2eqv    ! only for info 
     531      CALL lim_var_agg(1) 
     532 
     533      !------------------------------------------------------------------------------| 
     534      !  3) Add frazil ice growing in leads. 
     535      !------------------------------------------------------------------------------| 
     536      CALL lim_thd_lac 
     537      CALL lim_var_glo2eqv    ! only for info 
     538       
     539      IF(ln_ctl) THEN   ! Control print 
     540         CALL prt_ctl_info(' ') 
     541         CALL prt_ctl_info(' - Cell values : ') 
     542         CALL prt_ctl_info('   ~~~~~~~~~~~~~ ') 
     543         CALL prt_ctl(tab2d_1=area , clinfo1=' lim_itd_th  : cell area :') 
     544         CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_itd_th  : at_i      :') 
     545         CALL prt_ctl(tab2d_1=vt_i , clinfo1=' lim_itd_th  : vt_i      :') 
     546         CALL prt_ctl(tab2d_1=vt_s , clinfo1=' lim_itd_th  : vt_s      :') 
     547         DO jl = 1, jpl 
     548            CALL prt_ctl_info(' ') 
     549            CALL prt_ctl_info(' - Category : ', ivar1=jl) 
     550            CALL prt_ctl_info('   ~~~~~~~~~~') 
     551            CALL prt_ctl(tab2d_1=a_i   (:,:,jl)   , clinfo1= ' lim_itd_th  : a_i      : ') 
     552            CALL prt_ctl(tab2d_1=ht_i  (:,:,jl)   , clinfo1= ' lim_itd_th  : ht_i     : ') 
     553            CALL prt_ctl(tab2d_1=ht_s  (:,:,jl)   , clinfo1= ' lim_itd_th  : ht_s     : ') 
     554            CALL prt_ctl(tab2d_1=v_i   (:,:,jl)   , clinfo1= ' lim_itd_th  : v_i      : ') 
     555            CALL prt_ctl(tab2d_1=v_s   (:,:,jl)   , clinfo1= ' lim_itd_th  : v_s      : ') 
     556            CALL prt_ctl(tab2d_1=e_s   (:,:,1,jl) , clinfo1= ' lim_itd_th  : e_s      : ') 
     557            CALL prt_ctl(tab2d_1=t_su  (:,:,jl)   , clinfo1= ' lim_itd_th  : t_su     : ') 
     558            CALL prt_ctl(tab2d_1=t_s   (:,:,1,jl) , clinfo1= ' lim_itd_th  : t_snow   : ') 
     559            CALL prt_ctl(tab2d_1=sm_i  (:,:,jl)   , clinfo1= ' lim_itd_th  : sm_i     : ') 
     560            CALL prt_ctl(tab2d_1=smv_i (:,:,jl)   , clinfo1= ' lim_itd_th  : smv_i    : ') 
     561            DO jk = 1, nlay_i 
     562               CALL prt_ctl_info(' ') 
     563               CALL prt_ctl_info(' - Layer : ', ivar1=jk) 
     564               CALL prt_ctl_info('   ~~~~~~~') 
     565               CALL prt_ctl(tab2d_1=t_i(:,:,jk,jl) , clinfo1= ' lim_itd_th  : t_i      : ') 
     566               CALL prt_ctl(tab2d_1=e_i(:,:,jk,jl) , clinfo1= ' lim_itd_th  : e_i      : ') 
     567            END DO 
     568         END DO 
     569      ENDIF 
    524570      ! 
    525571      ! conservation test 
    526572      IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limthd', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
    527       ! 
     573 
    528574      IF( nn_timing == 1 )  CALL timing_stop('limthd') 
    529575 
     
    571617      INTEGER, INTENT(in) ::   kideb, kiut   ! bounds for the spatial loop 
    572618      INTEGER             ::   ji            ! dummy loop indices 
    573  
    574       WRITE(numout,*) ' Lateral melting ON ' 
     619      REAL(wp)            ::   zhi_bef       ! ice thickness before thermo 
     620      REAL(wp)            ::   zdh_mel       ! net melting 
     621 
    575622      DO ji = kideb, kiut 
    576          IF( ht_i_1d(ji) > epsi10 .AND. dh_i_melt_1d(ji) < 0._wp ) THEN      
    577             a_i_1d(ji) = MAX( 0._wp, a_i_1d(ji) + a_i_1d(ji) * dh_i_melt_1d(ji) / ( 2._wp * ht_i_1d(ji) ) )  
    578          END IF 
     623         zhi_bef = ht_i_1d(ji) - ( dh_i_surf(ji) + dh_i_bott(ji) + dh_snowice(ji) ) 
     624         zdh_mel = MIN( 0._wp,     dh_i_surf(ji) + dh_i_bott(ji) + dh_snowice(ji) ) 
     625         IF( zdh_mel < 0._wp )   & 
     626            &   a_i_1d(ji) = MAX( 0._wp, a_i_1d(ji) + a_i_1d(ji) * zdh_mel / ( 2._wp * MAX( zhi_bef, epsi10 ) ) )  
    579627      END DO 
    580628      at_i_1d(:) = a_i_1d(:) 
Note: See TracChangeset for help on using the changeset viewer.