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

Last change on this file since 3 was 3, checked in by opalod, 21 years ago

Initial revision

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