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 4672 – NEMO

Changeset 4672


Ignore:
Timestamp:
2014-06-17T17:06:59+02:00 (10 years ago)
Author:
clem
Message:
 
Location:
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3
Files:
8 edited

Legend:

Unmodified
Added
Removed

  • branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/ice.F90

    r4659 r4672  
    167167   !                                          !!** ice-dynamic namelist (namicedyn) ** 
    168168   INTEGER , PUBLIC ::   nevp   = 400          !: number of iterations for subcycling 
    169    INTEGER , PUBLIC ::   nlay_i = 5            !: number of layers in the ice 
    170169 
    171170   !                                          !!** ice-dynamic namelist (namicedyn) ** 

  • branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90

    r4659 r4672  
    9191      REAL(wp) :: zinda, zindb, zareamin  
    9292      REAL(wp) :: zfric_u, zqld, zqfr 
    93       REAL(wp), POINTER, DIMENSION(:) :: zdq, zq_ini, zhfx, zqfx 
    94       REAL(wp)                        :: zhfx_err, ztest 
    9593      ! 
    9694      REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b  
    9795      !!------------------------------------------------------------------- 
    9896      IF( nn_timing == 1 )  CALL timing_start('limthd') 
    99  
    100       CALL wrk_alloc( jpij, zdq, zq_ini, zhfx, zqfx ) 
    101     
    102       ! init debug 
    103       zdq(:) = 0._wp ; zq_ini(:) = 0._wp ; zhfx(:) = 0._wp ; zqfx(:) = 0._wp       
    10497 
    10598      ! conservation test 
    … …  
    333326            ! 4.3) Thermodynamic processes 
    334327            !-------------------------------- 
    335             ! --- diag error on heat diffusion - PART 1 --- ! 
    336             DO ji = 1, nbpb 
    337                zq_ini(ji) = ( SUM( q_i_b(ji,1:nlay_i) ) * ht_i_b(ji) / REAL( nlay_i ) +  & 
    338                   &           SUM( q_s_b(ji,1:nlay_s) ) * ht_s_b(ji) / REAL( nlay_s ) )  
    339             END DO 
    340328 
    341329            !---------------------------------! 
    342330            ! Ice/Snow Temperature profile    ! 
    343331            !---------------------------------! 
    344             CALL lim_thd_dif( 1, nbpb, jl ) 
    345  
    346             ! --- computes sea ice energy of melting compulsory for limthd_dh --- ! 
    347             CALL lim_thd_enmelt( 1, nbpb ) 
    348  
    349             DO ji = 1, nbpb 
    350               ! --- diag error on heat diffusion - PART 2 --- ! 
    351                zdq(ji)        = - zq_ini(ji) + ( SUM( q_i_b(ji,1:nlay_i) ) * ht_i_b(ji) / REAL( nlay_i ) +  & 
    352                   &                              SUM( q_s_b(ji,1:nlay_s) ) * ht_s_b(ji) / REAL( nlay_s ) ) 
    353                zhfx_err       = ( fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - ftr_ice_1d(ji) - fc_bo_i(ji) + zdq(ji) * r1_rdtice )  
    354                hfx_err_1d(ji) = hfx_err_1d(ji) + zhfx_err * a_i_b(ji) 
    355                ! --- correction of qns_ice and surface conduction flux --- ! 
    356                qns_ice_1d(ji) = qns_ice_1d(ji) - zhfx_err  
    357                fc_su     (ji) = fc_su     (ji) - zhfx_err  
    358                ! --- Heat flux at the ice surface in W.m-2 --- ! 
    359                ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1 
    360                hfx_in (ii,ij) = hfx_in (ii,ij) + a_i_b(ji) * ( qsr_ice_1d(ji) + qns_ice_1d(ji) ) 
    361             END DO 
     332            CALL lim_thd_dif( 1, nbpb ) 
    362333 
    363334            !---------------------------------! 
    364335            ! Ice/Snow thicnkess              ! 
    365336            !---------------------------------! 
    366             ! --- diag error on heat remapping - PART 1 --- ! 
    367             DO ji = 1, nbpb 
    368                zq_ini(ji) = ( SUM( q_i_b(ji,1:nlay_i) ) * ht_i_b(ji) / REAL( nlay_i ) + & 
    369                   &           SUM( q_s_b(ji,1:nlay_s) ) * ht_s_b(ji) / REAL( nlay_s ) )  
    370             END DO 
    371  
    372             CALL lim_thd_dh( 1, nbpb, jl )     
     337            CALL lim_thd_dh( 1, nbpb )     
    373338 
    374339            ! --- Ice enthalpy remapping --- ! 
    375             CALL lim_thd_ent( 1, nbpb, jl, q_i_b(1:nbpb,:) )  
    376             !                                 
    377             ! --- diag error on heat remapping - PART 2 --- ! 
    378             DO ji = 1, nbpb 
    379                zdq(ji)        = - ( zq_ini(ji) + dq_i(ji) + dq_s(ji) )   & 
    380                   &             + ( SUM( q_i_b(ji,1:nlay_i) ) * ht_i_b(ji) / REAL( nlay_i ) +  & 
    381                   &                 SUM( q_s_b(ji,1:nlay_s) ) * ht_s_b(ji) / REAL( nlay_s ) ) 
    382                hfx_err_rem_1d(ji) = hfx_err_rem_1d(ji) + zdq(ji) * a_i_b(ji) * r1_rdtice 
    383             END DO 
    384  
     340            CALL lim_thd_ent( 1, nbpb, q_i_b(1:nbpb,:) )  
     341                                             
    385342            !---------------------------------! 
    386343            ! --- Ice salinity --- ! 
    … …  
    528485      IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limthd', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
    529486      ! 
    530       CALL wrk_dealloc( jpij, zdq, zq_ini, zhfx, zqfx ) 
    531  
    532487      IF( nn_timing == 1 )  CALL timing_stop('limthd') 
    533    END SUBROUTINE lim_thd 
    534  
    535   
    536    SUBROUTINE lim_thd_enmelt( kideb, kiut ) 
    537       !!----------------------------------------------------------------------- 
    538       !!                   ***  ROUTINE lim_thd_enmelt ***  
    539       !!                  
    540       !! ** Purpose :   Computes sea ice energy of melting q_i (J.m-3) from temperature 
    541       !! 
    542       !! ** Method  :   Formula (Bitz and Lipscomb, 1999) 
    543       !!------------------------------------------------------------------- 
    544       INTEGER, INTENT(in) ::   kideb, kiut   ! bounds for the spatial loop 
    545       !! 
    546       INTEGER  ::   ji, jk   ! dummy loop indices 
    547       REAL(wp) ::   ztmelts, zindb  ! local scalar  
    548       !!------------------------------------------------------------------- 
    549       ! 
    550       DO jk = 1, nlay_i             ! Sea ice energy of melting 
    551          DO ji = kideb, kiut 
    552             ztmelts      = - tmut  * s_i_b(ji,jk) + rtt  
    553             zindb        = MAX( 0._wp , SIGN( 1._wp , -(t_i_b(ji,jk) - rtt) - epsi10 ) ) 
    554             q_i_b(ji,jk) = rhoic * ( cpic * ( ztmelts - t_i_b(ji,jk) )                                             & 
    555                &                   + lfus * ( 1.0 - zindb * ( ztmelts-rtt ) / MIN( t_i_b(ji,jk)-rtt, -epsi10 ) )   & 
    556                &                   - rcp  *                 ( ztmelts-rtt )  )  
    557          END DO 
    558       END DO 
    559       DO jk = 1, nlay_s             ! Snow energy of melting 
    560          DO ji = kideb, kiut 
    561             q_s_b(ji,jk) = rhosn * ( cpic * ( rtt - t_s_b(ji,jk) ) + lfus ) 
    562          END DO 
    563       END DO 
    564       ! 
    565    END SUBROUTINE lim_thd_enmelt 
     488   END SUBROUTINE lim_thd  
    566489 
    567490   SUBROUTINE lim_thd_temp( kideb, kiut ) 

  • branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90

    r4649 r4672  
    4343CONTAINS 
    4444 
    45    SUBROUTINE lim_thd_dh( kideb, kiut, jl ) 
     45   SUBROUTINE lim_thd_dh( kideb, kiut ) 
    4646      !!------------------------------------------------------------------ 
    4747      !!                ***  ROUTINE lim_thd_dh  *** 
    … …  
    6868      !!------------------------------------------------------------------ 
    6969      INTEGER , INTENT(in) ::   kideb, kiut   ! Start/End point on which the  the computation is applied 
    70       INTEGER , INTENT(in) ::   jl            ! Thickness cateogry number 
    7170      !!  
    7271      INTEGER  ::   ji , jk        ! dummy loop indices 
    … …  
    151150      zintermelt(:)   = 0._wp 
    152151      icount    (:)   = 0 
    153  
    154       ! debug 
    155       dq_i(:) = 0._wp 
    156       dq_s(:) = 0._wp 
    157152 
    158153      ! initialize layer thicknesses and enthalpies 
    … …  
    273268         zh_s  (ji) = ht_s_b(ji) / REAL( nlay_s ) 
    274269 
    275          ! clem debug: variation of enthalpy (J.m-2) 
    276          dq_s(ji) = dq_s(ji) + ( zdh_s_pre(ji) + zdh_s_mel(ji) ) * zqprec(ji)   
    277270         ENDIF 
    278271      END DO 
    … …  
    297290            ht_s_b(ji) = MAX( 0._wp , ht_s_b(ji) + zdeltah(ji,jk) ) 
    298291 
    299             ! clem debug: variation of enthalpy (J.m-2) 
    300             dq_s(ji) = dq_s(ji) + zdeltah(ji,jk) * q_s_b(ji,jk)   
    301292         END DO 
    302293      END DO 
    … …  
    325316            ! new snow thickness 
    326317            ht_s_b(ji)     =  MAX( 0._wp , ht_s_b(ji) + zdh_s_sub(ji) ) 
    327             ! clem debug: variation of enthalpy (J.m-2) 
    328             dq_s(ji) = dq_s(ji) + zdh_s_sub(ji) * q_s_b(ji,1)   
    329318         END DO 
    330319      ENDIF 
    … …  
    397386            icount(ji)  = icount(ji) + zindh 
    398387            zh_i(ji,jk) = MAX( 0._wp , zh_i(ji,jk) + zdeltah(ji,jk) ) 
    399  
    400             ! clem debug: variation of enthalpy (J.m-2) 
    401             dq_i(ji) = dq_i(ji) + zdeltah(ji,jk) * q_i_b(ji,jk)   
    402388 
    403389            ! update heat content (J.m-2) and layer thickness 
    … …  
    513499            wfx_bog_1d(ji) =  wfx_bog_1d(ji) - rhoic * a_i_b(ji) * dh_i_bott(ji) * r1_rdtice 
    514500 
    515             ! clem debug: variation of enthalpy (J.m-2) 
    516             dq_i(ji) = dq_i(ji) + dh_i_bott(ji) * q_i_b(ji,nlay_i+1)   
    517  
    518501            ! update heat content (J.m-2) and layer thickness 
    519502            qh_i_old(ji,nlay_i+1) = qh_i_old(ji,nlay_i+1) + dh_i_bott(ji) * q_i_b(ji,nlay_i+1) 
    … …  
    557540                  ! Contribution to mass flux 
    558541                  wfx_res_1d(ji) =  wfx_res_1d(ji) - rhoic * a_i_b(ji) * zdeltah(ji,jk) * r1_rdtice 
    559  
    560                   ! clem debug: variation of enthalpy (J.m-2) 
    561                   dq_i(ji) = dq_i(ji) + zdeltah(ji,jk) * q_i_b(ji,jk)   
    562542 
    563543                  ! update heat content (J.m-2) and layer thickness 
    … …  
    598578                  ! Contribution to mass flux 
    599579                  wfx_bom_1d(ji) =  wfx_bom_1d(ji) - rhoic * a_i_b(ji) * zdeltah(ji,jk) * r1_rdtice 
    600  
    601                   ! clem debug: variation of enthalpy (J.m-2) 
    602                   dq_i(ji) = dq_i(ji) + zdeltah(ji,jk) * q_i_b(ji,jk)   
    603580 
    604581                  ! update heat content (J.m-2) and layer thickness 
    … …  
    664641!         ! Contribution to mass flux 
    665642!         wfx_snw_1d(ji)  =  wfx_snw_1d(ji) - rhosn * a_i_b(ji) * zdeltah(ji,1) * r1_rdtice 
    666 !         ! clem debug: variation of enthalpy (J.m-2) 
    667 !         dq_s(ji) = dq_s(ji) + zdeltah(ji,1) * q_s_b(ji,1)   
    668643!     
    669644         ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1 
    … …  
    721696         wfx_snw_1d(ji) = wfx_snw_1d(ji) + a_i_b(ji) * dh_snowice(ji) * rhosn * r1_rdtice 
    722697 
    723          ! clem debug: variation of enthalpy (J.m-2) 
    724          dq_s(ji) = dq_s(ji) - dh_snowice(ji) * q_s_b(ji,1)   
    725          dq_i(ji) = dq_i(ji) + dh_snowice(ji) * q_s_b(ji,1) + zfmdt * zEw   
    726  
    727698         ! update heat content (J.m-2) and layer thickness 
    728699         qh_i_old(ji,0) = qh_i_old(ji,0) + dh_snowice(ji) * q_s_b(ji,1) + zfmdt * zEw 

  • branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90

    r4649 r4672  
    4040CONTAINS 
    4141 
    42    SUBROUTINE lim_thd_dif( kideb , kiut , jl ) 
     42   SUBROUTINE lim_thd_dif( kideb , kiut ) 
    4343      !!------------------------------------------------------------------ 
    4444      !!                ***  ROUTINE lim_thd_dif  *** 
    … …  
    9393      !!------------------------------------------------------------------ 
    9494      INTEGER , INTENT(in) ::   kideb, kiut   ! Start/End point on which the  the computation is applied 
    95       INTEGER , INTENT(in) ::   jl            ! Thickness cateogry number 
    9695 
    9796      !! * Local variables 
    … …  
    146145      REAL(wp), POINTER, DIMENSION(:,:) ::   zdiagbis 
    147146      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ztrid   ! tridiagonal system terms 
     147      ! diag errors on heat 
     148      REAL(wp), POINTER, DIMENSION(:) :: zdq, zq_ini 
     149      REAL(wp)                        :: zhfx_err 
    148150      !!------------------------------------------------------------------      
    149151      !  
    … …  
    155157      CALL wrk_alloc( jpij, jkmax+2, zindterm, zindtbis, zdiagbis  ) 
    156158      CALL wrk_alloc( jpij, jkmax+2, 3, ztrid ) 
     159 
     160      CALL wrk_alloc( jpij, zdq, zq_ini ) 
     161 
     162      ! --- diag error on heat diffusion - PART 1 --- ! 
     163      zdq(:) = 0._wp ; zq_ini(:) = 0._wp       
     164      DO ji = kideb, kiut 
     165         zq_ini(ji) = ( SUM( q_i_b(ji,1:nlay_i) ) * ht_i_b(ji) / REAL( nlay_i ) +  & 
     166            &           SUM( q_s_b(ji,1:nlay_s) ) * ht_s_b(ji) / REAL( nlay_s ) )  
     167      END DO 
    157168 
    158169      !------------------------------------------------------------------------------! 
    … …  
    671682         DO layer  =  1, nlay_s 
    672683            DO ji = kideb , kiut 
    673                ii = MOD( npb(ji) - 1, jpi ) + 1 
    674                ij = ( npb(ji) - 1 ) / jpi + 1 
    675684               t_s_b(ji,layer) = MAX(  MIN( t_s_b(ji,layer), rtt ), 190._wp  ) 
    676685               zerrit(ji)      = MAX(zerrit(ji),ABS(t_s_b(ji,layer) - ztstemp(ji,layer))) 
    … …  
    722731         IF( t_su_b(ji) < rtt ) THEN  ! case T_su < 0degC 
    723732            hfx_dif_1d(ji) = hfx_dif_1d(ji) + ( qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_b(ji) 
    724          ELSE                                    ! case T_su = 0degC 
     733         ELSE                         ! case T_su = 0degC 
    725734            hfx_dif_1d(ji) = hfx_dif_1d(ji) + ( fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_b(ji) 
    726735         ENDIF 
    727736      END DO 
    728737 
     738      ! --- computes sea ice energy of melting compulsory for limthd_dh --- ! 
     739      CALL lim_thd_enmelt( kideb, kiut ) 
     740 
     741      ! --- diag error on heat diffusion - PART 2 --- ! 
     742      DO ji = kideb, kiut 
     743         zdq(ji)        = - zq_ini(ji) + ( SUM( q_i_b(ji,1:nlay_i) ) * ht_i_b(ji) / REAL( nlay_i ) +  & 
     744            &                              SUM( q_s_b(ji,1:nlay_s) ) * ht_s_b(ji) / REAL( nlay_s ) ) 
     745         zhfx_err    = ( fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq(ji) * r1_rdtice )  
     746         hfx_err_1d(ji) = hfx_err_1d(ji) + zhfx_err * a_i_b(ji) 
     747         ! --- correction of qns_ice and surface conduction flux --- ! 
     748         qns_ice_1d(ji) = qns_ice_1d(ji) - zhfx_err  
     749         fc_su     (ji) = fc_su     (ji) - zhfx_err  
     750         ! --- Heat flux at the ice surface in W.m-2 --- ! 
     751         ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1 
     752         hfx_in (ii,ij) = hfx_in (ii,ij) + a_i_b(ji) * ( qsr_ice_1d(ji) + qns_ice_1d(ji) ) 
     753      END DO 
     754    
    729755      ! 
    730756      CALL wrk_dealloc( jpij, numeqmin, numeqmax, isnow ) 
    … …  
    735761      CALL wrk_dealloc( jpij, jkmax+2, zindterm, zindtbis, zdiagbis ) 
    736762      CALL wrk_dealloc( jpij, jkmax+2, 3, ztrid ) 
     763      CALL wrk_dealloc( jpij, zdq, zq_ini ) 
    737764 
    738765   END SUBROUTINE lim_thd_dif 
     766 
     767   SUBROUTINE lim_thd_enmelt( kideb, kiut ) 
     768      !!----------------------------------------------------------------------- 
     769      !!                   ***  ROUTINE lim_thd_enmelt ***  
     770      !!                  
     771      !! ** Purpose :   Computes sea ice energy of melting q_i (J.m-3) from temperature 
     772      !! 
     773      !! ** Method  :   Formula (Bitz and Lipscomb, 1999) 
     774      !!------------------------------------------------------------------- 
     775      INTEGER, INTENT(in) ::   kideb, kiut   ! bounds for the spatial loop 
     776      ! 
     777      INTEGER  ::   ji, jk   ! dummy loop indices 
     778      REAL(wp) ::   ztmelts, zindb  ! local scalar  
     779      !!------------------------------------------------------------------- 
     780      ! 
     781      DO jk = 1, nlay_i             ! Sea ice energy of melting 
     782         DO ji = kideb, kiut 
     783            ztmelts      = - tmut  * s_i_b(ji,jk) + rtt  
     784            zindb        = MAX( 0._wp , SIGN( 1._wp , -(t_i_b(ji,jk) - rtt) - epsi10 ) ) 
     785            q_i_b(ji,jk) = rhoic * ( cpic * ( ztmelts - t_i_b(ji,jk) )                                             & 
     786               &                   + lfus * ( 1.0 - zindb * ( ztmelts-rtt ) / MIN( t_i_b(ji,jk)-rtt, -epsi10 ) )   & 
     787               &                   - rcp  *                 ( ztmelts-rtt )  )  
     788         END DO 
     789      END DO 
     790      DO jk = 1, nlay_s             ! Snow energy of melting 
     791         DO ji = kideb, kiut 
     792            q_s_b(ji,jk) = rhosn * ( cpic * ( rtt - t_s_b(ji,jk) ) + lfus ) 
     793         END DO 
     794      END DO 
     795      ! 
     796   END SUBROUTINE lim_thd_enmelt 
    739797 
    740798#else 

  • branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_ent.F90

    r4649 r4672  
    3636   PRIVATE 
    3737 
    38    PUBLIC   lim_thd_ent         ! called by lim_thd 
     38   PUBLIC   lim_thd_ent         ! called by limthd and limthd_lac 
    3939 
    4040   REAL(wp) :: epsi20 = 1.e-20   ! constant values 
    … …  
    4848CONTAINS 
    4949  
    50    SUBROUTINE lim_thd_ent( kideb, kiut, jl, qnew ) 
     50   SUBROUTINE lim_thd_ent( kideb, kiut, qnew ) 
    5151      !!------------------------------------------------------------------- 
    5252      !!               ***   ROUTINE lim_thd_ent  *** 
    … …  
    7474      !!------------------------------------------------------------------- 
    7575      INTEGER , INTENT(in) ::   kideb, kiut   ! Start/End point on which the  the computation is applied 
    76       INTEGER , INTENT(in) ::   jl            ! Thickness cateogry number 
    7776 
    78       REAL(wp), INTENT(inout), DIMENSION(:,:) :: qnew          ! new enthlapies (remapped) 
     77      REAL(wp), INTENT(inout), DIMENSION(:,:) :: qnew          ! new enthlapies (J.m-3, remapped) 
    7978 
    80       INTEGER  :: ji,ii,ij   !  dummy loop indices 
     79      INTEGER  :: ji         !  dummy loop indices 
    8180      INTEGER  :: jk0, jk1   !  old/new layer indices 
    82       REAL(wp) :: ztmelts    ! temperature of melting 
    83       REAL(wp) :: zswitch, zaaa, zbbb, zccc, zdiscrim ! converting enthalpy to temperature 
     81      REAL(wp) :: zswitch 
    8482      ! 
    8583      REAL(wp), POINTER, DIMENSION(:,:) :: zqh_cum0, zh_cum0   ! old cumulative enthlapies and layers interfaces 
    … …  
    139137      DO jk1 = 1, nlay_i 
    140138         DO ji = kideb, kiut 
    141             zswitch      =  1._wp - MAX( 0._wp , SIGN( 1._wp , - zhnew(ji) + epsi10 ) )  
     139            zswitch      = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zhnew(ji) + epsi10 ) )  
    142140            qnew(ji,jk1) = zswitch * ( zqh_cum1(ji,jk1) - zqh_cum1(ji,jk1-1) ) / MAX( zhnew(ji), epsi10 ) 
    143141         ENDDO 
    144142      ENDDO 
     143 
     144      ! --- diag error on heat remapping --- ! 
     145      ! comment: if input h_i_old and qh_i_old are already multiplied by a_i (as in limthd_lac),  
     146      ! then we should not (* a_i) again but not important since this is just to check that remap error is ~0 
     147      DO ji = kideb, kiut 
     148         hfx_err_rem_1d(ji) = hfx_err_rem_1d(ji) + a_i_b(ji) * r1_rdtice *  & 
     149            &               ( SUM( qnew(ji,1:nlay_i) ) * zhnew(ji) - SUM( qh_i_old(ji,0:nlay_i+1) ) )  
     150      END DO 
     151       
    145152      ! 
    146153      CALL wrk_dealloc( jpij, nlay_i+3, zqh_cum0, zh_cum0, kjstart = 0 ) 

  • branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90

    r4659 r4672  
    329329         ! Salinity of new ice  
    330330         !---------------------- 
    331  
    332331         SELECT CASE ( num_sal ) 
    333332         CASE ( 1 )                    ! Sice = constant  
    … …  
    343342         END SELECT 
    344343 
    345  
    346344         !------------------------- 
    347345         ! Heat content of new ice 
    … …  
    354352               &                       - rcp  *         ( ztmelts - rtt )  ) 
    355353         END DO ! ji 
     354 
    356355         !---------------- 
    357356         ! Age of new ice 
    … …  
    395394         END DO 
    396395 
    397  
    398396         !----------------- 
    399397         ! Area of new ice 
    … …  
    407405         !------------------------------------------------------------------------------! 
    408406 
    409          !------------------------------------------- 
    410          ! Compute excessive new ice area and volume 
    411          !------------------------------------------- 
     407         !------------------------ 
     408         ! 6.1) lateral ice growth 
     409         !------------------------ 
    412410         ! If lateral ice growth gives an ice concentration gt 1, then 
    413411         ! we keep the excessive volume in memory and attribute it later to bottom accretion 
    … …  
    422420               zdv_res(ji) = 0._wp 
    423421            ENDIF 
    424          END DO ! ji 
    425  
    426          !------------------------------------------------ 
    427          ! Laterally redistribute new ice volume and area 
    428          !------------------------------------------------ 
     422         END DO 
     423 
     424         ! find which category to fill 
    429425         zat_i_1d(:) = 0._wp 
    430426         DO jl = 1, jpl 
    … …  
    433429                  za_i_1d (ji,jl) = za_i_1d (ji,jl) + za_newice(ji) 
    434430                  zv_i_1d (ji,jl) = zv_i_1d (ji,jl) + zv_newice(ji) 
    435                   jcat  (ji)    = jl 
     431                  jcat    (ji)    = jl 
    436432               ENDIF 
    437433               zat_i_1d(ji) = zat_i_1d(ji) + za_i_1d  (ji,jl) 
    … …  
    439435         END DO 
    440436 
    441          !---------------------------------- 
    442          ! Heat content - lateral accretion 
    443          !---------------------------------- 
    444          DO ji = 1, nbpac 
    445             jl = jcat(ji)                                                  ! categroy in which new ice is put 
     437         ! Heat content 
     438         DO ji = 1, nbpac 
     439            jl = jcat(ji)                                                    ! categroy in which new ice is put 
    446440            zswinew  (ji) = MAX( 0._wp , SIGN( 1._wp , - za_old(ji,jl) ) )   ! 0 if old ice 
    447441         END DO 
    … …  
    457451         END DO 
    458452 
    459          !----------------------------------------------- 
    460          ! Add excessive volume of new ice at the bottom 
    461          !----------------------------------------------- 
     453         !------------------------------------------------ 
     454         ! 6.2) bottom ice growth + ice enthalpy remapping 
     455         !------------------------------------------------ 
    462456         DO jl = 1, jpl 
     457 
     458            ! for remapping 
    463459            h_i_old (1:nbpac,0:nlay_i+1) = 0._wp 
    464460            qh_i_old(1:nbpac,0:nlay_i+1) = 0._wp 
    465  
    466461            DO jk = 1, nlay_i 
    467462               DO ji = 1, nbpac 
    … …  
    471466            END DO 
    472467 
     468            ! new volumes including lateral/bottom accretion + residual 
    473469            DO ji = 1, nbpac 
    474470               zinda          = MAX( 0._wp, SIGN( 1._wp , zat_i_1d(ji) - epsi20 ) ) 
    475471               zv_newfra      = zinda * ( zdv_res(ji) + zv_frazb(ji) ) * za_i_1d(ji,jl) / MAX( zat_i_1d(ji) , epsi20 ) 
    476472               za_i_1d(ji,jl) = zinda * za_i_1d(ji,jl)                
    477  
     473               zv_i_1d(ji,jl) = zv_i_1d(ji,jl) + zv_newfra 
     474 
     475               ! for remapping 
    478476               h_i_old (ji,nlay_i+1) = zv_newfra 
    479477               qh_i_old(ji,nlay_i+1) = ze_newice(ji) * zv_newfra 
    480  
    481                zv_i_1d(ji,jl) = zv_i_1d(ji,jl) + zv_newfra 
    482478            ENDDO 
    483479 
    484480            ! --- Ice enthalpy remapping --- ! 
    485             CALL lim_thd_ent( 1, nbpac, jl, ze_i_1d(1:nbpac,:,jl) )  
     481            IF( zv_newfra > 0._wp ) THEN 
     482               CALL lim_thd_ent( 1, nbpac, ze_i_1d(1:nbpac,:,jl) )  
     483            ENDIF 
    486484 
    487485         ENDDO 
    … …  
    500498         ! Update salinity 
    501499         !----------------- 
    502          !clem IF(  num_sal == 2  ) THEN 
    503             DO jl = 1, jpl 
    504                DO ji = 1, nbpac 
    505                   zdv   = zv_i_1d(ji,jl) - zv_old(ji,jl) 
    506                   zsmv_i_1d(ji,jl) = zsmv_i_1d(ji,jl) + zdv * zs_newice(ji) 
    507                END DO 
    508             END DO    
    509          !clem ENDIF 
     500         DO jl = 1, jpl 
     501            DO ji = 1, nbpac 
     502               zdv   = zv_i_1d(ji,jl) - zv_old(ji,jl) 
     503               zsmv_i_1d(ji,jl) = zsmv_i_1d(ji,jl) + zdv * zs_newice(ji) 
     504            END DO 
     505         END DO 
    510506 
    511507         !------------------------------------------------------------------------------! 
    512          ! 8) Change 2D vectors to 1D vectors  
     508         ! 7) Change 2D vectors to 1D vectors  
    513509         !------------------------------------------------------------------------------! 
    514510         DO jl = 1, jpl 
    … …  
    531527 
    532528      !------------------------------------------------------------------------------! 
    533       ! 9) Change units for e_i 
     529      ! 8) Change units for e_i 
    534530      !------------------------------------------------------------------------------!     
    535531      DO jl = 1, jpl 

  • branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/par_ice.F90

    r2528 r4672  
    1313   !                                             !!! ice thermodynamics 
    1414   INTEGER, PUBLIC, PARAMETER ::   jkmax    = 6   !: maximum number of ice layers 
     15   INTEGER, PUBLIC, PARAMETER ::   nlay_i   = 5   !: number of ice layers 
    1516   INTEGER, PUBLIC, PARAMETER ::   nlay_s   = 1   !: number of snow layers 
    1617 

  • branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90

    r4659 r4672  
    124124   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   q_s_b   !:    Snow enthalpy per unit volume 
    125125 
    126    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dq_i   !: variation of ice enthalpy (debug) 
    127    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dq_s   !: variation of snw enthalpy (debug) 
    128  
    129126   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   qh_i_old  !: ice heat content (q*h, J.m-2) 
    130127   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   h_i_old   !: ice thickness layer (m) 
    … …  
    171168         &      dh_s_tot  (jpij) , dh_i_surf(jpij) , dh_i_bott(jpij) ,    &     
    172169         &      dh_snowice(jpij) , sm_i_b   (jpij) , s_i_new  (jpij) ,    & 
    173          &      dq_i      (jpij) , dq_s     (jpij),  t_s_b(jpij,nlay_s),  & 
     170         &      t_s_b(jpij,nlay_s),                                       & 
    174171         &      t_i_b(jpij,jkmax), s_i_b(jpij,jkmax)                ,     &             
    175172         &      q_i_b(jpij,jkmax), q_s_b(jpij,jkmax)                ,     & 
Note: See TracChangeset for help on using the changeset viewer.