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 5167 for trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90 – NEMO

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

LIM3 bug fix. see ticket #1492 (forthcoming update) which also solve ticket #1497

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90

    r5146 r5167  
    9494      REAL(wp), POINTER, DIMENSION(:,:) ::  zqsr, zqns 
    9595      !!------------------------------------------------------------------- 
    96       CALL wrk_alloc( jpi, jpj, zqsr, zqns ) 
     96      CALL wrk_alloc( jpi,jpj, zqsr, zqns ) 
    9797 
    9898      IF( nn_timing == 1 )  CALL timing_start('limthd') 
     
    101101      IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limthd', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
    102102 
     103      CALL lim_var_glo2eqv 
    103104      !------------------------------------------------------------------------! 
    104105      ! 1) Initialization of some variables                                    ! 
     
    209210            ! Net heat flux on top of ice-ocean [W.m-2] 
    210211            ! ----------------------------------------- 
    211             !     First  step here      : heat flux at the ocean surface + precip 
    212             !     Second step below     : heat flux at the ice   surface (after limthd_dif)  
     212            !     heat flux at the ocean surface + precip 
     213            !   + heat flux at the ice   surface  
    213214            hfx_in(ji,jj) = hfx_in(ji,jj)                                                                                         &  
    214215               ! heat flux above the ocean 
     
    216217               ! latent heat of precip (note that precip is included in qns but not in qns_ice) 
    217218               &    +   ( 1._wp - pfrld(ji,jj) ) * sprecip(ji,jj) * ( cpic * ( MIN( tatm_ice(ji,jj), rt0_snow ) - rt0 ) - lfus )  & 
    218                &    +   ( 1._wp - pfrld(ji,jj) ) * ( tprecip(ji,jj) - sprecip(ji,jj) ) * rcp * ( tatm_ice(ji,jj) - rt0 ) 
     219               &    +   ( 1._wp - pfrld(ji,jj) ) * ( tprecip(ji,jj) - sprecip(ji,jj) ) * rcp * ( tatm_ice(ji,jj) - rt0 )          & 
     220               ! heat flux above the ice 
     221               &    +   SUM(    a_i_b(ji,jj,:)   * ( qns_ice(ji,jj,:) + qsr_ice(ji,jj,:) ) ) 
    219222 
    220223            ! ----------------------------------------------------------------------------- 
     
    226229            hfx_out(ji,jj) = hfx_out(ji,jj)                                                                                       &  
    227230               ! Non solar heat flux received by the ocean 
    228                &    +        pfrld(ji,jj) * qns(ji,jj)                                                                            & 
     231               &    +        pfrld(ji,jj) * zqns(ji,jj)                                                                            & 
    229232               ! latent heat of precip (note that precip is included in qns but not in qns_ice) 
    230233               &    +      ( pfrld(ji,jj)**rn_betas - pfrld(ji,jj) ) * sprecip(ji,jj)       & 
     
    311314            ! --- lateral melting if monocat --- ! 
    312315            !------------------------------------! 
    313             IF ( ( ( nn_monocat == 1 ) .OR. ( nn_monocat == 4 ) ) .AND. ( jpl == 1 ) ) THEN 
     316            IF ( ( nn_monocat == 1 .OR. nn_monocat == 4 ) .AND. jpl == 1 ) THEN 
    314317               CALL lim_thd_lam( 1, nbpb ) 
    315318            END IF 
     
    324327         ENDIF 
    325328         ! 
    326       END DO 
     329      END DO !jl 
    327330 
    328331      !------------------------------------------------------------------------------! 
     
    358361      ! Change thickness to volume 
    359362      !---------------------------------- 
    360       CALL lim_var_eqv2glo 
     363      v_i(:,:,:)   = ht_i(:,:,:) * a_i(:,:,:) 
     364      v_s(:,:,:)   = ht_s(:,:,:) * a_i(:,:,:) 
     365      smv_i(:,:,:) = sm_i(:,:,:) * v_i(:,:,:) 
    361366 
    362367      CALL lim_var_zapsmall 
     
    399404      ! 
    400405      ! 
    401       CALL wrk_dealloc( jpi, jpj, zqsr, zqns ) 
    402  
    403406      IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limthd', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
     407 
     408      CALL wrk_dealloc( jpi,jpj, zqsr, zqns ) 
     409 
    404410      !------------------------------------------------------------------------------| 
    405411      !  6) Transport of ice between thickness categories.                           | 
    406412      !------------------------------------------------------------------------------| 
     413      ! Given thermodynamic growth rates, transport ice between thickness categories. 
    407414      IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limitd_th_rem', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
    408415 
    409       ! Given thermodynamic growth rates, transport ice between thickness categories. 
    410       IF( jpl > 1 )   CALL lim_itd_th_rem( 1, jpl, kt ) 
    411       ! 
    412       CALL lim_var_glo2eqv    ! only for info 
    413       CALL lim_var_agg(1) 
     416      IF( jpl > 1 )      CALL lim_itd_th_rem( 1, jpl, kt ) 
    414417 
    415418      IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limitd_th_rem', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
     419 
    416420      !------------------------------------------------------------------------------| 
    417421      !  7) Add frazil ice growing in leads. 
    418422      !------------------------------------------------------------------------------| 
    419423      IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limthd_lac', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
     424 
    420425      CALL lim_thd_lac 
    421       CALL lim_var_glo2eqv    ! only for info 
    422426       
    423       ! conservation test 
    424427      IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limthd_lac', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
    425428 
    426       IF(ln_ctl) THEN   ! Control print 
     429      ! Control print 
     430      IF(ln_ctl) THEN 
     431         CALL lim_var_glo2eqv 
     432 
    427433         CALL prt_ctl_info(' ') 
    428434         CALL prt_ctl_info(' - Cell values : ') 
     
    503509      REAL(wp)            ::   zhi_bef            ! ice thickness before thermo 
    504510      REAL(wp)            ::   zdh_mel, zda_mel   ! net melting 
    505       REAL(wp)            ::   zv                 ! ice volume  
     511      REAL(wp)            ::   zvi, zvs           ! ice/snow volumes  
    506512 
    507513      DO ji = kideb, kiut 
    508514         zdh_mel = MIN( 0._wp, dh_i_surf(ji) + dh_i_bott(ji) + dh_snowice(ji) ) 
    509          IF( zdh_mel < 0._wp )  THEN 
    510             zv         = a_i_1d(ji) * ht_i_1d(ji) 
     515         IF( zdh_mel < 0._wp .AND. a_i_1d(ji) > 0._wp )  THEN 
     516            zvi          = a_i_1d(ji) * ht_i_1d(ji) 
     517            zvs          = a_i_1d(ji) * ht_s_1d(ji) 
    511518            ! lateral melting = concentration change 
    512519            zhi_bef     = ht_i_1d(ji) - zdh_mel 
    513             zda_mel     =  a_i_1d(ji) * zdh_mel / ( 2._wp * MAX( zhi_bef, epsi10 ) ) 
    514             a_i_1d(ji)  = MAX( 0._wp, a_i_1d(ji) + zda_mel )  
    515             ! adjust thickness 
    516             rswitch     = MAX( 0._wp , SIGN( 1._wp , a_i_1d(ji) - epsi20 ) ) 
    517             ht_i_1d(ji) = rswitch * zv / MAX( a_i_1d(ji), epsi20 ) 
     520            rswitch     = MAX( 0._wp , SIGN( 1._wp , zhi_bef - epsi20 ) ) 
     521            zda_mel     = rswitch * a_i_1d(ji) * zdh_mel / ( 2._wp * MAX( zhi_bef, epsi20 ) ) 
     522            a_i_1d(ji)  = MAX( epsi20, a_i_1d(ji) + zda_mel )  
     523             ! adjust thickness 
     524            ht_i_1d(ji) = zvi / a_i_1d(ji)             
     525            ht_s_1d(ji) = zvs / a_i_1d(ji)             
    518526            ! retrieve total concentration 
    519527            at_i_1d(ji) = a_i_1d(ji) 
     
    676684      INTEGER  ::   ios                 ! Local integer output status for namelist read 
    677685      NAMELIST/namicethd/ rn_hnewice, ln_frazil, rn_maxfrazb, rn_vfrazb, rn_Cfrazb,                       & 
    678          &                rn_himin, parsub, rn_betas, rn_kappa_i, nn_conv_dif, rn_terr_dif, nn_ice_thcon, & 
     686         &                rn_himin, rn_betas, rn_kappa_i, nn_conv_dif, rn_terr_dif, nn_ice_thcon, & 
    679687         &                nn_monocat, ln_it_qnsice 
    680688      !!------------------------------------------------------------------- 
     
    700708      ENDIF 
    701709 
    702       IF( lk_cpl .AND. parsub /= 0.0 )   CALL ctl_stop( 'In coupled mode, use parsub = 0. or send dqla' ) 
    703710      ! 
    704711      IF(lwp) THEN                          ! control print 
     
    712719         WRITE(numout,*)'      minimum ice thickness                                   rn_himin     = ', rn_himin  
    713720         WRITE(numout,*)'      numerical carac. of the scheme for diffusion in ice ' 
    714          WRITE(numout,*)'      switch for snow sublimation  (=1) or not (=0)           parsub       = ', parsub   
    715721         WRITE(numout,*)'      coefficient for ice-lead partition of snowfall          rn_betas     = ', rn_betas 
    716722         WRITE(numout,*)'      extinction radiation parameter in sea ice               rn_kappa_i   = ', rn_kappa_i 
Note: See TracChangeset for help on using the changeset viewer.