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

source: trunk/NEMO/OPA_SRC/TRA/tranpc.F90 @ 719

Last change on this file since 719 was 719, checked in by ctlod, 17 years ago

get back to the nemo_v2_3 version for trunk

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 10.9 KB
Line 
1MODULE tranpc
2   !!==============================================================================
3   !!                       ***  MODULE  tranpc  ***
4   !! Ocean active tracers:  non penetrative convection scheme
5   !!==============================================================================
6   !! History :  1.0  !  90-09  (G. Madec)  Original code
7   !!                 !  91-11  (G. Madec)
8   !!                 !  92-06  (M. Imbard)  periodic conditions on t and s
9   !!                 !  93-03  (M. Guyon)  symetrical conditions
10   !!                 !  96-01  (G. Madec)  statement function for e3
11   !!                                       suppression of common work arrays
12   !!            8.5  !  02-06  (G. Madec)  free form F90
13   !!----------------------------------------------------------------------
14
15   !!----------------------------------------------------------------------
16   !!   tra_npc      : apply the non penetrative convection scheme
17   !!   tra_npc_init : initialization and control of the scheme
18   !!----------------------------------------------------------------------
19   USE oce             ! ocean dynamics and active tracers
20   USE dom_oce         ! ocean space and time domain
21   USE trdmod          ! ocean active tracer trends
22   USE trdmod_oce      ! ocean variables trends
23   USE eosbn2          ! equation of state (eos routine)
24   USE lbclnk          ! lateral boundary conditions (or mpp link)
25   USE in_out_manager  ! I/O manager
26
27   IMPLICIT NONE
28   PRIVATE
29
30   PUBLIC   tra_npc    ! routine called by step.F90
31
32   !!* Namelist namnpc: non penetrative convection algorithm
33   INTEGER ::   nnpc1 =   1   ! nnpc1   non penetrative convective scheme frequency
34   INTEGER ::   nnpc2 =  15   ! nnpc2   non penetrative convective scheme print frequency
35
36   !! * Substitutions
37#  include "domzgr_substitute.h90"
38   !!----------------------------------------------------------------------
39   !!   OPA 9.0 , LOCEAN-IPSL (2005)
40   !! $Header$
41   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
42   !!----------------------------------------------------------------------
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 (now, after the swap)
52      !!      while conserving heat and salt contents.
53      !!
54      !! ** Method  :   The algorithm used converges in a maximium of jpk
55      !!      iterations. instabilities are treated when the vertical density
56      !!      gradient is less than 1.e-5.
57      !!      l_trdtra=T: the trend associated with this algorithm is saved.
58      !!
59      !! ** Action  : - (tn,sn) after the application od the npc scheme
60      !!              - save the associated trends (ttrd,strd) ('key_trdtra')
61      !!
62      !! References : Madec, et al., 1991, JPO, 21, 9, 1349-1371.
63      !!----------------------------------------------------------------------
64      USE oce, ONLY :    ztrdt => ua   ! use ua as 3D workspace   
65      USE oce, ONLY :    ztrds => va   ! use va as 3D workspace   
66      !!
67      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
68      !!
69      INTEGER  ::   ji, jj, jk   ! dummy loop indices
70      INTEGER  ::   inpcc        ! number of statically instable water column
71      INTEGER  ::   inpci        ! number of iteration for npc scheme
72      INTEGER  ::   jiter, jkdown, jkp        ! ???
73      INTEGER  ::   ikbot, ik, ikup, ikdown   ! ???
74      REAL(wp) ::   ze3tot, zta, zsa, zraua, ze3dwn
75      REAL(wp), DIMENSION(jpi,jpk)     ::   zwx, zwy, zwz   ! 2D arrays
76      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zrhop           ! 3D arrays
77      !!----------------------------------------------------------------------
78
79      IF( kt == nit000  )   CALL tra_npc_init   ! Initialisation
80
81      IF( MOD( kt, nnpc1 ) == 0 ) THEN
82
83         inpcc = 0
84         inpci = 0
85
86         CALL eos( tn, sn, rhd, zrhop )         ! Potential density
87
88
89         IF( l_trdtra )   THEN                  ! Save tn and sn trends
90            ztrdt(:,:,:) = tn(:,:,:) 
91            ztrds(:,:,:) = sn(:,:,:) 
92         ENDIF
93
94         !                                                ! ===============
95         DO jj = 1, jpj                                   !  Vertical slab
96            !                                             ! ===============
97            !  Static instability pointer
98            ! ----------------------------
99            DO jk = 1, jpkm1
100               DO ji = 1, jpi
101                  zwx(ji,jk) = ( zrhop(ji,jj,jk) - zrhop(ji,jj,jk+1) ) * tmask(ji,jj,jk+1)
102               END DO
103            END DO
104
105            ! 1.1 do not consider the boundary points
106
107            ! even if east-west cyclic b. c. do not considere ji=1 or jpi
108            DO jk = 1, jpkm1
109               zwx( 1 ,jk) = 0.e0
110               zwx(jpi,jk) = 0.e0
111            END DO
112            ! even if south-symmetric b. c. used, do not considere jj=1
113            IF( jj == 1 )   zwx(:,:) = 0.e0
114
115            DO jk = 1, jpkm1
116               DO ji = 1, jpi
117                  zwx(ji,jk) = 1.
118                  IF( zwx(ji,jk) < 1.e-5 ) zwx(ji,jk) = 0.e0
119               END DO
120            END DO
121
122            zwy(:,1) = 0.e0
123            DO ji = 1, jpi
124               DO jk = 1, jpkm1
125                  zwy(ji,1) = zwy(ji,1) + zwx(ji,jk)
126               END DO
127            END DO
128
129            zwz(1,1) = 0.e0
130            DO ji = 1, jpi
131               zwz(1,1) = zwz(1,1) + zwy(ji,1)
132            END DO
133
134            inpcc = inpcc + NINT( zwz(1,1) )
135
136
137            ! 2. Vertical mixing for each instable portion of the density profil
138            ! ------------------------------------------------------------------
139
140            IF( zwz(1,1) /= 0.e0 ) THEN         ! -->> the density profil is statically instable :
141               DO ji = 1, jpi
142                  IF( zwy(ji,1) /= 0.e0 ) THEN
143                     !
144                     ikbot = mbathy(ji,jj)      ! ikbot: ocean bottom level
145                     !
146                     DO jiter = 1, jpk          ! vertical iteration
147                        !
148                        ! search of ikup : the first static instability from the sea surface
149                        !
150                        ik = 0
151220                     CONTINUE
152                        ik = ik + 1
153                        IF( ik >= ikbot-1 ) GO TO 200
154                        zwx(ji,ik) = zrhop(ji,jj,ik) - zrhop(ji,jj,ik+1)
155                        IF( zwx(ji,ik) <= 0.e0 ) GO TO 220
156                        ikup = ik
157                        ! the density profil is instable below ikup
158                        ! ikdown : bottom of the instable portion of the density profil
159                        ! search of ikdown and vertical mixing from ikup to ikdown
160                        !
161                        ze3tot= fse3t(ji,jj,ikup)
162                        zta   = tn   (ji,jj,ikup)
163                        zsa   = sn   (ji,jj,ikup)
164                        zraua = zrhop(ji,jj,ikup)
165                        !
166                        DO jkdown = ikup+1, ikbot-1
167                           IF( zraua <= zrhop(ji,jj,jkdown) ) THEN
168                              ikdown = jkdown
169                              GO TO 240
170                           ENDIF
171                           ze3dwn =  fse3t(ji,jj,jkdown)
172                           ze3tot =  ze3tot + ze3dwn
173                           zta   = ( zta*(ze3tot-ze3dwn) + tn(ji,jj,jkdown)*ze3dwn )/ze3tot
174                           zsa   = ( zsa*(ze3tot-ze3dwn) + sn(ji,jj,jkdown)*ze3dwn )/ze3tot
175                           zraua = ( zraua*(ze3tot-ze3dwn) + zrhop(ji,jj,jkdown)*ze3dwn )/ze3tot
176                           inpci = inpci+1
177                        END DO
178                        ikdown = ikbot-1
179240                     CONTINUE
180                        !
181                        DO jkp = ikup, ikdown-1
182                           tn(ji,jj,jkp) = zta
183                           sn(ji,jj,jkp) = zsa
184                           zrhop(ji,jj,jkp) = zraua
185                        END DO
186                        IF (ikdown == ikbot-1 .AND. zraua >= zrhop(ji,jj,ikdown) ) THEN
187                           tn(ji,jj,ikdown) = zta
188                           sn(ji,jj,ikdown) = zsa
189                           zrhop(ji,jj,ikdown) = zraua
190                        ENDIF
191                     END DO
192                  ENDIF
193200               CONTINUE
194               END DO
195               ! <<-- no more static instability on slab jj
196            ENDIF
197            !                                             ! ===============
198         END DO                                           !   End of slab
199         !                                                ! ===============
200         !
201         IF( l_trdtra )   THEN         ! save the Non penetrative mixing trends for diagnostic
202            ztrdt(:,:,:) = tn(:,:,:) - ztrdt(:,:,:)
203            ztrds(:,:,:) = sn(:,:,:) - ztrds(:,:,:)
204            CALL trd_mod(ztrdt, ztrds, jptra_trd_npc, 'TRA', kt)
205         ENDIF
206     
207         ! Lateral boundary conditions on ( tn, sn )   ( Unchanged sign)
208         ! ------------------------------============
209         CALL lbc_lnk( tn, 'T', 1. )
210         CALL lbc_lnk( sn, 'T', 1. )
211     
212
213         !  2. non penetrative convective scheme statistics
214         !  -----------------------------------------------
215         IF( nnpc2 /= 0 .AND. MOD( kt, nnpc2 ) == 0 ) THEN
216            IF(lwp) WRITE(numout,*)' kt=',kt, ' number of statically instable',   &
217               &                   ' water column : ',inpcc, ' number of iteration : ',inpci
218         ENDIF
219         !
220      ENDIF
221      !
222   END SUBROUTINE tra_npc
223
224
225   SUBROUTINE tra_npc_init
226      !!----------------------------------------------------------------------
227      !!                  ***  ROUTINE tra_npc_init  ***
228      !!                   
229      !! ** Purpose :   initializations of the non-penetrative adjustment scheme
230      !!----------------------------------------------------------------------
231      NAMELIST/namnpc/ nnpc1, nnpc2
232      !
233      REWIND( numnam )           ! Namelist namzdf : vertical diffusion
234      READ  ( numnam, namnpc )
235      !
236      IF(lwp) THEN               ! Namelist print
237         WRITE(numout,*)
238         WRITE(numout,*) 'tra_npc_init : Non Penetrative Convection (npc) scheme'
239         WRITE(numout,*) '~~~~~~~~~~~~'
240         WRITE(numout,*) '       Namelist namnpc : set npc scheme parameters'
241         WRITE(numout,*) '          npc scheme frequency           nnpc1  = ', nnpc1
242         WRITE(numout,*) '          npc scheme print frequency     nnpc2  = ', nnpc2
243      ENDIF
244      !
245      IF ( nnpc1 == 0 ) THEN      ! Parameter controls
246          IF(lwp) WRITE(numout,cform_war)
247          IF(lwp) WRITE(numout,*) '             nnpc1 = ', nnpc1, ' is forced to 1'
248          nnpc1 = 1
249          nwarn = nwarn + 1
250      ENDIF
251      !
252   END SUBROUTINE tra_npc_init
253
254   !!======================================================================
255END MODULE tranpc
Note: See TracBrowser for help on using the repository browser.