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

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

2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps : interface changes to tra and trc routines for design compliance and consistency. Fully SETTE tested (non-AGRIF, only)

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