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 13605 for NEMO – NEMO

Changeset 13605 for NEMO


Ignore:
Timestamp:
2020-10-14T18:13:19+02:00 (3 years ago)
Author:
techene
Message:

#2385 add diags and diags at F-point

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/DIA/diawri.F90

    r13295 r13605  
    1919   !!            3.7  ! 2014-01  (G. Madec) remove eddy induced velocity from no-IOM output 
    2020   !!                 !                     change name of output variables in dia_wri_state 
     21   !!            4.0  ! 2020-10  (A. Nasser, S. Techene) add diagnostic for SWE 
    2122   !!---------------------------------------------------------------------- 
    2223 
     
    120121      REAL(wp)::   zztmp , zztmpx   ! local scalar 
    121122      REAL(wp)::   zztmp2, zztmpy   !   -      - 
     123      REAL(wp)::   ze3 
    122124      REAL(wp), DIMENSION(jpi,jpj)     ::   z2d   ! 2D workspace 
    123125      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   z3d   ! 3D workspace 
     
    136138      CALL iom_put("e3u_0", e3u_0(:,:,:) ) 
    137139      CALL iom_put("e3v_0", e3v_0(:,:,:) ) 
     140      CALL iom_put("e3f_0", e3f_0(:,:,:) ) 
    138141      ! 
    139142      IF ( iom_use("e3t") .OR. iom_use("e3tdef") ) THEN  ! time-varying e3t 
     
    162165         CALL iom_put( "e3w" , z3d(:,:,:) ) 
    163166      ENDIF 
     167      IF ( iom_use("e3f") ) THEN                         ! time-varying e3f caution here at Kaa 
     168          DO jk = 1, jpk 
     169            z3d(:,:,jk) =  e3f(:,:,jk) 
     170         END DO 
     171         CALL iom_put( "e3f" , z3d(:,:,:) ) 
     172      ENDIF 
    164173 
    165174      IF( ll_wd ) THEN                                   ! sea surface height (brought back to the reference used for wetting and drying) 
    166          CALL iom_put( "ssh" , (ssh(:,:,Kmm)+ssh_ref)*tmask(:,:,1) ) 
     175         CALL iom_put( "ssh" , (ssh(:,:,Kmm)+ssh_ref)*ssmask(:,:) ) 
    167176      ELSE 
    168177         CALL iom_put( "ssh" , ssh(:,:,Kmm) )              ! sea surface height 
    169178      ENDIF 
    170179 
    171       IF( iom_use("wetdep") )   &                  ! wet depth 
    172          CALL iom_put( "wetdep" , ht_0(:,:) + ssh(:,:,Kmm) ) 
     180      IF( iom_use("wetdep") )    CALL iom_put( "wetdep" , ht_0(:,:) + ssh(:,:,Kmm) )   ! wet depth 
     181          
     182#if defined key_qco 
     183      IF( iom_use("ht") )   CALL iom_put( "ht" , ht(:,:)     )   ! water column at t-point 
     184      IF( iom_use("hu") )   CALL iom_put( "hu" , hu(:,:,Kmm) )   ! water column at u-point 
     185      IF( iom_use("hv") )   CALL iom_put( "hv" , hv(:,:,Kmm) )   ! water column at v-point 
     186      IF( iom_use("hf") )   CALL iom_put( "hf" , hf_0(:,:)*( 1._wp + r3f(:,:) ) )   ! water column at f-point (caution here at Naa) 
     187#endif 
    173188       
    174189      CALL iom_put( "toce", ts(:,:,:,jp_tem,Kmm) )    ! 3D temperature 
     
    297312      ENDIF 
    298313      ! 
     314      IF ( iom_use("sKE") ) THEN                        ! surface kinetic energy at T point 
     315         z2d(:,:) = 0._wp 
     316         DO_2D( 0, 0, 0, 0 ) 
     317            z2d(ji,jj) = 0.25_wp * ( uu(ji  ,jj,1,Kmm) * uu(ji  ,jj,1,Kmm) * e1e2u(ji  ,jj) * e3u(ji  ,jj,1,Kmm)  & 
     318               &                   + uu(ji-1,jj,1,Kmm) * uu(ji-1,jj,1,Kmm) * e1e2u(ji-1,jj) * e3u(ji-1,jj,1,Kmm)  & 
     319               &                   + vv(ji,jj  ,1,Kmm) * vv(ji,jj  ,1,Kmm) * e1e2v(ji,jj  ) * e3v(ji,jj  ,1,Kmm)  &  
     320               &                   + vv(ji,jj-1,1,Kmm) * vv(ji,jj-1,1,Kmm) * e1e2v(ji,jj-1) * e3v(ji,jj-1,1,Kmm)  )  & 
     321               &                 * r1_e1e2t(ji,jj) / e3t(ji,jj,1,Kmm) * ssmask(ji,jj) 
     322         END_2D 
     323         CALL lbc_lnk( 'diawri', z2d, 'T', 1. ) 
     324         IF ( iom_use("sKE" ) )  CALL iom_put( "sKE" , z2d )    
     325      ENDIF 
     326      !     
     327      IF ( iom_use("sKEf") ) THEN                        ! surface kinetic energy at F point 
     328         z2d(:,:) = 0._wp                                ! CAUTION : only valid in SWE, not with bathymetry 
     329         DO_2D( 0, 0, 0, 0 ) 
     330            z2d(ji,jj) = 0.25_wp * ( uu(ji,jj  ,1,Kmm) * uu(ji,jj  ,1,Kmm) * e1e2u(ji,jj  ) * e3u(ji,jj  ,1,Kmm)  & 
     331               &                   + uu(ji,jj+1,1,Kmm) * uu(ji,jj+1,1,Kmm) * e1e2u(ji,jj+1) * e3u(ji,jj+1,1,Kmm)  & 
     332               &                   + vv(ji  ,jj,1,Kmm) * vv(ji,jj  ,1,Kmm) * e1e2v(ji  ,jj) * e3v(ji  ,jj,1,Kmm)  &  
     333               &                   + vv(ji+1,jj,1,Kmm) * vv(ji+1,jj,1,Kmm) * e1e2v(ji+1,jj) * e3v(ji+1,jj,1,Kmm)  )  & 
     334               &                 * r1_e1e2f(ji,jj) / e3f(ji,jj,1) * ssfmask(ji,jj) 
     335         END_2D 
     336         CALL lbc_lnk( 'diawri', z2d, 'F', 1. ) 
     337         CALL iom_put( "sKEf", z2d )                      
     338      ENDIF 
     339      ! 
    299340      CALL iom_put( "hdiv", hdiv )                  ! Horizontal divergence 
    300341      ! 
     
    376417       
    377418      IF (ln_dia25h)   CALL dia_25h( kt, Kmm )        ! 25h averaging 
     419       
     420      ! Output of vorticity terms 
     421      IF ( iom_use("relvor")    .OR. iom_use("plavor")    .OR.   & 
     422         & iom_use("relpotvor") .OR. iom_use("abspotvor") .OR.   & 
     423         & iom_use("Ens")                                        ) THEN 
     424         ! 
     425         z2d(:,:) = 0._wp  
     426         ze3 = 0._wp  
     427         DO_2D( 1, 0, 1, 0 ) 
     428            z2d(ji,jj) = (   e2v(ji+1,jj  ) * vv(ji+1,jj  ,1,Kmm) - e2v(ji,jj) * vv(ji,jj,1,Kmm)    & 
     429            &              - e1u(ji  ,jj+1) * uu(ji  ,jj+1,1,Kmm) + e1u(ji,jj) * uu(ji,jj,1,Kmm)  ) * r1_e1e2f(ji,jj) 
     430         END_2D 
     431         CALL lbc_lnk( 'diawri', z2d, 'F', 1. ) 
     432         CALL iom_put( "relvor", z2d )                  ! relative vorticity ( zeta )  
     433         ! 
     434         CALL iom_put( "plavor", ff_f )                 ! planetary vorticity ( f ) 
     435         ! 
     436         DO_2D( 1, 0, 1, 0 )   
     437            ze3 = (  e3t(ji,jj+1,1,Kmm) * e1e2t(ji,jj+1) + e3t(ji+1,jj+1,1,Kmm) * e1e2t(ji+1,jj+1)    & 
     438              &    + e3t(ji,jj  ,1,Kmm) * e1e2t(ji,jj  ) + e3t(ji+1,jj  ,1,Kmm) * e1e2t(ji+1,jj  )  ) * r1_e1e2f(ji,jj) 
     439            IF( ze3 /= 0._wp ) THEN   ;   ze3 = 4._wp / ze3 
     440            ELSE                      ;   ze3 = 0._wp 
     441            ENDIF 
     442            z2d(ji,jj) = ze3 * z2d(ji,jj)  
     443         END_2D 
     444         CALL lbc_lnk( 'diawri', z2d, 'F', 1. ) 
     445         CALL iom_put( "relpotvor", z2d )                  ! relative potential vorticity (zeta/h) 
     446         ! 
     447         DO_2D( 1, 0, 1, 0 ) 
     448            ze3 = (  e3t(ji,jj+1,1,Kmm) * e1e2t(ji,jj+1) + e3t(ji+1,jj+1,1,Kmm) * e1e2t(ji+1,jj+1)    & 
     449              &    + e3t(ji,jj  ,1,Kmm) * e1e2t(ji,jj  ) + e3t(ji+1,jj  ,1,Kmm) * e1e2t(ji+1,jj  )  ) * r1_e1e2f(ji,jj) 
     450            IF( ze3 /= 0._wp ) THEN   ;   ze3 = 4._wp / ze3 
     451            ELSE                      ;   ze3 = 0._wp 
     452            ENDIF 
     453            z2d(ji,jj) = ze3 * ff_f(ji,jj) + z2d(ji,jj)  
     454         END_2D 
     455         CALL lbc_lnk( 'diawri', z2d, 'F', 1. ) 
     456         CALL iom_put( "abspotvor", z2d )                  ! absolute potential vorticity ( q ) 
     457         ! 
     458         DO_2D( 1, 0, 1, 0 )   
     459            z2d(ji,jj) = 0.5_wp * z2d(ji,jj)  * z2d(ji,jj)  
     460         END_2D 
     461         CALL lbc_lnk( 'diawri', z2d, 'F', 1. ) 
     462         CALL iom_put( "Ens", z2d )                        ! potential enstrophy ( 1/2*q2 ) 
     463         ! 
     464      ENDIF 
    378465 
    379466      IF( ln_timing )   CALL timing_stop('dia_wri') 
Note: See TracChangeset for help on using the changeset viewer.