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

Last change on this file since 527 was 503, checked in by opalod, 18 years ago

nemo_v1_update_064 : CT : general trends update including the addition of mean windows analysis possibility in the mixed layer

  • 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   NAMELIST/namnpc/ nnpc1, nnpc2
36
37   !! * Substitutions
38#  include "domzgr_substitute.h90"
39   !!----------------------------------------------------------------------
40   !!   OPA 9.0 , LOCEAN-IPSL (2005)
41   !! $Header$
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 (now, after the swap)
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  : - (tn,sn) 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( tn, sn, rhd, zrhop )         ! Potential density
88
89
90         IF( l_trdtra )   THEN                  ! Save tn and sn trends
91            ztrdt(:,:,:) = tn(:,:,:) 
92            ztrds(:,:,:) = sn(:,:,:) 
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   = tn   (ji,jj,ikup)
164                        zsa   = sn   (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) + tn(ji,jj,jkdown)*ze3dwn )/ze3tot
175                           zsa   = ( zsa*(ze3tot-ze3dwn) + sn(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                           tn(ji,jj,jkp) = zta
184                           sn(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                           tn(ji,jj,ikdown) = zta
189                           sn(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(:,:,:) = tn(:,:,:) - ztrdt(:,:,:)
204            ztrds(:,:,:) = sn(:,:,:) - ztrds(:,:,:)
205            CALL trd_mod(ztrdt, ztrds, jptra_trd_npc, 'TRA', kt)
206         ENDIF
207     
208         ! Lateral boundary conditions on ( tn, sn )   ( Unchanged sign)
209         ! ------------------------------============
210         CALL lbc_lnk( tn, 'T', 1. )
211         CALL lbc_lnk( sn, '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      !
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.