# Changeset 8355

Ignore:
Timestamp:
2017-07-20T10:20:49+02:00 (3 years ago)
Message:

simplifications

Location:
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3
Files:
5 edited

Unmodified
Removed
• ## branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limcons.F90

 r8341 PRIVATE PUBLIC   lim_column_sum PUBLIC   lim_column_sum_energy PUBLIC   lim_cons_check PUBLIC   lim_cons_hsm PUBLIC   lim_cons_final !!---------------------------------------------------------------------- CONTAINS SUBROUTINE lim_column_sum( ksum, pin, pout ) !!------------------------------------------------------------------- !!               ***  ROUTINE lim_column_sum *** !! !! ** Purpose : Compute the sum of xin over nsum categories !! !! ** Method  : Arithmetics !! !! ** Action  : Gets xin(ji,jj,jl) and computes xout(ji,jj) !!--------------------------------------------------------------------- INTEGER                   , INTENT(in   ) ::   ksum   ! number of categories/layers REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pin    ! input field REAL(wp), DIMENSION(:,:)  , INTENT(  out) ::   pout   ! output field ! INTEGER ::   jl   ! dummy loop indices !!--------------------------------------------------------------------- ! pout(:,:) = pin(:,:,1) DO jl = 2, ksum pout(:,:) = pout(:,:) + pin(:,:,jl) END DO ! END SUBROUTINE lim_column_sum SUBROUTINE lim_column_sum_energy( ksum, klay, pin, pout) !!------------------------------------------------------------------- !!               ***  ROUTINE lim_column_sum_energy *** !! !! ** Purpose : Compute the sum of xin over nsum categories !!              and nlay layers !! !! ** Method  : Arithmetics !!--------------------------------------------------------------------- INTEGER                                , INTENT(in   ) ::   ksum   !: number of categories INTEGER                                , INTENT(in   ) ::   klay   !: number of vertical layers REAL(wp), DIMENSION(jpi,jpj,nlay_i,jpl), INTENT(in   ) ::   pin    !: input field REAL(wp), DIMENSION(jpi,jpj)           , INTENT(  out) ::   pout   !: output field ! INTEGER ::   jk, jl   ! dummy loop indices !!--------------------------------------------------------------------- ! pout(:,:) = 0._wp DO jl = 1, ksum DO jk = 2, klay pout(:,:) = pout(:,:) + pin(:,:,jk,jl) END DO END DO ! END SUBROUTINE lim_column_sum_energy SUBROUTINE lim_cons_check( px1, px2, pmax_err, cd_fieldid ) !!------------------------------------------------------------------- !!               ***  ROUTINE lim_cons_check *** !! !! ** Purpose : Test the conservation of a certain variable !!              For each physical grid cell, check that initial !!              and final values !!              of a conserved field are equal to within a small value. !! !! ** Method  : !!--------------------------------------------------------------------- REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   px1          !: initial field REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   px2          !: final field REAL(wp)                , INTENT(in   ) ::   pmax_err     !: max allowed error CHARACTER(len=15)       , INTENT(in   ) ::   cd_fieldid   !: field identifyer ! INTEGER  ::   ji, jj          ! dummy loop indices INTEGER  ::   inb_error       ! number of g.c where there is a cons. error LOGICAL  ::   llconserv_err   ! = .true. if conservation check failed REAL(wp) ::   zmean_error     ! mean error on error points !!--------------------------------------------------------------------- ! IF(lwp) WRITE(numout,*) ' lim_cons_check ' IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~~ ' llconserv_err = .FALSE. inb_error     = 0 zmean_error   = 0._wp IF( MAXVAL( px2(:,:) - px1(:,:) ) > pmax_err )   llconserv_err = .TRUE. IF( llconserv_err ) THEN DO jj = 1, jpj DO ji = 1, jpi IF( ABS( px2(ji,jj) - px1(ji,jj) ) > pmax_err ) THEN inb_error   = inb_error + 1 zmean_error = zmean_error + ABS( px2(ji,jj) - px1(ji,jj) ) ! IF(lwp) THEN WRITE (numout,*) ' ALERTE 99 ' WRITE (numout,*) ' Conservation error: ', cd_fieldid WRITE (numout,*) ' Point             : ', ji, jj WRITE (numout,*) ' lat, lon          : ', gphit(ji,jj), glamt(ji,jj) WRITE (numout,*) ' Initial value     : ', px1(ji,jj) WRITE (numout,*) ' Final value       : ', px2(ji,jj) WRITE (numout,*) ' Difference        : ', px2(ji,jj) - px1(ji,jj) ENDIF ENDIF END DO END DO ! ENDIF IF(lk_mpp)   CALL mpp_sum( inb_error   ) IF(lk_mpp)   CALL mpp_sum( zmean_error ) ! IF( inb_error > 0 .AND. lwp ) THEN zmean_error = zmean_error / REAL( inb_error, wp ) WRITE(numout,*) ' Conservation check for : ', cd_fieldid WRITE(numout,*) ' Number of error points : ', inb_error WRITE(numout,*) ' Mean error on these pts: ', zmean_error ENDIF ! END SUBROUTINE lim_cons_check SUBROUTINE lim_cons_hsm( icount, cd_routine, zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b )
• ## branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90

 r8233 CONTAINS SUBROUTINE lim_itd_th_rem( klbnd, kubnd, kt ) SUBROUTINE lim_itd_th_rem( kt ) !!------------------------------------------------------------------ !!                ***  ROUTINE lim_itd_th_rem *** !! References : W.H. Lipscomb, JGR 2001 !!------------------------------------------------------------------ INTEGER , INTENT (in) ::   klbnd   ! Start thickness category index point INTEGER , INTENT (in) ::   kubnd   ! End point on which the  the computation is applied INTEGER , INTENT (in) ::   kt      ! Ocean time step ! REAL(wp) ::   zx2, zwk2, zda0, zetamax           !   -      - REAL(wp) ::   zx3 CHARACTER (len = 15) :: fieldid INTEGER , POINTER, DIMENSION(:,:,:) ::   zdonor   ! donor category index REAL(wp), POINTER, DIMENSION(:,:,:) ::   hR          ! left boundary for the ITD for each thickness REAL(wp), POINTER, DIMENSION(:,:,:) ::   zht_i_b     ! old ice thickness REAL(wp), POINTER, DIMENSION(:,:,:) ::   dummy_es REAL(wp), POINTER, DIMENSION(:,:,:) ::   zdaice, zdvice          ! local increment of ice area and volume REAL(wp), POINTER, DIMENSION(:)     ::   zvetamin, zvetamax      ! maximum values for etas INTEGER , POINTER, DIMENSION(:)     ::   nind_i, nind_j          ! compressed indices for i/j directions INTEGER                             ::   nbrem                   ! number of cells with ice to transfer REAL(wp)                            ::   zslope                  ! used to compute local thermodynamic "speeds" REAL(wp), POINTER, DIMENSION(:,:)   ::   zhb0, zhb1              ! category boundaries for thinnes categories REAL(wp), POINTER, DIMENSION(:,:)   ::   vt_i_init, vt_i_final   !  ice volume summed over categories REAL(wp), POINTER, DIMENSION(:,:)   ::   vt_s_init, vt_s_final   !  snow volume summed over categories REAL(wp), POINTER, DIMENSION(:,:)   ::   et_i_init, et_i_final   !  ice energy summed over categories REAL(wp), POINTER, DIMENSION(:,:)   ::   et_s_init, et_s_final   !  snow energy summed over categories INTEGER , POINTER, DIMENSION(:,:)   ::   zremap_flag      ! compute remapping or not ???? REAL(wp), POINTER, DIMENSION(:,:,:) ::   zhbnew           ! new boundaries of ice categories CALL wrk_alloc( jpi,jpj, zremap_flag ) CALL wrk_alloc( jpi,jpj,jpl-1, zdonor ) CALL wrk_alloc( jpi,jpj,jpl, zdhice, g0, g1, hL, hR, zht_i_b, dummy_es ) CALL wrk_alloc( jpi,jpj,jpl, zdhice, g0, g1, hL, hR, zht_i_b ) CALL wrk_alloc( jpi,jpj,jpl-1, zdaice, zdvice ) CALL wrk_alloc( jpi,jpj,jpl+1, zhbnew, kkstart = 0 ) CALL wrk_alloc( (jpi+1)*(jpj+1), zvetamin, zvetamax ) CALL wrk_alloc( (jpi+1)*(jpj+1), nind_i, nind_j ) CALL wrk_alloc( jpi,jpj, zhb0,zhb1,vt_i_init,vt_i_final,vt_s_init,vt_s_final,et_i_init,et_i_final,et_s_init,et_s_final ) !!---------------------------------------------------------------------------------------------- !! 0) Conservation checkand changes in each ice category !!---------------------------------------------------------------------------------------------- IF( con_i ) THEN CALL lim_column_sum (jpl,   v_i, vt_i_init) CALL lim_column_sum (jpl,   v_s, vt_s_init) CALL lim_column_sum_energy (jpl, nlay_i,   e_i, et_i_init) dummy_es(:,:,:) = e_s(:,:,1,:) CALL lim_column_sum (jpl, dummy_es(:,:,:) , et_s_init) ENDIF CALL wrk_alloc( jpi,jpj, zhb0, zhb1 ) !!---------------------------------------------------------------------------------------------- WRITE(numout,*) 'lim_itd_th_rem  : Remapping the ice thickness distribution' WRITE(numout,*) '~~~~~~~~~~~~~~~' WRITE(numout,*) ' klbnd :       ', klbnd WRITE(numout,*) ' kubnd :       ', kubnd ENDIF zdhice(:,:,:) = 0._wp DO jl = klbnd, kubnd DO jl = 1, jpl DO jj = 1, jpj DO ji = 1, jpi END DO !----------------------------------------------------------------------------------------------- !  2) Compute fractional ice area in each grid cell !----------------------------------------------------------------------------------------------- at_i(:,:) = 0._wp DO jl = klbnd, kubnd at_i(:,:) = at_i(:,:) + a_i(:,:,jl) END DO !!      !----------------------------------------------------------------------------------------------- !!      !  2) Compute fractional ice area in each grid cell !!      !----------------------------------------------------------------------------------------------- !!      at_i(:,:) = 0._wp !!      DO jl = 1, jpl !!         at_i(:,:) = at_i(:,:) + a_i(:,:,jl) !!      END DO !----------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------- !  4) Compute new category boundaries !----------------------------------------------------------------------------------------------- !- 4.1 Compute category boundaries !  4) Compute new category boundaries: 1:jpl-1 !----------------------------------------------------------------------------------------------- zhbnew(:,:,:) = 0._wp DO jl = klbnd, kubnd - 1 DO ji = 1, nbrem ii = nind_i(ji) ij = nind_j(ji) ! zhbnew(ii,ij,jl) = hi_max(jl) IF    ( a_i_b(ii,ij,jl) > epsi10 .AND. a_i_b(ii,ij,jl+1) > epsi10 ) THEN !interpolate between adjacent category growth rates zslope           = ( zdhice(ii,ij,jl+1) - zdhice(ii,ij,jl) ) / ( zht_i_b(ii,ij,jl+1) - zht_i_b(ii,ij,jl) ) zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl) + zslope * ( hi_max(jl) - zht_i_b(ii,ij,jl) ) ELSEIF( a_i_b(ii,ij,jl) > epsi10) THEN zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl) ELSEIF( a_i_b(ii,ij,jl+1) > epsi10) THEN zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl+1) ENDIF END DO !- 4.2 Check that each zhbnew lies between adjacent values of ice thickness DO ji = 1, nbrem ii = nind_i(ji) ij = nind_j(ji) ! clem: we do not want ht_i to be too close to either HR or HL otherwise a division by nearly 0 is possible ! in lim_itd_fitline in the case (HR-HL) = 3(Hice - HL) or = 3(HR - Hice) IF    ( a_i(ii,ij,jl  ) > epsi10 .AND. ht_i(ii,ij,jl  ) > ( zhbnew(ii,ij,jl) - epsi10 ) ) THEN zremap_flag(ii,ij) = 0 ELSEIF( a_i(ii,ij,jl+1) > epsi10 .AND. ht_i(ii,ij,jl+1) < ( zhbnew(ii,ij,jl) + epsi10 ) ) THEN zremap_flag(ii,ij) = 0 ENDIF !- 4.3 Check that each zhbnew does not exceed maximal values hi_max IF( zhbnew(ii,ij,jl) < hi_max(jl-1) ) zremap_flag(ii,ij) = 0 IF( zhbnew(ii,ij,jl) > hi_max(jl+1) ) zremap_flag(ii,ij) = 0 ! clem bug: why is not the following instead? !!IF( zhbnew(ii,ij,jl) < hi_max(jl-1) ) zremap_flag(ii,ij) = 0 !!IF( zhbnew(ii,ij,jl) > hi_max(jl  ) ) zremap_flag(ii,ij) = 0 END DO END DO IF( nbrem > 0 ) THEN DO jl = 1, jpl - 1 DO ji = 1, nbrem ii = nind_i(ji) ij = nind_j(ji) ! ! --- New boundary categories Hn* = Hn + Fn*dt --- !!    Fn*dt = ( fn + (fn+1 - fn)/(hn+1 - hn) * (Hn - hn) ) * dt = zdhice + zslope * (Hmax - ht_i_b) IF    ( a_i_b(ii,ij,jl) >  epsi10 .AND. a_i_b(ii,ij,jl+1) >  epsi10 ) THEN !! a_i_b(jl+1) & a_i_b(jl) /= 0 zslope           = ( zdhice(ii,ij,jl+1) - zdhice(ii,ij,jl) ) / ( zht_i_b(ii,ij,jl+1) - zht_i_b(ii,ij,jl) ) zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl) + zslope * ( hi_max(jl) - zht_i_b(ii,ij,jl) ) !!               ELSEIF( a_i_b(ii,ij,jl) >  epsi10 .AND. a_i_b(ii,jj,jl+1) <= epsi10 ) THEN !! a_i_b(jl+1)=0 => Hn* = Hn + fn*dt ELSEIF( a_i_b(ii,ij,jl) >  epsi10 ) THEN !! a_i_b(jl+1)=0 => Hn* = Hn + fn*dt zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl) !!               ELSEIF( a_i_b(ii,ij,jl) <= epsi10 .AND. a_i_b(ii,ij,jl+1) >  epsi10 ) THEN !! a_i_b(jl)=0 => Hn* = Hn + fn+1*dt ELSEIF( a_i_b(ii,ij,jl+1) > epsi10 ) THEN !! a_i_b(jl)=0 => Hn* = Hn + fn+1*dt zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl+1) ELSE                                                                       !! a_i_b(jl+1) & a_i_b(jl) = 0 zhbnew(ii,ij,jl) = hi_max(jl) ENDIF ! --- 2 conditions for remapping --- ! ! 1) hn(t+1)+espi < Hn* < hn+1(t+1)-epsi !    Note: hn(t+1) must not be too close to either HR or HL otherwise a division by nearly 0 is possible !          in lim_itd_fitline in the case (HR-HL) = 3(Hice - HL) or = 3(HR - Hice) IF( a_i(ii,ij,jl  ) > epsi10 .AND. ht_i(ii,ij,jl  ) > ( zhbnew(ii,ij,jl) - epsi10 ) )   zremap_flag(ii,ij) = 0 IF( a_i(ii,ij,jl+1) > epsi10 .AND. ht_i(ii,ij,jl+1) < ( zhbnew(ii,ij,jl) + epsi10 ) )   zremap_flag(ii,ij) = 0 ! 2) Hn-1 < Hn* < Hn+1 IF( zhbnew(ii,ij,jl) < hi_max(jl-1) )   zremap_flag(ii,ij) = 0 IF( zhbnew(ii,ij,jl) > hi_max(jl+1) )   zremap_flag(ii,ij) = 0 END DO END DO ENDIF !----------------------------------------------------------------------------------------------- !  5) Identify cells where ITD is to be remapped !----------------------------------------------------------------------------------------------- !  6) Fill arrays with lowermost / uppermost boundaries of 'new' categories !  6) Compute new category boundaries: jpl !----------------------------------------------------------------------------------------------- DO jj = 1, jpj zhb1(ji,jj) = hi_max(1) IF( a_i(ji,jj,kubnd) > epsi10 ) THEN zhbnew(ji,jj,kubnd) = MAX( hi_max(kubnd-1), 3._wp * ht_i(ji,jj,kubnd) - 2._wp * zhbnew(ji,jj,kubnd-1) ) ! boundary for jpl IF( a_i(ji,jj,jpl) > epsi10 ) THEN zhbnew(ji,jj,jpl) = MAX( hi_max(jpl-1), 3._wp * ht_i(ji,jj,jpl) - 2._wp * zhbnew(ji,jj,jpl-1) ) ELSE !clem bug               zhbnew(ji,jj,kubnd) = hi_max(kubnd) zhbnew(ji,jj,kubnd) = hi_max(kubnd-1) ! not used anyway ENDIF ! clem: we do not want ht_i_b to be too close to either HR or HL otherwise a division by nearly 0 is possible ! in lim_itd_fitline in the case (HR-HL) = 3(Hice - HL) or = 3(HR - Hice) IF    ( zht_i_b(ji,jj,klbnd) < ( zhb0(ji,jj) + epsi10 ) )  THEN zremap_flag(ji,jj) = 0 ELSEIF( zht_i_b(ji,jj,klbnd) > ( zhb1(ji,jj) - epsi10 ) )  THEN zremap_flag(ji,jj) = 0 ENDIF !clem bug               zhbnew(ji,jj,jpl) = hi_max(jpl) zhbnew(ji,jj,jpl) = hi_max(jpl-1) ! not used anyway ENDIF ! 1 condition for remapping for the 1st category ! H0+epsi < h1(t) < H1-epsi !    h1(t) must not be too close to either HR or HL otherwise a division by nearly 0 is possible !    in lim_itd_fitline in the case (HR-HL) = 3(Hice - HL) or = 3(HR - Hice) IF( zht_i_b(ji,jj,1) < ( zhb0(ji,jj) + epsi10 ) )   zremap_flag(ji,jj) = 0 IF( zht_i_b(ji,jj,1) > ( zhb1(ji,jj) - epsi10 ) )   zremap_flag(ji,jj) = 0 END DO !----------------------------------------------------------------------------------------------- !- 7.1 g(h) for category 1 at start of time step CALL lim_itd_fitline( klbnd, zhb0, zhb1, zht_i_b(:,:,klbnd), g0(:,:,klbnd), g1(:,:,klbnd), hL(:,:,klbnd),   & &                  hR(:,:,klbnd), zremap_flag ) !- 7.2 Area lost due to melting of thin ice (first category,  klbnd) CALL lim_itd_fitline( 1, zhb0, zhb1, zht_i_b(:,:,1), g0(:,:,1), g1(:,:,1), hL(:,:,1),   & &                  hR(:,:,1), zremap_flag ) !- 7.2 Area lost due to melting of thin ice (first category,  1) DO ji = 1, nbrem ii = nind_i(ji) ij = nind_j(ji) IF( a_i(ii,ij,klbnd) > epsi10 ) THEN zdh0 = zdhice(ii,ij,klbnd) !decrease of ice thickness in the lower category IF( a_i(ii,ij,1) > epsi10 ) THEN zdh0 = zdhice(ii,ij,1) !decrease of ice thickness in the lower category IF( zdh0 < 0.0 ) THEN      !remove area from category 1 zdh0 = MIN( -zdh0, hi_max(klbnd) ) zdh0 = MIN( -zdh0, hi_max(1) ) !Integrate g(1) from 0 to dh0 to estimate area melted zetamax = MIN( zdh0, hR(ii,ij,klbnd) ) - hL(ii,ij,klbnd) zetamax = MIN( zdh0, hR(ii,ij,1) ) - hL(ii,ij,1) IF( zetamax > 0.0 ) THEN zx1    = zetamax zx2    = 0.5 * zetamax * zetamax zda0   = g1(ii,ij,klbnd) * zx2 + g0(ii,ij,klbnd) * zx1                        ! ice area removed zdamax = a_i(ii,ij,klbnd) * (1.0 - ht_i(ii,ij,klbnd) / zht_i_b(ii,ij,klbnd) ) ! Constrain new thickness <= ht_i zda0   = g1(ii,ij,1) * zx2 + g0(ii,ij,1) * zx1                        ! ice area removed zdamax = a_i(ii,ij,1) * (1.0 - ht_i(ii,ij,1) / zht_i_b(ii,ij,1) ) ! Constrain new thickness <= ht_i zda0   = MIN( zda0, zdamax )                                                  ! ice area lost due to melting !     of thin ice (zdamax > 0) ! Remove area, conserving volume ht_i(ii,ij,klbnd) = ht_i(ii,ij,klbnd) * a_i(ii,ij,klbnd) / ( a_i(ii,ij,klbnd) - zda0 ) a_i(ii,ij,klbnd)  = a_i(ii,ij,klbnd) - zda0 v_i(ii,ij,klbnd)  = a_i(ii,ij,klbnd) * ht_i(ii,ij,klbnd) ! clem-useless ? ht_i(ii,ij,1) = ht_i(ii,ij,1) * a_i(ii,ij,1) / ( a_i(ii,ij,1) - zda0 ) a_i(ii,ij,1)  = a_i(ii,ij,1) - zda0 v_i(ii,ij,1)  = a_i(ii,ij,1) * ht_i(ii,ij,1) ! clem-useless ? ENDIF ELSE ! if ice accretion zdh0 > 0 ! zhbnew was 0, and is shifted to the right to account for thin ice growth in openwater (F0 = f1) zhbnew(ii,ij,klbnd-1) = MIN( zdh0, hi_max(klbnd) ) zhbnew(ii,ij,0) = MIN( zdh0, hi_max(1) ) ENDIF !- 7.3 g(h) for each thickness category DO jl = klbnd, kubnd DO jl = 1, jpl CALL lim_itd_fitline( jl, zhbnew(:,:,jl-1), zhbnew(:,:,jl), ht_i(:,:,jl),  & &                  g0(:,:,jl), g1(:,:,jl), hL(:,:,jl), hR(:,:,jl), zremap_flag ) !----------------------------------------------------------------------------------------------- !  8) Compute area and volume to be shifted across each boundary !----------------------------------------------------------------------------------------------- DO jl = klbnd, kubnd - 1 !  8) Compute area and volume to be shifted across each boundary (Eq. 18) !----------------------------------------------------------------------------------------------- DO jl = 1, jpl - 1 DO jj = 1, jpj DO ji = 1, jpi ij = nind_j(ji) IF (zhbnew(ii,ij,jl) > hi_max(jl)) THEN ! transfer from jl to jl+1 ! left and right integration limits in eta space zvetamin(ji) = MAX( hi_max(jl), hL(ii,ij,jl) ) - hL(ii,ij,jl) zvetamax(ji) = MIN( zhbnew(ii,ij,jl), hR(ii,ij,jl) ) - hL(ii,ij,jl) ! left and right integration limits in eta space IF (zhbnew(ii,ij,jl) > hi_max(jl)) THEN ! Hn* > Hn => transfer from jl to jl+1 zetamin = MAX( hi_max(jl), hL(ii,ij,jl) ) - hL(ii,ij,jl)       ! hi_max(jl) - hL zetamax = MIN( zhbnew(ii,ij,jl), hR(ii,ij,jl) ) - hL(ii,ij,jl) ! hR - hL zdonor(ii,ij,jl) = jl ELSE                                    ! zhbnew(jl) <= hi_max(jl) ; transfer from jl+1 to jl ! left and right integration limits in eta space zvetamin(ji) = 0.0 zvetamax(ji) = MIN( hi_max(jl), hR(ii,ij,jl+1) ) - hL(ii,ij,jl+1) ELSE                                    ! Hn* <= Hn => transfer from jl+1 to jl zetamin = 0.0 zetamax = MIN( hi_max(jl), hR(ii,ij,jl+1) ) - hL(ii,ij,jl+1)   ! hi_max(jl) - hL zdonor(ii,ij,jl) = jl + 1 ENDIF zetamax = MAX( zvetamax(ji), zvetamin(ji) ) ! no transfer if etamax < etamin zetamin = zvetamin(ji) ENDIF zetamax = MAX( zetamax, zetamin ) ! no transfer if etamax < etamin zx1  = zetamax - zetamin !! 9) Shift ice between categories !!---------------------------------------------------------------------------------------------- CALL lim_itd_shiftice ( klbnd, kubnd, zdonor, zdaice, zdvice ) CALL lim_itd_shiftice ( zdonor, zdaice, zdvice ) !!---------------------------------------------------------------------------------------------- END DO !!---------------------------------------------------------------------------------------------- !! 11) Conservation check !!---------------------------------------------------------------------------------------------- IF ( con_i ) THEN CALL lim_column_sum (jpl,   v_i, vt_i_final) fieldid = ' v_i : limitd_th ' CALL lim_cons_check (vt_i_init, vt_i_final, 1.0e-6, fieldid) CALL lim_column_sum_energy (jpl, nlay_i,  e_i, et_i_final) fieldid = ' e_i : limitd_th ' CALL lim_cons_check (et_i_init, et_i_final, 1.0e-3, fieldid) CALL lim_column_sum (jpl,   v_s, vt_s_final) fieldid = ' v_s : limitd_th ' CALL lim_cons_check (vt_s_init, vt_s_final, 1.0e-6, fieldid) dummy_es(:,:,:) = e_s(:,:,1,:) CALL lim_column_sum (jpl, dummy_es(:,:,:) , et_s_final) fieldid = ' e_s : limitd_th ' CALL lim_cons_check (et_s_init, et_s_final, 1.0e-3, fieldid) ENDIF CALL wrk_dealloc( jpi,jpj, zremap_flag ) CALL wrk_dealloc( jpi,jpj,jpl-1, zdonor ) CALL wrk_dealloc( jpi,jpj,jpl, zdhice, g0, g1, hL, hR, zht_i_b, dummy_es ) CALL wrk_dealloc( jpi,jpj,jpl, zdhice, g0, g1, hL, hR, zht_i_b ) CALL wrk_dealloc( jpi,jpj,jpl-1, zdaice, zdvice ) CALL wrk_dealloc( jpi,jpj,jpl+1, zhbnew, kkstart = 0 ) CALL wrk_dealloc( (jpi+1)*(jpj+1), zvetamin, zvetamax ) CALL wrk_dealloc( (jpi+1)*(jpj+1), nind_i, nind_j ) CALL wrk_dealloc( jpi,jpj, zhb0,zhb1,vt_i_init,vt_i_final,vt_s_init,vt_s_final,et_i_init,et_i_final,et_s_init,et_s_final ) CALL wrk_dealloc( jpi,jpj, zhb0, zhb1 ) END SUBROUTINE lim_itd_th_rem !!                ***  ROUTINE lim_itd_fitline *** !! !! ** Purpose :   fit g(h) with a line using area, volume constraints !! !! ** Method  :   Fit g(h) with a line, satisfying area and volume constraints. !!              To reduce roundoff errors caused by large values of g0 and g1, !!              we actually compute g(eta), where eta = h - hL, and hL is the !!              left boundary. !! ** Purpose :   build g(h) satisfying area and volume constraints (Eq. 6 and 9) !! !! ** Method  :   g(h) is linear and written as: g(eta) = g1(eta) + g0 !!                with eta = h - HL !! !!------------------------------------------------------------------ INTEGER                     , INTENT(in   ) ::   num_cat      ! category index hR(ji,jj) = HbR(ji,jj) ! Change hL or hR if hice falls outside central third of range zh13 = 1.0 / 3.0 * ( 2.0 * hL(ji,jj) + hR(ji,jj) ) zh23 = 1.0 / 3.0 * ( hL(ji,jj) + 2.0 * hR(ji,jj) ) IF    ( hice(ji,jj) < zh13 ) THEN   ;   hR(ji,jj) = 3._wp * hice(ji,jj) - 2._wp * hL(ji,jj) ELSEIF( hice(ji,jj) > zh23 ) THEN   ;   hL(ji,jj) = 3._wp * hice(ji,jj) - 2._wp * hR(ji,jj) ! Change hL or hR if hice falls outside central third of range, ! so that hice is in the central third of the range [HL HR] zh13 = 1.0 / 3.0 * ( 2.0 * hL(ji,jj) +       hR(ji,jj) ) zh23 = 1.0 / 3.0 * (       hL(ji,jj) + 2.0 * hR(ji,jj) ) IF    ( hice(ji,jj) < zh13 ) THEN   ;   hR(ji,jj) = 3._wp * hice(ji,jj) - 2._wp * hL(ji,jj) ! move HR to the left ELSEIF( hice(ji,jj) > zh23 ) THEN   ;   hL(ji,jj) = 3._wp * hice(ji,jj) - 2._wp * hR(ji,jj) ! move HL to the right ENDIF zwk1 = 6._wp * a_i(ji,jj,num_cat) * zdhr zwk2 = ( hice(ji,jj) - hL(ji,jj) ) * zdhr g0(ji,jj) = zwk1 * ( 2._wp / 3._wp - zwk2 ) g1(ji,jj) = 2._wp * zdhr * zwk1 * ( zwk2 - 0.5 ) g0(ji,jj) = zwk1 * ( 2._wp / 3._wp - zwk2 )      ! Eq. 14 g1(ji,jj) = 2._wp * zdhr * zwk1 * ( zwk2 - 0.5 ) ! Eq. 14 ! ELSE  ! remap_flag = .false. or a_i < epsi10 SUBROUTINE lim_itd_shiftice( klbnd, kubnd, zdonor, zdaice, zdvice ) SUBROUTINE lim_itd_shiftice( zdonor, zdaice, zdvice ) !!------------------------------------------------------------------ !!                ***  ROUTINE lim_itd_shiftice *** !! ** Method  : !!------------------------------------------------------------------ INTEGER                           , INTENT(in   ) ::   klbnd    ! Start thickness category index point INTEGER                           , INTENT(in   ) ::   kubnd    ! End point on which the  the computation is applied INTEGER , DIMENSION(jpi,jpj,jpl-1), INTENT(in   ) ::   zdonor   ! donor category index REAL(wp), DIMENSION(jpi,jpj,jpl-1), INTENT(inout) ::   zdaice   ! ice area transferred across boundary !---------------------------------------------------------------------------------------------- DO jl = klbnd, kubnd DO jl = 1, jpl zaTsfn(:,:,jl) = a_i(:,:,jl) * t_su(:,:,jl) END DO !------------------------------------------------------------------------------- DO jl = klbnd, kubnd - 1 DO jl = 1, jpl - 1 nbrem = 0 DO jj = 1, jpj ENDIF !-------------- ! Ice areas !-------------- a_i(ii,ij,jl1) = a_i(ii,ij,jl1) - zdaice(ii,ij,jl) a_i(ii,ij,jl2) = a_i(ii,ij,jl2) + zdaice(ii,ij,jl) !-------------- ! Ice volumes !-------------- v_i(ii,ij,jl1) = v_i(ii,ij,jl1) - zdvice(ii,ij,jl) v_i(ii,ij,jl2) = v_i(ii,ij,jl2) + zdvice(ii,ij,jl) !-------------- ! Snow volumes !-------------- zdvsnow        = v_s(ii,ij,jl1) * zworka(ii,ij) v_s(ii,ij,jl1) = v_s(ii,ij,jl1) - zdvsnow v_s(ii,ij,jl2) = v_s(ii,ij,jl2) + zdvsnow !-------------------- ! Snow heat content !-------------------- zdesnow            = e_s(ii,ij,1,jl1) * zworka(ii,ij) e_s(ii,ij,1,jl1)   = e_s(ii,ij,1,jl1) - zdesnow e_s(ii,ij,1,jl2)   = e_s(ii,ij,1,jl2) + zdesnow !-------------- ! Ice age !-------------- zdo_aice           = oa_i(ii,ij,jl1) * zdaice(ii,ij,jl) oa_i(ii,ij,jl1)    = oa_i(ii,ij,jl1) - zdo_aice oa_i(ii,ij,jl2)    = oa_i(ii,ij,jl2) + zdo_aice !-------------- ! Ice salinity !-------------- zdsm_vice          = smv_i(ii,ij,jl1) * zworka(ii,ij) smv_i(ii,ij,jl1)   = smv_i(ii,ij,jl1) - zdsm_vice smv_i(ii,ij,jl2)   = smv_i(ii,ij,jl2) + zdsm_vice !--------------------- ! Surface temperature !--------------------- zdaTsf             = t_su(ii,ij,jl1) * zdaice(ii,ij,jl) zaTsfn(ii,ij,jl1)  = zaTsfn(ii,ij,jl1) - zdaTsf ! MV MP 2016 IF ( nn_pnd_scheme > 0 ) THEN !--------------------- ! Pond fraction !--------------------- zdapnd             = a_ip(ii,ij,jl1) * zdaice(ii,ij,jl) a_ip(ii,ij,jl1)    = a_ip(ii,ij,jl1) - zdapnd a_ip(ii,ij,jl2)    = a_ip(ii,ij,jl2) + zdapnd !--------------------- ! Pond volume !--------------------- zdvpnd             = v_ip(ii,ij,jl1) * zdaice(ii,ij,jl) v_ip(ii,ij,jl1)    = v_ip(ii,ij,jl1) - zdvpnd END DO !------------------ ! Ice heat content !------------------ DO jk = 1, nlay_i DO ji = 1, nbrem !----------------------------------------------------------------- DO jl = klbnd, kubnd DO jl = 1, jpl DO jj = 1, jpj DO ji = 1, jpi SUBROUTINE lim_itd_th_reb( klbnd, kubnd ) SUBROUTINE lim_itd_th_reb !!------------------------------------------------------------------ !!                ***  ROUTINE lim_itd_th_reb *** !! ** Method  : !!------------------------------------------------------------------ INTEGER , INTENT (in) ::   klbnd   ! Start thickness category index point INTEGER , INTENT (in) ::   kubnd   ! End point on which the  the computation is applied ! INTEGER ::   ji,jj, jl   ! dummy loop indices INTEGER ::   zshiftflag          ! = .true. if ice must be shifted CHARACTER (len = 15) :: fieldid INTEGER , POINTER, DIMENSION(:,:,:) ::   zdonor           ! donor category index REAL(wp), POINTER, DIMENSION(:,:,:) ::   zdaice, zdvice   ! ice area and volume transferred REAL(wp), POINTER, DIMENSION(:,:) ::   vt_i_init, vt_i_final   ! ice volume summed over categories REAL(wp), POINTER, DIMENSION(:,:) ::   vt_s_init, vt_s_final   ! snow volume summed over categories !!------------------------------------------------------------------ CALL wrk_alloc( jpi,jpj,jpl, zdonor )   ! interger CALL wrk_alloc( jpi,jpj,jpl, zdaice, zdvice ) CALL wrk_alloc( jpi,jpj, vt_i_init, vt_i_final, vt_s_init, vt_s_final ) ! IF( con_i ) THEN                 ! conservation check CALL lim_column_sum (jpl,   v_i, vt_i_init) CALL lim_column_sum (jpl,   v_s, vt_s_init) ENDIF ! !------------------------------------------------------------------------------ ! 1) Compute ice thickness. !------------------------------------------------------------------------------ DO jl = klbnd, kubnd DO jl = 1, jpl DO jj = 1, jpj DO ji = 1, jpi ! Initialize shift arrays !------------------------- DO jl = klbnd, kubnd DO jl = 1, jpl zdonor(:,:,jl) = 0 zdaice(:,:,jl) = 0._wp !------------------------- DO jl = klbnd, kubnd - 1  ! loop over category boundaries DO jl = 1, jpl - 1  ! loop over category boundaries !--------------------------------------- IF( zshiftflag == 1 ) THEN            ! Shift ice between categories CALL lim_itd_shiftice( klbnd, kubnd, zdonor, zdaice, zdvice ) CALL lim_itd_shiftice( zdonor, zdaice, zdvice ) ! Reset shift parameters zdonor(:,:,jl) = 0 !---------------------------- DO jl = kubnd - 1, 1, -1       ! loop over category boundaries DO jl = jpl - 1, 1, -1       ! loop over category boundaries !----------------------------------------- IF( zshiftflag == 1 ) THEN            ! Shift ice between categories CALL lim_itd_shiftice( klbnd, kubnd, zdonor, zdaice, zdvice ) CALL lim_itd_shiftice( zdonor, zdaice, zdvice ) ! Reset shift parameters zdonor(:,:,jl) = 0 END DO !------------------------------------------------------------------------------ ! 3) Conservation check !------------------------------------------------------------------------------ IF( con_i ) THEN CALL lim_column_sum (jpl,   v_i, vt_i_final) fieldid = ' v_i : limitd_reb ' CALL lim_cons_check (vt_i_init, vt_i_final, 1.0e-6, fieldid) CALL lim_column_sum (jpl,   v_s, vt_s_final) fieldid = ' v_s : limitd_reb ' CALL lim_cons_check (vt_s_init, vt_s_final, 1.0e-6, fieldid) ENDIF ! CALL wrk_dealloc( jpi,jpj,jpl, zdonor ) CALL wrk_dealloc( jpi,jpj,jpl, zdaice, zdvice ) CALL wrk_dealloc( jpi,jpj, vt_i_init, vt_i_final, vt_s_init, vt_s_final ) END SUBROUTINE lim_itd_th_reb
• ## branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90

 r8342 IF( ln_limdiachk ) CALL lim_cons_hsm(0, 'limitd_th_rem', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) IF( jpl > 1 )      CALL lim_itd_th_rem( 1, jpl, kt ) IF( jpl > 1 )      CALL lim_itd_th_rem( kt ) IF( ln_limdiachk ) CALL lim_cons_hsm(1, 'limitd_th_rem', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
• ## branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90

 r8316 ! Rebin categories with thickness out of bounds !---------------------------------------------------- IF ( jpl > 1 ) CALL lim_itd_th_reb(1, jpl) IF ( jpl > 1 ) CALL lim_itd_th_reb !-----------------
• ## branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90

 r8316 ! Rebin categories with thickness out of bounds !---------------------------------------------------- IF ( jpl > 1 ) CALL lim_itd_th_reb( 1, jpl ) IF ( jpl > 1 ) CALL lim_itd_th_reb !-----------------
Note: See TracChangeset for help on using the changeset viewer.