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 branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/TRA/tranpc.F90 @ 4431

Last change on this file since 4431 was 3211, checked in by spickles2, 12 years ago

Stephen Pickles, 11 Dec 2011

Commit to bring the rest of the DCSE NEMO development branch
in line with the latest development version. This includes
array index re-ordering of all OPA_SRC/.

  • Property svn:keywords set to Id
File size: 10.3 KB
Line 
1MODULE tranpc
2   !!==============================================================================
3   !!                       ***  MODULE  tranpc  ***
4   !! Ocean active tracers:  non penetrative convection scheme
5   !!==============================================================================
6   !! History :  1.0  ! 1990-09  (G. Madec)  Original code
7   !!                 ! 1996-01  (G. Madec)  statement function for e3
8   !!   NEMO     1.0  ! 2002-06  (G. Madec)  free form F90
9   !!            3.0  ! 2008-06  (G. Madec)  applied on ta, sa and called before tranxt in step.F90
10   !!            3.3  ! 2010-05  (C. Ethe, G. Madec)  merge TRC-TRA
11   !!----------------------------------------------------------------------
12
13   !!----------------------------------------------------------------------
14   !!   tra_npc : apply the non penetrative convection scheme
15   !!----------------------------------------------------------------------
16   USE oce             ! ocean dynamics and active tracers
17   USE dom_oce         ! ocean space and time domain
18   USE zdf_oce         ! ocean vertical physics
19   USE trdmod_oce      ! ocean active tracer trends
20   USE trdtra          ! ocean active tracer trends
21   USE eosbn2          ! equation of state (eos routine)
22   USE lbclnk          ! lateral boundary conditions (or mpp link)
23   USE in_out_manager  ! I/O manager
24   USE lib_mpp         ! MPP library
25
26   IMPLICIT NONE
27   PRIVATE
28
29   PUBLIC   tra_npc       ! routine called by step.F90
30
31   !! * Control permutation of array indices
32#  include "oce_ftrans.h90"
33#  include "dom_oce_ftrans.h90"
34#  include "zdf_oce_ftrans.h90"
35
36   !! * Substitutions
37#  include "domzgr_substitute.h90"
38   !!----------------------------------------------------------------------
39   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
40   !! $Id$
41   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
42   !!----------------------------------------------------------------------
43CONTAINS
44
45   SUBROUTINE tra_npc( kt )
46      !!----------------------------------------------------------------------
47      !!                  ***  ROUTINE tranpc  ***
48      !!
49      !! ** Purpose :   Non penetrative convective adjustment scheme. solve
50      !!      the static instability of the water column on after fields
51      !!      while conserving heat and salt contents.
52      !!
53      !! ** Method  :   The algorithm used converges in a maximium of jpk
54      !!      iterations. instabilities are treated when the vertical density
55      !!      gradient is less than 1.e-5.
56      !!      l_trdtra=T: the trend associated with this algorithm is saved.
57      !!
58      !! ** Action  : - (ta,sa) after the application od the npc scheme
59      !!              - save the associated trends (ttrd,strd) ('key_trdtra')
60      !!
61      !! References : Madec, et al., 1991, JPO, 21, 9, 1349-1371.
62      !!----------------------------------------------------------------------
63      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released, wrk_in_use_xz, wrk_not_released_xz
64      USE wrk_nemo, ONLY:   ztrdt => wrk_3d_1 , ztrds => wrk_3d_2 , zrhop => wrk_3d_3
65      USE wrk_nemo, ONLY:   zwx   => wrk_xz_1 , zwy   => wrk_xz_2 , zwz   => wrk_xz_3
66
67      !! DCSE_NEMO: need additional directives for renamed module variables
68!FTRANS ztrdt ztrds zrhop :I :I :z
69
70      !
71      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
72      !
73      INTEGER  ::   ji, jj, jk   ! dummy loop indices
74      INTEGER  ::   inpcc        ! number of statically instable water column
75      INTEGER  ::   inpci        ! number of iteration for npc scheme
76      INTEGER  ::   jiter, jkdown, jkp        ! ???
77      INTEGER  ::   ikbot, ik, ikup, ikdown   ! ???
78      REAL(wp) ::   ze3tot, zta, zsa, zraua, ze3dwn
79      !!----------------------------------------------------------------------
80
81      ! Strictly 1 and 2 3D workspaces only needed if(l_trdtra) but it doesn't
82      ! cost us anything and makes code simpler.
83      IF( wrk_in_use(3, 1,2,3) .OR. wrk_in_use_xz(1,2,3) ) THEN
84         CALL ctl_stop('tra_npc: requested workspace arrays unavailable')   ;   RETURN
85      ENDIF
86
87      IF( MOD( kt, nn_npc ) == 0 ) THEN
88
89         inpcc = 0
90         inpci = 0
91
92         CALL eos( tsa, rhd, zrhop )         ! Potential density
93
94         IF( l_trdtra )   THEN                    !* Save ta and sa trends
95            ztrdt(:,:,:) = tsa(:,:,:,jp_tem) 
96            ztrds(:,:,:) = tsa(:,:,:,jp_sal)
97         ENDIF
98
99         !                                                ! ===============
100         DO jj = 1, jpj                                   !  Vertical slab
101            !                                             ! ===============
102            !  Static instability pointer
103            ! ----------------------------
104#if defined key_z_first
105            DO ji = 1, jpi
106               DO jk = 1, jpkm1
107#else
108            DO jk = 1, jpkm1
109               DO ji = 1, jpi
110#endif
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.e0
129               END DO
130            END DO
131
132            zwy(:,1) = 0.e0
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.e0
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
147            ! 2. Vertical mixing for each instable portion of the density profil
148            ! ------------------------------------------------------------------
149
150            IF( zwz(1,1) /= 0.e0 ) THEN         ! -->> the density profil is statically instable :
151               DO ji = 1, jpi
152                  IF( zwy(ji,1) /= 0.e0 ) THEN
153                     !
154                     ikbot = mbkt(ji,jj)        ! ikbot: ocean bottom T-level
155                     !
156                     DO jiter = 1, jpk          ! vertical iteration
157                        !
158                        ! search of ikup : the first static instability from the sea surface
159                        !
160                        ik = 0
161220                     CONTINUE
162                        ik = ik + 1
163                        IF( ik >= ikbot ) GO TO 200
164                        zwx(ji,ik) = zrhop(ji,jj,ik) - zrhop(ji,jj,ik+1)
165                        IF( zwx(ji,ik) <= 0.e0 ) GO TO 220
166                        ikup = ik
167                        ! the density profil is instable below ikup
168                        ! ikdown : bottom of the instable portion of the density profil
169                        ! search of ikdown and vertical mixing from ikup to ikdown
170                        !
171                        ze3tot= fse3t(ji,jj,ikup)
172                        zta   = tsa  (ji,jj,ikup,jp_tem)
173                        zsa   = tsa  (ji,jj,ikup,jp_sal)
174                        zraua = zrhop(ji,jj,ikup)
175                        !
176                        DO jkdown = ikup+1, ikbot-1
177                           IF( zraua <= zrhop(ji,jj,jkdown) ) THEN
178                              ikdown = jkdown
179                              GO TO 240
180                           ENDIF
181                           ze3dwn =  fse3t(ji,jj,jkdown)
182                           ze3tot =  ze3tot + ze3dwn
183                           zta   = ( zta*(ze3tot-ze3dwn) + tsa(ji,jj,jkdown,jp_tem)*ze3dwn )/ze3tot
184                           zsa   = ( zsa*(ze3tot-ze3dwn) + tsa(ji,jj,jkdown,jp_sal)*ze3dwn )/ze3tot
185                           zraua = ( zraua*(ze3tot-ze3dwn) + zrhop(ji,jj,jkdown)*ze3dwn )/ze3tot
186                           inpci = inpci+1
187                        END DO
188                        ikdown = ikbot-1
189240                     CONTINUE
190                        !
191                        DO jkp = ikup, ikdown-1
192                           tsa  (ji,jj,jkp,jp_tem) = zta
193                           tsa  (ji,jj,jkp,jp_sal) = zsa
194                           zrhop(ji,jj,jkp       ) = zraua
195                        END DO
196                        IF (ikdown == ikbot-1 .AND. zraua >= zrhop(ji,jj,ikdown) ) THEN
197                           tsa  (ji,jj,jkp,jp_tem) = zta
198                           tsa  (ji,jj,jkp,jp_sal) = zsa
199                           zrhop(ji,jj,ikdown    ) = zraua
200                        ENDIF
201                     END DO
202                  ENDIF
203200               CONTINUE
204               END DO
205               ! <<-- no more static instability on slab jj
206            ENDIF
207            !                                             ! ===============
208         END DO                                           !   End of slab
209         !                                                ! ===============
210         !
211         IF( l_trdtra )   THEN         ! save the Non penetrative mixing trends for diagnostic
212            ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)
213            ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:)
214            CALL trd_tra( kt, 'TRA', jp_tem, jptra_trd_npc, ztrdt )
215            CALL trd_tra( kt, 'TRA', jp_sal, jptra_trd_npc, ztrds )
216         ENDIF
217     
218         ! Lateral boundary conditions on ( ta, sa )   ( Unchanged sign)
219         ! ------------------------------============
220         CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1. )   ;   CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1. )
221     
222
223         !  2. non penetrative convective scheme statistics
224         !  -----------------------------------------------
225         IF( nn_npcp /= 0 .AND. MOD( kt, nn_npcp ) == 0 ) THEN
226            IF(lwp) WRITE(numout,*)' kt=',kt, ' number of statically instable',   &
227               &                   ' water column : ',inpcc, ' number of iteration : ',inpci
228         ENDIF
229         !
230      ENDIF
231      !
232      IF( wrk_not_released(3, 1,2,3) .OR.   &
233          wrk_not_released_xz(1,2,3) )   CALL ctl_stop('tra_npc: failed to release workspace arrays')
234      !
235   END SUBROUTINE tra_npc
236
237   !!======================================================================
238END MODULE tranpc
Note: See TracBrowser for help on using the repository browser.