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/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/TRA – NEMO

source: NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/TRA/tranpc.F90 @ 13151

Last change on this file since 13151 was 13151, checked in by gm, 4 years ago

result from merge with qco r12983

  • Property svn:keywords set to Id
File size: 16.5 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
[12377]36#  include "do_loop_substitute.h90"
[13151]37#  include "domzgr_substitute.h90"
[3]38   !!----------------------------------------------------------------------
[9598]39   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
[4990]40   !! $Id$
[10068]41   !! Software governed by the CeCILL license (see ./LICENSE)
[3]42   !!----------------------------------------------------------------------
43CONTAINS
44
[12377]45   SUBROUTINE tra_npc( kt, Kmm, Krhs, pts, Kaa )
[3]46      !!----------------------------------------------------------------------
47      !!                  ***  ROUTINE tranpc  ***
48      !!
[4990]49      !! ** Purpose : Non-penetrative convective adjustment scheme. solve
[1111]50      !!      the static instability of the water column on after fields
[3]51      !!      while conserving heat and salt contents.
52      !!
[4990]53      !! ** Method  : updated algorithm able to deal with non-linear equation of state
54      !!              (i.e. static stability computed locally)
[3]55      !!
[6140]56      !! ** Action  : - tsa: after tracers with the application of the npc scheme
[4990]57      !!              - send the associated trends for on-line diagnostics (l_trdtra=T)
[3]58      !!
[4990]59      !! References :     Madec, et al., 1991, JPO, 21, 9, 1349-1371.
[503]60      !!----------------------------------------------------------------------
[12377]61      INTEGER,                                   INTENT(in   ) :: kt              ! ocean time-step index
62      INTEGER,                                   INTENT(in   ) :: Kmm, Krhs, Kaa  ! time level indices
63      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts             ! active tracers and RHS of tracer equation
[2715]64      !
[503]65      INTEGER  ::   ji, jj, jk   ! dummy loop indices
66      INTEGER  ::   inpcc        ! number of statically instable water column
[5386]67      INTEGER  ::   jiter, ikbot, ikp, ikup, ikdown, ilayer, ik_low   ! local integers
[4990]68      LOGICAL  ::   l_bottom_reached, l_column_treated
69      REAL(wp) ::   zta, zalfa, zsum_temp, zsum_alfa, zaw, zdz, zsum_z
[12489]70      REAL(wp) ::   zsa, zbeta, zsum_sali, zsum_beta, zbw, zrw, z1_rDt
[12377]71      REAL(wp), PARAMETER ::   zn2_zero = 1.e-14_wp             ! acceptance criteria for neutrality (N2==0)
72      REAL(wp), DIMENSION(        jpk     )   ::   zvn2         ! vertical profile of N2 at 1 given point...
73      REAL(wp), DIMENSION(        jpk,jpts)   ::   zvts, zvab   ! vertical profile of T & S , and  alpha & betaat 1 given point
[13151]74      REAL(wp), DIMENSION(jpi,jpj,jpk     )   ::   zn2          ! N^2
[12377]75      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts)   ::   zab          ! alpha and beta
76      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrdt, ztrds ! 3D workspace
[4990]77      !
[5386]78      LOGICAL, PARAMETER :: l_LB_debug = .FALSE. ! set to true if you want to follow what is
79      INTEGER :: ilc1, jlc1, klc1, nncpu         ! actually happening in a water column at point "ilc1, jlc1"
80      LOGICAL :: lp_monitor_point = .FALSE.      ! in CPU domain "nncpu"
[3]81      !!----------------------------------------------------------------------
[3294]82      !
[9019]83      IF( ln_timing )   CALL timing_start('tra_npc')
[3294]84      !
[1537]85      IF( MOD( kt, nn_npc ) == 0 ) THEN
[4990]86         !
87         IF( l_trdtra )   THEN                    !* Save initial after fields
[9019]88            ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) )
[13151]89            ztrdt(:,:,:) = pts(:,:,:,jp_tem,Kaa)
[12377]90            ztrds(:,:,:) = pts(:,:,:,jp_sal,Kaa)
[216]91         ENDIF
[9019]92         !
[4990]93         IF( l_LB_debug ) THEN
[5386]94            ! Location of 1 known convection site to follow what's happening in the water column
[13151]95            ilc1 = 45 ;  jlc1 = 3 ; !  ORCA2 4x4, Antarctic coast, more than 2 unstable portions in the  water column...
[5386]96            nncpu = 1  ;            ! the CPU domain contains the convection spot
[13151]97            klc1 =  mbkt(ilc1,jlc1) ! bottom of the ocean for debug point...
[4990]98         ENDIF
[9019]99         !
[12377]100         CALL eos_rab( pts(:,:,:,:,Kaa), zab, Kmm )         ! after alpha and beta (given on T-points)
101         CALL bn2    ( pts(:,:,:,:,Kaa), zab, zn2, Kmm )    ! after Brunt-Vaisala  (given on W-points)
[9019]102         !
[4990]103         inpcc = 0
[9019]104         !
[12377]105         DO_2D_00_00
106            !
107            IF( tmask(ji,jj,2) == 1 ) THEN      ! At least 2 ocean points
[13151]108               !                                     ! consider one ocean column
[12377]109               zvts(:,jp_tem) = pts(ji,jj,:,jp_tem,Kaa)      ! temperature
110               zvts(:,jp_sal) = pts(ji,jj,:,jp_sal,Kaa)      ! salinity
[4990]111               !
[13151]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
[12377]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:
[13151]120                  lp_monitor_point = (narea == nncpu).AND.lp_monitor_point
[12377]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 )
[6140]130                  !
[12377]131                  jiter = jiter + 1
[13151]132                  !
[12377]133                  IF( jiter >= 400 ) EXIT
[6140]134                  !
[12377]135                  l_bottom_reached = .FALSE.
[6140]136                  !
[12377]137                  DO WHILE ( .NOT. l_bottom_reached )
[4990]138                     !
[12377]139                     ikp = ikp + 1
[6140]140                     !
[12377]141                     !! Testing level ikp for instability
142                     !! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
143                     IF( zvn2(ikp) <  -zn2_zero ) THEN ! Instability found!
[6140]144                        !
[12377]145                        ilayer = ilayer + 1    ! yet another instable portion of the water column found....
[6140]146                        !
[13151]147                        IF( lp_monitor_point ) THEN
[12377]148                           WRITE(numout,*)
149                           IF( ilayer == 1 .AND. jiter == 1 ) THEN   ! first time a column is spoted with an instability
[5386]150                              WRITE(numout,*)
[12377]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                        !
[13151]162                        IF( jiter == 1 )   inpcc = inpcc + 1
[12377]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
[5386]176                              ENDIF
[12377]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
[13151]187
[12377]188                        DO jk = ikup, ikbot      ! Inside the instable (and overlying neutral) portion of the column
[6140]189                           !
[12377]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
[13151]196                           !
[12377]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
[13151]201
[12377]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')
[5386]204
[12377]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
[4990]210
[12377]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
[4990]219
[12377]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
[13151]227
228
[12377]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!
[13151]232
[12377]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!
[13151]235
[12377]236                        DO jk = ikup, ik_low              ! we must go 1 point deeper than ikdown!
[5386]237
[12377]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
[5386]243
[12377]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)
[5386]248
[12377]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) )  )
[5386]253
[12377]254                        END DO
[13151]255
[12377]256                        ikp = MIN(ikdown+1,ikbot)
[5386]257
[13151]258
[12377]259                     ENDIF  !IF( zvn2(ikp) < 0. )
[5386]260
261
[12377]262                     IF( ikp == ikbot ) l_bottom_reached = .TRUE.
263                     !
264                  END DO ! DO WHILE ( .NOT. l_bottom_reached )
[4990]265
[12377]266                  IF( ikp /= ikbot )   CALL ctl_stop( 'tra_npc :  PROBLEM #3')
[13151]267
[12377]268                  ! ******* At this stage ikp == ikbot ! *******
[13151]269
[12377]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,*)
[4990]280                     ENDIF
281                     !
[12377]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 )
[4990]289
[12377]290               !! Updating pts:
291               pts(ji,jj,:,jp_tem,Kaa) = zvts(:,jp_tem)
292               pts(ji,jj,:,jp_sal,Kaa) = zvts(:,jp_sal)
[4990]293
[12377]294               !! LB:  Potentially some other global variable beside theta and S can be treated here
295               !!      like BGC tracers.
[4990]296
[12377]297               IF( lp_monitor_point )   WRITE(numout,*)
[4990]298
[12377]299            ENDIF ! IF( tmask(ji,jj,3) == 1 ) THEN
[4990]300
[12377]301         END_2D
[4990]302         !
303         IF( l_trdtra ) THEN         ! send the Non penetrative mixing trends for diagnostic
[12489]304            z1_rDt = 1._wp / (2._wp * rn_Dt)
305            ztrdt(:,:,:) = ( pts(:,:,:,jp_tem,Kaa) - ztrdt(:,:,:) ) * z1_rDt
306            ztrds(:,:,:) = ( pts(:,:,:,jp_sal,Kaa) - ztrds(:,:,:) ) * z1_rDt
[12377]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 )
[9019]309            DEALLOCATE( ztrdt, ztrds )
[216]310         ENDIF
[4990]311         !
[12377]312         CALL lbc_lnk_multi( 'tranpc', pts(:,:,:,jp_tem,Kaa), 'T', 1., pts(:,:,:,jp_sal,Kaa), 'T', 1. )
[4990]313         !
[5386]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,*)
[3]317         ENDIF
[503]318         !
[4990]319      ENDIF   ! IF( MOD( kt, nn_npc ) == 0 ) THEN
[503]320      !
[9019]321      IF( ln_timing )   CALL timing_stop('tra_npc')
[3294]322      !
[3]323   END SUBROUTINE tra_npc
324
325   !!======================================================================
326END MODULE tranpc
Note: See TracBrowser for help on using the repository browser.