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

Last change on this file since 1291 was 1146, checked in by rblod, 16 years ago

Add svn Id (first try), see ticket #210

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 11.0 KB
RevLine 
[3]1MODULE tranpc
2   !!==============================================================================
3   !!                       ***  MODULE  tranpc  ***
4   !! Ocean active tracers:  non penetrative convection scheme
5   !!==============================================================================
[503]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
[1111]13   !!            9.0  !  08-06  (G. Madec)  applied on ta, sa and called before tranxt in step.F90
[503]14   !!----------------------------------------------------------------------
[3]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
[216]22   USE trdmod          ! ocean active tracer trends
23   USE trdmod_oce      ! ocean variables trends
[3]24   USE eosbn2          ! equation of state (eos routine)
25   USE lbclnk          ! lateral boundary conditions (or mpp link)
[216]26   USE in_out_manager  ! I/O manager
[3]27
28   IMPLICIT NONE
29   PRIVATE
30
[503]31   PUBLIC   tra_npc    ! routine called by step.F90
[3]32
[503]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
[3]36
37   !! * Substitutions
38#  include "domzgr_substitute.h90"
39   !!----------------------------------------------------------------------
[247]40   !!   OPA 9.0 , LOCEAN-IPSL (2005)
[1146]41   !! $Id$
[503]42   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
[3]43   !!----------------------------------------------------------------------
44
45CONTAINS
46
47   SUBROUTINE tra_npc( kt )
48      !!----------------------------------------------------------------------
49      !!                  ***  ROUTINE tranpc  ***
50      !!
51      !! ** Purpose :   Non penetrative convective adjustment scheme. solve
[1111]52      !!      the static instability of the water column on after fields
[3]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.
[503]58      !!      l_trdtra=T: the trend associated with this algorithm is saved.
[3]59      !!
[1111]60      !! ** Action  : - (ta,sa) after the application od the npc scheme
[3]61      !!              - save the associated trends (ttrd,strd) ('key_trdtra')
62      !!
[503]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
[3]69      !!
[503]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
[3]78      !!----------------------------------------------------------------------
[216]79
[503]80      IF( kt == nit000  )   CALL tra_npc_init   ! Initialisation
[3]81
82      IF( MOD( kt, nnpc1 ) == 0 ) THEN
83
84         inpcc = 0
85         inpci = 0
86
[1111]87         CALL eos( ta, sa, rhd, zrhop )         ! Potential density
[3]88
89
[1111]90         IF( l_trdtra )   THEN                  ! Save ta and sa trends
91            ztrdt(:,:,:) = ta(:,:,:) 
92            ztrds(:,:,:) = sa(:,:,:) 
[216]93         ENDIF
94
[3]95         !                                                ! ===============
96         DO jj = 1, jpj                                   !  Vertical slab
97            !                                             ! ===============
[503]98            !  Static instability pointer
99            ! ----------------------------
[3]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
[503]114            IF( jj == 1 )   zwx(:,:) = 0.e0
[3]115
116            DO jk = 1, jpkm1
117               DO ji = 1, jpi
118                  zwx(ji,jk) = 1.
[503]119                  IF( zwx(ji,jk) < 1.e-5 ) zwx(ji,jk) = 0.e0
[3]120               END DO
121            END DO
122
[503]123            zwy(:,1) = 0.e0
[3]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
[503]130            zwz(1,1) = 0.e0
[3]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
[503]141            IF( zwz(1,1) /= 0.e0 ) THEN         ! -->> the density profil is statically instable :
[3]142               DO ji = 1, jpi
[503]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                        !
[3]149                        ! search of ikup : the first static instability from the sea surface
[503]150                        !
[3]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)
[503]156                        IF( zwx(ji,ik) <= 0.e0 ) GO TO 220
[3]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
[503]161                        !
[3]162                        ze3tot= fse3t(ji,jj,ikup)
[1111]163                        zta   = ta   (ji,jj,ikup)
164                        zsa   = sa   (ji,jj,ikup)
[3]165                        zraua = zrhop(ji,jj,ikup)
[503]166                        !
[3]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
[1111]174                           zta   = ( zta*(ze3tot-ze3dwn) + ta(ji,jj,jkdown)*ze3dwn )/ze3tot
175                           zsa   = ( zsa*(ze3tot-ze3dwn) + sa(ji,jj,jkdown)*ze3dwn )/ze3tot
[3]176                           zraua = ( zraua*(ze3tot-ze3dwn) + zrhop(ji,jj,jkdown)*ze3dwn )/ze3tot
177                           inpci = inpci+1
178                        END DO
179                        ikdown = ikbot-1
180240                     CONTINUE
[503]181                        !
[3]182                        DO jkp = ikup, ikdown-1
[1111]183                           ta(ji,jj,jkp) = zta
184                           sa(ji,jj,jkp) = zsa
[3]185                           zrhop(ji,jj,jkp) = zraua
186                        END DO
187                        IF (ikdown == ikbot-1 .AND. zraua >= zrhop(ji,jj,ikdown) ) THEN
[1111]188                           ta(ji,jj,ikdown) = zta
189                           sa(ji,jj,ikdown) = zsa
[3]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         !                                                ! ===============
[503]201         !
202         IF( l_trdtra )   THEN         ! save the Non penetrative mixing trends for diagnostic
[1111]203            ztrdt(:,:,:) = ta(:,:,:) - ztrdt(:,:,:)
204            ztrds(:,:,:) = sa(:,:,:) - ztrds(:,:,:)
[503]205            CALL trd_mod(ztrdt, ztrds, jptra_trd_npc, 'TRA', kt)
[216]206         ENDIF
[3]207     
[1111]208         ! Lateral boundary conditions on ( ta, sa )   ( Unchanged sign)
[3]209         ! ------------------------------============
[1111]210         CALL lbc_lnk( ta, 'T', 1. )
211         CALL lbc_lnk( sa, 'T', 1. )
[3]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',   &
[503]218               &                   ' water column : ',inpcc, ' number of iteration : ',inpci
[3]219         ENDIF
[503]220         !
[3]221      ENDIF
[503]222      !
[3]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      !!----------------------------------------------------------------------
[541]232      NAMELIST/namnpc/ nnpc1, nnpc2
[503]233      !
234      REWIND( numnam )           ! Namelist namzdf : vertical diffusion
[3]235      READ  ( numnam, namnpc )
[503]236      !
237      IF(lwp) THEN               ! Namelist print
[3]238         WRITE(numout,*)
239         WRITE(numout,*) 'tra_npc_init : Non Penetrative Convection (npc) scheme'
240         WRITE(numout,*) '~~~~~~~~~~~~'
[503]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
[3]244      ENDIF
[503]245      !
246      IF ( nnpc1 == 0 ) THEN      ! Parameter controls
[3]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
[503]252      !
[3]253   END SUBROUTINE tra_npc_init
254
255   !!======================================================================
256END MODULE tranpc
Note: See TracBrowser for help on using the repository browser.