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

Last change on this file since 4424 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
RevLine 
[3]1MODULE tranpc
2   !!==============================================================================
3   !!                       ***  MODULE  tranpc  ***
4   !! Ocean active tracers:  non penetrative convection scheme
5   !!==============================================================================
[1537]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
[2528]10   !!            3.3  ! 2010-05  (C. Ethe, G. Madec)  merge TRC-TRA
[503]11   !!----------------------------------------------------------------------
[3]12
13   !!----------------------------------------------------------------------
[1537]14   !!   tra_npc : apply the non penetrative convection scheme
[3]15   !!----------------------------------------------------------------------
16   USE oce             ! ocean dynamics and active tracers
17   USE dom_oce         ! ocean space and time domain
[1537]18   USE zdf_oce         ! ocean vertical physics
[2528]19   USE trdmod_oce      ! ocean active tracer trends
[2715]20   USE trdtra          ! ocean active tracer trends
[3]21   USE eosbn2          ! equation of state (eos routine)
22   USE lbclnk          ! lateral boundary conditions (or mpp link)
[216]23   USE in_out_manager  ! I/O manager
[2715]24   USE lib_mpp         ! MPP library
[3]25
26   IMPLICIT NONE
27   PRIVATE
28
[2715]29   PUBLIC   tra_npc       ! routine called by step.F90
[3]30
[3211]31   !! * Control permutation of array indices
32#  include "oce_ftrans.h90"
33#  include "dom_oce_ftrans.h90"
34#  include "zdf_oce_ftrans.h90"
35
[3]36   !! * Substitutions
37#  include "domzgr_substitute.h90"
38   !!----------------------------------------------------------------------
[2528]39   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
[1146]40   !! $Id$
[2528]41   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
[3]42   !!----------------------------------------------------------------------
43CONTAINS
44
45   SUBROUTINE tra_npc( kt )
46      !!----------------------------------------------------------------------
47      !!                  ***  ROUTINE tranpc  ***
48      !!
49      !! ** Purpose :   Non penetrative convective adjustment scheme. solve
[1111]50      !!      the static instability of the water column on after fields
[3]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.
[503]56      !!      l_trdtra=T: the trend associated with this algorithm is saved.
[3]57      !!
[1111]58      !! ** Action  : - (ta,sa) after the application od the npc scheme
[3]59      !!              - save the associated trends (ttrd,strd) ('key_trdtra')
60      !!
[503]61      !! References : Madec, et al., 1991, JPO, 21, 9, 1349-1371.
62      !!----------------------------------------------------------------------
[2715]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
[3211]66
67      !! DCSE_NEMO: need additional directives for renamed module variables
68!FTRANS ztrdt ztrds zrhop :I :I :z
69
[2715]70      !
[503]71      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
[2715]72      !
[503]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
[3]79      !!----------------------------------------------------------------------
[216]80
[2715]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
[1537]87      IF( MOD( kt, nn_npc ) == 0 ) THEN
[3]88
89         inpcc = 0
90         inpci = 0
91
[2528]92         CALL eos( tsa, rhd, zrhop )         ! Potential density
[3]93
[2528]94         IF( l_trdtra )   THEN                    !* Save ta and sa trends
[2715]95            ztrdt(:,:,:) = tsa(:,:,:,jp_tem) 
96            ztrds(:,:,:) = tsa(:,:,:,jp_sal)
[216]97         ENDIF
98
[3]99         !                                                ! ===============
100         DO jj = 1, jpj                                   !  Vertical slab
101            !                                             ! ===============
[503]102            !  Static instability pointer
103            ! ----------------------------
[3211]104#if defined key_z_first
105            DO ji = 1, jpi
106               DO jk = 1, jpkm1
107#else
[3]108            DO jk = 1, jpkm1
109               DO ji = 1, jpi
[3211]110#endif
[3]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
[503]123            IF( jj == 1 )   zwx(:,:) = 0.e0
[3]124
125            DO jk = 1, jpkm1
126               DO ji = 1, jpi
127                  zwx(ji,jk) = 1.
[503]128                  IF( zwx(ji,jk) < 1.e-5 ) zwx(ji,jk) = 0.e0
[3]129               END DO
130            END DO
131
[503]132            zwy(:,1) = 0.e0
[3]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
[503]139            zwz(1,1) = 0.e0
[3]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
[503]150            IF( zwz(1,1) /= 0.e0 ) THEN         ! -->> the density profil is statically instable :
[3]151               DO ji = 1, jpi
[503]152                  IF( zwy(ji,1) /= 0.e0 ) THEN
153                     !
[2528]154                     ikbot = mbkt(ji,jj)        ! ikbot: ocean bottom T-level
[503]155                     !
156                     DO jiter = 1, jpk          ! vertical iteration
157                        !
[3]158                        ! search of ikup : the first static instability from the sea surface
[503]159                        !
[3]160                        ik = 0
161220                     CONTINUE
162                        ik = ik + 1
[2528]163                        IF( ik >= ikbot ) GO TO 200
[3]164                        zwx(ji,ik) = zrhop(ji,jj,ik) - zrhop(ji,jj,ik+1)
[503]165                        IF( zwx(ji,ik) <= 0.e0 ) GO TO 220
[3]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
[503]170                        !
[3]171                        ze3tot= fse3t(ji,jj,ikup)
[2528]172                        zta   = tsa  (ji,jj,ikup,jp_tem)
173                        zsa   = tsa  (ji,jj,ikup,jp_sal)
[3]174                        zraua = zrhop(ji,jj,ikup)
[503]175                        !
[3]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
[2528]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
[3]185                           zraua = ( zraua*(ze3tot-ze3dwn) + zrhop(ji,jj,jkdown)*ze3dwn )/ze3tot
186                           inpci = inpci+1
187                        END DO
188                        ikdown = ikbot-1
189240                     CONTINUE
[503]190                        !
[3]191                        DO jkp = ikup, ikdown-1
[2528]192                           tsa  (ji,jj,jkp,jp_tem) = zta
193                           tsa  (ji,jj,jkp,jp_sal) = zsa
194                           zrhop(ji,jj,jkp       ) = zraua
[3]195                        END DO
196                        IF (ikdown == ikbot-1 .AND. zraua >= zrhop(ji,jj,ikdown) ) THEN
[2528]197                           tsa  (ji,jj,jkp,jp_tem) = zta
198                           tsa  (ji,jj,jkp,jp_sal) = zsa
199                           zrhop(ji,jj,ikdown    ) = zraua
[3]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         !                                                ! ===============
[503]210         !
211         IF( l_trdtra )   THEN         ! save the Non penetrative mixing trends for diagnostic
[2528]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]216         ENDIF
[3]217     
[1111]218         ! Lateral boundary conditions on ( ta, sa )   ( Unchanged sign)
[3]219         ! ------------------------------============
[2715]220         CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1. )   ;   CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1. )
[3]221     
222
223         !  2. non penetrative convective scheme statistics
224         !  -----------------------------------------------
[1537]225         IF( nn_npcp /= 0 .AND. MOD( kt, nn_npcp ) == 0 ) THEN
[3]226            IF(lwp) WRITE(numout,*)' kt=',kt, ' number of statically instable',   &
[503]227               &                   ' water column : ',inpcc, ' number of iteration : ',inpci
[3]228         ENDIF
[503]229         !
[3]230      ENDIF
[503]231      !
[2715]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      !
[3]235   END SUBROUTINE tra_npc
236
237   !!======================================================================
238END MODULE tranpc
Note: See TracBrowser for help on using the repository browser.