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 11081 for NEMO/branches/UKMO/NEMO_4.0_mirror/src/ICE/icedyn_rdgrft.F90 – NEMO

Ignore:
Timestamp:
2019-06-06T16:11:54+02:00 (5 years ago)
Author:
davestorkey
Message:

UKMO/NEMO_4.0_mirror : update to be a copy of rev 11079 of release-4.0.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/UKMO/NEMO_4.0_mirror/src/ICE/icedyn_rdgrft.F90

    r10888 r11081  
    142142      INTEGER, PARAMETER ::   jp_itermax = 20     
    143143      !!------------------------------------------------------------------- 
    144       ! clem: The redistribution of ice between categories can lead to small negative values (as for the remapping in ice_itd_rem) 
    145       !       likely due to truncation error ( i.e. 1. - 1. /= 0 ) 
    146       !       I do not think it should be a concern since small areas and volumes are erased (in ice_var_zapsmall.F90) 
    147        
    148144      ! controls 
    149145      IF( ln_timing    )   CALL timing_start('icedyn_rdgrft')                                                             ! timing 
     
    156152      ENDIF       
    157153 
    158       CALL ice_var_zapsmall   ! Zero out categories with very small areas 
    159  
    160154      !-------------------------------- 
    161155      ! 0) Identify grid cells with ice 
    162156      !-------------------------------- 
     157      at_i(:,:) = SUM( a_i, dim=3 ) 
     158      ! 
    163159      npti = 0   ;   nptidx(:) = 0 
    164160      ipti = 0   ;   iptidx(:) = 0 
    165161      DO jj = 1, jpj 
    166162         DO ji = 1, jpi 
    167             IF ( at_i(ji,jj) > 0._wp ) THEN 
     163            IF ( at_i(ji,jj) > epsi10 ) THEN 
    168164               npti           = npti + 1 
    169165               nptidx( npti ) = (jj - 1) * jpi + ji 
     
    178174         
    179175         ! just needed here 
    180          CALL tab_2d_1d( npti, nptidx(1:npti), zdivu(1:npti), divu_i(:,:) ) 
    181          CALL tab_2d_1d( npti, nptidx(1:npti), zdelt(1:npti), delta_i(:,:) ) 
     176         CALL tab_2d_1d( npti, nptidx(1:npti), zdivu   (1:npti)      , divu_i ) 
     177         CALL tab_2d_1d( npti, nptidx(1:npti), zdelt   (1:npti)      , delta_i ) 
    182178         ! needed here and in the iteration loop 
    183          CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d  (1:npti,1:jpl), a_i(:,:,:) ) 
    184          CALL tab_3d_2d( npti, nptidx(1:npti), v_i_2d  (1:npti,1:jpl), v_i(:,:,:) ) 
    185          CALL tab_2d_1d( npti, nptidx(1:npti), ato_i_1d(1:npti)      , ato_i(:,:) ) 
     179         CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d  (1:npti,1:jpl), a_i  ) 
     180         CALL tab_3d_2d( npti, nptidx(1:npti), v_i_2d  (1:npti,1:jpl), v_i  ) 
     181         CALL tab_2d_1d( npti, nptidx(1:npti), ato_i_1d(1:npti)      , ato_i ) 
    186182 
    187183         DO ji = 1, npti 
     
    310306 
    311307      !                       ! Ice thickness needed for rafting 
    312       WHERE( pa_i(1:npti,:) > 0._wp )   ;   zhi(1:npti,:) = pv_i(1:npti,:) / pa_i(1:npti,:) 
    313       ELSEWHERE                         ;   zhi(1:npti,:) = 0._wp 
     308      WHERE( pa_i(1:npti,:) > epsi20 )   ;   zhi(1:npti,:) = pv_i(1:npti,:) / pa_i(1:npti,:) 
     309      ELSEWHERE                          ;   zhi(1:npti,:) = 0._wp 
    314310      END WHERE 
    315311 
     
    329325      zasum(1:npti) = pato_i(1:npti) + SUM( pa_i(1:npti,:), dim=2 ) 
    330326      ! 
    331       WHERE( zasum(1:npti) > 0._wp )   ;   z1_asum(1:npti) = 1._wp / zasum(1:npti) 
    332       ELSEWHERE                        ;   z1_asum(1:npti) = 0._wp 
     327      WHERE( zasum(1:npti) > epsi20 )   ;   z1_asum(1:npti) = 1._wp / zasum(1:npti) 
     328      ELSEWHERE                         ;   z1_asum(1:npti) = 0._wp 
    333329      END WHERE 
    334330      ! 
     
    455451      ! Based on the ITD of ridging and ridged ice, convert the net closing rate to a gross closing rate.   
    456452      ! NOTE: 0 < aksum <= 1 
    457       WHERE( zaksum(1:npti) > 0._wp )   ;   closing_gross(1:npti) = pclosing_net(1:npti) / zaksum(1:npti) 
    458       ELSEWHERE                         ;   closing_gross(1:npti) = 0._wp 
     453      WHERE( zaksum(1:npti) > epsi20 )   ;   closing_gross(1:npti) = pclosing_net(1:npti) / zaksum(1:npti) 
     454      ELSEWHERE                          ;   closing_gross(1:npti) = 0._wp 
    459455      END WHERE 
    460456       
     
    466462         DO ji = 1, npti 
    467463            zfac = apartf(ji,jl) * closing_gross(ji) * rdt_ice 
    468             IF( zfac > pa_i(ji,jl) ) THEN 
     464            IF( zfac > pa_i(ji,jl) .AND. apartf(ji,jl) /= 0._wp ) THEN 
    469465               closing_gross(ji) = pa_i(ji,jl) / apartf(ji,jl) * r1_rdtice 
    470466            ENDIF 
     
    510506      REAL(wp), DIMENSION(jpij) ::   zswitch, fvol    ! new ridge volume going to jl2 
    511507      REAL(wp), DIMENSION(jpij) ::   z1_ai            ! 1 / a 
     508      REAL(wp), DIMENSION(jpij) ::   zvti             ! sum(v_i) 
    512509      ! 
    513510      REAL(wp), DIMENSION(jpij,nlay_s) ::   esrft     ! snow energy of rafting ice 
     
    518515      INTEGER , DIMENSION(jpij) ::   itest_rdg, itest_rft   ! test for conservation 
    519516      !!------------------------------------------------------------------- 
    520  
     517      ! 
     518      zvti(1:npti) = SUM( v_i_2d(1:npti,:), dim=2 )   ! total ice volume 
     519      ! 
    521520      ! 1) Change in open water area due to closing and opening 
    522521      !-------------------------------------------------------- 
     
    535534            IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN   ! only if ice is ridging 
    536535 
    537                z1_ai(ji) = 1._wp / a_i_2d(ji,jl1) 
    538  
     536               IF( a_i_2d(ji,jl1) > epsi20 ) THEN   ;   z1_ai(ji) = 1._wp / a_i_2d(ji,jl1) 
     537               ELSE                                 ;   z1_ai(ji) = 0._wp 
     538               ENDIF 
     539                
    539540               ! area of ridging / rafting ice (airdg1) and of new ridge (airdg2) 
    540541               airdg1 = aridge(ji,jl1) * closing_gross(ji) * rdt_ice 
     
    549550 
    550551               ! volume and enthalpy (J/m2, >0) of seawater trapped into ridges 
    551                vsw = v_i_2d(ji,jl1) * afrdg * rn_porordg 
     552               IF    ( zvti(ji) <= 10. ) THEN ; vsw = v_i_2d(ji,jl1) * afrdg * rn_porordg                                           ! v <= 10m then porosity = rn_porordg 
     553               ELSEIF( zvti(ji) >= 20. ) THEN ; vsw = 0._wp                                                                         ! v >= 20m then porosity = 0 
     554               ELSE                           ; vsw = v_i_2d(ji,jl1) * afrdg * rn_porordg * MAX( 0._wp, 2._wp - 0.1_wp * zvti(ji) ) ! v > 10m and v < 20m then porosity = linear transition to 0 
     555               ENDIF 
    552556               ersw(ji) = -rhoi * vsw * rcp * sst_1d(ji)   ! clem: if sst>0, then ersw <0 (is that possible?) 
    553557 
    554558               ! volume etc of ridging / rafting ice and new ridges (vi, vs, sm, oi, es, ei) 
    555559               virdg1     = v_i_2d (ji,jl1)   * afrdg 
    556                virdg2(ji) = v_i_2d (ji,jl1)   * afrdg * ( 1. + rn_porordg ) 
     560               virdg2(ji) = v_i_2d (ji,jl1)   * afrdg + vsw 
    557561               vsrdg(ji)  = v_s_2d (ji,jl1)   * afrdg 
    558562               sirdg1     = sv_i_2d(ji,jl1)   * afrdg 
     
    726730      END DO ! jl1 
    727731      ! 
     732      ! roundoff errors 
     733      !---------------- 
    728734      ! In case ridging/rafting lead to very small negative values (sometimes it happens) 
    729       WHERE( a_i_2d(1:npti,:) < 0._wp )   a_i_2d(1:npti,:) = 0._wp 
    730       WHERE( v_i_2d(1:npti,:) < 0._wp )   v_i_2d(1:npti,:) = 0._wp 
     735      CALL ice_var_roundoff( a_i_2d, v_i_2d, v_s_2d, sv_i_2d, oa_i_2d, a_ip_2d, v_ip_2d, ze_s_2d, ze_i_2d ) 
    731736      ! 
    732737   END SUBROUTINE rdgrft_shift 
     
    854859         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_snw_dyn_1d(1:npti), wfx_snw_dyn(:,:) ) 
    855860         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d    (1:npti), wfx_pnd    (:,:) ) 
    856  
     861         ! 
    857862         !                 !---------------------! 
    858863      CASE( 2 )            !==  from 1D to 2D  ==! 
     
    945950         CALL ctl_stop( 'ice_dyn_rdgrft_init: choose one and only one participation function (ln_partf_lin or ln_partf_exp)' ) 
    946951      ENDIF 
    947       !                              ! allocate tke arrays 
     952      ! 
     953      IF( .NOT. ln_icethd ) THEN 
     954         rn_porordg = 0._wp 
     955         rn_fsnwrdg = 1._wp ; rn_fsnwrft = 1._wp 
     956         rn_fpndrdg = 1._wp ; rn_fpndrft = 1._wp 
     957         IF( lwp ) THEN 
     958            WRITE(numout,*) '      ==> only ice dynamics is activated, thus some parameters must be changed' 
     959            WRITE(numout,*) '            rn_porordg   = ', rn_porordg 
     960            WRITE(numout,*) '            rn_fsnwrdg   = ', rn_fsnwrdg  
     961            WRITE(numout,*) '            rn_fpndrdg   = ', rn_fpndrdg  
     962            WRITE(numout,*) '            rn_fsnwrft   = ', rn_fsnwrft  
     963            WRITE(numout,*) '            rn_fpndrft   = ', rn_fpndrft  
     964         ENDIF 
     965      ENDIF 
     966      !                              ! allocate arrays 
    948967      IF( ice_dyn_rdgrft_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'ice_dyn_rdgrft_init: unable to allocate arrays' ) 
    949968      ! 
Note: See TracChangeset for help on using the changeset viewer.