source: NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/OCE/TRA/tranpc.F90 @ 13517

Last change on this file since 13517 was 13517, checked in by hadcv, 5 months ago

Tiling for modules after tra_ldf

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