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

Changeset 13646


Ignore:
Timestamp:
2020-10-20T17:33:01+02:00 (4 years ago)
Author:
clem
Message:

4.0-HEAD: change diagnostics of the ellipse for ice rheology (already done in the trunk). And avoid annoying prints from bbl

Location:
NEMO/releases/r4.0/r4.0-HEAD/src
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • NEMO/releases/r4.0/r4.0-HEAD/src/ICE/icedyn_rhg_evp.F90

    r13589 r13646  
    137137      ! 
    138138      REAL(wp), DIMENSION(jpi,jpj) ::   zdelta, zp_delt                 ! delta and P/delta at T points 
     139      REAL(wp), DIMENSION(jpi,jpj) ::   zten_i                          ! tension 
    139140      REAL(wp), DIMENSION(jpi,jpj) ::   zbeta                           ! beta coef from Kimmritz 2017 
    140141      ! 
    … …  
    168169      REAL(wp), DIMENSION(jpi,jpj) ::   zu_ice, zv_ice 
    169170      !! --- diags 
    170       REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zsig1, zsig2, zsig3 
     171      REAL(wp) ::   zsig1, zsig2, zsig12, zfac, z1_strength 
     172      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zsig_I, zsig_II, zsig1_p, zsig2_p          
    171173      !! --- SIMIP diags 
    172174      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xmtrp_ice ! X-component of ice mass transport (kg/s) 
    … …  
    764766            zdt2 = zdt * zdt 
    765767             
     768            zten_i(ji,jj) = zdt 
     769 
    766770            ! shear**2 at T points (doc eq. A16) 
    767771            zds2 = ( zds(ji,jj  ) * zds(ji,jj  ) * e1e2f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e1e2f(ji-1,jj  )  & 
    … …  
    778782             
    779783            ! delta at T points 
    780             zdelta(ji,jj)   = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 )   
    781             rswitch         = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zdelta(ji,jj) ) ) ! 0 if delta=0 
    782             pdelta_i(ji,jj) = zdelta(ji,jj) + rn_creepl * rswitch 
     784            zfac            = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 ) ! delta   
     785            rswitch         = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zfac ) ) ! 0 if delta=0 
     786            pdelta_i(ji,jj) = zfac + rn_creepl * rswitch ! delta+creepl 
    783787 
    784788         END DO 
    785789      END DO 
    786       CALL lbc_lnk_multi( 'icedyn_rhg_evp', pshear_i, 'T', 1., pdivu_i, 'T', 1., pdelta_i, 'T', 1. ) 
     790      CALL lbc_lnk_multi( 'icedyn_rhg_evp', pshear_i, 'T', 1., pdivu_i, 'T', 1., pdelta_i, 'T', 1., zten_i, 'T', 1., & 
     791         &                                  zs1     , 'T', 1., zs2    , 'T', 1., zs12    , 'F', 1. ) 
    787792       
    788793      ! --- Store the stress tensor for the next time step --- ! 
    789       CALL lbc_lnk_multi( 'icedyn_rhg_evp', zs1, 'T', 1., zs2, 'T', 1., zs12, 'F', 1. ) 
    790794      pstress1_i (:,:) = zs1 (:,:) 
    791795      pstress2_i (:,:) = zs2 (:,:) 
    … …  
    816820      IF( iom_use('icestr') )   CALL iom_put( 'icestr' , strength * zmsk00 )   ! strength 
    817821 
    818       ! --- stress tensor --- ! 
    819       IF( iom_use('isig1') .OR. iom_use('isig2') .OR. iom_use('isig3') .OR. iom_use('normstr') .OR. iom_use('sheastr') ) THEN 
    820          ! 
    821          ALLOCATE( zsig1(jpi,jpj) , zsig2(jpi,jpj) , zsig3(jpi,jpj) ) 
     822      ! --- Stress tensor invariants (SIMIP diags) --- ! 
     823      IF( iom_use('normstr') .OR. iom_use('sheastr') ) THEN 
     824         ! 
     825         ALLOCATE( zsig_I(jpi,jpj) , zsig_II(jpi,jpj) ) 
    822826         !          
    823          DO jj = 2, jpjm1 
    824             DO ji = 2, jpim1 
    825                zdum1 = ( zmsk00(ji-1,jj) * pstress12_i(ji-1,jj) + zmsk00(ji  ,jj-1) * pstress12_i(ji  ,jj-1) +  &  ! stress12_i at T-point 
    826                   &      zmsk00(ji  ,jj) * pstress12_i(ji  ,jj) + zmsk00(ji-1,jj-1) * pstress12_i(ji-1,jj-1) )  & 
    827                   &    / MAX( 1._wp, zmsk00(ji-1,jj) + zmsk00(ji,jj-1) + zmsk00(ji,jj) + zmsk00(ji-1,jj-1) ) 
    828  
    829                zshear = SQRT( pstress2_i(ji,jj) * pstress2_i(ji,jj) + 4._wp * zdum1 * zdum1 ) ! shear stress   
    830  
    831                zdum2 = zmsk00(ji,jj) / MAX( 1._wp, strength(ji,jj) ) 
    832  
    833 !!               zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) + zshear ) ! principal stress (y-direction, see Hunke & Dukowicz 2002) 
    834 !!               zsig2(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) - zshear ) ! principal stress (x-direction, see Hunke & Dukowicz 2002) 
    835 !!               zsig3(ji,jj) = zdum2**2 * ( ( pstress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 ) ! quadratic relation linking compressive stress to shear stress 
    836 !!                                                                                                               ! (scheme converges if this value is ~1, see Bouillon et al 2009 (eq. 11)) 
    837                zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) )          ! compressive stress, see Bouillon et al. 2015 
    838                zsig2(ji,jj) = 0.5_wp * zdum2 * ( zshear )                     ! shear stress 
    839                zsig3(ji,jj) = zdum2**2 * ( ( pstress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 ) 
    840             END DO 
    841          END DO 
    842          CALL lbc_lnk_multi( 'icedyn_rhg_evp', zsig1, 'T', 1., zsig2, 'T', 1., zsig3, 'T', 1. ) 
    843          ! 
    844          CALL iom_put( 'isig1' , zsig1 ) 
    845          CALL iom_put( 'isig2' , zsig2 ) 
    846          CALL iom_put( 'isig3' , zsig3 ) 
    847          ! 
    848          ! Stress tensor invariants (normal and shear stress N/m) 
    849          IF( iom_use('normstr') )   CALL iom_put( 'normstr' ,       ( zs1(:,:) + zs2(:,:) )                       * zmsk00(:,:) ) ! Normal stress 
    850          IF( iom_use('sheastr') )   CALL iom_put( 'sheastr' , SQRT( ( zs1(:,:) - zs2(:,:) )**2 + 4*zs12(:,:)**2 ) * zmsk00(:,:) ) ! Shear stress 
    851  
    852          DEALLOCATE( zsig1 , zsig2 , zsig3 ) 
     827         DO jj = 1, jpj 
     828            DO ji = 1, jpi 
     829             
     830               ! Ice stresses 
     831               ! sigma1, sigma2, sigma12 are some useful recombination of the stresses (Hunke and Dukowicz MWR 2002, Bouillon et al., OM2013) 
     832               ! These are NOT stress tensor components, neither stress invariants, neither stress principal components 
     833               ! I know, this can be confusing... 
     834               zfac             =   strength(ji,jj) / ( pdelta_i(ji,jj) + rn_creepl )  
     835               zsig1            =   zfac * ( pdivu_i(ji,jj) - pdelta_i(ji,jj) ) 
     836               zsig2            =   zfac * z1_ecc2 * zten_i(ji,jj) 
     837               zsig12           =   zfac * z1_ecc2 * pshear_i(ji,jj) 
     838                
     839               ! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008) 
     840               zsig_I (ji,jj)   =   zsig1 * 0.5_wp                                           ! 1st stress invariant, aka average normal stress, aka negative pressure 
     841               zsig_II(ji,jj)   =   SQRT ( MAX( 0._wp, zsig2 * zsig2 * 0.25_wp + zsig12 ) )  ! 2nd  ''       '', aka maximum shear stress 
     842                
     843            END DO 
     844         END DO          
     845         ! 
     846         ! Stress tensor invariants (normal and shear stress N/m) - SIMIP diags - definitions following Coon (1974) and Feltham (2008) 
     847         IF( iom_use('normstr') )   CALL iom_put( 'normstr', zsig_I (:,:) * zmsk00(:,:) ) ! Normal stress 
     848         IF( iom_use('sheastr') )   CALL iom_put( 'sheastr', zsig_II(:,:) * zmsk00(:,:) ) ! Maximum shear stress 
     849          
     850         DEALLOCATE ( zsig_I, zsig_II ) 
     851          
    853852      ENDIF 
     853 
     854      ! --- Normalized stress tensor principal components --- ! 
     855      ! This are used to plot the normalized yield curve, see Lemieux & Dupont, 2020 
     856      ! Recommendation 1 : we use ice strength, not replacement pressure 
     857      ! Recommendation 2 : need to use deformations at PREVIOUS iterate for viscosities 
     858!!$      IF( iom_use('sig1_pnorm') .OR. iom_use('sig2_pnorm') ) THEN 
     859!!$         ! 
     860!!$         ALLOCATE( zsig1_p(jpi,jpj) , zsig2_p(jpi,jpj) , zsig_I(jpi,jpj) , zsig_II(jpi,jpj) )          
     861!!$         !          
     862!!$         DO jj = 1, jpj 
     863!!$            DO ji = 1, jpi 
     864!!$             
     865!!$               ! Ice stresses computed with **viscosities** (delta, p/delta) at **previous** iterates  
     866!!$               !                        and **deformations** at current iterates 
     867!!$               !                        following Lemieux & Dupont (2020) 
     868!!$               zfac             =   zp_delt(ji,jj) 
     869!!$               zsig1            =   zfac * ( pdivu_i(ji,jj) - ( zdelta(ji,jj) + rn_creepl ) ) 
     870!!$               zsig2            =   zfac * z1_ecc2 * zten_i(ji,jj) 
     871!!$               zsig12           =   zfac * z1_ecc2 * pshear_i(ji,jj) 
     872!!$                
     873!!$               ! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008), T-point 
     874!!$               zsig_I(ji,jj)    =   zsig1 * 0.5_wp                                           ! 1st stress invariant, aka average normal stress, aka negative pressure 
     875!!$               zsig_II(ji,jj)   =   SQRT ( MAX( 0._wp, zsig2 * zsig2 * 0.25_wp + zsig12 ) )  ! 2nd  ''       '', aka maximum shear stress 
     876!!$       
     877!!$               ! Normalized  principal stresses (used to display the ellipse) 
     878!!$               z1_strength      =   1._wp / MAX( 1._wp, strength(ji,jj) ) 
     879!!$               zsig1_p(ji,jj)   =   ( zsig_I(ji,jj) + zsig_II(ji,jj) ) * z1_strength 
     880!!$               zsig2_p(ji,jj)   =   ( zsig_I(ji,jj) - zsig_II(ji,jj) ) * z1_strength 
     881!!$            END DO 
     882!!$         END DO                
     883!!$         ! 
     884!!$         CALL iom_put( 'sig1_pnorm' , zsig1_p )  
     885!!$         CALL iom_put( 'sig2_pnorm' , zsig2_p )  
     886!!$       
     887!!$         DEALLOCATE( zsig1_p , zsig2_p , zsig_I, zsig_II ) 
     888!!$          
     889!!$      ENDIF 
    854890 
    855891      ! --- SIMIP --- ! 

  • NEMO/releases/r4.0/r4.0-HEAD/src/OCE/TRA/trabbl.F90

    r11536 r13646  
    513513      IF( tra_bbl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'tra_bbl_init : unable to allocate arrays' ) 
    514514      ! 
    515       IF( nn_bbl_adv == 1 )    WRITE(numout,*) '       * Advective BBL using upper velocity' 
    516       IF( nn_bbl_adv == 2 )    WRITE(numout,*) '       * Advective BBL using velocity = F( delta rho)' 
     515      IF(lwp) THEN 
     516         IF( nn_bbl_adv == 1 )    WRITE(numout,*) '       * Advective BBL using upper velocity' 
     517         IF( nn_bbl_adv == 2 )    WRITE(numout,*) '       * Advective BBL using velocity = F( delta rho)' 
     518      ENDIF 
    517519      ! 
    518520      !                             !* vertical index of  "deep" bottom u- and v-points 
Note: See TracChangeset for help on using the changeset viewer.