source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/tranpc.F90 @ 10954

Last change on this file since 10954 was 10954, checked in by acc, 18 months ago

2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps : Convert TRA modules and all knock on effects of these conversions. SETTE tested

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