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.
trcldf.F90 in trunk/NEMOGCM/NEMO/TOP_SRC/TRP – NEMO

source: trunk/NEMOGCM/NEMO/TOP_SRC/TRP/trcldf.F90 @ 7698

Last change on this file since 7698 was 7698, checked in by mocavero, 7 years ago

update trunk with OpenMP parallelization
  • Property svn:keywords set to Id
File size: 16.4 KB
Line 
1MODULE trcldf
2   !!======================================================================
3   !!                       ***  MODULE  trcldf  ***
4   !! Ocean Passive tracers : lateral diffusive trends
5   !!=====================================================================
6   !! History :  1.0  ! 2005-11  (G. Madec)  Original code
7   !!            3.0  ! 2008-01  (C. Ethe, G. Madec)  merge TRC-TRA
8   !!            3.7  ! 2014-03  (G. Madec)  LDF simplification
9   !!----------------------------------------------------------------------
10#if defined key_top
11   !!----------------------------------------------------------------------
12   !!   'key_top'                                                TOP models
13   !!----------------------------------------------------------------------
14   !!   trc_ldf       : update the tracer trend with the lateral diffusion
15   !!   trc_ldf_ini   : initialization, namelist read, and parameters control
16   !!----------------------------------------------------------------------
17   USE trc            ! ocean passive tracers variables
18   USE oce_trc        ! ocean dynamics and active tracers
19   USE ldfslp         ! lateral diffusion: iso-neutral slope
20   USE traldf_lap_blp ! lateral diffusion: lap/bilaplacian iso-level      operator  (tra_ldf_lap/_blp   routine)
21   USE traldf_iso     ! lateral diffusion: laplacian iso-neutral standard operator  (tra_ldf_iso        routine)
22   USE traldf_triad   ! lateral diffusion: laplacian iso-neutral triad    operator  (tra_ldf_     triad routine)
23   USE trd_oce        ! trends: ocean variables
24   USE trdtra         ! trends manager: tracers
25   !
26   USE prtctl_trc     ! Print control
27
28   IMPLICIT NONE
29   PRIVATE
30
31   PUBLIC   trc_ldf   
32   PUBLIC   trc_ldf_ini   
33   !
34   LOGICAL , PUBLIC ::   ln_trcldf_lap       !:   laplacian operator
35   LOGICAL , PUBLIC ::   ln_trcldf_blp       !: bilaplacian operator
36   LOGICAL , PUBLIC ::   ln_trcldf_lev       !: iso-level   direction
37   LOGICAL , PUBLIC ::   ln_trcldf_hor       !: horizontal  direction (rotation to geopotential)
38   LOGICAL , PUBLIC ::   ln_trcldf_iso       !: iso-neutral direction (standard)
39   LOGICAL , PUBLIC ::   ln_trcldf_triad     !: iso-neutral direction (triad)
40   REAL(wp), PUBLIC ::   rn_ahtrc_0          !:   laplacian diffusivity coefficient for passive tracer [m2/s]
41   REAL(wp), PUBLIC ::   rn_bhtrc_0          !: bilaplacian      -          --     -       -   [m4/s]
42   REAL(wp), PUBLIC ::   rn_fact_lap         !: Enhanced zonal diffusivity coefficent in the equatorial domain
43   !
44   !                      !!: ** lateral mixing namelist (nam_trcldf) **
45   REAL(wp) ::  rldf       ! ratio between active and passive tracers diffusive coefficient
46   
47   INTEGER  ::  nldf = 0   ! type of lateral diffusion used defined from ln_trcldf_... namlist logicals)
48   
49   !! * Substitutions
50#  include "vectopt_loop_substitute.h90"
51   !!----------------------------------------------------------------------
52   !! NEMO/TOP 3.7 , NEMO Consortium (2014)
53   !! $Id$
54   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
55   !!----------------------------------------------------------------------
56CONTAINS
57
58   SUBROUTINE trc_ldf( kt )
59      !!----------------------------------------------------------------------
60      !!                  ***  ROUTINE tra_ldf  ***
61      !!
62      !! ** Purpose :   compute the lateral ocean tracer physics.
63      !!
64      !!----------------------------------------------------------------------
65      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
66      !
67      INTEGER            :: ji, jj, jk, jn
68      REAL(wp)           :: zdep
69      CHARACTER (len=22) :: charout
70      REAL(wp), POINTER, DIMENSION(:,:,:)   ::   zahu, zahv
71      REAL(wp), POINTER, DIMENSION(:,:,:,:) ::   ztrtrd
72      !!----------------------------------------------------------------------
73      !
74      IF( nn_timing == 1 )   CALL timing_start('trc_ldf')
75      !
76      IF( l_trdtrc )  THEN
77         CALL wrk_alloc( jpi,jpj,jpk,jptra,   ztrtrd )
78!$OMP PARALLEL DO schedule(static) private(jn,jk,jj,ji)
79         DO jn = 1, jptra
80            DO jk = 1, jpk
81               DO jj = 1, jpj
82                  DO ji = 1, jpi
83                     ztrtrd(ji,jj,jk,jn)  = tra(ji,jj,jk,jn)
84                  END DO
85               END DO
86            END DO
87         END DO
88      ENDIF
89      !                                  !* set the lateral diffusivity coef. for passive tracer     
90      CALL wrk_alloc( jpi,jpj,jpk,   zahu, zahv )
91!$OMP PARALLEL
92!$OMP DO schedule(static) private(jk,jj,ji)
93      DO jk = 1, jpk
94         DO jj = 1, jpj
95            DO ji = 1, jpi
96               zahu(ji,jj,jk) = rldf * ahtu(ji,jj,jk) 
97               zahv(ji,jj,jk) = rldf * ahtv(ji,jj,jk)
98            END DO
99         END DO
100      END DO
101      !                                  !* Enhanced zonal diffusivity coefficent in the equatorial domain
102!$OMP DO schedule(static) private(jk,jj,ji,zdep)
103      DO jk= 1, jpk
104         DO jj = 1, jpj
105            DO ji = 1, jpi
106               IF( gdept_n(ji,jj,jk) > 200. .AND. gphit(ji,jj) < 5. .AND. gphit(ji,jj) > -5. ) THEN
107                  zdep = MAX( gdept_n(ji,jj,jk) - 1000., 0. ) / 1000.
108                  zahu(ji,jj,jk) = zahu(ji,jj,jk) * MAX( 1., rn_fact_lap * EXP( -zdep ) )
109               ENDIF
110            END DO
111         END DO
112      END DO
113!$OMP END DO NOWAIT
114!$OMP END PARALLEL
115      !
116      SELECT CASE ( nldf )                     !* compute lateral mixing trend and add it to the general trend
117      !
118      CASE ( np_lap   )                               ! iso-level laplacian
119         CALL tra_ldf_lap  ( kt, nittrc000,'TRC', zahu, zahv, gtru, gtrv, gtrui, gtrvi, trb,      tra, jptra,  1   )
120         !
121      CASE ( np_lap_i )                               ! laplacian : standard iso-neutral operator (Madec)
122         CALL tra_ldf_iso  ( kt, nittrc000,'TRC', zahu, zahv, gtru, gtrv, gtrui, gtrvi, trb, trb, tra, jptra,  1   )
123         !
124      CASE ( np_lap_it )                              ! laplacian : triad iso-neutral operator (griffies)
125         CALL tra_ldf_triad( kt, nittrc000,'TRC', zahu, zahv, gtru, gtrv, gtrui, gtrvi, trb, trb, tra, jptra,  1   )
126         !
127      CASE ( np_blp , np_blp_i , np_blp_it )          ! bilaplacian: all operator (iso-level, -neutral)
128         CALL tra_ldf_blp  ( kt, nittrc000,'TRC', zahu, zahv, gtru, gtrv, gtrui, gtrvi, trb     , tra, jptra, nldf )
129         !
130      END SELECT
131      !
132      IF( l_trdtrc )   THEN                    ! send the trends for further diagnostics
133        DO jn = 1, jptra
134!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
135           DO jk = 1, jpk
136              DO jj = 1, jpj
137                 DO ji = 1, jpi
138                    ztrtrd(ji,jj,jk,jn) = tra(ji,jj,jk,jn) - ztrtrd(ji,jj,jk,jn)
139                 END DO
140              END DO
141           END DO
142           CALL trd_tra( kt, 'TRC', jn, jptra_ldf, ztrtrd(:,:,:,jn) )
143        END DO
144        CALL wrk_dealloc( jpi, jpj, jpk, jptra, ztrtrd )
145      ENDIF
146      !               
147      IF( ln_ctl ) THEN                        ! print mean trends (used for debugging)
148         WRITE(charout, FMT="('ldf ')")
149         CALL prt_ctl_trc_info(charout)
150         CALL prt_ctl_trc( tab4d=tra, mask=tmask, clinfo=ctrcnm, clinfo2='trd' )
151      ENDIF
152      !
153      CALL wrk_dealloc( jpi,jpj,jpk,   zahu, zahv )
154      !
155      IF( nn_timing == 1 )   CALL timing_stop('trc_ldf')
156      !
157   END SUBROUTINE trc_ldf
158
159
160   SUBROUTINE trc_ldf_ini
161      !!----------------------------------------------------------------------
162      !!                  ***  ROUTINE ldf_ctl  ***
163      !!
164      !! ** Purpose :   Define the operator for the lateral diffusion
165      !!
166      !! ** Method  :   set nldf from the namtra_ldf logicals
167      !!      nldf ==  0   laplacian operator
168      !!      nldf ==  1   Rotated laplacian operator
169      !!      nldf ==  2   bilaplacian operator
170      !!      nldf ==  3   Rotated bilaplacian
171      !!----------------------------------------------------------------------
172      INTEGER ::   ioptio, ierr   ! temporary integers
173      INTEGER ::   ios            ! Local integer output status for namelist read
174      !!
175      NAMELIST/namtrc_ldf/ ln_trcldf_lap, ln_trcldf_blp,                                  &
176         &                 ln_trcldf_lev, ln_trcldf_hor, ln_trcldf_iso, ln_trcldf_triad,  &
177         &                 rn_ahtrc_0   , rn_bhtrc_0, rn_fact_lap 
178      !!----------------------------------------------------------------------
179      !
180      REWIND( numnat_ref )             !  namtrc_ldf in reference namelist
181      READ  ( numnat_ref, namtrc_ldf, IOSTAT = ios, ERR = 903)
182903   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namtrc_ldf in reference namelist', lwp )
183      !
184      REWIND( numnat_cfg )             !  namtrc_ldf in configuration namelist
185      READ  ( numnat_cfg, namtrc_ldf, IOSTAT = ios, ERR = 904 )
186904   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namtrc_ldf in configuration namelist', lwp )
187      IF(lwm) WRITE ( numont, namtrc_ldf )
188      !
189      IF(lwp) THEN                     ! Namelist print
190         WRITE(numout,*)
191         WRITE(numout,*) 'trc_ldf_ini : lateral tracer diffusive operator'
192         WRITE(numout,*) '~~~~~~~~~~~'
193         WRITE(numout,*) '   Namelist namtrc_ldf : set lateral mixing parameters (type, direction, coefficients)'
194         WRITE(numout,*) '      operator'
195         WRITE(numout,*) '           laplacian                 ln_trcldf_lap   = ', ln_trcldf_lap
196         WRITE(numout,*) '         bilaplacian                 ln_trcldf_blp   = ', ln_trcldf_blp
197         WRITE(numout,*) '      direction of action'
198         WRITE(numout,*) '         iso-level                   ln_trcldf_lev   = ', ln_trcldf_lev
199         WRITE(numout,*) '         horizontal (geopotential)   ln_trcldf_hor   = ', ln_trcldf_hor
200         WRITE(numout,*) '         iso-neutral (standard)      ln_trcldf_iso   = ', ln_trcldf_iso
201         WRITE(numout,*) '         iso-neutral (triad)         ln_trcldf_triad = ', ln_trcldf_triad
202         WRITE(numout,*) '      diffusivity coefficient'
203         WRITE(numout,*) '           laplacian                 rn_ahtrc_0      = ', rn_ahtrc_0
204         WRITE(numout,*) '         bilaplacian                 rn_bhtrc_0      = ', rn_bhtrc_0
205         WRITE(numout,*) '      enhanced zonal diffusivity     rn_fact_lap     = ', rn_fact_lap
206
207      ENDIF
208      !     
209      !                                ! control the namelist parameters
210      ioptio = 0
211      IF( ln_trcldf_lap )   ioptio = ioptio + 1
212      IF( ln_trcldf_blp )   ioptio = ioptio + 1
213      IF( ioptio >  1   )   CALL ctl_stop( 'trc_ldf_ctl: use ONE or NONE of the 2 lap/bilap operator type on tracer' )
214      IF( ioptio == 0   )   nldf = np_no_ldf   ! No lateral diffusion
215     
216      IF( ln_trcldf_lap .AND. ln_trcldf_blp )   CALL ctl_stop( 'trc_ldf_ctl: bilaplacian should be used on both TRC and TRA' )
217      IF( ln_trcldf_blp .AND. ln_trcldf_lap )   CALL ctl_stop( 'trc_ldf_ctl:   laplacian should be used on both TRC and TRA' )
218      !
219      ioptio = 0
220      IF( ln_trcldf_lev )   ioptio = ioptio + 1
221      IF( ln_trcldf_hor )   ioptio = ioptio + 1
222      IF( ln_trcldf_iso )   ioptio = ioptio + 1
223      IF( ioptio /= 1   )   CALL ctl_stop( 'trc_ldf_ctl: use only ONE direction (level/hor/iso)' )
224      !
225      ! defined the type of lateral diffusion from ln_trcldf_... logicals
226      ! CAUTION : nldf = 1 is used in trazdf_imp, change it carefully
227      ierr = 0
228      IF( ln_trcldf_lap ) THEN      !==  laplacian operator  ==!
229         IF ( ln_zco ) THEN                ! z-coordinate
230            IF ( ln_trcldf_lev   )   nldf = np_lap     ! iso-level = horizontal (no rotation)
231            IF ( ln_trcldf_hor   )   nldf = np_lap     ! iso-level = horizontal (no rotation)
232            IF ( ln_trcldf_iso   )   nldf = np_lap_i   ! iso-neutral: standard  (   rotation)
233            IF ( ln_trcldf_triad )   nldf = np_lap_it  ! iso-neutral: triad     (   rotation)
234         ENDIF
235         IF ( ln_zps ) THEN             ! z-coordinate with partial step
236            IF ( ln_trcldf_lev   )   ierr = 1         ! iso-level not allowed
237            IF ( ln_trcldf_hor   )   nldf = np_lap     ! horizontal (no rotation)
238            IF ( ln_trcldf_iso   )   nldf = np_lap_i   ! iso-neutral: standard (rotation)
239            IF ( ln_trcldf_triad )   nldf = np_lap_it  ! iso-neutral: triad    (rotation)
240         ENDIF
241         IF ( ln_sco ) THEN             ! s-coordinate
242            IF ( ln_trcldf_lev   )   nldf = np_lap     ! iso-level  (no rotation)
243            IF ( ln_trcldf_hor   )   nldf = np_lap_it  ! horizontal (   rotation)       !!gm   a checker....
244            IF ( ln_trcldf_iso   )   nldf = np_lap_i   ! iso-neutral: standard (rotation)
245            IF ( ln_trcldf_triad )   nldf = np_lap_it  ! iso-neutral: triad    (rotation)
246         ENDIF
247         !                                ! diffusivity ratio: passive / active tracers
248         IF( ABS(rn_aht_0) < 2._wp*TINY(1._wp) ) THEN
249            IF( ABS(rn_ahtrc_0) < 2._wp*TINY(1._wp) ) THEN
250               rldf = 1.0_wp
251            ELSE
252               CALL ctl_stop( 'trc_ldf_ctl : cannot define rldf, rn_aht_0==0, rn_ahtrc_0 /=0' )
253            ENDIF
254         ELSE
255            rldf = rn_ahtrc_0 / rn_aht_0
256         ENDIF
257      ENDIF
258      !
259      IF( ln_trcldf_blp ) THEN      !==  bilaplacian operator  ==!
260         IF ( ln_zco ) THEN                ! z-coordinate
261            IF ( ln_trcldf_lev   )   nldf = np_blp     ! iso-level = horizontal (no rotation)
262            IF ( ln_trcldf_hor   )   nldf = np_blp     ! iso-level = horizontal (no rotation)
263            IF ( ln_trcldf_iso   )   nldf = np_blp_i   ! iso-neutral: standard (rotation)
264            IF ( ln_trcldf_triad )   nldf = np_blp_it  ! iso-neutral: triad    (rotation)
265         ENDIF
266         IF ( ln_zps ) THEN             ! z-coordinate with partial step
267            IF ( ln_trcldf_lev   )   ierr = 1         ! iso-level not allowed
268            IF ( ln_trcldf_hor   )   nldf = np_blp     ! horizontal (no rotation)
269            IF ( ln_trcldf_iso   )   nldf = np_blp_i   ! iso-neutral: standard (rotation)
270            IF ( ln_trcldf_triad )   nldf = np_blp_it  ! iso-neutral: triad    (rotation)
271         ENDIF
272         IF ( ln_sco ) THEN             ! s-coordinate
273            IF ( ln_trcldf_lev   )   nldf = np_blp     ! iso-level  (no rotation)
274            IF ( ln_trcldf_hor   )   nldf = np_blp_it  ! horizontal (   rotation)       !!gm   a checker....
275            IF ( ln_trcldf_iso   )   nldf = np_blp_i   ! iso-neutral: standard (rotation)
276            IF ( ln_trcldf_triad )   nldf = np_blp_it  ! iso-neutral: triad    (rotation)
277         ENDIF
278         !                                ! diffusivity ratio: passive / active tracers
279         IF( ABS(rn_bht_0) < 2._wp*TINY(1._wp) ) THEN
280            IF( ABS(rn_bhtrc_0) < 2._wp*TINY(1._wp) ) THEN
281               rldf = 1.0_wp
282            ELSE
283               CALL ctl_stop( 'trc_ldf_ctl : cannot define rldf, rn_aht_0==0, rn_ahtrc_0 /=0' )
284            ENDIF
285         ELSE
286            rldf = SQRT(  ABS( rn_bhtrc_0 / rn_bht_0 )  )
287         ENDIF
288      ENDIF
289      !
290      IF( ierr == 1 )   CALL ctl_stop( 'trc_ldf_ctl: iso-level in z-partial step, not allowed' )
291      IF( ln_ldfeiv .AND. .NOT.ln_trcldf_iso )   CALL ctl_stop( 'trc_ldf_ctl: eiv requires isopycnal laplacian diffusion' )
292      IF( nldf == 1 .OR. nldf == 3 )   l_ldfslp = .TRUE.    ! slope of neutral surfaces required
293      !
294      IF(lwp) THEN
295         WRITE(numout,*)
296         SELECT CASE( nldf )
297         CASE( np_no_ldf )   ;   WRITE(numout,*) '          NO lateral diffusion'
298         CASE( np_lap    )   ;   WRITE(numout,*) '          laplacian iso-level operator'
299         CASE( np_lap_i  )   ;   WRITE(numout,*) '          Rotated laplacian operator (standard)'
300         CASE( np_lap_it )   ;   WRITE(numout,*) '          Rotated laplacian operator (triad)'
301         CASE( np_blp    )   ;   WRITE(numout,*) '          bilaplacian iso-level operator'
302         CASE( np_blp_i  )   ;   WRITE(numout,*) '          Rotated bilaplacian operator (standard)'
303         CASE( np_blp_it )   ;   WRITE(numout,*) '          Rotated bilaplacian operator (triad)'
304         END SELECT
305      ENDIF
306      !
307   END SUBROUTINE trc_ldf_ini
308#else
309   !!----------------------------------------------------------------------
310   !!   Default option                                         Empty module
311   !!----------------------------------------------------------------------
312CONTAINS
313   SUBROUTINE trc_ldf( kt )
314      INTEGER, INTENT(in) :: kt
315      WRITE(*,*) 'trc_ldf: You should not have seen this print! error?', kt
316   END SUBROUTINE trc_ldf
317#endif
318   !!======================================================================
319END MODULE trcldf
Note: See TracBrowser for help on using the repository browser.