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

Changeset 15515


Ignore:
Timestamp:
2021-11-16T14:15:05+01:00 (6 months ago)
Author:
vancop
Message:

some bugfixes on diagnostics

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2020/SI3-03_VP_rheology/src/ICE/icedyn_rhg_vp.F90

    r15499 r15515  
    164164      ! 
    165165      REAL(wp), DIMENSION(jpi,jpj) ::   zdeltat, zdeltastar_t           ! Delta & Delta* at T-points 
    166       REAL(wp), DIMENSION(jpi,jpj) ::   zten_i                          ! Tension 
     166      REAL(wp), DIMENSION(jpi,jpj) ::   ztens, zshear                   ! Tension, shear 
    167167      REAL(wp), DIMENSION(jpi,jpj) ::   zp_deltastar_t                  ! P/delta* at T points 
    168168      REAL(wp), DIMENSION(jpi,jpj) ::   zzt, zet                        ! Viscosity pre-factors at T points 
     
    538538                
    539539               ! Temporary zef factor at F-point 
    540                zef(ji,jj)      = zp_deltastar_f * r1_e1e2f(ji,jj) * z1_ecc2 * zfmask(ji,jj) 
     540               zef(ji,jj)      = zp_deltastar_f * r1_e1e2f(ji,jj) * z1_ecc2 * zfmask(ji,jj) * 0.5_wp ! BUG missing un facteur 2 Nov 12 
    541541 
    542542            END DO 
     
    12191219            zdt2 = zdt * zdt 
    12201220             
    1221             zten_i(ji,jj) = zdt 
     1221            ztens(ji,jj)    = zdt 
    12221222             
    12231223            ! shear**2 at T points (doc eq. A16) 
     
    12261226               &   ) * 0.25_wp * r1_e1e2t(ji,jj) 
    12271227             
    1228             ! shear at T points 
    1229             pshear_i(ji,jj) = SQRT( zdt2 + zds2 ) 
     1228            ! maximum shear rate at T points (includees tension, output only) 
     1229            pshear_i(ji,jj) = SQRT( zdt2 + zds2 ) ! i think this is maximum shear rate and not actual shear --- i'm not totally sure here 
     1230 
     1231            ! shear at T-points 
     1232            zshear(ji,jj)   = SQRT( zds2 ) 
    12301233 
    12311234            ! divergence at T points 
     
    12331236               &             + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   & 
    12341237               &             ) * r1_e1e2t(ji,jj) 
     1238 
     1239            zdiv2          =  pdivu_i(ji,jj) *  pdivu_i(ji,jj) 
    12351240             
    12361241            ! delta at T points 
    1237             zdelta         = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 )   
     1242            zdelta         = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 ) 
    12381243            rswitch        = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zdelta ) ) ! 0 if delta=0 
    12391244             
     
    12611266      ! 
    12621267      ! ---- Sea ice stresses at T-points 
    1263       IF ( iom_use('normstr') .OR. iom_use('sheastr') .OR. iom_use('intstrx') .OR. iom_use('intstry') ) THEN 
     1268      IF ( iom_use('normstr')    .OR. iom_use('sheastr')    .OR. & 
     1269     &     iom_use('intstrx')    .OR. iom_use('intstry')    .OR. & 
     1270     &     iom_use('sig1_pnorm') .OR. iom_use('sig2_pnorm') ) THEN 
    12641271       
    12651272         DO jj = 2, jpj - 1 
     
    12681275               zfac                    =   zp_deltastar_t(ji,jj)  
    12691276               zs1(ji,jj)              =   zfac * ( pdivu_i(ji,jj) - pdelta_i(ji,jj) ) 
    1270                zs2(ji,jj)              =   zfac * z1_ecc2 * zten_i(ji,jj) 
    1271                zs12(ji,jj)             =   zfac * z1_ecc2 * pshear_i(ji,jj) 
     1277               zs2(ji,jj)              =   zfac * z1_ecc2 * ztens(ji,jj) 
     1278               zs12(ji,jj)             =   zfac * z1_ecc2 * zshear(ji,jj) * 0.5_wp ! Bug 12 nov 
    12721279            END DO 
    12731280         END DO 
     
    13471354      ! --- Divergence, shear and strength --- ! 
    13481355      IF( iom_use('icediv') )   CALL iom_put( 'icediv' , pdivu_i  * zmsk00 )   ! divergence 
    1349       IF( iom_use('iceshe') )   CALL iom_put( 'iceshe' , pshear_i * zmsk00 )   ! shear 
     1356      IF( iom_use('iceshe') )   CALL iom_put( 'iceshe' , pshear_i * zmsk00 )   ! maximum shear rate 
    13501357      IF( iom_use('icedlt') )   CALL iom_put( 'icedlt' , pdelta_i * zmsk00 )   ! delta 
    13511358      IF( iom_use('icestr') )   CALL iom_put( 'icestr' , strength * zmsk00 )   ! strength 
     
    13671374            DO ji = 2, jpi - 1 
    13681375               ! Stress invariants 
    1369                zsig_I(ji,jj)    =   zs1(ji,jj) * 0.5_wp                                                   ! 1st invariant, aka average normal stress aka negative pressure 
    1370                zsig_II(ji,jj)   =   SQRT ( zs2(ji,jj) * zs2(ji,jj) * 0.25_wp + zs12(ji,jj)*zs12(ji,jj) )  ! 2nd invariant, aka maximum shear stress                
     1376               zsig_I(ji,jj)    =   zs1(ji,jj) * 0.5_wp   ! 1st invariant, aka average normal stress aka negative pressure 
     1377               zsig_II(ji,jj)   =   0.5_wp * SQRT ( zs2(ji,jj) * zs2(ji,jj) + 4. * zs12(ji,jj) * zs12(ji,jj) )   ! 2nd invariant, aka maximum shear stress                
    13711378            END DO 
    13721379         END DO 
     
    14001407               !                        following Lemieux & Dupont (2020) 
    14011408               zfac             =   zp_deltastar_t(ji,jj) 
    1402                zsig1            =   zfac * ( pdivu_i(ji,jj) - zdeltastar_t(ji,jj) ) 
    1403                zsig2            =   zfac * z1_ecc2 * zten_i(ji,jj) 
    1404                zsig12           =   zfac * z1_ecc2 * pshear_i(ji,jj) 
     1409               zsig1            =   zfac * ( pdivu_i(ji,jj) - zdeltat(ji,jj) ) 
     1410               zsig2            =   zfac * z1_ecc2 * ztens(ji,jj) 
     1411               zsig12           =   zfac * z1_ecc2 * zshear(ji,jj) * 0.5_wp ! Bugfix 12 Nov 
    14051412                
    14061413               ! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008), T-point 
    1407                zsig_I(ji,jj)    =   zsig1 * 0.5_wp                              ! 1st invariant 
    1408                zsig_II(ji,jj)   =   SQRT ( zsig2 * zsig2 * 0.25_wp + zsig12 * zsig12 )   ! 2nd invariant 
     1414               zsig_I(ji,jj)    =   zsig1 * 0.5_wp                                       ! 1st invariant 
     1415               zsig_II(ji,jj)   =   0.5_wp * SQRT ( zsig2 * zsig2 + 4. *zsig12 * zsig12 )   ! 2nd invariant 
    14091416 
    14101417               ! Normalized  principal stresses (used to display the ellipse) 
     
    18481855      !                                                ! ==================== ! 
    18491856 
    1850       it_inn_file =  ( kt - 1 ) * kitinntotmax + kitinntot ! time step in the file 
    1851       it_out_file =  ( kt - 1 ) * kitoutmax  + kitout 
     1857      it_inn_file =  ( kt - nit000 ) * kitinntotmax + kitinntot ! time step in the file 
     1858      it_out_file =  ( kt - nit000 ) * kitoutmax    + kitout 
    18521859 
    18531860      ! write variables 
     
    18871894 
    18881895         ! close file 
    1889          IF( kt == nitend )   istatus = NF90_CLOSE( ncvgid ) 
     1896         ! IF( kt == nitend )   istatus = NF90_CLOSE( ncvgid ) ! BUG Ed blockley 
     1897         IF( kt == nitend - nn_fsbc + 1 .AND. kitinntot == kitinntotmax )    istatus = NF90_CLOSE( ncvgid ) 
    18901898      ENDIF 
    18911899       
Note: See TracChangeset for help on using the changeset viewer.