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

Changeset 3325


Ignore:
Timestamp:
2012-03-12T15:44:43+01:00 (12 years ago)
Author:
gm
Message:

Ediag branche: #927 add Kinetic Energy trend diagnostics (trdken.F90)

Location:
branches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/NEMO/OPA_SRC/TRD
Files:
5 edited
1 copied

Legend:

Unmodified
Added
Removed

  • branches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/NEMO/OPA_SRC/TRD/trddyn.F90

    r3318 r3325  
    55   !!===================================================================== 
    66   !! History :  3.5  !  2012-02  (G. Madec) creation from trdmod: split DYN and TRA trends 
    7    !!                           and manage  3D trends output for U, V, and KE 
     7   !!                                        and manage  3D trends output for U, V, and KE 
    88   !!---------------------------------------------------------------------- 
    99 
    … …  
    2121   USE sbc_oce        ! surface boundary condition: ocean 
    2222   USE phycst         ! physical constants 
    23    USE trdvor         ! ocean vorticity trends  
    24    USE trdglo         ! trends:global domain averaged 
    25    USE trdmld         ! ocean active mixed layer tracers trends  
     23   USE trdken         ! trends: Kinetic ENergy  
     24   USE trdglo         ! trends: global domain averaged 
     25   USE trdvor         ! trends: vertical averaged vorticity  
     26   USE trdmld         ! trends: mixed layer averaged  
    2627   USE in_out_manager ! I/O manager 
    2728   USE iom            ! I/O manager library 
    … …  
    3132   IMPLICIT NONE 
    3233   PRIVATE 
    33  
    34    REAL(wp) ::   r2dt    ! time-step, = 2 rdttra except at nit000 (=rdttra) if neuler=0 
    3534 
    3635   PUBLIC trd_dyn        ! called by all dynXX modules 
    … …  
    5756      INTEGER                   , INTENT(in   ) ::   ktrd           ! trend index 
    5857      INTEGER                   , INTENT(in   ) ::   kt             ! time step 
    59       !! 
    60       INTEGER ::   ji, jj   ! dummy loop indices 
    61       REAL(wp), POINTER, DIMENSION(:,:) ::   ztswu, ztswv    ! 2D workspace  
    6258      !!---------------------------------------------------------------------- 
     59      ! 
     60      putrd(:,:,:) = putrd(:,:,:) * umask(:,:,:)                       ! mask the trends 
     61      pvtrd(:,:,:) = pvtrd(:,:,:) * vmask(:,:,:) 
     62      ! 
    6363 
    64       CALL wrk_alloc( jpi, jpj, ztswu, ztswv ) 
    65  
    66       IF( neuler == 0 .AND. kt == nit000    ) THEN   ;   r2dt =      rdt      ! = rdtra (restart with Euler time stepping) 
    67       ELSEIF(               kt <= nit000 + 1) THEN   ;   r2dt = 2. * rdt      ! = 2 rdttra (leapfrog) 
    68       ENDIF 
     64!!gm NB : here a lbc_lnk should probably be added 
    6965 
    7066      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    … …  
    8177      !  Kinetic Energy trends 
    8278      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    83 !!gm      IF( ln_KE_trd  )   CALL trd_KE( putrd, pvtrd, ktrd, kt ) 
     79      IF( ln_KE_trd  )   CALL trd_ken( putrd, pvtrd, ktrd, kt ) 
    8480 
    8581      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    … …  
    9288      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    9389!!gm      IF( ln_dyn_mld )   CALL trd_mld_dyn    
    94       ! 
    95       CALL wrk_dealloc( jpi, jpj, ztswu, ztswv ) 
    9690      ! 
    9791   END SUBROUTINE trd_dyn 
    … …  
    113107      REAL(wp), POINTER, DIMENSION(:,:,:) ::   z3dx, z3dy                 ! 3D workspace  
    114108      !!---------------------------------------------------------------------- 
    115       ! 
    116       putrd(:,:,:) = putrd(:,:,:) * umask(:,:,:)                       ! mask the trends 
    117       pvtrd(:,:,:) = pvtrd(:,:,:) * vmask(:,:,:) 
    118  
    119 !!gm NB : here a lbc_lnk should probably be added 
    120  
    121109      ! 
    122110      SELECT CASE( ktrd ) 
    … …  
    151139                              CALL iom_put( "vtrd_zdf", pvtrd ) 
    152140                              !                                    ! wind stress trends 
    153                               z2dx(:,:) = ( utau_b(ji,jj) + utau(ji,jj) ) / ( fse3u(:,:,1) * rau0 ) 
    154                               z2dy(:,:) = ( vtau_b(ji,jj) + vtau(ji,jj) ) / ( fse3v(:,:,1) * rau0 ) 
     141                              z2dx(:,:) = ( utau_b(:,:) + utau(:,:) ) / ( fse3u(:,:,1) * rau0 ) 
     142                              z2dy(:,:) = ( vtau_b(:,:) + vtau(:,:) ) / ( fse3v(:,:,1) * rau0 ) 
    155143                              CALL iom_put( "utrd_tau", z2dx ) 
    156144                              CALL iom_put( "vtrd_tau", z2dy ) 

  • branches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/NEMO/OPA_SRC/TRD/trdglo.F90

    r3318 r3325  
    184184 
    185185 
    186    SUBROUTINE trd_glo_init 
    187       !!--------------------------------------------------------------------- 
    188       !!                  ***  ROUTINE trd_glo_init  *** 
    189       !!  
    190       !! ** Purpose :   Read the namtrd namelist 
    191       !!---------------------------------------------------------------------- 
    192       INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    193       !!---------------------------------------------------------------------- 
    194  
    195       IF(lwp) THEN 
    196          WRITE(numout,*) 
    197          WRITE(numout,*) 'trd_glo_init : integral constraints properties trends' 
    198          WRITE(numout,*) '~~~~~~~~~~~~~' 
    199       ENDIF 
    200  
    201       ! Total volume at t-points: 
    202       tvolt = 0._wp 
    203       DO jk = 1, jpkm1 
    204          tvolt = SUM( e1e2t(:,:) * fse3t(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:) ) 
    205       END DO 
    206       IF( lk_mpp )   CALL mpp_sum( tvolt )   ! sum over the global domain 
    207  
    208       IF(lwp) WRITE(numout,*) '                total ocean volume at T-point   tvolt = ',tvolt 
    209  
    210       ! Initialization of potential to kinetic energy conversion 
    211       rpktrd = 0._wp 
    212  
    213       ! Total volume at u-, v- points: 
    214       tvolu = 0._wp 
    215       tvolv = 0._wp 
    216  
    217       DO jk = 1, jpk 
    218          DO jj = 2, jpjm1 
    219             DO ji = fs_2, fs_jpim1   ! vector opt. 
    220                tvolu = tvolu + e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) * tmask_i(ji+1,jj  ) * tmask_i(ji,jj) * umask(ji,jj,jk) 
    221                tvolv = tvolv + e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) * tmask_i(ji  ,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,jk) 
    222             END DO 
    223          END DO 
    224       END DO 
    225       IF( lk_mpp )   CALL mpp_sum( tvolu )   ! sums over the global domain 
    226       IF( lk_mpp )   CALL mpp_sum( tvolv ) 
    227  
    228       IF(lwp) THEN 
    229          WRITE(numout,*) '                total ocean volume at U-point   tvolu = ',tvolu 
    230          WRITE(numout,*) '                total ocean volume at V-point   tvolv = ',tvolv 
    231       ENDIF 
    232       ! 
    233    END SUBROUTINE trd_glo_init 
    234  
    235  
    236186   SUBROUTINE glo_dyn_wri( kt ) 
    237187      !!--------------------------------------------------------------------- 
    … …  
    556506   END SUBROUTINE glo_tra_wri 
    557507 
     508 
     509   SUBROUTINE trd_glo_init 
     510      !!--------------------------------------------------------------------- 
     511      !!                  ***  ROUTINE trd_glo_init  *** 
     512      !!  
     513      !! ** Purpose :   Read the namtrd namelist 
     514      !!---------------------------------------------------------------------- 
     515      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
     516      !!---------------------------------------------------------------------- 
     517 
     518      IF(lwp) THEN 
     519         WRITE(numout,*) 
     520         WRITE(numout,*) 'trd_glo_init : integral constraints properties trends' 
     521         WRITE(numout,*) '~~~~~~~~~~~~~' 
     522      ENDIF 
     523 
     524      ! Total volume at t-points: 
     525      tvolt = 0._wp 
     526      DO jk = 1, jpkm1 
     527         tvolt = SUM( e1e2t(:,:) * fse3t(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:) ) 
     528      END DO 
     529      IF( lk_mpp )   CALL mpp_sum( tvolt )   ! sum over the global domain 
     530 
     531      IF(lwp) WRITE(numout,*) '                total ocean volume at T-point   tvolt = ',tvolt 
     532 
     533      ! Initialization of potential to kinetic energy conversion 
     534      rpktrd = 0._wp 
     535 
     536      ! Total volume at u-, v- points: 
     537!!gm :  bug?  je suis quasi sur que le produit des tmask_i ne correspond pas exactement au umask_i et vmask_i ! 
     538      tvolu = 0._wp 
     539      tvolv = 0._wp 
     540 
     541      DO jk = 1, jpk 
     542         DO jj = 2, jpjm1 
     543            DO ji = fs_2, fs_jpim1   ! vector opt. 
     544               tvolu = tvolu + e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) * tmask_i(ji+1,jj  ) * tmask_i(ji,jj) * umask(ji,jj,jk) 
     545               tvolv = tvolv + e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) * tmask_i(ji  ,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,jk) 
     546            END DO 
     547         END DO 
     548      END DO 
     549      IF( lk_mpp )   CALL mpp_sum( tvolu )   ! sums over the global domain 
     550      IF( lk_mpp )   CALL mpp_sum( tvolv ) 
     551 
     552      IF(lwp) THEN 
     553         WRITE(numout,*) '                total ocean volume at U-point   tvolu = ',tvolu 
     554         WRITE(numout,*) '                total ocean volume at V-point   tvolv = ',tvolv 
     555      ENDIF 
     556      ! 
     557   END SUBROUTINE trd_glo_init 
     558 
    558559   !!====================================================================== 
    559560END MODULE trdglo 

  • branches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/NEMO/OPA_SRC/TRD/trdini.F90

    r3318 r3325  
    1111   !!---------------------------------------------------------------------- 
    1212   USE trd_oce        ! trends: ocean variables 
    13 !   USE ldftra_oce     ! ocean active tracers lateral physics 
    14    USE trdglo         ! ocean bassin integral constraints properties 
     13   USE trdken         ! trends: 3D kinetic energy 
     14   USE trdglo         ! trends: global domain integral constraints properties 
    1515   USE trdmld         ! ocean active mixed layer tracers trends  
    1616   USE trdvor         ! ocean vorticity trends  
    … …  
    7979      ! 
    8080       
    81       IF( ln_KE_trd .OR. ln_PE_trd .OR. ln_dyn_mld )   & 
    82          CALL ctl_stop( 'KE, PE, aur ML on momentum are not yet coded we stop' ) 
     81      IF( ln_PE_trd .OR. ln_dyn_mld )   & 
     82         CALL ctl_stop( 'PE, or ML on momentum are not yet coded we stop' ) 
    8383 
    8484      ! 
    85       IF( ln_glo_trd )    CALL trd_glo_init       ! integral constraints trends 
    86       IF( ln_tra_mld )    CALL trd_mld_init       ! mixed-layer trends (active  tracers)   
    87       IF( ln_vor_trd )    CALL trd_vor_init       ! vorticity trends         
     85      IF( ln_glo_trd )   CALL trd_glo_init      ! integral constraints trends 
     86      IF( ln_tra_mld )   CALL trd_mld_init      ! mixed-layer trends (active  tracers)   
     87      IF( ln_vor_trd )   CALL trd_vor_init      ! vorticity trends 
     88      IF( ln_KE_trd  )   CALL trd_ken_init      ! 3D Kinetic energy diagnostics 
     89 
    8890      ! 
    8991!!gm  : Potential BUG : 3D output only for vector invariant form!  add a ctl_stop or code the flux form case 

  • branches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/NEMO/OPA_SRC/TRD/trdken.F90

    r3318 r3325  
    1 MODULE trddyn 
     1MODULE trdken 
    22   !!====================================================================== 
    3    !!                       ***  MODULE  trddyn  *** 
    4    !! Ocean diagnostics:  ocean dynamic trends 
     3   !!                       ***  MODULE  trdken  *** 
     4   !! Ocean diagnostics:  compute and output 3D kinetic energy trends 
    55   !!===================================================================== 
    6    !! History :  3.5  !  2012-02  (G. Madec) creation from trdmod: split DYN and TRA trends 
    7    !!                           and manage  3D trends output for U, V, and KE 
    8    !!---------------------------------------------------------------------- 
    9  
    10    !!---------------------------------------------------------------------- 
    11    !!   trd_dyn       : manage the type of momentum trend diagnostics (3D I/O, domain averaged, KE) 
    12    !!   trd_dyn_iom   : output 3D momentum and/or tracer trends using IOM 
    13    !!   trd_dyn_init  : initialization step 
     6   !! History :  3.5  !  2012-02  (G. Madec) original code 
     7   !!---------------------------------------------------------------------- 
     8 
     9   !!---------------------------------------------------------------------- 
     10   !!   trd_ken       : compute and output 3D Kinetic energy trends using IOM 
     11   !!   trd_ken_init  : initialisation 
    1412   !!---------------------------------------------------------------------- 
    1513   USE oce            ! ocean dynamics and tracers variables 
    … …  
    3230   PRIVATE 
    3331 
    34    REAL(wp) ::   r2dt    ! time-step, = 2 rdttra except at nit000 (=rdttra) if neuler=0 
    35  
    36    PUBLIC trd_dyn        ! called by all dynXX modules 
     32   PUBLIC   trd_ken       ! called by trddyn module 
     33   PUBLIC   trd_ken_init  ! called by trdini module 
     34 
     35   INTEGER ::   nkstp       ! current time step  
     36   REAL(wp)::   r1_2_rau0   ! = 0.5 / rau0  
     37 
     38   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   bu, bv   ! volume of u- and v-boxes 
     39   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   r1_bt    ! inverse of t-box volume 
    3740 
    3841   !! * Substitutions 
    … …  
    4649CONTAINS 
    4750 
    48    SUBROUTINE trd_dyn( putrd, pvtrd, ktrd, kt ) 
    49       !!--------------------------------------------------------------------- 
    50       !!                  ***  ROUTINE trd_mod  *** 
     51   INTEGER FUNCTION trd_ken_alloc() 
     52      !!--------------------------------------------------------------------- 
     53      !!                  ***  FUNCTION trd_ken_alloc  *** 
     54      !!--------------------------------------------------------------------- 
     55      ALLOCATE( bu(jpi,jpj,jpk) , bv(jpi,jpj,jpk) , r1_bt(jpi,jpj,jpk) , STAT= trd_ken_alloc ) 
     56      ! 
     57      IF( lk_mpp             )   CALL mpp_sum ( trd_ken_alloc ) 
     58      IF( trd_ken_alloc /= 0 )   CALL ctl_warn('trd_ken_alloc: failed to allocate arrays') 
     59   END FUNCTION trd_ken_alloc 
     60 
     61 
     62   SUBROUTINE trd_ken( putrd, pvtrd, ktrd, kt ) 
     63      !!--------------------------------------------------------------------- 
     64      !!                  ***  ROUTINE trd_ken  *** 
    5165      !!  
    52       !! ** Purpose :   Dispatch momentum trend computation, e.g. 3D output,  
    53       !!              integral constraints, barotropic vorticity, kinetic enrgy,  
    54       !!              and/or mixed layer budget. 
    55       !!---------------------------------------------------------------------- 
    56       REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   putrd, pvtrd   ! U and V trends  
     66      !! ** Purpose :   output 3D Kinetic Energy trends using IOM 
     67      !! 
     68      !! ** Method  : - apply lbc to the input masked velocity trends  
     69      !!              - compute the associated KE trend: 
     70      !!          zke = 0.5 * (  mi-1[ un * putrd * bu ] + mj-1[ vn * pvtrd * bv]  ) / bt 
     71      !!      where bu, bv, bt are the volume of u-, v- and t-boxes.  
     72      !!              - vertical diffusion case (jpdyn_zdf):  
     73      !!          diagnose separately the KE trend associated with wind stress 
     74      !!              - bottom friction case (jpdyn_bfr): 
     75      !!          explicit case (ln_bfrimp=F): bottom trend put in the 1st level  
     76      !!                                       of putrd, pvtrd 
     77      ! 
     78      ! 
     79      !!---------------------------------------------------------------------- 
     80      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   putrd, pvtrd   ! U and V masked trends 
    5781      INTEGER                   , INTENT(in   ) ::   ktrd           ! trend index 
    5882      INTEGER                   , INTENT(in   ) ::   kt             ! time step 
    59       !! 
    60       INTEGER ::   ji, jj   ! dummy loop indices 
    61       REAL(wp), POINTER, DIMENSION(:,:) ::   ztswu, ztswv    ! 2D workspace  
    62       !!---------------------------------------------------------------------- 
    63  
    64       CALL wrk_alloc( jpi, jpj, ztswu, ztswv ) 
    65  
    66       IF( neuler == 0 .AND. kt == nit000    ) THEN   ;   r2dt =      rdt      ! = rdtra (restart with Euler time stepping) 
    67       ELSEIF(               kt <= nit000 + 1) THEN   ;   r2dt = 2. * rdt      ! = 2 rdttra (leapfrog) 
     83      ! 
     84      INTEGER ::   ji, jj, jk       ! dummy loop indices 
     85      INTEGER ::   ikbu  , ikbv     ! local integers 
     86      INTEGER ::   ikbum1, ikbvm1   !   -       - 
     87      REAL(wp), POINTER, DIMENSION(:,:)   ::   z2dx, z2dy, zke2d   ! 2D workspace  
     88      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zke                 ! 3D workspace  
     89      !!---------------------------------------------------------------------- 
     90      ! 
     91      CALL lbc_lnk( putrd, 'U', -1. )   ;   CALL lbc_lnk( pvtrd, 'V', -1. )      ! lateral boundary conditions 
     92      ! 
     93      IF ( lk_vvl .AND. kt /= nkstp ) THEN   ! Variable volume: set box volume at the 1st call of kt time step 
     94         nkstp = kt 
     95         DO jk = 1, jpkm1 
     96            bu   (:,:,jk) =  e1u(:,:) * e2u(:,:) * fse3u_n(:,:,jk) 
     97            bv   (:,:,jk) =  e1v(:,:) * e2v(:,:) * fse3v_n(:,:,jk) 
     98            r1_bt(:,:,jk) = 1._wp / ( e1e2t(:,:) * fse3t_n(:,:,jk) ) * tmask(:,:,jk) 
     99         END DO 
    68100      ENDIF 
    69  
    70       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    71       !   3D output of momentum and/or tracers trends using IOM interface 
    72       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    73       IF( ln_dyn_trd )   CALL trd_dyn_iom( putrd, pvtrd, ktrd, kt ) 
    74           
    75       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    76       !  Integral Constraints Properties for momentum and/or tracers trends 
    77       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    78       IF( ln_glo_trd  )   CALL trd_glo( putrd, pvtrd, ktrd, 'DYN', kt ) 
    79  
    80       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    81       !  Kinetic Energy trends 
    82       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    83 !!gm      IF( ln_KE_trd  )   CALL trd_KE( putrd, pvtrd, ktrd, kt ) 
    84  
    85       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    86       !  Vorticity trends 
    87       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    88       IF( ln_vor_trd  )   CALL trd_vor( putrd, pvtrd, ktrd, kt ) 
    89  
    90       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    91       !  Mixed layer trends for active tracers 
    92       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    93 !!gm      IF( ln_dyn_mld )   CALL trd_mld_dyn    
    94       ! 
    95       CALL wrk_dealloc( jpi, jpj, ztswu, ztswv ) 
    96       ! 
    97    END SUBROUTINE trd_dyn 
    98  
    99  
    100    SUBROUTINE trd_dyn_iom( putrd, pvtrd, ktrd, kt ) 
    101       !!--------------------------------------------------------------------- 
    102       !!                  ***  ROUTINE trd_dyn_iom  *** 
    103       !!  
    104       !! ** Purpose :   output 3D trends using IOM 
    105       !!---------------------------------------------------------------------- 
    106       REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   putrd, pvtrd   ! U and V trends 
    107       INTEGER                   , INTENT(in   ) ::   ktrd           ! trend index 
    108       INTEGER                   , INTENT(in   ) ::   kt             ! time step 
    109       ! 
    110       INTEGER ::   ji, jj, jk   ! dummy loop indices 
    111       INTEGER ::   ikbu, ikbv   ! local integers 
    112       REAL(wp), POINTER, DIMENSION(:,:)   ::   z2dx, z2dy, ztswu, ztswv   ! 2D workspace  
    113       REAL(wp), POINTER, DIMENSION(:,:,:) ::   z3dx, z3dy                 ! 3D workspace  
    114       !!---------------------------------------------------------------------- 
    115       ! 
    116       putrd(:,:,:) = putrd(:,:,:) * umask(:,:,:)                       ! mask the trends 
    117       pvtrd(:,:,:) = pvtrd(:,:,:) * vmask(:,:,:) 
    118  
    119 !!gm NB : here a lbc_lnk should probably be added 
    120  
    121       ! 
    122       SELECT CASE( ktrd ) 
    123       CASE( jpdyn_hpg )   ;   CALL iom_put( "utrd_hpg", putrd )    ! hydrostatic pressure gradient 
    124                               CALL iom_put( "vtrd_hpg", pvtrd ) 
    125       CASE( jpdyn_spg )   ;   CALL iom_put( "utrd_spg", putrd )    ! surface pressure gradient 
    126                               CALL iom_put( "vtrd_spg", pvtrd ) 
    127       CASE( jpdyn_pvo )   ;   CALL iom_put( "utrd_pvo", putrd )    ! planetary vorticity 
    128                               CALL iom_put( "vtrd_pvo", pvtrd ) 
    129       CASE( jpdyn_rvo )   ;   CALL iom_put( "utrd_rvo", putrd )    ! relative  vorticity     (or metric term) 
    130                               CALL iom_put( "vtrd_rvo", pvtrd ) 
    131       CASE( jpdyn_keg )   ;   CALL iom_put( "utrd_keg", putrd )    ! Kinetic Energy gradient (or had) 
    132                               CALL iom_put( "vtrd_keg", pvtrd ) 
    133          z3dx(:,:,:) = 0._wp                                           ! U.dxU & V.dyV (approximation) 
    134          z3dy(:,:,:) = 0._wp 
    135          DO jk = 1, jpkm1                                                  ! no mask as un,vn are masked 
    136             DO jj = 2, jpjm1 
    137                 DO ji = 2, jpim1 
    138                   z3dx(ji,jj,jk) = un(ji,jj,jk) * ( un(ji+1,jj,jk) - un(ji-1,jj,jk) ) / ( 2._wp * e1u(ji,jj) ) 
    139                   z3dy(ji,jj,jk) = vn(ji,jj,jk) * ( vn(ji,jj+1,jk) - vn(ji,jj-1,jk) ) / ( 2._wp * e2v(ji,jj) ) 
    140                END DO 
     101      ! 
     102      zke(:,:,jpk) = 0._wp 
     103      zke(1,:, : ) = 0._wp 
     104      zke(:,1, : ) = 0._wp 
     105      DO jk = 1, jpkm1 
     106         DO jj = 2, jpj 
     107            DO ji = 2, jpj 
     108               zke(ji,jj,jk) = 0.5_wp * (   un(ji  ,jj,jk) * putrd(ji  ,jj,jk) * bu(ji  ,jj,jk)  & 
     109                  &                       + un(ji-1,jj,jk) * putrd(ji-1,jj,jk) * bu(ji-1,jj,jk)  & 
     110                  &                       + vn(ji,jj  ,jk) * pvtrd(ji,jj  ,jk) * bv(ji,jj  ,jk)  & 
     111                  &                       + vn(ji,jj-1,jk) * pvtrd(ji,jj-1,jk) * bv(ji,jj-1,jk)  ) * r1_bt(ji,jj,jk) 
    141112            END DO 
    142113         END DO 
    143          CALL lbc_lnk( z3dx, 'U', -1. )   ;    CALL lbc_lnk( z3dy, 'V', -1. ) 
    144                               CALL iom_put( "utrd_udx", z3dx  )  
    145                               CALL iom_put( "vtrd_vdy", z3dy  ) 
    146       CASE( jpdyn_zad )   ;   CALL iom_put( "utrd_zad", putrd )    ! vertical   advection 
    147                               CALL iom_put( "vtrd_zad", pvtrd ) 
    148       CASE( jpdyn_ldf )   ;   CALL iom_put( "utrd_ldf", putrd )    ! lateral diffusion 
    149                               CALL iom_put( "vtrd_ldf", pvtrd ) 
    150       CASE( jpdyn_zdf )   ;   CALL iom_put( "utrd_zdf", putrd )    ! vertical diffusion  
    151                               CALL iom_put( "vtrd_zdf", pvtrd ) 
    152                               !                                    ! wind stress trends 
    153                               z2dx(:,:) = ( utau_b(ji,jj) + utau(ji,jj) ) / ( fse3u(:,:,1) * rau0 ) 
    154                               z2dy(:,:) = ( vtau_b(ji,jj) + vtau(ji,jj) ) / ( fse3v(:,:,1) * rau0 ) 
    155                               CALL iom_put( "utrd_tau", z2dx ) 
    156                               CALL iom_put( "vtrd_tau", z2dy ) 
    157       CASE( jpdyn_bfr ) 
    158          IF( .NOT.ln_bfrimp )     CALL iom_put( "utrd_bfr", putrd )    ! bottom friction (explicit case) 
    159          IF( .NOT.ln_bfrimp )     CALL iom_put( "vtrd_bfr", pvtrd ) 
     114      END DO 
     115      ! 
     116      SELECT CASE( ktrd ) 
     117      CASE( jpdyn_hpg )   ;   CALL iom_put( "ketrd_hpg", zke )    ! hydrostatic pressure gradient 
     118      CASE( jpdyn_spg )   ;   CALL iom_put( "ketrd_spg", zke )    ! surface pressure gradient 
     119      CASE( jpdyn_pvo )   ;   CALL iom_put( "ketrd_pvo", zke )    ! planetary vorticity 
     120      CASE( jpdyn_rvo )   ;   CALL iom_put( "ketrd_rvo", zke )    ! relative  vorticity     (or metric term) 
     121      CASE( jpdyn_keg )   ;   CALL iom_put( "ketrd_keg", zke )    ! Kinetic Energy gradient (or had) 
     122      CASE( jpdyn_zad )   ;   CALL iom_put( "ketrd_zad", zke )    ! vertical   advection 
     123      CASE( jpdyn_ldf )   ;   CALL iom_put( "ketrd_ldf", zke )    ! lateral diffusion 
     124      CASE( jpdyn_zdf )   ;   CALL iom_put( "ketrd_zdf", zke )    ! vertical diffusion  
     125                              !                                   ! wind stress trends 
     126         z2dx(:,:) = un(:,:,1) * ( utau_b(:,:) + utau(:,:) ) * e1u(:,:) * e2u(:,:) * umask(:,:,1) 
     127         z2dy(:,:) = vn(:,:,1) * ( vtau_b(:,:) + vtau(:,:) ) * e1v(:,:) * e2v(:,:) * vmask(:,:,1) 
     128         zke2d(1,:) = 0._wp   ;   zke2d(:,1) = 0._wp 
     129         DO jj = 2, jpj 
     130            DO ji = 2, jpj 
     131               zke2d(ji,jj) = r1_2_rau0 * (   z2dx(ji,jj) + z2dx(ji-1,jj)   & 
     132                  &                         + z2dy(ji,jj) + z2dy(ji,jj-1)   ) * r1_bt(ji,jj,1) 
     133            END DO            ! 
     134         END DO               ! 
     135                              CALL iom_put( "ketrd_tau", zke2d ) 
     136      CASE( jpdyn_bfr )   ;   CALL iom_put( "ketrd_bfr", zke )    ! bottom friction (explicit)  
     137!!gm TO BE DONE properly 
     138!         IF( .NOT.ln_bfrimp ) THEN 
     139!            DO jj = 1, jpj    !    
     140!               DO ji = 1, jpi 
     141!                  ikbu = mbku(ji,jj)         ! deepest ocean u- & v-levels 
     142!                  ikbv = mbkv(ji,jj)    
     143!                  z2dx(ji,jj) = un(ji,jj,ikbu) * bfrua(ji,jj) * un(ji,jj,ikbu) 
     144!                  z2dy(ji,jj) = vn(ji,jj,ikbu) * bfrva(ji,jj) * vn(ji,jj,ikbv) 
     145!               END DO 
     146!            END DO 
     147!            zke2d(1,:) = 0._wp   ;   zke2d(:,1) = 0._wp 
     148!            DO jj = 2, jpj 
     149!               DO ji = 2, jpj 
     150!                  zke2d(ji,jj) = 0.5_wp * (   z2dx(ji,jj) + z2dx(ji-1,jj)   & 
     151!                     &                      + z2dy(ji,jj) + z2dy(ji,jj-1)   ) * r1_bt(ji,jj,  BEURK!!! 
     152!               END DO 
     153!            END DO 
     154!                              CALL iom_put( "ketrd_bfr", zke2d )    ! bottom friction (explicit case) 
    160155!!gm only valid if ln_bfrimp=T otherwise the bottom stress as to be recomputed at the end of the compuation.... 
    161  
    162       CASE( jpdyn_atf )   ;   CALL iom_put( "utrd_atf", putrd )    ! asselin filter trends  
    163                               CALL iom_put( "vtrd_atf", pvtrd ) 
    164          IF( ln_bfrimp ) THEN                                          ! bottom friction (implicit case) 
    165             z3dx(:,:,:) = 0._wp   ;   z3dy(:,:,:) = 0._wp                 ! after velocity known (now filed at this stage) 
    166             DO jk = 1, jpkm1 
    167                DO jj = 2, jpjm1 
    168                   DO ji = 2, jpim1 
    169                      ikbu = mbku(ji,jj)          ! deepest ocean u- & v-levels 
    170                      ikbv = mbkv(ji,jj) 
    171                      z3dx(ji,jj,jk) = bfrua(ji,jj) * un(ji,jj,ikbu) / fse3u(ji,jj,ikbu) 
    172                      z3dy(ji,jj,jk) = bfrva(ji,jj) * vn(ji,jj,ikbv) / fse3v(ji,jj,ikbv) 
    173                   END DO 
    174                END DO 
    175             END DO 
    176                               CALL iom_put( "utrd_bfr", z3dx )    ! bottom friction (implicit) 
    177                               CALL iom_put( "vtrd_bfr", z3dy ) 
    178          ENDIF 
     156!         ENDIF 
     157!!gm end 
     158      CASE( jpdyn_atf )   ;   CALL iom_put( "ketrd_atf", zke )    ! asselin filter trends  
     159 
     160 
     161!! a faire !!!!  idee changer dynnxt pour avoir un appel a jpdynbfr avant le swap !!! 
     162!! reflechir à une possible sauvegarde du "vrai" un,vn pour le calcul de atf.... 
     163! 
     164!         IF( ln_bfrimp ) THEN                                          ! bottom friction (implicit case) 
     165!            DO jj = 1, jpj                                                  ! after velocity known (now filed at this stage) 
     166!               DO ji = 1, jpi 
     167!                  ikbu = mbku(ji,jj)          ! deepest ocean u- & v-levels 
     168!                  ikbv = mbkv(ji,jj) 
     169!                  z2dx(ji,jj) = un(ji,jj,ikbu) * bfrua(ji,jj) * un(ji,jj,ikbu) / fse3u(ji,jj,ikbu) 
     170!                  z2dy(ji,jj) = un(ji,jj,ikbu) * bfrva(ji,jj) * vn(ji,jj,ikbv) / fse3v(ji,jj,ikbv) 
     171!               END DO 
     172!            END DO 
     173!            zke2d(1,:) = 0._wp   ;   zke2d(:,1) = 0._wp 
     174!            DO jj = 2, jpj 
     175!               DO ji = 2, jpj 
     176!                  zke2d(ji,jj) = 0.5_wp * (   z2dx(ji,jj) + z2dx(ji-1,jj)   & 
     177!                     &                      + z2dy(ji,jj) + z2dy(ji,jj-1)   ) 
     178!               END DO 
     179!            END DO 
     180!                              CALL iom_put( "ketrd_bfr", zke2d ) 
     181!         ENDIF 
    179182         ! 
    180183      END SELECT 
    181184      ! 
    182       CALL wrk_dealloc( jpi, jpj     , z2dx, z2dy, ztswu, ztswv ) 
    183       CALL wrk_dealloc( jpi, jpj, jpk, z3dx, z3dy ) 
    184       ! 
    185    END SUBROUTINE trd_dyn_iom 
     185      CALL wrk_dealloc( jpi, jpj     , z2dx, z2dy, zke2d ) 
     186      CALL wrk_dealloc( jpi, jpj, jpk, zke ) 
     187      ! 
     188   END SUBROUTINE trd_ken 
     189 
     190 
     191   SUBROUTINE trd_ken_init 
     192      !!--------------------------------------------------------------------- 
     193      !!                  ***  ROUTINE trd_ken_init  *** 
     194      !!  
     195      !! ** Purpose :   initialisation of 3D Kinetic Energy trend diagnostic 
     196      !!---------------------------------------------------------------------- 
     197      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
     198      !!---------------------------------------------------------------------- 
     199      ! 
     200      IF(lwp) THEN 
     201         WRITE(numout,*) 
     202         WRITE(numout,*) 'trd_ken_init : 3D Kinetic Energy trends' 
     203         WRITE(numout,*) '~~~~~~~~~~~~~' 
     204      ENDIF 
     205      !                           ! allocate box volume arrays 
     206      IF ( trd_ken_alloc() /= 0 )   CALL ctl_stop('trd_ken_alloc: failed to allocate arrays') 
     207      ! 
     208      IF ( .NOT.lk_vvl ) THEN     ! constant volume: bu, bv, 1/bt computed one for all 
     209         DO jk = 1, jpkm1 
     210            bu   (:,:,jk) =  e1u(:,:) * e2u(:,:) * fse3u_n(:,:,jk) 
     211            bv   (:,:,jk) =  e1v(:,:) * e2v(:,:) * fse3v_n(:,:,jk) 
     212            r1_bt(:,:,jk) = 1._wp / ( e1e2t(:,:) * fse3t_n(:,:,jk) ) 
     213         END DO 
     214      ENDIF 
     215      ! 
     216      nkstp     = 0 
     217      r1_2_rau0 = 0.5_wp / rau0 
     218      ! 
     219   END SUBROUTINE trd_ken_init 
    186220 
    187221   !!====================================================================== 
    188 END MODULE trddyn 
     222END MODULE trdken 

  • branches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/NEMO/OPA_SRC/TRD/trdvor.F90

    r3318 r3325  
    116116         CALL trd_vor_zint( putrd, pvtrd, jpvor_zdf )                             ! zdf trend including surf./bot. stresses  
    117117         CALL trd_vor_zint( ztswu, ztswv, jpvor_swf )                             ! surface wind stress  
    118       CASE ( jpdyn_bfr ) 
     118      CASE( jpdyn_bfr ) 
    119119         CALL trd_vor_zint( putrd, pvtrd, jpvor_bfr )                             ! Bottom stress 
    120120         ! 
    … …  
    176176      !  ===================================== 
    177177 
    178       SELECT CASE (ktrd)  
    179       ! 
    180       CASE (jpvor_bfr)        ! bottom friction 
     178      SELECT CASE( ktrd )  
     179      ! 
     180      CASE( jpvor_bfr )        ! bottom friction 
    181181         DO jj = 2, jpjm1 
    182182            DO ji = fs_2, fs_jpim1  
    … …  
    188188         END DO 
    189189         ! 
    190       CASE (jpvor_swf)        ! wind stress 
     190      CASE( jpvor_swf )        ! wind stress 
    191191         zudpvor(:,:) = putrdvor(:,:) * fse3u(:,:,1) * e1u(:,:) * umask(:,:,1) 
    192192         zvdpvor(:,:) = pvtrdvor(:,:) * fse3v(:,:,1) * e2v(:,:) * vmask(:,:,1) 
    … …  
    199199    
    200200      ! Curl 
    201       DO ji=1,jpim1 
    202          DO jj=1,jpjm1 
     201      DO ji = 1, jpim1 
     202         DO jj = 1, jpjm1 
    203203            vortrd(ji,jj,ktrd) = (    zvdpvor(ji+1,jj) - zvdpvor(ji,jj)       & 
    204204                 &                - ( zudpvor(ji,jj+1) - zudpvor(ji,jj) )   ) / ( e1f(ji,jj) * e2f(ji,jj) ) 

  • branches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/NEMO/OPA_SRC/TRD/trdvor_oce.F90

    r3318 r3325  
    44   !! Ocean trends :   set vorticity trend variables 
    55   !!====================================================================== 
    6    !! History :  9.0  !  04-2006  (L. Brunier, A-M. Treguier) Original code  
     6   !! History :  1.0  !  04-2006  (L. Brunier, A-M. Treguier) Original code  
    77   !!---------------------------------------------------------------------- 
    8  
    9    !!---------------------------------------------------------------------- 
     8    
    109   USE par_oce      ! ocean parameters 
    1110 
Note: See TracChangeset for help on using the changeset viewer.