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 @ 1111

Last change on this file since 1111 was 1111, checked in by ctlod, 16 years ago

trunk: remove a small bug related to the non-penetrative convective adjustment scheme, see ticket: #203

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