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.
traldf_lap.F90 in branches/2014/dev_CNRS0_NOC1_LDF/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

source: branches/2014/dev_CNRS0_NOC1_LDF/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_lap.F90 @ 4596

Last change on this file since 4596 was 4596, checked in by gm, 10 years ago

#1260: LDF simplification + bilap iso-neutral for TRA and GYRE
  • Property svn:keywords set to Id
File size: 11.7 KB
Line 
1MODULE traldf_lap
2   !!==============================================================================
3   !!                       ***  MODULE  traldf_lap  ***
4   !! Ocean tracers:  lateral diffusivity trend  (laplacian and bilaplacian)
5   !!==============================================================================
6   !! History :  OPA  ! 1987-06  (P. Andrich, D. L Hostis)  Original code
7   !!                 ! 1991-11  (G. Madec)
8   !!                 ! 1995-11  (G. Madec)  suppress volumetric scale factors
9   !!                 ! 1996-01  (G. Madec)  statement function for e3
10   !!            NEMO ! 2002-06  (G. Madec)  F90: Free form and module
11   !!            1.0  ! 2004-08  (C. Talandier) New trends organization
12   !!                 ! 2005-11  (G. Madec)  add zps case
13   !!            3.0  ! 2010-06  (C. Ethe, G. Madec) Merge TRA-TRC
14   !!            3.7  ! 2014-01  (G. Madec) re-entrant laplacian
15   !!----------------------------------------------------------------------
16
17   !!----------------------------------------------------------------------
18   !!   tra_ldf_lap   : update the tracer trend with the lateral diffusion : iso-level laplacian operator
19   !!   tra_ldf_bilap : update the tracer trend with the lateral diffusion : iso-level bilaplacian operator
20   !!----------------------------------------------------------------------
21   USE oce             ! ocean dynamics and active tracers
22   USE dom_oce         ! ocean space and time domain
23   USE ldftra          ! lateral physics: eddy diffusivity
24   USE diaptr          ! poleward transport diagnostics
25   USE trc_oce         ! share passive tracers/Ocean variables
26   USE zpshde          ! partial step: hor. derivative     (zps_hde routine)
27   !
28   USE in_out_manager  ! I/O manager
29   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
30   USE lib_mpp         ! distribued memory computing library
31   USE timing          ! Timing
32   USE wrk_nemo        ! Memory allocation
33
34   IMPLICIT NONE
35   PRIVATE
36
37   PUBLIC   tra_ldf_lap   ! routine called by step.F90
38   PUBLIC   tra_ldf_blp   ! routine called by step.F90
39
40   REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:) ::   e1ur, e2vr   ! scale factor coefficients
41
42   !! * Substitutions
43#  include "domzgr_substitute.h90"
44#  include "vectopt_loop_substitute.h90"
45   !!----------------------------------------------------------------------
46   !! NEMO/OPA 3.7 , NEMO Consortium (2014)
47   !! $Id$
48   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
49   !!----------------------------------------------------------------------
50CONTAINS
51
52   SUBROUTINE tra_ldf_lap( kt, kit000, cdtype, pgu, pgv, pahu, pahv,   &
53      &                                        ptb, pta, kjpt, kpass ) 
54      !!----------------------------------------------------------------------
55      !!                  ***  ROUTINE tra_ldf_lap  ***
56      !!                   
57      !! ** Purpose :   Compute the before horizontal tracer (t & s) diffusive
58      !!      trend and add it to the general trend of tracer equation.
59      !!
60      !! ** Method  :   Second order diffusive operator evaluated using before
61      !!      fields (forward time scheme). The horizontal diffusive trends of
62      !!      the tracer is given by:
63      !!          difft = 1/(e1t*e2t*e3t) {  di-1[ pahu e2u*e3u/e1u di(tb) ]
64      !!                                   + dj-1[ pahv e1v*e3v/e2v dj(tb) ] }
65      !!      Add this trend to the general tracer trend pta :
66      !!          pta = pta + difft
67      !!
68      !! ** Action  : - Update pta arrays with the before iso-level
69      !!                harmonic mixing trend.
70      !!----------------------------------------------------------------------
71      INTEGER                              , INTENT(in   ) ::   kt         ! ocean time-step index
72      INTEGER                              , INTENT(in   ) ::   kit000     ! first time step index
73      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator)
74      INTEGER                              , INTENT(in   ) ::   kjpt       ! number of tracers
75      INTEGER                              , INTENT(in   ) ::   kpass      ! =1/2 first or second passage
76      REAL(wp), DIMENSION(jpi,jpj    ,kjpt), INTENT(in   ) ::   pgu, pgv   ! tracer gradient at pstep levels
77      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(in   ) ::   pahu, pahv ! eddy diffusivity at u- and v-points  [m2/s]
78      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb        ! before and now tracer fields
79      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta        ! tracer trend
80      !
81      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
82      INTEGER  ::   iku, ikv         ! local integers
83      REAL(wp) ::   zsign            ! local scalars
84      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztu, ztv, zaheeu, zaheev
85      !!----------------------------------------------------------------------
86      !
87      IF( kt == nit000 .AND. lwp )  THEN
88         WRITE(numout,*)
89         WRITE(numout,*) 'tra_ldf_lap : iso-level laplacian diffusion on ', cdtype, ', pass=', kpass
90         WRITE(numout,*) '~~~~~~~~~~~ '
91      ENDIF
92      !
93      IF( nn_timing == 1 )  CALL timing_start('tra_ldf_lap')
94      !
95      CALL wrk_alloc( jpi, jpj, jpk, ztu, ztv, zaheeu, zaheev ) 
96      !
97      DO jk = 1, jpkm1           !==  Initialization of metric arrays used for all tracers  ==!
98         DO jj = 1, jpjm1
99            DO ji = 1, fs_jpim1   ! vector opt.
100               zaheeu(ji,jj,jk) = pahu(ji,jj,jk) * e2u(ji,jj) * fse3u(ji,jj,jk) / e1u(ji,jj)   !!gm   * umask(ji,jj,jk)
101               zaheev(ji,jj,jk) = pahv(ji,jj,jk) * e1v(ji,jj) * fse3v(ji,jj,jk) / e2v(ji,jj)   !!gm   * vmask(ji,jj,jk)
102            END DO
103         END DO
104      END DO
105      IF( kpass == 1 ) THEN   ;   zsign =  1._wp      ! bilaplacian operator require a minus sign (eddy diffusivity >0)
106      ELSE                    ;   zsign = -1._wp
107      ENDIF
108      !
109      !                             ! =========== !
110      DO jn = 1, kjpt               ! tracer loop !
111         !                          ! =========== !   
112         !                               
113         DO jk = 1, jpkm1              !== First derivative (gradient)  ==!
114            DO jj = 1, jpjm1
115               DO ji = 1, fs_jpim1
116                  ztu(ji,jj,jk) = zaheeu(ji,jj,jk) * ( ptb(ji+1,jj  ,jk,jn) - ptb(ji,jj,jk,jn) )
117                  ztv(ji,jj,jk) = zaheev(ji,jj,jk) * ( ptb(ji  ,jj+1,jk,jn) - ptb(ji,jj,jk,jn) )
118               END DO
119            END DO
120         END DO 
121         IF( ln_zps ) THEN                ! set gradient at partial step level
122            DO jj = 1, jpjm1
123               DO ji = 1, fs_jpim1
124                  iku = mbku(ji,jj)             ! last level
125                  ikv = mbkv(ji,jj)
126                  ztu(ji,jj,iku) = zaheeu(ji,jj,iku) * pgu(ji,jj,jn)
127                  ztv(ji,jj,ikv) = zaheev(ji,jj,ikv) * pgv(ji,jj,jn)
128               END DO
129            END DO 
130         ENDIF
131         !
132         DO jk = 1, jpkm1              !== Second derivative (divergence) added to the general tracer trends  ==!
133            DO jj = 2, jpjm1
134               DO ji = fs_2, fs_jpim1
135                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zsign * ( ztu(ji,jj,jk) - ztu(ji-1,jj,jk)     &
136                     &                                          + ztv(ji,jj,jk) - ztv(ji,jj-1,jk) )   &
137                     &                                        / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) )
138               END DO
139            END DO
140         END DO 
141         !
142         !                             !== "Poleward" diffusive heat or salt transports  ==!
143         IF( ( kpass == 1 .AND. .NOT.ln_traldf_blp ) .OR.  &     !==  first pass only (  laplacian)  ==!
144             ( kpass == 2 .AND.      ln_traldf_blp ) ) THEN      !==  2nd   pass only (bilaplacian)  ==!
145            IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN
146               IF( jn  == jp_tem)   htr_ldf(:) = ptr_vj( ztv(:,:,:) )
147               IF( jn  == jp_sal)   str_ldf(:) = ptr_vj( ztv(:,:,:) )
148            ENDIF
149         ENDIF
150         !                          ! ==================
151      END DO                        ! end of tracer loop
152      !                             ! ==================
153      !
154      CALL wrk_dealloc( jpi, jpj, jpk, ztu, ztv, zaheeu, zaheev ) 
155      !
156      IF( nn_timing == 1 )  CALL timing_stop('tra_ldf_lap')
157      !
158   END SUBROUTINE tra_ldf_lap
159   
160
161   SUBROUTINE tra_ldf_blp( kt, kit000, cdtype, pgu, pgv, pahu, pahv,   &
162      &                                        ptb, pta, kjpt        )
163      !!----------------------------------------------------------------------
164      !!                 ***  ROUTINE tra_ldf_blp  ***
165      !!                   
166      !! ** Purpose :   Compute the before horizontal tracer diffusive
167      !!      trend and add it to the general trend of tracer equation.
168      !!
169      !! ** Method  :   The lateral diffusive trends is provided by a bilaplacian
170      !!      operator applied to before field (forward in time).
171      !!      It is computed by two successive calls to tra_ldf_lap routine
172      !!
173      !! ** Action :   pta   updated with the before rotated bilaplacian diffusion
174      !!----------------------------------------------------------------------
175      INTEGER                              , INTENT(in   ) ::   kt         ! ocean time-step index
176      INTEGER                              , INTENT(in   ) ::   kit000     ! first time step index
177      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator)
178      INTEGER                              , INTENT(in   ) ::   kjpt       ! number of tracers
179      REAL(wp), DIMENSION(jpi,jpj    ,kjpt), INTENT(in   ) ::   pgu, pgv   ! tracer gradient at pstep levels
180      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(in   ) ::   pahu, pahv ! eddy diffusivity at u- and v-points  [m2/s]
181      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb        ! before and now tracer fields
182      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta        ! tracer trend
183      !
184      INTEGER ::   ji, jj, jk, jn   ! dummy loop indices
185      REAL(wp), POINTER, DIMENSION(:,:,:,:) :: zlap         ! laplacian at t-point
186      REAL(wp), POINTER, DIMENSION(:,:,:)   :: zglu, zglv   ! gradient of the laplacian at partial step level (u- and v-points)
187      !!----------------------------------------------------------------------
188      !
189      IF( nn_timing == 1 )  CALL timing_start('tra_ldf_iso_blp')
190      !
191      CALL wrk_alloc( jpi, jpj, jpk, kjpt, zlap ) 
192      CALL wrk_alloc( jpi, jpj     , kjpt, zglu, zglv ) 
193      !
194      IF( kt == nit000 )  THEN
195         IF(lwp) WRITE(numout,*)
196         IF(lwp) WRITE(numout,*) 'tra_ldf_iso_blp : iso-neutral biharmonic operator on ', cdtype
197         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~'
198      ENDIF
199
200      zlap(:,:,:,:) = 0._wp
201      CALL tra_ldf_lap( kt, kit000, cdtype, pgu, pgv, pahu, pahv, ptb, zlap, kjpt, 1 )   ! rotated laplacian applied to ptb (output in zlap)
202      !
203      DO jn = 1, kjpt
204         CALL lbc_lnk( zlap(:,:,:,jn) , 'T', 1. )                     ! Lateral boundary conditions (unchanged sign)
205      END DO
206      !
207      IF( ln_zps )   CALL zps_hde( kt, jpts, zlap, zglu, zglv )       ! Partial steps: hor. gradient of laplacian at the partial step level     
208      !
209      CALL tra_ldf_lap( kt, kit000, cdtype, pgu, pgv, pahu, pahv, zlap, pta, kjpt, 2 )   ! rotated laplacian applied to zlap (output in pta)
210      !
211      CALL wrk_dealloc( jpi, jpj, jpk, kjpt, zlap ) 
212      CALL wrk_dealloc( jpi, jpj     , kjpt, zglu, zglv ) 
213      !
214      IF( nn_timing == 1 )  CALL timing_stop('tra_ldf_iso_blp')
215      !
216   END SUBROUTINE tra_ldf_blp
217
218   !!==============================================================================
219END MODULE traldf_lap
Note: See TracBrowser for help on using the repository browser.