Ignore:
Timestamp:
2017-09-05T19:53:41+02:00 (3 years ago)
Author:
clem
Message:

changes in style - part2 -

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerhg_evp.F90

    r8486 r8498  
    3131   USE prtctl         ! Print control 
    3232   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
     33   USE iom            ! 
    3334#if defined key_agrif 
    3435   USE agrif_lim3_interp 
     
    111112      INTEGER ::   ji, jj       ! dummy loop indices 
    112113      INTEGER ::   jter         ! local integers 
    113       CHARACTER (len=50) ::   charout 
    114114 
    115115      REAL(wp) ::   zrhoco                                                   ! rau0 * rn_cio 
     
    121121      REAL(wp) ::   zTauO, zTauB, zTauE, zvel                                ! temporary scalars 
    122122 
    123       REAL(wp) ::   zsig1, zsig2                                             ! internal ice stress 
    124123      REAL(wp) ::   zresm                                                    ! Maximal error on ice velocity 
    125124      REAL(wp) ::   zintb, zintn                                             ! dummy argument 
    126125      REAL(wp) ::   zfac_x, zfac_y 
     126      REAL(wp) ::   zshear, zdum1, zdum2 
    127127       
    128128      REAL(wp), DIMENSION(jpi,jpj) ::   z1_e1t0, z1_e2t0                ! scale factors 
     
    152152      REAL(wp), PARAMETER          ::   zepsi  = 1.0e-20_wp             ! tolerance parameter 
    153153      REAL(wp), PARAMETER          ::   zmmin  = 1._wp                  ! ice mass (kg/m2) below which ice velocity equals ocean velocity 
     154      !! --- diags 
     155      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zswi, zsig1, zsig2, zsig3 
     156      !! --- SIMIP diags 
     157      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_sig1      ! Average normal stress in sea ice    
     158      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_sig2      ! Maximum shear stress in sea ice 
     159      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_dssh_dx   ! X-direction sea-surface tilt term (N/m2) 
     160      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_dssh_dy   ! X-direction sea-surface tilt term (N/m2) 
     161      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_corstrx   ! X-direction coriolis stress (N/m2) 
     162      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_corstry   ! Y-direction coriolis stress (N/m2) 
     163      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_intstrx   ! X-direction internal stress (N/m2) 
     164      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_intstry   ! Y-direction internal stress (N/m2) 
     165      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_utau_oi   ! X-direction ocean-ice stress 
     166      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_vtau_oi   ! Y-direction ocean-ice stress   
     167      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xmtrp_ice ! X-component of ice mass transport (kg/s) 
     168      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_ymtrp_ice ! Y-component of ice mass transport (kg/s) 
     169      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xmtrp_snw ! X-component of snow mass transport (kg/s) 
     170      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_ymtrp_snw ! Y-component of snow mass transport (kg/s) 
     171      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xatrp     ! X-component of area transport (m2/s) 
     172      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_yatrp     ! Y-component of area transport (m2/s)       
    154173      !!------------------------------------------------------------------- 
    155174 
     
    671690 
    672691      !------------------------------------------------------------------------------! 
    673       ! 5) SIMIP diagnostics 
    674       !------------------------------------------------------------------------------! 
    675  
    676 !!gm  encapsulate with a flag (iom_use of the variable or better a flag defined one for all testing if one of the 
    677 !!    diag is output.  NB the diag_...  are should only be ALLOCATED if the flag is true ! 
    678  
    679       DO jj = 2, jpjm1 
    680          DO ji = 2, jpim1 
    681              rswitch  = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice 
    682  
    683              ! Stress tensor invariants (normal and shear stress N/m) 
    684              diag_sig1(ji,jj) = ( zs1(ji,jj) + zs2(ji,jj) ) * rswitch                                 ! normal stress 
    685              diag_sig2(ji,jj) = SQRT( ( zs1(ji,jj) - zs2(ji,jj) )**2 + 4*zs12(ji,jj)**2 ) * rswitch   ! shear stress 
    686  
    687              ! Stress terms of the momentum equation (N/m2) 
    688              diag_dssh_dx(ji,jj) = zspgU(ji,jj) * rswitch    ! sea surface slope stress term 
    689              diag_dssh_dy(ji,jj) = zspgV(ji,jj) * rswitch 
    690  
    691              diag_corstrx(ji,jj) = zCorx(ji,jj) * rswitch    ! Coriolis stress term 
    692              diag_corstry(ji,jj) = zCory(ji,jj) * rswitch 
    693  
    694              diag_intstrx(ji,jj) = zfU(ji,jj)   * rswitch    ! internal stress term 
    695              diag_intstry(ji,jj) = zfV(ji,jj)   * rswitch 
    696             
    697              diag_utau_oi(ji,jj) = ztaux_oi(ji,jj) * rswitch  ! oceanic stress 
    698              diag_vtau_oi(ji,jj) = ztauy_oi(ji,jj) * rswitch 
    699  
    700              ! 2D ice mass, snow mass, area transport arrays (X, Y) 
    701              zfac_x = 0.5 * u_ice(ji,jj) * e2u(ji,jj) * rswitch 
    702              zfac_y = 0.5 * v_ice(ji,jj) * e1v(ji,jj) * rswitch 
    703  
    704              diag_xmtrp_ice(ji,jj) = rhoic * zfac_x * ( vt_i(ji+1,jj) + vt_i(ji,jj) ) ! ice mass transport, X-component 
    705              diag_ymtrp_ice(ji,jj) = rhoic * zfac_y * ( vt_i(ji,jj+1) + vt_i(ji,jj) ) !        ''           Y-   '' 
    706  
    707              diag_xmtrp_snw(ji,jj) = rhosn * zfac_x * ( vt_s(ji+1,jj) + vt_s(ji,jj) ) ! snow mass transport, X-component 
    708              diag_ymtrp_snw(ji,jj) = rhosn * zfac_y * ( vt_s(ji,jj+1) + vt_s(ji,jj) ) !          ''          Y-   '' 
    709  
    710              diag_xatrp(ji,jj)     = zfac_x * ( at_i(ji+1,jj) + at_i(ji,jj) )         ! area transport,      X-component 
    711              diag_yatrp(ji,jj)     = zfac_y * ( at_i(ji,jj+1) + at_i(ji,jj) )         !        ''            Y-   '' 
    712  
    713          END DO 
    714       END DO 
    715  
    716       CALL lbc_lnk_multi(   diag_sig1   , 'T',  1., diag_sig2   , 'T',  1.,   & 
    717          &                  diag_dssh_dx, 'U', -1., diag_dssh_dy, 'V', -1.,   & 
    718          &                  diag_corstrx, 'U', -1., diag_corstry, 'V', -1.,   &  
    719          &                  diag_intstrx, 'U', -1., diag_intstry, 'V', -1.    ) 
    720  
    721       CALL lbc_lnk_multi(   diag_utau_oi, 'U', -1., diag_vtau_oi, 'V', -1.    ) 
    722  
    723       CALL lbc_lnk_multi(   diag_xmtrp_ice, 'U', -1., diag_xmtrp_snw, 'U', -1.,   & 
    724          &                  diag_xatrp    , 'U', -1., diag_ymtrp_ice, 'V', -1.,   & 
    725          &                  diag_ymtrp_snw, 'V', -1., diag_yatrp    , 'V', -1.    ) 
    726  
    727       ! 
    728       !------------------------------------------------------------------------------! 
    729       ! 6) Control prints of residual and charge ellipse 
    730       !------------------------------------------------------------------------------! 
    731       ! 
    732       ! print the residual for convergence 
    733       IF(ln_ctl) THEN 
    734          WRITE(charout,FMT="('ice_rhg_evp  : res =',D23.16, ' iter =',I4)") zresm, jter 
    735          CALL prt_ctl_info(charout) 
    736          CALL prt_ctl(tab2d_1=u_ice, clinfo1=' ice_rhg_evp  : u_ice :', tab2d_2=v_ice, clinfo2=' v_ice :') 
     692      ! 5) diagnostics 
     693      !------------------------------------------------------------------------------! 
     694      ! --- ellipse --- ! 
     695      IF( iom_use('isig1') .OR. iom_use('isig2') .OR. iom_use('isig3') ) THEN 
    737696         ! 
    738          ! print charge ellipse 
    739          ! This can be desactivated once the user is sure that the stress state 
    740       ! lie on the charge ellipse. See Bouillon et al. (2008) for more details 
    741          IF( MOD(kt_ice+nn_fsbc-1,nwrite) == 0 ) THEN 
    742             WRITE(charout,FMT="('ice_rhg_evp  :', I4, I6, I1, I1, A10)") 1000, kt_ice, 0, 0, ' ch. ell. ' 
    743             CALL prt_ctl_info(charout) 
    744             DO jj = 2, jpjm1 
    745                DO ji = 2, jpim1 
    746                   IF (strength(ji,jj) > 1.0) THEN 
    747                      zsig1 = ( zs1(ji,jj) + SQRT(zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 ) ) / ( 2*strength(ji,jj) )  
    748                      zsig2 = ( zs1(ji,jj) - SQRT(zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 ) ) / ( 2*strength(ji,jj) ) 
    749                      WRITE(charout,FMT="('ice_rhg_evp  :', I4, I4, D23.16, D23.16, D23.16, D23.16, A10)") 
    750                      CALL prt_ctl_info(charout) 
    751                   ENDIF 
    752                END DO 
    753             END DO 
    754             WRITE(charout,FMT="('ice_rhg_evp  :', I4, I6, I1, I1, A10)") 2000, kt_ice, 0, 0, ' ch. ell. ' 
    755             CALL prt_ctl_info(charout) 
    756          ENDIF 
     697         ALLOCATE( zswi(jpi,jpj) , zsig1(jpi,jpj) , zsig2(jpi,jpj) , zsig3(jpi,jpj) ) 
     698         ! 
     699         DO jj = 1, jpj 
     700            DO ji = 1, jpi 
     701               zswi(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice 
     702            END DO 
     703         END DO 
     704          
     705         DO jj = 2, jpjm1 
     706            DO ji = 2, jpim1 
     707               zdum1 = ( zswi(ji-1,jj) * stress12_i(ji-1,jj) + zswi(ji  ,jj-1) * stress12_i(ji  ,jj-1) +  &  ! stress12_i at T-point 
     708                  &     zswi(ji  ,jj) * stress12_i(ji  ,jj) + zswi(ji-1,jj-1) * stress12_i(ji-1,jj-1) )  & 
     709                  &   / MAX( 1._wp, zswi(ji-1,jj) + zswi(ji,jj-1) + zswi(ji,jj) + zswi(ji-1,jj-1) ) 
     710 
     711               zshear = SQRT( stress2_i(ji,jj) * stress2_i(ji,jj) + 4._wp * zdum1 * zdum1 ) ! shear stress   
     712 
     713               zdum2 = zswi(ji,jj) / MAX( 1._wp, strength(ji,jj) ) 
     714 
     715!!               zsig1(ji,jj) = 0.5_wp * zdum2 * ( stress1_i(ji,jj) + zshear ) ! principal stress (y-direction, see Hunke & Dukowicz 2002) 
     716!!               zsig2(ji,jj) = 0.5_wp * zdum2 * ( stress1_i(ji,jj) - zshear ) ! principal stress (x-direction, see Hunke & Dukowicz 2002) 
     717!!               zsig3(ji,jj) = zdum2**2 * ( ( stress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 ) ! quadratic relation linking compressive stress to shear stress 
     718!!                                                                                                               ! (scheme converges if this value is ~1, see Bouillon et al 2009 (eq. 11)) 
     719               zsig1(ji,jj) = 0.5_wp * zdum2 * ( stress1_i(ji,jj) )          ! compressive stress, see Bouillon et al. 2015 
     720               zsig2(ji,jj) = 0.5_wp * zdum2 * ( zshear )                    ! shear stress 
     721               zsig3(ji,jj) = zdum2**2 * ( ( stress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 ) 
     722            END DO 
     723         END DO 
     724         CALL lbc_lnk_multi( zsig1, 'T', 1., zsig2, 'T', 1., zsig3, 'T', 1. ) 
     725         ! 
     726         IF( iom_use('isig1') )   CALL iom_put( "isig1" , zsig1 ) 
     727         IF( iom_use('isig2') )   CALL iom_put( "isig2" , zsig2 ) 
     728         IF( iom_use('isig3') )   CALL iom_put( "isig3" , zsig3 ) 
     729         ! 
     730         DEALLOCATE( zswi , zsig1 , zsig2 , zsig3 ) 
     731      ENDIF 
     732       
     733      ! --- SIMIP --- ! 
     734      IF ( iom_use( 'normstr'  ) .OR. iom_use( 'sheastr'  ) .OR. iom_use( 'dssh_dx'  ) .OR. iom_use( 'dssh_dy'  ) .OR. & 
     735         & iom_use( 'corstrx'  ) .OR. iom_use( 'corstry'  ) .OR. iom_use( 'intstrx'  ) .OR. iom_use( 'intstry'  ) .OR. & 
     736         & iom_use( 'utau_oi'  ) .OR. iom_use( 'vtau_oi'  ) .OR. iom_use( 'xmtrpice' ) .OR. iom_use( 'ymtrpice' ) .OR. & 
     737         & iom_use( 'xmtrpsnw' ) .OR. iom_use( 'ymtrpsnw' ) .OR. iom_use( 'xatrp'    ) .OR. iom_use( 'yatrp'    ) ) THEN 
     738 
     739         ALLOCATE( zdiag_sig1     (jpi,jpj) , zdiag_sig2     (jpi,jpj) , zdiag_dssh_dx  (jpi,jpj) , zdiag_dssh_dy  (jpi,jpj) ,  & 
     740            &      zdiag_corstrx  (jpi,jpj) , zdiag_corstry  (jpi,jpj) , zdiag_intstrx  (jpi,jpj) , zdiag_intstry  (jpi,jpj) ,  & 
     741            &      zdiag_utau_oi  (jpi,jpj) , zdiag_vtau_oi  (jpi,jpj) , zdiag_xmtrp_ice(jpi,jpj) , zdiag_ymtrp_ice(jpi,jpj) ,  & 
     742            &      zdiag_xmtrp_snw(jpi,jpj) , zdiag_ymtrp_snw(jpi,jpj) , zdiag_xatrp    (jpi,jpj) , zdiag_yatrp    (jpi,jpj) ) 
     743          
     744         DO jj = 2, jpjm1 
     745            DO ji = 2, jpim1 
     746               rswitch  = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice 
     747                
     748               ! Stress tensor invariants (normal and shear stress N/m) 
     749               zdiag_sig1(ji,jj) = ( zs1(ji,jj) + zs2(ji,jj) ) * rswitch                                 ! normal stress 
     750               zdiag_sig2(ji,jj) = SQRT( ( zs1(ji,jj) - zs2(ji,jj) )**2 + 4*zs12(ji,jj)**2 ) * rswitch   ! shear stress 
     751                
     752               ! Stress terms of the momentum equation (N/m2) 
     753               zdiag_dssh_dx(ji,jj) = zspgU(ji,jj) * rswitch    ! sea surface slope stress term 
     754               zdiag_dssh_dy(ji,jj) = zspgV(ji,jj) * rswitch 
     755                
     756               zdiag_corstrx(ji,jj) = zCorx(ji,jj) * rswitch    ! Coriolis stress term 
     757               zdiag_corstry(ji,jj) = zCory(ji,jj) * rswitch 
     758                
     759               zdiag_intstrx(ji,jj) = zfU(ji,jj)   * rswitch    ! internal stress term 
     760               zdiag_intstry(ji,jj) = zfV(ji,jj)   * rswitch 
     761                
     762               zdiag_utau_oi(ji,jj) = ztaux_oi(ji,jj) * rswitch  ! oceanic stress 
     763               zdiag_vtau_oi(ji,jj) = ztauy_oi(ji,jj) * rswitch 
     764                
     765               ! 2D ice mass, snow mass, area transport arrays (X, Y) 
     766               zfac_x = 0.5 * u_ice(ji,jj) * e2u(ji,jj) * rswitch 
     767               zfac_y = 0.5 * v_ice(ji,jj) * e1v(ji,jj) * rswitch 
     768                
     769               zdiag_xmtrp_ice(ji,jj) = rhoic * zfac_x * ( vt_i(ji+1,jj) + vt_i(ji,jj) ) ! ice mass transport, X-component 
     770               zdiag_ymtrp_ice(ji,jj) = rhoic * zfac_y * ( vt_i(ji,jj+1) + vt_i(ji,jj) ) !        ''           Y-   '' 
     771                
     772               zdiag_xmtrp_snw(ji,jj) = rhosn * zfac_x * ( vt_s(ji+1,jj) + vt_s(ji,jj) ) ! snow mass transport, X-component 
     773               zdiag_ymtrp_snw(ji,jj) = rhosn * zfac_y * ( vt_s(ji,jj+1) + vt_s(ji,jj) ) !          ''          Y-   '' 
     774                
     775               zdiag_xatrp(ji,jj)     = zfac_x * ( at_i(ji+1,jj) + at_i(ji,jj) )         ! area transport,      X-component 
     776               zdiag_yatrp(ji,jj)     = zfac_y * ( at_i(ji,jj+1) + at_i(ji,jj) )         !        ''            Y-   '' 
     777                
     778            END DO 
     779         END DO 
     780          
     781         CALL lbc_lnk_multi(   zdiag_sig1   , 'T',  1., zdiag_sig2   , 'T',  1.,   & 
     782            &                  zdiag_dssh_dx, 'U', -1., zdiag_dssh_dy, 'V', -1.,   & 
     783            &                  zdiag_corstrx, 'U', -1., zdiag_corstry, 'V', -1.,   &  
     784            &                  zdiag_intstrx, 'U', -1., zdiag_intstry, 'V', -1.    ) 
     785          
     786         CALL lbc_lnk_multi(   zdiag_utau_oi, 'U', -1., zdiag_vtau_oi, 'V', -1.    ) 
     787          
     788         CALL lbc_lnk_multi(   zdiag_xmtrp_ice, 'U', -1., zdiag_xmtrp_snw, 'U', -1.,   & 
     789            &                  zdiag_xatrp    , 'U', -1., zdiag_ymtrp_ice, 'V', -1.,   & 
     790            &                  zdiag_ymtrp_snw, 'V', -1., zdiag_yatrp    , 'V', -1.    ) 
     791          
     792         IF ( iom_use( 'normstr'  ) )   CALL iom_put( 'normstr'  ,  zdiag_sig1(:,:)      )   ! Normal stress 
     793         IF ( iom_use( 'sheastr'  ) )   CALL iom_put( 'sheastr'  ,  zdiag_sig2(:,:)      )   ! Shear stress 
     794         IF ( iom_use( 'dssh_dx'  ) )   CALL iom_put( 'dssh_dx'  ,  zdiag_dssh_dx(:,:)   )   ! Sea-surface tilt term in force balance (x) 
     795         IF ( iom_use( 'dssh_dy'  ) )   CALL iom_put( 'dssh_dy'  ,  zdiag_dssh_dy(:,:)   )   ! Sea-surface tilt term in force balance (y) 
     796         IF ( iom_use( 'corstrx'  ) )   CALL iom_put( 'corstrx'  ,  zdiag_corstrx(:,:)   )   ! Coriolis force term in force balance (x) 
     797         IF ( iom_use( 'corstry'  ) )   CALL iom_put( 'corstry'  ,  zdiag_corstry(:,:)   )   ! Coriolis force term in force balance (y) 
     798         IF ( iom_use( 'intstrx'  ) )   CALL iom_put( 'intstrx'  ,  zdiag_intstrx(:,:)   )   ! Internal force term in force balance (x) 
     799         IF ( iom_use( 'intstry'  ) )   CALL iom_put( 'intstry'  ,  zdiag_intstry(:,:)   )   ! Internal force term in force balance (y) 
     800         IF ( iom_use( 'utau_oi'  ) )   CALL iom_put( 'utau_oi'  ,  zdiag_utau_oi(:,:)   )   ! Ocean stress term in force balance (x) 
     801         IF ( iom_use( 'vtau_oi'  ) )   CALL iom_put( 'vtau_oi'  ,  zdiag_vtau_oi(:,:)   )   ! Ocean stress term in force balance (y) 
     802         IF ( iom_use( 'xmtrpice' ) )   CALL iom_put( 'xmtrpice' ,  zdiag_xmtrp_ice(:,:) )   ! X-component of sea-ice mass transport (kg/s) 
     803         IF ( iom_use( 'ymtrpice' ) )   CALL iom_put( 'ymtrpice' ,  zdiag_ymtrp_ice(:,:) )   ! Y-component of sea-ice mass transport  
     804         IF ( iom_use( 'xmtrpsnw' ) )   CALL iom_put( 'xmtrpsnw' ,  zdiag_xmtrp_snw(:,:) )   ! X-component of snow mass transport (kg/s) 
     805         IF ( iom_use( 'ymtrpsnw' ) )   CALL iom_put( 'ymtrpsnw' ,  zdiag_ymtrp_snw(:,:) )   ! Y-component of snow mass transport 
     806         IF ( iom_use( 'xatrp'    ) )   CALL iom_put( 'xatrp'    ,  zdiag_xatrp(:,:)     )   ! X-component of ice area transport 
     807         IF ( iom_use( 'yatrp'    ) )   CALL iom_put( 'yatrp'    ,  zdiag_yatrp(:,:)     )   ! Y-component of ice area transport 
     808 
     809         DEALLOCATE( zdiag_sig1      , zdiag_sig2      , zdiag_dssh_dx   , zdiag_dssh_dy   ,  & 
     810            &        zdiag_corstrx   , zdiag_corstry   , zdiag_intstrx   , zdiag_intstry   ,  & 
     811            &        zdiag_utau_oi   , zdiag_vtau_oi   , zdiag_xmtrp_ice , zdiag_ymtrp_ice ,  & 
     812            &        zdiag_xmtrp_snw , zdiag_ymtrp_snw , zdiag_xatrp     , zdiag_yatrp     ) 
     813 
    757814      ENDIF 
    758815      ! 
Note: See TracChangeset for help on using the changeset viewer.