Ignore:
Timestamp:
2019-09-30T11:07:57+02:00 (14 months ago)
Author:
davestorkey
Message:

UKMO/NEMO_4.0_momentum_trends branch : first set of code changes.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/UKMO/NEMO_4.0_momentum_trends/src/OCE/TRD/trddyn.F90

    r10888 r11613  
    1818   USE sbc_oce        ! surface boundary condition: ocean 
    1919   USE zdf_oce        ! ocean vertical physics: variables 
    20 !!gm   USE zdfdrg         ! ocean vertical physics: bottom friction 
    2120   USE trd_oce        ! trends: ocean variables 
    2221   USE trdken         ! trends: Kinetic ENergy  
     
    3534   PUBLIC trd_dyn        ! called by all dynXXX modules 
    3635 
     36   INTERFACE trd_dyn 
     37      module procedure trd_dyn_3d, trd_dyn_2d 
     38   END INTERFACE 
     39 
     40   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), SAVE :: zutrd_hpg, zvtrd_hpg 
     41   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), SAVE :: zutrd_pvo, zvtrd_pvo 
     42 
    3743   !! * Substitutions 
    3844#  include "vectopt_loop_substitute.h90" 
     
    4450CONTAINS 
    4551 
    46    SUBROUTINE trd_dyn( putrd, pvtrd, ktrd, kt ) 
     52   SUBROUTINE trd_dyn_3d( putrd, pvtrd, ktrd, kt ) 
    4753      !!--------------------------------------------------------------------- 
    48       !!                  ***  ROUTINE trd_mod  *** 
     54      !!                  ***  ROUTINE trd_dyn_3d  *** 
    4955      !!  
    5056      !! ** Purpose :   Dispatch momentum trend computation, e.g. 3D output,  
     
    6369!!gm NB : here a lbc_lnk should probably be added 
    6470 
    65       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    66       !   3D output of momentum and/or tracers trends using IOM interface 
    67       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    68       IF( ln_dyn_trd )   CALL trd_dyn_iom( putrd, pvtrd, ktrd, kt ) 
     71      SELECT CASE( ktrd ) 
     72      CASE( jpdyn_hpg_save )  
     73         ! 
     74         ! save 3D HPG trends to possibly have barotropic part corrected later before writing out 
     75         ALLOCATE( zutrd_hpg(jpi,jpj,jpk), zvtrd_hpg(jpi,jpj,jpk) ) 
     76         zutrd_hpg(:,:,:) = putrd(:,:,:) 
     77         zvtrd_hpg(:,:,:) = pvtrd(:,:,:) 
     78 
     79      CASE( jpdyn_pvo_save )  
     80         ! 
     81         ! save 3D coriolis trends to possibly have barotropic part corrected later before writing out 
     82         ALLOCATE( zutrd_pvo(jpi,jpj,jpk), zvtrd_pvo(jpi,jpj,jpk) ) 
     83         zutrd_pvo(:,:,:) = putrd(:,:,:) 
     84         zvtrd_pvo(:,:,:) = pvtrd(:,:,:) 
     85 
     86      CASE DEFAULT 
     87 
     88         IF( ktrd == jpdyn_spg ) THEN 
     89            ! For explicit scheme SPG trends come here as 3D fields 
     90            ! Add SPG trend to 3D HPG trend and also output as 2D diagnostic in own right. 
     91            CALL trd_dyn_iom_2d( putrd(:,:,1), pvtrd(:,:,1), jpdyn_spg, kt )  
     92            putrd(:,:,:) = putrd(:,:,:) + zutrd_hpg(:,:,:)  
     93            pvtrd(:,:,:) = pvtrd(:,:,:) + zvtrd_hpg(:,:,:)  
     94            DEALLOCATE( zutrd_hpg, zvtrd_hpg ) 
     95         ENDIF 
     96 
     97         !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     98         !   3D output of momentum and/or tracers trends using IOM interface 
     99         !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     100         IF( ln_dyn_trd )   CALL trd_dyn_iom_3d( putrd, pvtrd, ktrd, kt ) 
     101 
     102         !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     103         !  Integral Constraints Properties for momentum and/or tracers trends 
     104         !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     105         IF( ln_glo_trd )   CALL trd_glo( putrd, pvtrd, ktrd, 'DYN', kt ) 
     106 
     107         !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     108         !  Kinetic Energy trends 
     109         !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     110         IF( ln_KE_trd  )   CALL trd_ken( putrd, pvtrd, ktrd, kt ) 
     111 
     112         !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     113         !  Vorticity trends 
     114         !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     115         IF( ln_vor_trd )   CALL trd_vor( putrd, pvtrd, ktrd, kt ) 
     116 
     117         !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     118         !  Mixed layer trends for active tracers 
     119         !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     120         !!gm      IF( ln_dyn_mxl )   CALL trd_mxl_dyn    
     121         ! 
     122      END SELECT 
     123      ! 
     124   END SUBROUTINE trd_dyn_3d 
     125 
     126 
     127   SUBROUTINE trd_dyn_2d( putrd, pvtrd, ktrd, kt ) 
     128      !!--------------------------------------------------------------------- 
     129      !!                  ***  ROUTINE trd_mod  *** 
     130      !!  
     131      !! ** Purpose :   Dispatch momentum trend computation, e.g. 2D output,  
     132      !!              integral constraints, barotropic vorticity, kinetic enrgy,  
     133      !!              and/or mixed layer budget. 
     134      !!---------------------------------------------------------------------- 
     135      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   putrd, pvtrd   ! U and V trends  
     136      INTEGER                 , INTENT(in   ) ::   ktrd           ! trend index 
     137      INTEGER                 , INTENT(in   ) ::   kt             ! time step 
     138      INTEGER                                 ::   jk 
     139      !!---------------------------------------------------------------------- 
     140      ! 
     141      putrd(:,:) = putrd(:,:) * umask(:,:,1)                       ! mask the trends 
     142      pvtrd(:,:) = pvtrd(:,:) * vmask(:,:,1) 
     143      ! 
     144 
     145!!gm NB : here a lbc_lnk should probably be added 
     146 
     147      SELECT CASE(ktrd) 
     148 
     149      CASE ( jpdyn_hpg_corr ) 
     150         ! 
     151         ! Remove "first-guess" SPG trend from 3D HPG trend.  
     152         DO jk = 1, jpkm1 
     153            zutrd_hpg(:,:,jk) = zutrd_hpg(:,:,jk) - putrd(:,:) 
     154            zvtrd_hpg(:,:,jk) = zvtrd_hpg(:,:,jk) - pvtrd(:,:) 
     155         ENDDO 
     156         CALL trd_dyn_iom_2d( putrd, pvtrd, jpdyn_hpg_corr, kt )  
     157 
     158      CASE( jpdyn_pvo_corr ) 
     159         ! 
     160         ! Remove "first-guess" barotropic coriolis trend from 3D PVO trend.  
     161         DO jk = 1, jpkm1 
     162            zutrd_pvo(:,:,jk) = zutrd_pvo(:,:,jk) - putrd(:,:) 
     163            zvtrd_pvo(:,:,jk) = zvtrd_pvo(:,:,jk) - pvtrd(:,:) 
     164         ENDDO 
     165         CALL trd_dyn_iom_2d( putrd, pvtrd, jpdyn_pvo_corr, kt )  
     166 
     167      CASE( jpdyn_spg ) 
     168          ! 
     169          ! For split-explicit scheme SPG trends come here as 2D fields 
     170          ! Add SPG trend to 3D HPG trend and also output as 2D diagnostic in own right. 
     171          DO jk = 1, jpkm1 
     172             zutrd_hpg(:,:,jk) = zutrd_hpg(:,:,jk) + putrd(:,:) 
     173             zvtrd_hpg(:,:,jk) = zvtrd_hpg(:,:,jk) + pvtrd(:,:) 
     174          ENDDO 
     175          CALL trd_dyn_iom_2d( putrd(:,:), pvtrd(:,:), jpdyn_spg, kt )  
     176          CALL trd_dyn_3d( zutrd_hpg, zvtrd_hpg, jpdyn_hpg, kt ) 
     177          DEALLOCATE( zutrd_hpg, zvtrd_hpg ) 
     178 
     179      CASE( jpdyn_pvo ) 
     180          ! 
     181          ! Add 2D PVO trend to 3D PVO trend and also output as diagnostic in own right. 
     182          DO jk = 1, jpkm1 
     183             zutrd_pvo(:,:,jk) = zutrd_pvo(:,:,jk) + putrd(:,:) 
     184             zvtrd_pvo(:,:,jk) = zvtrd_pvo(:,:,jk) + pvtrd(:,:) 
     185          ENDDO 
     186          CALL trd_dyn_iom_2d( putrd, pvtrd, jpdyn_pvo, kt )  
     187          CALL trd_dyn_3d( zutrd_pvo, zvtrd_pvo, jpdyn_pvo, kt ) 
     188          DEALLOCATE( zutrd_pvo, zvtrd_pvo ) 
     189 
     190      CASE DEFAULT  
     191 
     192         !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     193         !   2D output of momentum and/or tracers trends using IOM interface 
     194         !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     195         IF( ln_dyn_trd )   CALL trd_dyn_iom_2d( putrd, pvtrd, ktrd, kt ) 
    69196          
    70       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    71       !  Integral Constraints Properties for momentum and/or tracers trends 
    72       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    73       IF( ln_glo_trd )   CALL trd_glo( putrd, pvtrd, ktrd, 'DYN', kt ) 
    74  
    75       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    76       !  Kinetic Energy trends 
    77       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    78       IF( ln_KE_trd  )   CALL trd_ken( putrd, pvtrd, ktrd, kt ) 
    79  
    80       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    81       !  Vorticity trends 
    82       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    83       IF( ln_vor_trd )   CALL trd_vor( putrd, pvtrd, ktrd, kt ) 
     197      END SELECT 
     198 
     199!!$   CALLS TO THESE ROUTINES FOR 2D DIAGOSTICS NOT CODED YET 
     200!!$      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     201!!$      !  Integral Constraints Properties for momentum and/or tracers trends 
     202!!$      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     203!!$      IF( ln_glo_trd )   CALL trd_glo( putrd, pvtrd, ktrd, 'DYN', kt ) 
     204!!$ 
     205!!$      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     206!!$      !  Kinetic Energy trends 
     207!!$      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     208!!$      IF( ln_KE_trd  )   CALL trd_ken( putrd, pvtrd, ktrd, kt ) 
     209!!$ 
     210!!$      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     211!!$      !  Vorticity trends 
     212!!$      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     213!!$      IF( ln_vor_trd )   CALL trd_vor( putrd, pvtrd, ktrd, kt ) 
    84214 
    85215      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     
    88218!!gm      IF( ln_dyn_mxl )   CALL trd_mxl_dyn    
    89219      ! 
    90    END SUBROUTINE trd_dyn 
    91  
    92  
    93    SUBROUTINE trd_dyn_iom( putrd, pvtrd, ktrd, kt ) 
     220   END SUBROUTINE trd_dyn_2d 
     221 
     222 
     223   SUBROUTINE trd_dyn_iom_3d( putrd, pvtrd, ktrd, kt ) 
    94224      !!--------------------------------------------------------------------- 
    95225      !!                  ***  ROUTINE trd_dyn_iom  *** 
     
    110240      CASE( jpdyn_hpg )   ;   CALL iom_put( "utrd_hpg", putrd )    ! hydrostatic pressure gradient 
    111241                              CALL iom_put( "vtrd_hpg", pvtrd ) 
    112       CASE( jpdyn_spg )   ;   CALL iom_put( "utrd_spg", putrd )    ! surface pressure gradient 
    113                               CALL iom_put( "vtrd_spg", pvtrd ) 
    114242      CASE( jpdyn_pvo )   ;   CALL iom_put( "utrd_pvo", putrd )    ! planetary vorticity 
    115243                              CALL iom_put( "vtrd_pvo", pvtrd ) 
     
    147275                              CALL iom_put( "vtrd_tau", z2dy ) 
    148276                              DEALLOCATE( z2dx , z2dy ) 
    149 !!gm  to be changed : computation should be done in dynzdf.F90 
    150 !!gm                + missing the top friction  
    151 !                              !                                    ! bottom stress tends (implicit case) 
    152 !                              IF( ln_drgimp ) THEN 
    153 !                                 ALLOCATE( z3dx(jpi,jpj,jpk) , z3dy(jpi,jpj,jpk) ) 
    154 !                             z3dx(:,:,:) = 0._wp   ;   z3dy(:,:,:) = 0._wp  ! after velocity known (now filed at this stage) 
    155 !                            DO jk = 1, jpkm1 
    156 !                                    DO jj = 2, jpjm1 
    157 !                                       DO ji = 2, jpim1 
    158 !                                      ikbu = mbku(ji,jj)          ! deepest ocean u- & v-levels 
    159 !                                          ikbv = mbkv(ji,jj) 
    160 !                                          z3dx(ji,jj,jk) = 0.5 * ( rCdU_bot(ji+1,jj) + rCdU_bot(ji,jj) ) &  
    161 !                                               &         * un(ji,jj,ikbu) / e3u_n(ji,jj,ikbu) 
    162 !                                          z3dy(ji,jj,jk) = 0.5 * ( rCdU_bot(ji,jj+1) + rCdU_bot(ji,jj) ) & 
    163 !                                               &         * vn(ji,jj,ikbv) / e3v_n(ji,jj,ikbv) 
    164 !                                    END DO 
    165 !                                 END DO 
    166 !                              END DO 
    167 !                              CALL lbc_lnk_multi( 'trddyn', z3dx, 'U', -1., z3dy, 'V', -1. ) 
    168 !                              CALL iom_put( "utrd_bfr", z3dx ) 
    169 !                              CALL iom_put( "vtrd_bfr", z3dy ) 
    170 !                                 DEALLOCATE( z3dx , z3dy ) 
    171 !                              ENDIF 
    172 !!gm end 
    173       CASE( jpdyn_bfr )       ! called if ln_drgimp=F 
    174                               CALL iom_put( "utrd_bfr", putrd )    ! bottom friction (explicit case) 
     277      CASE( jpdyn_bfr )   ;   CALL iom_put( "utrd_bfr", putrd )    ! bottom friction (explicit case) 
    175278                              CALL iom_put( "vtrd_bfr", pvtrd ) 
     279      CASE( jpdyn_bfri)   ;   CALL iom_put( "utrd_bfri", putrd )    ! bottom friction (implicit case) 
     280                              CALL iom_put( "vtrd_bfri", pvtrd ) 
    176281      CASE( jpdyn_atf )   ;   CALL iom_put( "utrd_atf", putrd )        ! asselin filter trends  
    177282                              CALL iom_put( "vtrd_atf", pvtrd ) 
    178283      END SELECT 
    179284      ! 
    180    END SUBROUTINE trd_dyn_iom 
     285   END SUBROUTINE trd_dyn_iom_3d 
     286 
     287 
     288   SUBROUTINE trd_dyn_iom_2d( putrd, pvtrd, ktrd, kt ) 
     289      !!--------------------------------------------------------------------- 
     290      !!                  ***  ROUTINE trd_dyn_iom  *** 
     291      !!  
     292      !! ** Purpose :   output 2D trends using IOM 
     293      !!---------------------------------------------------------------------- 
     294      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   putrd, pvtrd   ! U and V trends 
     295      INTEGER                 , INTENT(in   ) ::   ktrd           ! trend index 
     296      INTEGER                 , INTENT(in   ) ::   kt             ! time step 
     297      ! 
     298      INTEGER ::   ji, jj, jk   ! dummy loop indices 
     299      INTEGER ::   ikbu, ikbv   ! local integers 
     300      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   z2dx, z2dy   ! 2D workspace  
     301      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   z3dx, z3dy   ! 3D workspace  
     302      !!---------------------------------------------------------------------- 
     303      ! 
     304      SELECT CASE( ktrd ) 
     305      CASE( jpdyn_spg )      ;   CALL iom_put( "utrd_spg2d", putrd )      ! surface pressure gradient 
     306                                 CALL iom_put( "vtrd_spg2d", pvtrd ) 
     307      CASE( jpdyn_pvo )      ;   CALL iom_put( "utrd_pvo2d", putrd )      ! planetary vorticity (barotropic part) 
     308                                 CALL iom_put( "vtrd_pvo2d", pvtrd ) 
     309      CASE( jpdyn_hpg_corr ) ;   CALL iom_put( "utrd_hpg_corr", putrd )   ! horizontal pressure gradient correction 
     310                                 CALL iom_put( "vtrd_hpg_corr", pvtrd ) 
     311      CASE( jpdyn_pvo_corr ) ;   CALL iom_put( "utrd_pvo_corr", putrd )   ! planetary vorticity correction 
     312                                 CALL iom_put( "vtrd_pvo_corr", pvtrd ) 
     313      CASE( jpdyn_bfr )      ;   CALL iom_put( "utrd_bfr2d", putrd )      ! bottom friction due to barotropic currents 
     314                                 CALL iom_put( "vtrd_bfr2d", pvtrd ) 
     315      END SELECT 
     316      ! 
     317   END SUBROUTINE trd_dyn_iom_2d 
    181318 
    182319   !!====================================================================== 
Note: See TracChangeset for help on using the changeset viewer.