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 11822 for NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/ICE/icedyn_rdgrft.F90 – NEMO

Ignore:
Timestamp:
2019-10-29T11:41:36+01:00 (4 years ago)
Author:
acc
Message:

Branch 2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps. Sette tested updates to branch to align with trunk changes between 10721 and 11740. Sette tests are passing but results differ from branch before these changes (except for GYRE_PISCES and VORTEX) and branch results already differed from trunk because of algorithmic fixes. Will need more checks to confirm correctness.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/ICE/icedyn_rdgrft.F90

    r10531 r11822  
    8686      !!                ***  ROUTINE ice_dyn_rdgrft_alloc *** 
    8787      !!------------------------------------------------------------------- 
    88       ALLOCATE( closing_net(jpij), opning(jpij)   , closing_gross(jpij),   & 
    89          &      apartf(jpij,0:jpl), hrmin(jpij,jpl), hraft(jpij,jpl)    , aridge(jpij,jpl),  & 
    90          &      hrmax(jpij,jpl), hi_hrdg(jpij,jpl)  , araft (jpij,jpl),  & 
     88      ALLOCATE( closing_net(jpij)  , opning(jpij)      , closing_gross(jpij) ,               & 
     89         &      apartf(jpij,0:jpl) , hrmin  (jpij,jpl) , hraft(jpij,jpl) , aridge(jpij,jpl), & 
     90         &      hrmax (jpij,jpl)   , hi_hrdg(jpij,jpl) , araft(jpij,jpl) ,                   & 
    9191         &      ze_i_2d(jpij,nlay_i,jpl), ze_s_2d(jpij,nlay_s,jpl), STAT=ice_dyn_rdgrft_alloc ) 
    9292 
     
    137137      REAL(wp) ::   zfac                       ! local scalar 
    138138      INTEGER , DIMENSION(jpij) ::   iptidx        ! compute ridge/raft or not 
    139       REAL(wp), DIMENSION(jpij) ::   zdivu_adv     ! divu as implied by transport scheme  (1/s) 
    140139      REAL(wp), DIMENSION(jpij) ::   zdivu, zdelt  ! 1D divu_i & delta_i 
    141140      ! 
    142141      INTEGER, PARAMETER ::   jp_itermax = 20     
    143142      !!------------------------------------------------------------------- 
    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        
    148143      ! controls 
    149144      IF( ln_timing    )   CALL timing_start('icedyn_rdgrft')                                                             ! timing 
    150145      IF( ln_icediachk )   CALL ice_cons_hsm(0, 'icedyn_rdgrft', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 
     146      IF( ln_icediachk )   CALL ice_cons2D  (0, 'icedyn_rdgrft',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft) ! conservation 
    151147 
    152148      IF( kt == nit000 ) THEN 
     
    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), zdelt   (1:npti)      , delta_i ) 
    182177         ! 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(:,:) ) 
     178         CALL tab_2d_1d( npti, nptidx(1:npti), zdivu   (1:npti)      , divu_i) ! zdivu is used as a work array here (no change in divu_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 
     
    190186            closing_net(ji) = rn_csrdg * 0.5_wp * ( zdelt(ji) - ABS( zdivu(ji) ) ) - MIN( zdivu(ji), 0._wp ) 
    191187            ! 
    192             ! divergence given by the advection scheme 
    193             !   (which may not be equal to divu as computed from the velocity field) 
    194             IF    ( ln_adv_Pra ) THEN 
    195                zdivu_adv(ji) = ( 1._wp - ato_i_1d(ji) - SUM( a_i_2d(ji,:) ) ) * r1_rdtice 
    196             ELSEIF( ln_adv_UMx ) THEN 
    197                zdivu_adv(ji) = zdivu(ji) 
    198             ENDIF 
    199             ! 
    200             IF( zdivu_adv(ji) < 0._wp )   closing_net(ji) = MAX( closing_net(ji), -zdivu_adv(ji) )   ! make sure the closing rate is large enough 
    201             !                                                                                        ! to give asum = 1.0 after ridging 
     188            IF( zdivu(ji) < 0._wp )   closing_net(ji) = MAX( closing_net(ji), -zdivu(ji) )   ! make sure the closing rate is large enough 
     189            !                                                                                ! to give asum = 1.0 after ridging 
    202190            ! Opening rate (non-negative) that will give asum = 1.0 after ridging. 
    203             opning(ji) = closing_net(ji) + zdivu_adv(ji) 
     191            opning(ji) = closing_net(ji) + zdivu(ji) 
    204192         END DO 
    205193         ! 
     
    218206               ato_i_1d   (ipti)   = ato_i_1d   (ji) 
    219207               closing_net(ipti)   = closing_net(ji) 
    220                zdivu_adv  (ipti)   = zdivu_adv  (ji) 
     208               zdivu      (ipti)   = zdivu      (ji) 
    221209               opning     (ipti)   = opning     (ji) 
    222210            ENDIF 
     
    262250               ELSE 
    263251                  iterate_ridging  = 1 
    264                   zdivu_adv  (ji) = zfac * r1_rdtice 
    265                   closing_net(ji) = MAX( 0._wp, -zdivu_adv(ji) ) 
    266                   opning     (ji) = MAX( 0._wp,  zdivu_adv(ji) ) 
     252                  zdivu      (ji) = zfac * r1_rdtice 
     253                  closing_net(ji) = MAX( 0._wp, -zdivu(ji) ) 
     254                  opning     (ji) = MAX( 0._wp,  zdivu(ji) ) 
    267255               ENDIF 
    268256            END DO 
     
    280268 
    281269      ! controls 
     270      IF( ln_ctl       )   CALL ice_prt3D   ('icedyn_rdgrft')                                                             ! prints 
     271      IF( ln_icectl    )   CALL ice_prt     (kt, iiceprt, jiceprt,-1, ' - ice dyn rdgrft - ')                             ! prints 
    282272      IF( ln_icediachk )   CALL ice_cons_hsm(1, 'icedyn_rdgrft', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 
    283       IF( ln_ctl       )   CALL ice_prt3D   ('icedyn_rdgrft')                                                             ! prints 
     273      IF( ln_icediachk )   CALL ice_cons2D  (1, 'icedyn_rdgrft',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft) ! conservation 
    284274      IF( ln_timing    )   CALL timing_stop ('icedyn_rdgrft')                                                             ! timing 
    285275      ! 
     
    310300 
    311301      !                       ! 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 
     302      WHERE( pa_i(1:npti,:) > epsi10 )   ;   zhi(1:npti,:) = pv_i(1:npti,:) / pa_i(1:npti,:) 
     303      ELSEWHERE                          ;   zhi(1:npti,:) = 0._wp 
    314304      END WHERE 
    315305 
     
    329319      zasum(1:npti) = pato_i(1:npti) + SUM( pa_i(1:npti,:), dim=2 ) 
    330320      ! 
    331       WHERE( zasum(1:npti) > 0._wp )   ;   z1_asum(1:npti) = 1._wp / zasum(1:npti) 
    332       ELSEWHERE                        ;   z1_asum(1:npti) = 0._wp 
     321      WHERE( zasum(1:npti) > epsi10 )   ;   z1_asum(1:npti) = 1._wp / zasum(1:npti) 
     322      ELSEWHERE                         ;   z1_asum(1:npti) = 0._wp 
    333323      END WHERE 
    334324      ! 
     
    455445      ! Based on the ITD of ridging and ridged ice, convert the net closing rate to a gross closing rate.   
    456446      ! 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 
     447      WHERE( zaksum(1:npti) > epsi10 )   ;   closing_gross(1:npti) = pclosing_net(1:npti) / zaksum(1:npti) 
     448      ELSEWHERE                          ;   closing_gross(1:npti) = 0._wp 
    459449      END WHERE 
    460450       
     
    466456         DO ji = 1, npti 
    467457            zfac = apartf(ji,jl) * closing_gross(ji) * rdt_ice 
    468             IF( zfac > pa_i(ji,jl) ) THEN 
     458            IF( zfac > pa_i(ji,jl) .AND. apartf(ji,jl) /= 0._wp ) THEN 
    469459               closing_gross(ji) = pa_i(ji,jl) / apartf(ji,jl) * r1_rdtice 
    470460            ENDIF 
     
    510500      REAL(wp), DIMENSION(jpij) ::   zswitch, fvol    ! new ridge volume going to jl2 
    511501      REAL(wp), DIMENSION(jpij) ::   z1_ai            ! 1 / a 
     502      REAL(wp), DIMENSION(jpij) ::   zvti             ! sum(v_i) 
    512503      ! 
    513504      REAL(wp), DIMENSION(jpij,nlay_s) ::   esrft     ! snow energy of rafting ice 
     
    518509      INTEGER , DIMENSION(jpij) ::   itest_rdg, itest_rft   ! test for conservation 
    519510      !!------------------------------------------------------------------- 
    520  
     511      ! 
     512      zvti(1:npti) = SUM( v_i_2d(1:npti,:), dim=2 )   ! total ice volume 
     513      ! 
    521514      ! 1) Change in open water area due to closing and opening 
    522515      !-------------------------------------------------------- 
     
    535528            IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN   ! only if ice is ridging 
    536529 
    537                z1_ai(ji) = 1._wp / a_i_2d(ji,jl1) 
    538  
     530               IF( a_i_2d(ji,jl1) > epsi10 ) THEN   ;   z1_ai(ji) = 1._wp / a_i_2d(ji,jl1) 
     531               ELSE                                 ;   z1_ai(ji) = 0._wp 
     532               ENDIF 
     533                
    539534               ! area of ridging / rafting ice (airdg1) and of new ridge (airdg2) 
    540535               airdg1 = aridge(ji,jl1) * closing_gross(ji) * rdt_ice 
     
    549544 
    550545               ! volume and enthalpy (J/m2, >0) of seawater trapped into ridges 
    551                vsw = v_i_2d(ji,jl1) * afrdg * rn_porordg 
     546               IF    ( zvti(ji) <= 10. ) THEN ; vsw = v_i_2d(ji,jl1) * afrdg * rn_porordg                                           ! v <= 10m then porosity = rn_porordg 
     547               ELSEIF( zvti(ji) >= 20. ) THEN ; vsw = 0._wp                                                                         ! v >= 20m then porosity = 0 
     548               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 
     549               ENDIF 
    552550               ersw(ji) = -rhoi * vsw * rcp * sst_1d(ji)   ! clem: if sst>0, then ersw <0 (is that possible?) 
    553551 
    554552               ! volume etc of ridging / rafting ice and new ridges (vi, vs, sm, oi, es, ei) 
    555553               virdg1     = v_i_2d (ji,jl1)   * afrdg 
    556                virdg2(ji) = v_i_2d (ji,jl1)   * afrdg * ( 1. + rn_porordg ) 
     554               virdg2(ji) = v_i_2d (ji,jl1)   * afrdg + vsw 
    557555               vsrdg(ji)  = v_s_2d (ji,jl1)   * afrdg 
    558556               sirdg1     = sv_i_2d(ji,jl1)   * afrdg 
     
    588586               ! virtual salt flux to keep salinity constant 
    589587               IF( nn_icesal /= 2 )  THEN 
    590                   sirdg2(ji)     = sirdg2(ji)     - vsw * ( sss_1d(ji) - s_i_1d(ji) )        ! ridge salinity = s_i 
     588                  sirdg2(ji)     = sirdg2(ji)     - vsw * ( sss_1d(ji) - s_i_1d(ji) )       ! ridge salinity = s_i 
    591589                  sfx_bri_1d(ji) = sfx_bri_1d(ji) + sss_1d(ji) * vsw * rhoi * r1_rdtice  &  ! put back sss_m into the ocean 
    592590                     &                            - s_i_1d(ji) * vsw * rhoi * r1_rdtice     ! and get  s_i  from the ocean  
     
    726724      END DO ! jl1 
    727725      ! 
     726      ! roundoff errors 
     727      !---------------- 
    728728      ! 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 
     729      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 ) 
    731730      ! 
    732731   END SUBROUTINE rdgrft_shift 
     
    854853         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_snw_dyn_1d(1:npti), wfx_snw_dyn(:,:) ) 
    855854         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d    (1:npti), wfx_pnd    (:,:) ) 
    856  
     855         ! 
    857856         !                 !---------------------! 
    858857      CASE( 2 )            !==  from 1D to 2D  ==! 
     
    911910      REWIND( numnam_ice_ref )              ! Namelist namicetdme in reference namelist : Ice mechanical ice redistribution 
    912911      READ  ( numnam_ice_ref, namdyn_rdgrft, IOSTAT = ios, ERR = 901) 
    913 901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namdyn_rdgrft in reference namelist', lwp ) 
     912901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namdyn_rdgrft in reference namelist' ) 
    914913      REWIND( numnam_ice_cfg )              ! Namelist namdyn_rdgrft in configuration namelist : Ice mechanical ice redistribution 
    915914      READ  ( numnam_ice_cfg, namdyn_rdgrft, IOSTAT = ios, ERR = 902 ) 
    916 902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namdyn_rdgrft in configuration namelist', lwp ) 
     915902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namdyn_rdgrft in configuration namelist' ) 
    917916      IF(lwm) WRITE ( numoni, namdyn_rdgrft ) 
    918917      ! 
     
    945944         CALL ctl_stop( 'ice_dyn_rdgrft_init: choose one and only one participation function (ln_partf_lin or ln_partf_exp)' ) 
    946945      ENDIF 
    947       !                              ! allocate tke arrays 
     946      ! 
     947      IF( .NOT. ln_icethd ) THEN 
     948         rn_porordg = 0._wp 
     949         rn_fsnwrdg = 1._wp ; rn_fsnwrft = 1._wp 
     950         rn_fpndrdg = 1._wp ; rn_fpndrft = 1._wp 
     951         IF( lwp ) THEN 
     952            WRITE(numout,*) '      ==> only ice dynamics is activated, thus some parameters must be changed' 
     953            WRITE(numout,*) '            rn_porordg   = ', rn_porordg 
     954            WRITE(numout,*) '            rn_fsnwrdg   = ', rn_fsnwrdg  
     955            WRITE(numout,*) '            rn_fpndrdg   = ', rn_fpndrdg  
     956            WRITE(numout,*) '            rn_fsnwrft   = ', rn_fsnwrft  
     957            WRITE(numout,*) '            rn_fpndrft   = ', rn_fpndrft  
     958         ENDIF 
     959      ENDIF 
     960      !                              ! allocate arrays 
    948961      IF( ice_dyn_rdgrft_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'ice_dyn_rdgrft_init: unable to allocate arrays' ) 
    949962      ! 
Note: See TracChangeset for help on using the changeset viewer.