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/2014/dev_r4650_UKMO3_masked_damping/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

source: branches/2014/dev_r4650_UKMO3_masked_damping/NEMOGCM/NEMO/OPA_SRC/TRA/tranpc.F90 @ 5086

Last change on this file since 5086 was 5086, checked in by timgraham, 9 years ago

Merged head of trunk into branch in preparation for putting code back onto the trunk
In working copy ran the command:
svn merge svn+sshtimgraham@…/ipsl/forge/projets/nemo/svn/trunk

Also recompiled NEMO_book.pdf with merged input files

  • Property svn:keywords set to Id
File size: 16.9 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.7  ! 2014-06  (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 wrk_nemo        ! Memory Allocation
29   USE timing          ! Timing
30
31   IMPLICIT NONE
32   PRIVATE
33
34   PUBLIC   tra_npc    ! routine called by step.F90
35
36   !! * Substitutions
37#  include "domzgr_substitute.h90"
38#  include "vectopt_loop_substitute.h90"
39   !!----------------------------------------------------------------------
40   !! NEMO/OPA 3.6 , NEMO Consortium (2014)
41   !! $Id$
42   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
43   !!----------------------------------------------------------------------
44CONTAINS
45
46   SUBROUTINE tra_npc( kt )
47      !!----------------------------------------------------------------------
48      !!                  ***  ROUTINE tranpc  ***
49      !!
50      !! ** Purpose : Non-penetrative convective adjustment scheme. solve
51      !!      the static instability of the water column on after fields
52      !!      while conserving heat and salt contents.
53      !!
54      !! ** Method  : updated algorithm able to deal with non-linear equation of state
55      !!              (i.e. static stability computed locally)
56      !!
57      !! ** Action  : - (ta,sa) after the application od the npc scheme
58      !!              - send the associated trends for on-line diagnostics (l_trdtra=T)
59      !!
60      !! References :     Madec, et al., 1991, JPO, 21, 9, 1349-1371.
61      !!----------------------------------------------------------------------
62      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
63      !
64      INTEGER  ::   ji, jj, jk   ! dummy loop indices
65      INTEGER  ::   inpcc        ! number of statically instable water column
66      INTEGER  ::   jiter, ikbot, ik, ikup, ikdown, ilayer, ikm   ! 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), POINTER, DIMENSION(:)       ::   zvn2   ! vertical profile of N2 at 1 given point...
71      REAL(wp), POINTER, DIMENSION(:,:)     ::   zvts   ! vertical profile of T and S at 1 given point...
72      REAL(wp), POINTER, DIMENSION(:,:)     ::   zvab   ! vertical profile of alpha and beta
73      REAL(wp), POINTER, DIMENSION(:,:,:)   ::   zn2    ! N^2
74      REAL(wp), POINTER, DIMENSION(:,:,:,:) ::   zab    ! alpha and beta
75      REAL(wp), POINTER, DIMENSION(:,:,:)   ::   ztrdt, ztrds   ! 3D workspace
76      !
77      !!LB debug:
78      LOGICAL, PARAMETER :: l_LB_debug = .FALSE.
79      INTEGER :: ilc1, jlc1, klc1, nncpu
80      LOGICAL :: lp_monitor_point = .FALSE.
81      !!LB debug.
82      !!----------------------------------------------------------------------
83      !
84      IF( nn_timing == 1 )  CALL timing_start('tra_npc')
85      !
86      IF( MOD( kt, nn_npc ) == 0 ) THEN
87         !
88         CALL wrk_alloc( jpi, jpj, jpk, zn2 )    ! N2
89         CALL wrk_alloc( jpi, jpj, jpk, 2, zab ) ! Alpha and Beta
90         CALL wrk_alloc( jpk, 2, zvts, zvab )    ! 1D column vector at point ji,jj
91         CALL wrk_alloc( jpk, zvn2 )             ! 1D column vector at point ji,jj
92
93         IF( l_trdtra )   THEN                    !* Save initial after fields
94            CALL wrk_alloc( jpi, jpj, jpk, ztrdt, ztrds )
95            ztrdt(:,:,:) = tsa(:,:,:,jp_tem) 
96            ztrds(:,:,:) = tsa(:,:,:,jp_sal)
97         ENDIF
98
99         !LB debug:
100         IF( lwp .AND. l_LB_debug ) THEN
101            WRITE(numout,*)
102            WRITE(numout,*) 'LOLO: entering tra_npc, kt, narea =', kt, narea
103         ENDIF
104         !LBdebug: Monitoring of 1 column subject to convection...
105         IF( l_LB_debug ) THEN
106            ! Location of 1 known convection spot to follow what's happening in the water column
107            ilc1 = 54 ;  jlc1 = 15 ; ! Labrador ORCA1 4x4 cpus:
108            nncpu = 15  ; ! the CPU domain contains the convection spot
109            !ilc1 = 14 ;  jlc1 = 13 ; ! Labrador ORCA1 8x8 cpus:
110            !nncpu = 54  ; ! the CPU domain contains the convection spot
111            klc1 =  mbkt(ilc1,jlc1) ! bottom of the ocean for debug point...         
112         ENDIF
113         !LBdebug.
114
115         CALL eos_rab( tsa, zab )         ! after alpha and beta
116         CALL bn2    ( tsa, zab, zn2 )    ! after Brunt-Vaisala
117       
118         inpcc = 0
119
120         DO jj = 2, jpjm1                 ! interior column only
121            DO ji = fs_2, fs_jpim1
122               !
123               IF( tmask(ji,jj,2) == 1 ) THEN      ! At least 2 ocean points
124                  !                                     ! consider one ocean column
125                  zvts(:,jp_tem) = tsa(ji,jj,:,jp_tem)      ! temperature
126                  zvts(:,jp_sal) = tsa(ji,jj,:,jp_sal)      ! salinity
127
128                  zvab(:,jp_tem)  = zab(ji,jj,:,jp_tem)     ! Alpha
129                  zvab(:,jp_sal)  = zab(ji,jj,:,jp_sal)     ! Beta 
130                  zvn2(:)         = zn2(ji,jj,:)            ! N^2
131                 
132                  IF( l_LB_debug ) THEN                  !LB debug:
133                     lp_monitor_point = .FALSE.
134                     IF( ( ji == ilc1 ).AND.( jj == jlc1 ) ) lp_monitor_point = .TRUE.
135                     ! writing only if on CPU domain where conv region is:
136                     lp_monitor_point = (narea == nncpu).AND.lp_monitor_point 
137                     
138                     IF(lp_monitor_point) THEN
139                        WRITE(numout,*) '' ;WRITE(numout,*) '' ;
140                       WRITE(numout,'("Time step = ",i6.6," !!!")')  kt
141                        WRITE(numout,'(" *** BEFORE anything, N^2 for point ",i3,",",i3,":" )') ji,jj
142                        DO jk = 1, klc1
143                           WRITE(numout,*) jk, zvn2(jk)
144                        END DO
145                        WRITE(numout,*) ' '
146                     ENDIF
147                  ENDIF                                  !LB debug  end
148
149                  ikbot = mbkt(ji,jj)   ! ikbot: ocean bottom T-level
150                  ik = 1                ! because N2 is irrelevant at the surface level (will start at ik=2)
151                  ilayer = 0
152                  jiter  = 0
153                  l_column_treated = .FALSE.
154                 
155                  DO WHILE ( .NOT. l_column_treated )
156                     !
157                     jiter = jiter + 1
158                   
159                     IF( jiter >= 400 ) EXIT
160                   
161                     l_bottom_reached = .FALSE.
162
163                     DO WHILE ( .NOT. l_bottom_reached )
164
165                        ik = ik + 1
166                       
167                        !! Checking level ik for instability
168                        !! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
169
170                        IF( zvn2(ik) < 0. ) THEN ! Instability found!
171
172                           ikm  = ik              ! first level whith negative N2
173                           ilayer = ilayer + 1    ! yet another layer found....
174                           IF(jiter == 1) inpcc = inpcc + 1
175
176                           IF(l_LB_debug .AND. lp_monitor_point) &
177                              & WRITE(numout,*) 'Negative N2 at ik =', ikm, ' layer nb.', ilayer, &
178                              & ' inpcc =', inpcc
179
180                           !! Case we mix with upper regions where N2==0:
181                           !! All the points above ikup where N2 == 0 must also be mixed => we go
182                           !! upward to find a new ikup, where the layer doesn't have N2==0
183                           ikup = ikm
184                           DO jk = ikm, 2, -1
185                              ikup = ikup - 1
186                              IF( (zvn2(jk-1) > 0.).OR.(ikup == 1) ) EXIT
187                           END DO
188                         
189                           ! adjusting ikup if the upper part of the unstable column was neutral (N2=0)
190                           IF((zvn2(ikup+1) == 0.).AND.(ikup /= 1)) ikup = ikup+1 ;
191
192                         
193                           IF( lp_monitor_point )   WRITE(numout,*) ' => ikup is =', ikup, ' layer nb.', ilayer
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+1      ! Inside the instable (and overlying neutral) portion of the column
202                              !
203                              IF(l_LB_debug .AND. lp_monitor_point) WRITE(numout,*) '     -> summing for jk =', jk
204                              !
205                              zdz       = fse3t(ji,jj,jk)
206                              zsum_temp = zsum_temp + zvts(jk,jp_tem)*zdz
207                              zsum_sali = zsum_sali + zvts(jk,jp_sal)*zdz
208                              zsum_alfa = zsum_alfa + zvab(jk,jp_tem)*zdz
209                              zsum_beta = zsum_beta + zvab(jk,jp_sal)*zdz
210                              zsum_z    = zsum_z    + zdz
211                              !
212                              !! EXIT if we found the bottom of the unstable portion of the water column   
213                              IF( (zvn2(jk+1) > 0.).OR.(jk == ikbot ).OR.((jk==ikm).AND.(zvn2(jk+1) == 0.)) )   EXIT
214                           END DO
215                         
216                           !ik     = jk !LB remove?
217                           ikdown = jk ! for the current unstable layer, ikdown is the deepest point with a negative N2
218                         
219                           IF(l_LB_debug .AND. lp_monitor_point) &
220                              &    WRITE(numout,*) '  => ikdown =', ikdown, '  layer nb.', ilayer
221                         
222                           ! Mixing Temperature and salinity between ikup and ikdown:
223                           zta   = zsum_temp/zsum_z
224                           zsa   = zsum_sali/zsum_z
225                           zalfa = zsum_alfa/zsum_z
226                           zbeta = zsum_beta/zsum_z
227
228                           IF(l_LB_debug .AND. lp_monitor_point) THEN
229                              WRITE(numout,*) '  => Mean temp. in that portion =', zta
230                              WRITE(numout,*) '  => Mean sali. in that portion =', zsa
231                              WRITE(numout,*) '  => Mean Alpha in that portion =', zalfa
232                              WRITE(numout,*) '  => Mean Beta  in that portion =', zbeta
233                           ENDIF
234
235                           !! Homogenaizing the temperature, salinity, alpha and beta in this portion of the column
236                           DO jk = ikup, ikdown
237                              zvts(jk,jp_tem) = zta
238                              zvts(jk,jp_sal) = zsa
239                              zvab(jk,jp_tem) = zalfa
240                              zvab(jk,jp_sal) = zbeta
241                           END DO
242                           !
243                           !! Before updating N2, it is possible that another unstable
244                           !! layer exists underneath the one we just homogeneized!
245                           ik = ikdown
246                           !
247                        ENDIF  ! IF( zvn2(ik+1) < 0. ) THEN
248                        !
249                        IF( ik == ikbot ) l_bottom_reached = .TRUE.
250                        !
251                     END DO ! DO WHILE ( .NOT. l_bottom_reached )
252
253                     IF( ik /= ikbot )   STOP 'ERROR: tranpc.F90 => PROBLEM #1'
254                   
255                     ! ******* At this stage ik == ikbot ! *******
256                   
257                     IF( ilayer > 0 ) THEN
258                        !! least an unstable layer has been found
259                        !! Temperature, Salinity, Alpha and Beta have been homogenized in the unstable portion
260                        !! => Need to re-compute N2! will use Alpha and Beta!
261                        !
262                        DO jk = ikup+1, ikdown+1   ! we must go 1 point deeper than ikdown!     
263                           !! Doing exactly as in eosbn2.F90:
264                           !! * Except that we only are interested in the sign of N2 !!!
265                           !!   => just considering the vertical gradient of density
266                           zrw =   (fsdepw(ji,jj,jk  ) - fsdept(ji,jj,jk)) &
267                               & / (fsdept(ji,jj,jk-1) - fsdept(ji,jj,jk))
268                           zaw = zvab(jk,jp_tem) * (1._wp - zrw) + zvab(jk-1,jp_tem) * zrw
269                           zbw = zvab(jk,jp_sal) * (1._wp - zrw) + zvab(jk-1,jp_sal) * zrw
270                         
271                           !zvn2(jk) = grav*( zaw * ( zvts(jk-1,jp_tem) - zvts(jk,jp_tem) )     &
272                           !     &           - zbw * ( zvts(jk-1,jp_sal) - zvts(jk,jp_sal) )  )  &
273                           !     &       / fse3w(ji,jj,jk) * tmask(ji,jj,jk)
274                           zvn2(jk) = (  zaw * ( zvts(jk-1,jp_tem) - zvts(jk,jp_tem) )     &
275                                &      - zbw * ( zvts(jk-1,jp_sal) - zvts(jk,jp_sal) )  ) 
276                        END DO
277
278                        IF(l_LB_debug .AND. lp_monitor_point) THEN
279                           WRITE(numout, '(" *** After iteration #",i3.3,", N^2 for point ",i3,",",i3,":" )') &
280                              & jiter, ji,jj
281                           DO jk = 1, klc1
282                              WRITE(numout,*) jk, zvn2(jk)
283                           END DO
284                           WRITE(numout,*) ' '
285                        ENDIF
286
287                        ik     = 1  ! starting again at the surface for the next iteration
288                        ilayer = 0
289                     ENDIF
290                     !
291                     IF( ik >= ikbot ) THEN
292                        IF(l_LB_debug .AND. lp_monitor_point) WRITE(numout,*) '    --- exiting jiter loop ---'
293                        l_column_treated = .TRUE.
294                     ENDIF
295                     !
296                  END DO ! DO WHILE ( .NOT. l_column_treated )
297
298                  !! Updating tsa:
299                  tsa(ji,jj,:,jp_tem) = zvts(:,jp_tem)
300                  tsa(ji,jj,:,jp_sal) = zvts(:,jp_sal)
301
302                  !! lolo:  Should we update something else????
303                  !! => like alpha and beta?
304
305                  IF(l_LB_debug .AND. lp_monitor_point) WRITE(numout,*) ''
306
307               ENDIF ! IF( tmask(ji,jj,3) == 1 ) THEN
308
309            END DO ! ji
310         END DO ! jj
311         !
312         IF( l_trdtra ) THEN         ! send the Non penetrative mixing trends for diagnostic
313            z1_r2dt = 1._wp / (2._wp * rdt)
314            ztrdt(:,:,:) = ( tsa(:,:,:,jp_tem) - ztrdt(:,:,:) ) * z1_r2dt
315            ztrds(:,:,:) = ( tsa(:,:,:,jp_sal) - ztrds(:,:,:) ) * z1_r2dt
316            CALL trd_tra( kt, 'TRA', jp_tem, jptra_npc, ztrdt )
317            CALL trd_tra( kt, 'TRA', jp_sal, jptra_npc, ztrds )
318            CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds )
319         ENDIF
320         !
321         CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1. )   ;   CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1. )
322         !
323         IF(lwp) THEN
324            WRITE(numout,*) 'LOLO: exiting tra_npc, kt =', kt
325            WRITE(numout,*)' => number of statically instable water column : ',inpcc
326            WRITE(numout,*) '' ; WRITE(numout,*) ''
327         ENDIF
328         !
329         CALL wrk_dealloc(jpi, jpj, jpk, zn2 )
330         CALL wrk_dealloc(jpi, jpj, jpk, 2, zab )
331         CALL wrk_dealloc(jpk, zvn2 )
332         CALL wrk_dealloc(jpk, 2, zvts, zvab )
333         !
334      ENDIF   ! IF( MOD( kt, nn_npc ) == 0 ) THEN
335      !
336      IF( nn_timing == 1 )  CALL timing_stop('tra_npc')
337      !
338   END SUBROUTINE tra_npc
339
340   !!======================================================================
341END MODULE tranpc
Note: See TracBrowser for help on using the repository browser.