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 NEMO/trunk/src/OCE/TRA – NEMO

source: NEMO/trunk/src/OCE/TRA/tranpc.F90 @ 12165

Last change on this file since 12165 was 10425, checked in by smasson, 5 years ago

trunk: merge back dev_r10164_HPC09_ESIWACE_PREP_MERGE@10424 into the trunk

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