source: branches/UKMO/dev_r5518_GO6_package_text_diagnostics/NEMOGCM/NEMO/OPA_SRC/TRA/traldf.F90 @ 10774

Last change on this file since 10774 was 10774, checked in by andmirek, 20 months ago

GMED 450 add flush after prints

File size: 18.1 KB
Line 
1MODULE traldf
2   !!======================================================================
3   !!                       ***  MODULE  traldf  ***
4   !! Ocean Active tracers : lateral diffusive trends
5   !!=====================================================================
6   !! History :  9.0  ! 2005-11 (G. Madec)  Original code
7   !!       NEMO 3.0  ! 2008-01  (C. Ethe, G. Madec)  merge TRC-TRA
8   !!----------------------------------------------------------------------
9
10   !!----------------------------------------------------------------------
11   !!   tra_ldf      : update the tracer trend with the lateral diffusion
12   !!   tra_ldf_init : initialization, namelist read, and parameters control
13   !!       ldf_ano  : compute lateral diffusion for constant T-S profiles
14   !!----------------------------------------------------------------------
15   USE oce             ! ocean dynamics and tracers
16   USE dom_oce         ! ocean space and time domain
17   USE phycst          ! physical constants
18   USE ldftra_oce      ! ocean tracer   lateral physics
19   USE ldfslp          ! ???
20   USE traldf_bilapg   ! lateral mixing            (tra_ldf_bilapg routine)
21   USE traldf_bilap    ! lateral mixing             (tra_ldf_bilap routine)
22   USE traldf_iso      ! lateral mixing               (tra_ldf_iso routine)
23   USE traldf_iso_grif ! lateral mixing          (tra_ldf_iso_grif routine)
24   USE traldf_lap      ! lateral mixing               (tra_ldf_lap routine)
25   USE trd_oce         ! trends: ocean variables
26   USE trdtra          ! trends manager: tracers
27   !
28   USE prtctl          ! Print control
29   USE in_out_manager  ! I/O manager
30   USE lib_mpp         ! distribued memory computing library
31   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
32   USE wrk_nemo        ! Memory allocation
33   USE timing          ! Timing
34
35   IMPLICIT NONE
36   PRIVATE
37
38   PUBLIC   tra_ldf        ! called by step.F90
39   PUBLIC   tra_ldf_init   ! called by opa.F90
40   !
41   INTEGER ::   nldf = 0   ! type of lateral diffusion used defined from ln_traldf_... namlist logicals)
42
43   REAL, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   t0_ldf, s0_ldf   !: lateral diffusion trends of T & S for a cst profile
44   !                                                               !  (key_traldf_ano only)
45
46   !! * Substitutions
47#  include "domzgr_substitute.h90"
48#  include "vectopt_loop_substitute.h90"
49   !!----------------------------------------------------------------------
50   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
51   !! $Id$
52   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
53   !!----------------------------------------------------------------------
54CONTAINS
55
56   SUBROUTINE tra_ldf( kt )
57      !!----------------------------------------------------------------------
58      !!                  ***  ROUTINE tra_ldf  ***
59      !!
60      !! ** Purpose :   compute the lateral ocean tracer physics.
61      !!----------------------------------------------------------------------
62      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
63      !!
64      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdt, ztrds
65      !!----------------------------------------------------------------------
66      !
67      IF( nn_timing == 1 )  CALL timing_start('tra_ldf')
68      !
69      rldf = 1     ! For active tracers the
70      r_fact_lap(:,:,:) = 1.0
71
72      IF( l_trdtra )   THEN                    !* Save ta and sa trends
73         CALL wrk_alloc( jpi, jpj, jpk, ztrdt, ztrds ) 
74         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) 
75         ztrds(:,:,:) = tsa(:,:,:,jp_sal)
76      ENDIF
77
78      SELECT CASE ( nldf )                       ! compute lateral mixing trend and add it to the general trend
79      CASE ( 0 )   ;   CALL tra_ldf_lap     ( kt, nit000, 'TRA', gtsu, gtsv, gtui, gtvi,        &
80                               &                                   tsb, tsa, jpts        )  ! iso-level laplacian
81      CASE ( 1 )                                                                              ! rotated laplacian
82         IF( ln_traldf_grif ) THEN                                                         
83                       CALL tra_ldf_iso_grif( kt, nit000,'TRA', gtsu, gtsv, tsb, tsa, jpts, ahtb0 )      ! Griffies operator
84         ELSE                                                                               
85                       CALL tra_ldf_iso     ( kt, nit000, 'TRA', gtsu, gtsv, gtui, gtvi,        &
86                               &                                  tsb, tsa, jpts, ahtb0 )      ! Madec operator
87         ENDIF
88      CASE ( 2 )   ;   CALL tra_ldf_bilap   ( kt, nit000, 'TRA', gtsu, gtsv, gtui, gtvi,        &
89                               &                                   tsb, tsa, jpts        )  ! iso-level bilaplacian
90      CASE ( 3 )   ;   CALL tra_ldf_bilapg  ( kt, nit000, 'TRA',             tsb, tsa, jpts        )  ! s-coord. geopot. bilap.
91         !
92      CASE ( -1 )                                ! esopa: test all possibility with control print
93         CALL tra_ldf_lap   ( kt, nit000, 'TRA', gtsu, gtsv, gtui, gtvi,        &
94         &                                       tsb, tsa, jpts        ) 
95         CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' ldf0 - Ta: ', mask1=tmask,               &
96         &             tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
97         IF( ln_traldf_grif ) THEN
98            CALL tra_ldf_iso_grif( kt, nit000, 'TRA', gtsu, gtsv, tsb, tsa, jpts, ahtb0 )
99         ELSE
100            CALL tra_ldf_iso     ( kt, nit000, 'TRA', gtsu, gtsv, gtui, gtvi,        &
101            &                                               tsb, tsa, jpts, ahtb0 ) 
102         ENDIF
103         CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' ldf1 - Ta: ', mask1=tmask,               &
104         &             tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
105         CALL tra_ldf_bilap ( kt, nit000, 'TRA', gtsu, gtsv, gtui, gtvi,        &
106         &                                       tsb, tsa, jpts        ) 
107         CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' ldf2 - Ta: ', mask1=tmask,               &
108         &             tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
109         CALL tra_ldf_bilapg( kt, nit000, 'TRA',             tsb, tsa, jpts        ) 
110         CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' ldf3 - Ta: ', mask1=tmask,               &
111         &             tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
112      END SELECT
113
114#if defined key_traldf_ano
115      tsa(:,:,:,jp_tem) = tsa(:,:,:,jp_tem) - t0_ldf(:,:,:)      ! anomaly: substract the reference diffusivity
116      tsa(:,:,:,jp_sal) = tsa(:,:,:,jp_sal) - s0_ldf(:,:,:)
117#endif
118
119      IF( l_trdtra )   THEN                      ! save the horizontal diffusive trends for further diagnostics
120         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)
121         ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:)
122         CALL trd_tra( kt, 'TRA', jp_tem, jptra_ldf, ztrdt )
123         CALL trd_tra( kt, 'TRA', jp_sal, jptra_ldf, ztrds )
124         CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds ) 
125      ENDIF
126      !                                          ! print mean trends (used for debugging)
127      IF(ln_ctl)   CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' ldf  - Ta: ', mask1=tmask,               &
128         &                       tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
129      !
130      IF( nn_timing == 1 )  CALL timing_stop('tra_ldf')
131      !
132   END SUBROUTINE tra_ldf
133
134
135   SUBROUTINE tra_ldf_init
136      !!----------------------------------------------------------------------
137      !!                  ***  ROUTINE tra_ldf_init  ***
138      !!
139      !! ** Purpose :   Choice of the operator for the lateral tracer diffusion
140      !!
141      !! ** Method  :   set nldf from the namtra_ldf logicals
142      !!      nldf == -1   ESOPA test: ALL operators are used
143      !!      nldf ==  0   laplacian operator
144      !!      nldf ==  1   Rotated laplacian operator
145      !!      nldf ==  2   bilaplacian operator
146      !!      nldf ==  3   Rotated bilaplacian
147      !!----------------------------------------------------------------------
148      INTEGER ::   ioptio, ierr         ! temporary integers
149      !!----------------------------------------------------------------------
150
151      !  Define the lateral mixing oparator for tracers
152      ! ===============================================
153   
154      IF(lwp) THEN                    ! Namelist print
155         WRITE(numout,*)
156         WRITE(numout,*) 'tra_ldf_init : lateral tracer diffusive operator'
157         WRITE(numout,*) '~~~~~~~~~~~'
158         WRITE(numout,*) '   Namelist namtra_ldf already read in ldftra module'
159         WRITE(numout,*) '   see ldf_tra_init report for lateral mixing parameters'
160         WRITE(numout,*)
161         IF(lflush) CALL flush(numout)
162      ENDIF
163
164      !                               ! control the input
165      ioptio = 0
166      IF( ln_traldf_lap   )   ioptio = ioptio + 1
167      IF( ln_traldf_bilap )   ioptio = ioptio + 1
168      IF( ioptio >  1 )   CALL ctl_stop( '          use ONE or NONE of the 2 lap/bilap operator type on tracer' )
169      IF( ioptio == 0 )   nldf = -2   ! No lateral diffusion
170      ioptio = 0
171      IF( ln_traldf_level )   ioptio = ioptio + 1
172      IF( ln_traldf_hor   )   ioptio = ioptio + 1
173      IF( ln_traldf_iso   )   ioptio = ioptio + 1
174      IF( ioptio >  1 )   CALL ctl_stop( '          use only ONE direction (level/hor/iso)' )
175
176      ! defined the type of lateral diffusion from ln_traldf_... logicals
177      ! CAUTION : nldf = 1 is used in trazdf_imp, change it carefully
178      ierr = 0
179      IF( ln_traldf_lap ) THEN       ! laplacian operator
180         IF ( ln_zco ) THEN                ! z-coordinate
181            IF ( ln_traldf_level )   nldf = 0      ! iso-level  (no rotation)
182            IF ( ln_traldf_hor   )   nldf = 0      ! horizontal (no rotation)
183            IF ( ln_traldf_iso   )   nldf = 1      ! isoneutral (   rotation)
184         ENDIF
185         IF ( ln_zps ) THEN             ! zps-coordinate
186            IF ( ln_traldf_level )   ierr = 1      ! iso-level not allowed
187            IF ( ln_traldf_hor   )   nldf = 0      ! horizontal (no rotation)
188            IF ( ln_traldf_iso   )   nldf = 1      ! isoneutral (   rotation)
189         ENDIF
190         IF ( ln_sco ) THEN             ! s-coordinate
191            IF ( ln_traldf_level )   nldf = 0      ! iso-level  (no rotation)
192            IF ( ln_traldf_hor   )   nldf = 1      ! horizontal (   rotation)
193            IF ( ln_traldf_iso   )   nldf = 1      ! isoneutral (   rotation)
194         ENDIF
195      ENDIF
196
197      IF( ln_traldf_bilap ) THEN      ! bilaplacian operator
198         IF ( ln_zco ) THEN                ! z-coordinate
199            IF ( ln_traldf_level )   nldf = 2      ! iso-level  (no rotation)
200            IF ( ln_traldf_hor   )   nldf = 2      ! horizontal (no rotation)
201            IF ( ln_traldf_iso   )   ierr = 2      ! isoneutral (   rotation)
202         ENDIF
203         IF ( ln_zps ) THEN             ! zps-coordinate
204            IF ( ln_traldf_level )   ierr = 1      ! iso-level not allowed
205            IF ( ln_traldf_hor   )   nldf = 2      ! horizontal (no rotation)
206            IF ( ln_traldf_iso   )   ierr = 2      ! isoneutral (   rotation)
207         ENDIF
208         IF ( ln_sco ) THEN             ! s-coordinate
209            IF ( ln_traldf_level )   nldf = 2      ! iso-level  (no rotation)
210            IF ( ln_traldf_hor   )   nldf = 3      ! horizontal (   rotation)
211            IF ( ln_traldf_iso   )   ierr = 2      ! isoneutral (   rotation)
212         ENDIF
213      ENDIF
214
215      IF( nldf == 3 )   CALL ctl_warn( 'geopotential bilaplacian tracer diffusion in s-coords not thoroughly tested' )
216      IF( ierr == 1 )   CALL ctl_stop( ' iso-level in z-coordinate - partial step, not allowed' )
217      IF( ierr == 2 )   CALL ctl_stop( ' isoneutral bilaplacian operator does not exist' )
218      IF( ln_traldf_grif .AND. ln_isfcav         )   &
219           CALL ctl_stop( ' ice shelf and traldf_grif not tested')
220      IF( lk_traldf_eiv .AND. .NOT.ln_traldf_iso )   &
221           CALL ctl_stop( '          eddy induced velocity on tracers',   &
222           &              ' the eddy induced velocity on tracers requires isopycnal laplacian diffusion' )
223      IF( nldf == 1 .OR. nldf == 3 ) THEN      ! rotation
224         IF( .NOT.lk_ldfslp )   CALL ctl_stop( '          the rotation of the diffusive tensor require key_ldfslp' )
225         l_traldf_rot = .TRUE.                 ! needed for trazdf_imp
226      ENDIF
227
228      IF( lk_esopa ) THEN
229         IF(lwp) WRITE(numout,*) '          esopa control: use all lateral physics options'
230         IF(lwp .AND. lflush) CALL flush(numout)
231         nldf = -1
232      ENDIF
233
234      IF(lwp) THEN
235         WRITE(numout,*)
236         IF( nldf == -2 )   WRITE(numout,*) '          NO lateral diffusion'
237         IF( nldf == -1 )   WRITE(numout,*) '          ESOPA test All scheme used'
238         IF( nldf ==  0 )   WRITE(numout,*) '          laplacian operator'
239         IF( nldf ==  1 )   WRITE(numout,*) '          Rotated laplacian operator'
240         IF( nldf ==  2 )   WRITE(numout,*) '          bilaplacian operator'
241         IF( nldf ==  3 )   WRITE(numout,*) '          Rotated bilaplacian'
242         IF(lflush) CALL flush(numout)
243      ENDIF
244
245      ! Reference T & S diffusivity (if necessary)
246      ! ===========================
247      CALL ldf_ano
248      !
249   END SUBROUTINE tra_ldf_init
250
251#if defined key_traldf_ano
252   !!----------------------------------------------------------------------
253   !!   'key_traldf_ano'               T & S lateral diffusion on anomalies
254   !!----------------------------------------------------------------------
255
256   SUBROUTINE ldf_ano
257      !!----------------------------------------------------------------------
258      !!                  ***  ROUTINE ldf_ano  ***
259      !!
260      !! ** Purpose :   initializations of
261      !!----------------------------------------------------------------------
262      !
263      USE zdf_oce         ! vertical mixing
264      USE trazdf          ! vertical mixing: double diffusion
265      USE zdfddm          ! vertical mixing: double diffusion
266      !
267      INTEGER  ::   jk              ! Dummy loop indice
268      INTEGER  ::   ierr            ! local integer
269      LOGICAL  ::   llsave          ! local logical
270      REAL(wp) ::   zt0, zs0, z12   ! local scalar
271      REAL(wp), POINTER, DIMENSION(:,:,:) :: zt_ref, zs_ref, ztb, zsb, zavt     
272      !!----------------------------------------------------------------------
273      !
274      IF( nn_timing == 1 )  CALL timing_start('ldf_ano')
275      !
276      CALL wrk_alloc( jpi, jpj, jpk, zt_ref, zs_ref, ztb, zsb, zavt ) 
277      !
278
279      IF(lwp) THEN
280         WRITE(numout,*)
281         WRITE(numout,*) 'tra:ldf_ano : lateral diffusion acting on anomalies'
282         WRITE(numout,*) '~~~~~~~~~~~'
283         IF(lflush) CALL flush(numout)
284      ENDIF
285
286      !                              ! allocate trabbl arrays
287      ALLOCATE( t0_ldf(jpi,jpj,jpk) , s0_ldf(jpi,jpj,jpk) , STAT=ierr )
288      IF( lk_mpp    )   CALL mpp_sum( ierr )
289      IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'ldf_ano: unable to allocate arrays' )
290
291      ! defined the T & S reference profiles
292      ! ------------------------------------
293      zt0 =10.e0                               ! homogeneous ocean
294      zs0 =35.e0
295      zt_ref(:,:,:) = 10.0 * tmask(:,:,:)
296      zs_ref(:,:,:) = 35.0 * tmask(:,:,:)
297      IF(lwp) WRITE(numout,*) '              homogeneous ocean T = ', zt0, ' S = ',zs0
298      IF(lwp .AND. lflush) CALL flush(numout)
299
300      ! Initialisation of gtui/gtvi in case of no cavity
301      IF ( .NOT. ln_isfcav ) THEN
302         gtui(:,:,:) = 0.0_wp
303         gtvi(:,:,:) = 0.0_wp
304      END IF
305      !                                        ! T & S profile (to be coded +namelist parameter
306
307      ! prepare the ldf computation
308      ! ---------------------------
309      llsave = l_trdtra
310      l_trdtra = .false.      ! desactivate trend computation
311      t0_ldf(:,:,:) = 0.e0
312      s0_ldf(:,:,:) = 0.e0
313      ztb   (:,:,:) = tsb (:,:,:,jp_tem)
314      zsb   (:,:,:) = tsb (:,:,:,jp_sal)
315      ua    (:,:,:) = tsa (:,:,:,jp_tem)
316      va    (:,:,:) = tsa (:,:,:,jp_sal)
317      zavt  (:,:,:) = avt(:,:,:)
318      IF( lk_zdfddm ) THEN CALL ctl_stop( ' key_traldf_ano with key_zdfddm not implemented' )
319      ! set tb, sb to reference values and avr to zero
320      tsb (:,:,:,jp_tem) = zt_ref(:,:,:)
321      tsb (:,:,:,jp_sal) = zs_ref(:,:,:)
322      tsa (:,:,:,jp_tem) = 0.e0
323      tsa (:,:,:,jp_sal) = 0.e0
324      avt(:,:,:)         = 0.e0
325
326      ! Compute the ldf trends
327      ! ----------------------
328      CALL tra_ldf( nit000 + 1 )      ! horizontal components (+1: no more init)
329      CALL tra_zdf( nit000     )      ! vertical component (if necessary nit000 to performed the init)
330
331      ! finalise the computation and recover all arrays
332      ! -----------------------------------------------
333      l_trdtra = llsave
334      z12 = 2.e0
335      IF( neuler == 1)   z12 = 1.e0
336      IF( ln_zdfexp ) THEN      ! ta,sa are the trends
337         t0_ldf(:,:,:) = tsa(:,:,:,jp_tem)
338         s0_ldf(:,:,:) = tsa(:,:,:,jp_sal)
339      ELSE
340         DO jk = 1, jpkm1
341            t0_ldf(:,:,jk) = ( tsa(:,:,jk,jp_tem) - tsb(:,:,jk,jp_tem) ) / ( z12 *rdttra(jk) )
342            s0_ldf(:,:,jk) = ( tsa(:,:,jk,jp_sal) - tsb(:,:,jk,jp_sal) ) / ( z12 *rdttra(jk) )
343         END DO
344      ENDIF
345      tsb(:,:,:,jp_tem) = ztb (:,:,:)
346      tsb(:,:,:,jp_sal) = zsb (:,:,:)
347      tsa(:,:,:,jp_tem) = ua  (:,:,:)
348      tsa(:,:,:,jp_sal) = va  (:,:,:)
349      avt(:,:,:)        = zavt(:,:,:)
350      !
351      CALL wrk_dealloc( jpi, jpj, jpk, zt_ref, zs_ref, ztb, zsb, zavt ) 
352      !
353      IF( nn_timing == 1 )  CALL timing_stop('ldf_ano')
354      !
355   END SUBROUTINE ldf_ano
356
357#else
358   !!----------------------------------------------------------------------
359   !!   default option :   Dummy code   NO T & S background profiles
360   !!----------------------------------------------------------------------
361   SUBROUTINE ldf_ano
362      IF(lwp) THEN
363         WRITE(numout,*)
364         WRITE(numout,*) 'tra:ldf_ano : lateral diffusion acting on the full fields'
365         WRITE(numout,*) '~~~~~~~~~~~'
366         IF(lflush) CALL flush(numout)
367      ENDIF
368   END SUBROUTINE ldf_ano
369#endif
370
371   !!======================================================================
372END MODULE traldf
Note: See TracBrowser for help on using the repository browser.