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.
tranpc.F90 in branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

source: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/TRA/tranpc.F90 @ 13320

Last change on this file since 13320 was 11101, checked in by frrh, 5 years ago

Merge changes from Met Office GMED ticket 450 to reduce unnecessary
text output from NEMO.
This output, which is typically not switchable, is rarely of interest
in normal (non-debugging) runs and simply redunantley consumes extra
file space.
Further, the presence of this text output has been shown to
significantly degrade performance of models which are run during
Met Office HPC RAID (disk) checks.
The new code introduces switches which are configurable via the
changes made in the associated Met Office MOCI ticket 399.

File size: 17.7 KB
RevLine 
[3]1MODULE tranpc
2   !!==============================================================================
3   !!                       ***  MODULE  tranpc  ***
[4990]4   !! Ocean active tracers:  non penetrative convective adjustment scheme
[3]5   !!==============================================================================
[1537]6   !! History :  1.0  ! 1990-09  (G. Madec)  Original code
7   !!                 ! 1996-01  (G. Madec)  statement function for e3
8   !!   NEMO     1.0  ! 2002-06  (G. Madec)  free form F90
9   !!            3.0  ! 2008-06  (G. Madec)  applied on ta, sa and called before tranxt in step.F90
[2528]10   !!            3.3  ! 2010-05  (C. Ethe, G. Madec)  merge TRC-TRA
[5386]11   !!            3.6  ! 2015-05  (L. Brodeau) new algorithm based on local Brunt-Vaisala freq.
[503]12   !!----------------------------------------------------------------------
[3]13
14   !!----------------------------------------------------------------------
[1537]15   !!   tra_npc : apply the non penetrative convection scheme
[3]16   !!----------------------------------------------------------------------
[4990]17   USE oce             ! ocean dynamics and active tracers
[3]18   USE dom_oce         ! ocean space and time domain
[4990]19   USE phycst          ! physical constants
[1537]20   USE zdf_oce         ! ocean vertical physics
[4990]21   USE trd_oce         ! ocean active tracer trends
[2715]22   USE trdtra          ! ocean active tracer trends
[4990]23   USE eosbn2          ! equation of state (eos routine)
24   !
[3]25   USE lbclnk          ! lateral boundary conditions (or mpp link)
[216]26   USE in_out_manager  ! I/O manager
[2715]27   USE lib_mpp         ! MPP library
[3294]28   USE wrk_nemo        ! Memory Allocation
29   USE timing          ! Timing
[3]30
31   IMPLICIT NONE
32   PRIVATE
33
[4990]34   PUBLIC   tra_npc    ! routine called by step.F90
[3]35
36   !! * Substitutions
37#  include "domzgr_substitute.h90"
[4990]38#  include "vectopt_loop_substitute.h90"
[3]39   !!----------------------------------------------------------------------
[4990]40   !! NEMO/OPA 3.6 , NEMO Consortium (2014)
41   !! $Id$
[2528]42   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
[3]43   !!----------------------------------------------------------------------
44CONTAINS
45
46   SUBROUTINE tra_npc( kt )
47      !!----------------------------------------------------------------------
48      !!                  ***  ROUTINE tranpc  ***
49      !!
[4990]50      !! ** Purpose : Non-penetrative convective adjustment scheme. solve
[1111]51      !!      the static instability of the water column on after fields
[3]52      !!      while conserving heat and salt contents.
53      !!
[4990]54      !! ** Method  : updated algorithm able to deal with non-linear equation of state
55      !!              (i.e. static stability computed locally)
[3]56      !!
[1111]57      !! ** Action  : - (ta,sa) after the application od the npc scheme
[4990]58      !!              - send the associated trends for on-line diagnostics (l_trdtra=T)
[3]59      !!
[4990]60      !! References :     Madec, et al., 1991, JPO, 21, 9, 1349-1371.
[503]61      !!----------------------------------------------------------------------
62      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
[2715]63      !
[503]64      INTEGER  ::   ji, jj, jk   ! dummy loop indices
65      INTEGER  ::   inpcc        ! number of statically instable water column
[5386]66      INTEGER  ::   jiter, ikbot, ikp, ikup, ikdown, ilayer, ik_low   ! local integers
[4990]67      LOGICAL  ::   l_bottom_reached, l_column_treated
68      REAL(wp) ::   zta, zalfa, zsum_temp, zsum_alfa, zaw, zdz, zsum_z
69      REAL(wp) ::   zsa, zbeta, zsum_sali, zsum_beta, zbw, zrw, z1_r2dt
[5386]70      REAL(wp), PARAMETER :: zn2_zero = 1.e-14_wp       ! acceptance criteria for neutrality (N2==0)
[4990]71      REAL(wp), POINTER, DIMENSION(:)       ::   zvn2   ! vertical profile of N2 at 1 given point...
72      REAL(wp), POINTER, DIMENSION(:,:)     ::   zvts   ! vertical profile of T and S at 1 given point...
73      REAL(wp), POINTER, DIMENSION(:,:)     ::   zvab   ! vertical profile of alpha and beta
74      REAL(wp), POINTER, DIMENSION(:,:,:)   ::   zn2    ! N^2
75      REAL(wp), POINTER, DIMENSION(:,:,:,:) ::   zab    ! alpha and beta
76      REAL(wp), POINTER, DIMENSION(:,:,:)   ::   ztrdt, ztrds   ! 3D workspace
77      !
[5386]78      LOGICAL, PARAMETER :: l_LB_debug = .FALSE. ! set to true if you want to follow what is
79      INTEGER :: ilc1, jlc1, klc1, nncpu         ! actually happening in a water column at point "ilc1, jlc1"
80      LOGICAL :: lp_monitor_point = .FALSE.      ! in CPU domain "nncpu"
[3]81      !!----------------------------------------------------------------------
[3294]82      !
83      IF( nn_timing == 1 )  CALL timing_start('tra_npc')
84      !
[1537]85      IF( MOD( kt, nn_npc ) == 0 ) THEN
[4990]86         !
87         CALL wrk_alloc( jpi, jpj, jpk, zn2 )    ! N2
88         CALL wrk_alloc( jpi, jpj, jpk, 2, zab ) ! Alpha and Beta
89         CALL wrk_alloc( jpk, 2, zvts, zvab )    ! 1D column vector at point ji,jj
90         CALL wrk_alloc( jpk, zvn2 )             ! 1D column vector at point ji,jj
[3]91
[4990]92         IF( l_trdtra )   THEN                    !* Save initial after fields
[3294]93            CALL wrk_alloc( jpi, jpj, jpk, ztrdt, ztrds )
[2715]94            ztrdt(:,:,:) = tsa(:,:,:,jp_tem) 
95            ztrds(:,:,:) = tsa(:,:,:,jp_sal)
[216]96         ENDIF
97
[4990]98         IF( l_LB_debug ) THEN
[5386]99            ! Location of 1 known convection site to follow what's happening in the water column
100            ilc1 = 45 ;  jlc1 = 3 ; !  ORCA2 4x4, Antarctic coast, more than 2 unstable portions in the  water column...           
101            nncpu = 1  ;            ! the CPU domain contains the convection spot
[4990]102            klc1 =  mbkt(ilc1,jlc1) ! bottom of the ocean for debug point...         
103         ENDIF
[5386]104         
105         CALL eos_rab( tsa, zab )         ! after alpha and beta (given on T-points)
106         CALL bn2    ( tsa, zab, zn2 )    ! after Brunt-Vaisala  (given on W-points)
[4990]107       
108         inpcc = 0
[3]109
[4990]110         DO jj = 2, jpjm1                 ! interior column only
111            DO ji = fs_2, fs_jpim1
112               !
113               IF( tmask(ji,jj,2) == 1 ) THEN      ! At least 2 ocean points
114                  !                                     ! consider one ocean column
115                  zvts(:,jp_tem) = tsa(ji,jj,:,jp_tem)      ! temperature
116                  zvts(:,jp_sal) = tsa(ji,jj,:,jp_sal)      ! salinity
[3]117
[4990]118                  zvab(:,jp_tem)  = zab(ji,jj,:,jp_tem)     ! Alpha
119                  zvab(:,jp_sal)  = zab(ji,jj,:,jp_sal)     ! Beta 
120                  zvn2(:)         = zn2(ji,jj,:)            ! N^2
121                 
122                  IF( l_LB_debug ) THEN                  !LB debug:
123                     lp_monitor_point = .FALSE.
124                     IF( ( ji == ilc1 ).AND.( jj == jlc1 ) ) lp_monitor_point = .TRUE.
125                     ! writing only if on CPU domain where conv region is:
[5386]126                     lp_monitor_point = (narea == nncpu).AND.lp_monitor_point                     
[4990]127                  ENDIF                                  !LB debug  end
[3]128
[4990]129                  ikbot = mbkt(ji,jj)   ! ikbot: ocean bottom T-level
[5386]130                  ikp = 1                  ! because N2 is irrelevant at the surface level (will start at ikp=2)
[4990]131                  ilayer = 0
132                  jiter  = 0
133                  l_column_treated = .FALSE.
134                 
135                  DO WHILE ( .NOT. l_column_treated )
136                     !
137                     jiter = jiter + 1
138                   
139                     IF( jiter >= 400 ) EXIT
140                   
141                     l_bottom_reached = .FALSE.
[3]142
[4990]143                     DO WHILE ( .NOT. l_bottom_reached )
[3]144
[5386]145                        ikp = ikp + 1
[4990]146                       
[5386]147                        !! Testing level ikp for instability
[4990]148                        !! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[5386]149                        IF( zvn2(ikp) <  -zn2_zero ) THEN ! Instability found!
[3]150
[5386]151                           ilayer = ilayer + 1    ! yet another instable portion of the water column found....
[3]152
[5386]153                           IF( lp_monitor_point ) THEN
154                              WRITE(numout,*)
155                              IF( ilayer == 1 .AND. jiter == 1 ) THEN   ! first time a column is spoted with an instability
156                                 WRITE(numout,*)
157                                 WRITE(numout,*) 'Time step = ',kt,' !!!'
158                              ENDIF
159                              WRITE(numout,*)  ' * Iteration #',jiter,': found instable portion #',ilayer,   &
160                                 &                                    ' in column! Starting at ikp =', ikp
161                              WRITE(numout,*)  ' *** N2 for point (i,j) = ',ji,' , ',jj
162                              DO jk = 1, klc1
163                                 WRITE(numout,*) jk, zvn2(jk)
164                              END DO
165                              WRITE(numout,*)
[11101]166                              IF(lflush) CALL flush(numout)
[5386]167                           ENDIF
168                           
[3]169
[5386]170                           IF( jiter == 1 )   inpcc = inpcc + 1 
[4990]171
[11101]172                           IF( lp_monitor_point )   THEN
173                              WRITE(numout, *) 'Negative N2 at ikp =',ikp,' for layer #', ilayer
174                              IF(lflush) CALL flush(numout)
175                           ENDIF
[4990]176
[5386]177                           !! ikup is the uppermost point where mixing will start:
178                           ikup = ikp - 1 ! ikup is always "at most at ikp-1", less if neutral levels overlying
179                           
180                           !! If the points above ikp-1 have N2 == 0 they must also be mixed:
181                           IF( ikp > 2 ) THEN
182                              DO jk = ikp-1, 2, -1
183                                 IF( ABS(zvn2(jk)) < zn2_zero ) THEN
184                                    ikup = ikup - 1  ! 1 more upper level has N2=0 and must be added for the mixing
185                                 ELSE
186                                    EXIT
187                                 ENDIF
188                              END DO
189                           ENDIF
190                           
191                           IF( ikup < 1 )   CALL ctl_stop( 'tra_npc :  PROBLEM #1')
192
[4990]193                           zsum_temp = 0._wp
194                           zsum_sali = 0._wp
195                           zsum_alfa = 0._wp
196                           zsum_beta = 0._wp
197                           zsum_z    = 0._wp
198                                                   
[5386]199                           DO jk = ikup, ikbot      ! Inside the instable (and overlying neutral) portion of the column
[4990]200                              !
201                              zdz       = fse3t(ji,jj,jk)
202                              zsum_temp = zsum_temp + zvts(jk,jp_tem)*zdz
203                              zsum_sali = zsum_sali + zvts(jk,jp_sal)*zdz
204                              zsum_alfa = zsum_alfa + zvab(jk,jp_tem)*zdz
205                              zsum_beta = zsum_beta + zvab(jk,jp_sal)*zdz
206                              zsum_z    = zsum_z    + zdz
[5386]207                              !                             
208                              IF( jk == ikbot ) EXIT ! avoid array-index overshoot in case ikbot = jpk, cause we're calling jk+1 next line
209                              !! EXIT when we have reached the last layer that is instable (N2<0) or neutral (N2=0):
210                              IF( zvn2(jk+1) > zn2_zero ) EXIT
[4990]211                           END DO
212                         
[5386]213                           ikdown = jk ! for the current unstable layer, ikdown is the deepest point with a negative or neutral N2
214                           IF( ikup == ikdown )   CALL ctl_stop( 'tra_npc :  PROBLEM #2')
215
216                           ! Mixing Temperature, salinity, alpha and beta from ikup to ikdown included:
[4990]217                           zta   = zsum_temp/zsum_z
218                           zsa   = zsum_sali/zsum_z
219                           zalfa = zsum_alfa/zsum_z
220                           zbeta = zsum_beta/zsum_z
221
[5386]222                           IF( lp_monitor_point ) THEN
223                              WRITE(numout,*) 'MIXED T, S, alfa and beta between ikup =',ikup,   &
224                                 &            ' and ikdown =',ikdown,', in layer #',ilayer
[4990]225                              WRITE(numout,*) '  => Mean temp. in that portion =', zta
226                              WRITE(numout,*) '  => Mean sali. in that portion =', zsa
[5386]227                              WRITE(numout,*) '  => Mean Alfa  in that portion =', zalfa
[4990]228                              WRITE(numout,*) '  => Mean Beta  in that portion =', zbeta
[11101]229                              IF(lflush) CALL flush(numout)
[4990]230                           ENDIF
231
232                           !! Homogenaizing the temperature, salinity, alpha and beta in this portion of the column
233                           DO jk = ikup, ikdown
234                              zvts(jk,jp_tem) = zta
235                              zvts(jk,jp_sal) = zsa
236                              zvab(jk,jp_tem) = zalfa
237                              zvab(jk,jp_sal) = zbeta
238                           END DO
[5386]239                           
240                           
241                           !! Updating N2 in the relvant portion of the water column
242                           !! Temperature, Salinity, Alpha and Beta have been homogenized in the unstable portion
243                           !! => Need to re-compute N2! will use Alpha and Beta!
244                           
245                           ikup   = MAX(2,ikup)         ! ikup can never be 1 !
246                           ik_low = MIN(ikdown+1,ikbot) ! we must go 1 point deeper than ikdown!
247                           
248                           DO jk = ikup, ik_low              ! we must go 1 point deeper than ikdown!
249
250                              !! Interpolating alfa and beta at W point:
251                              zrw =  (fsdepw(ji,jj,jk  ) - fsdept(ji,jj,jk)) &
252                                 & / (fsdept(ji,jj,jk-1) - fsdept(ji,jj,jk))
253                              zaw = zvab(jk,jp_tem) * (1._wp - zrw) + zvab(jk-1,jp_tem) * zrw
254                              zbw = zvab(jk,jp_sal) * (1._wp - zrw) + zvab(jk-1,jp_sal) * zrw
255
256                              !! N2 at W point, doing exactly as in eosbn2.F90:
257                              zvn2(jk) = grav*( zaw * ( zvts(jk-1,jp_tem) - zvts(jk,jp_tem) )     &
258                                 &            - zbw * ( zvts(jk-1,jp_sal) - zvts(jk,jp_sal) )  )  &
259                                 &       / fse3w(ji,jj,jk) * tmask(ji,jj,jk)
260
261                              !! OR, faster  => just considering the vertical gradient of density
262                              !! as only the signa maters...
263                              !zvn2(jk) = (  zaw * ( zvts(jk-1,jp_tem) - zvts(jk,jp_tem) )     &
264                              !     &      - zbw * ( zvts(jk-1,jp_sal) - zvts(jk,jp_sal) )  )
265
266                           END DO
267                       
268                           ikp = MIN(ikdown+1,ikbot)
269                           
270
271                        ENDIF  !IF( zvn2(ikp) < 0. )
272
273
274                        IF( ikp == ikbot ) l_bottom_reached = .TRUE.
[503]275                        !
[4990]276                     END DO ! DO WHILE ( .NOT. l_bottom_reached )
277
[5386]278                     IF( ikp /= ikbot )   CALL ctl_stop( 'tra_npc :  PROBLEM #3')
[4990]279                   
[5386]280                     ! ******* At this stage ikp == ikbot ! *******
[4990]281                   
[5386]282                     IF( ilayer > 0 ) THEN      !! least an unstable layer has been found
[503]283                        !
[5386]284                        IF( lp_monitor_point ) THEN
285                           WRITE(numout,*)
286                           WRITE(numout,*) 'After ',jiter,' iteration(s), we neutralized ',ilayer,' instable layer(s)'
287                           WRITE(numout,*) '   ==> N2 at i,j=',ji,',',jj,' now looks like this:'
[4990]288                           DO jk = 1, klc1
289                              WRITE(numout,*) jk, zvn2(jk)
290                           END DO
[5386]291                           WRITE(numout,*)
[11101]292                           IF(lflush) CALL flush(numout)
[3]293                        ENDIF
[5386]294                        !
295                        ikp    = 1     ! starting again at the surface for the next iteration
[4990]296                        ilayer = 0
297                     ENDIF
298                     !
[5386]299                     IF( ikp >= ikbot )   l_column_treated = .TRUE.
[4990]300                     !
301                  END DO ! DO WHILE ( .NOT. l_column_treated )
302
303                  !! Updating tsa:
304                  tsa(ji,jj,:,jp_tem) = zvts(:,jp_tem)
305                  tsa(ji,jj,:,jp_sal) = zvts(:,jp_sal)
306
[5386]307                  !! LB:  Potentially some other global variable beside theta and S can be treated here
308                  !!      like BGC tracers.
[4990]309
[5386]310                  IF( lp_monitor_point )   WRITE(numout,*)
[4990]311
312               ENDIF ! IF( tmask(ji,jj,3) == 1 ) THEN
313
314            END DO ! ji
315         END DO ! jj
316         !
317         IF( l_trdtra ) THEN         ! send the Non penetrative mixing trends for diagnostic
318            z1_r2dt = 1._wp / (2._wp * rdt)
319            ztrdt(:,:,:) = ( tsa(:,:,:,jp_tem) - ztrdt(:,:,:) ) * z1_r2dt
320            ztrds(:,:,:) = ( tsa(:,:,:,jp_sal) - ztrds(:,:,:) ) * z1_r2dt
321            CALL trd_tra( kt, 'TRA', jp_tem, jptra_npc, ztrdt )
322            CALL trd_tra( kt, 'TRA', jp_sal, jptra_npc, ztrds )
[3294]323            CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds )
[216]324         ENDIF
[4990]325         !
[2715]326         CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1. )   ;   CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1. )
[4990]327         !
[5386]328         IF( lwp .AND. l_LB_debug ) THEN
329            WRITE(numout,*) 'Exiting tra_npc , kt = ',kt,', => numb. of statically instable water-columns: ', inpcc
330            WRITE(numout,*)
[11101]331            IF(lflush) CALL flush(numout)
[3]332         ENDIF
[503]333         !
[4990]334         CALL wrk_dealloc(jpi, jpj, jpk, zn2 )
335         CALL wrk_dealloc(jpi, jpj, jpk, 2, zab )
336         CALL wrk_dealloc(jpk, zvn2 )
337         CALL wrk_dealloc(jpk, 2, zvts, zvab )
338         !
339      ENDIF   ! IF( MOD( kt, nn_npc ) == 0 ) THEN
[503]340      !
[3294]341      IF( nn_timing == 1 )  CALL timing_stop('tra_npc')
342      !
[3]343   END SUBROUTINE tra_npc
344
345   !!======================================================================
346END MODULE tranpc
Note: See TracBrowser for help on using the repository browser.