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.
trcdta.F90 in branches/UKMO/dev_r5518_GO6_under_ice_relax_dr_hook/NEMOGCM/NEMO/TOP_SRC – NEMO

source: branches/UKMO/dev_r5518_GO6_under_ice_relax_dr_hook/NEMOGCM/NEMO/TOP_SRC/trcdta.F90

Last change on this file was 11738, checked in by marc, 5 years ago

The Dr Hook changes from my perl code.
File size: 15.0 KB
Line 
1MODULE trcdta
2   !!======================================================================
3   !!                     ***  MODULE  trcdta  ***
4   !! TOP :  reads passive tracer data
5   !!=====================================================================
6   !! History :   1.0  !  2002-04  (O. Aumont)  original code
7   !!              -   !  2004-03  (C. Ethe)  module
8   !!              -   !  2005-03  (O. Aumont, A. El Moussaoui) F90
9   !!            3.4   !  2010-11  (C. Ethe, G. Madec)  use of fldread + dynamical allocation
10   !!            3.5   !  2013-08  (M. Vichi)  generalization for other BGC models
11   !!----------------------------------------------------------------------
12#if  defined key_top 
13   !!----------------------------------------------------------------------
14   !!   'key_top'                                                TOP model
15   !!----------------------------------------------------------------------
16   !!   trc_dta    : read and time interpolated passive tracer data
17   !!----------------------------------------------------------------------
18   USE par_trc       !  passive tracers parameters
19   USE oce_trc       !  shared variables between ocean and passive tracers
20   USE trc           !  passive tracers common variables
21   USE iom           !  I/O manager
22   USE lib_mpp       !  MPP library
23   USE fldread       !  read input fields
24
25   USE yomhook, ONLY: lhook, dr_hook
26   USE parkind1, ONLY: jprb, jpim
27
28   IMPLICIT NONE
29   PRIVATE
30
31   PUBLIC   trc_dta         ! called in trcini.F90 and trcdmp.F90
32   PUBLIC   trc_dta_init    ! called in trcini.F90
33
34   INTEGER  , SAVE, PUBLIC                             :: nb_trcdta   ! number of tracers to be initialised with data
35   INTEGER  , SAVE, PUBLIC, ALLOCATABLE, DIMENSION(:)  :: n_trc_index ! indice of tracer which is initialised with data
36   INTEGER  , SAVE, PUBLIC                             :: ntra        ! MAX( 1, nb_trcdta ) to avoid compilation error with bounds checking
37   REAL(wp) , SAVE, PUBLIC, ALLOCATABLE, DIMENSION(:)  :: rf_trfac    ! multiplicative factor for tracer values
38!$AGRIF_DO_NOT_TREAT
39   TYPE(FLD), SAVE, PUBLIC, ALLOCATABLE, DIMENSION(:)  :: sf_trcdta   ! structure of input SST (file informations, fields read)
40!$AGRIF_END_DO_NOT_TREAT
41   !! * Substitutions
42#  include "domzgr_substitute.h90"
43   !!----------------------------------------------------------------------
44   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
45   !! $Id$
46   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
47   !!----------------------------------------------------------------------
48CONTAINS
49
50   SUBROUTINE trc_dta_init(ntrc)
51      !!----------------------------------------------------------------------
52      !!                   ***  ROUTINE trc_dta_init  ***
53      !!                   
54      !! ** Purpose :   initialisation of passive tracer input data
55      !!
56      !! ** Method  : - Read namtsd namelist
57      !!              - allocates passive tracer data structure
58      !!----------------------------------------------------------------------
59      !
60      INTEGER,INTENT(IN) :: ntrc                             ! number of tracers
61      INTEGER            :: jl, jn                           ! dummy loop indices
62      INTEGER            :: ierr0, ierr1, ierr2, ierr3       ! temporary integers
63      INTEGER            :: ios                              ! Local integer output status for namelist read
64      CHARACTER(len=100) :: clndta, clntrc
65      REAL(wp)           :: zfact
66      !
67      CHARACTER(len=100)            :: cn_dir
68      TYPE(FLD_N), ALLOCATABLE, DIMENSION(:) :: slf_i ! array of namelist informations on the fields to read
69      TYPE(FLD_N), DIMENSION(jpmaxtrc) :: sn_trcdta
70      REAL(wp)   , DIMENSION(jpmaxtrc) :: rn_trfac    ! multiplicative factor for tracer values
71      !!
72      NAMELIST/namtrc_dta/ sn_trcdta, cn_dir, rn_trfac 
73      INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0
74      INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
75      REAL(KIND=jprb)               :: zhook_handle
76
77      CHARACTER(LEN=*), PARAMETER :: RoutineName='TRC_DTA_INIT'
78
79      IF (lhook) CALL dr_hook(RoutineName,zhook_in,zhook_handle)
80
81      !!----------------------------------------------------------------------
82      !
83      IF( nn_timing == 1 )  CALL timing_start('trc_dta_init')
84      !
85      !  Initialisation
86      ierr0 = 0  ;  ierr1 = 0  ;  ierr2 = 0  ;  ierr3 = 0 
87      ! Compute the number of tracers to be initialised with data
88      ALLOCATE( n_trc_index(ntrc), slf_i(ntrc), STAT=ierr0 )
89      IF( ierr0 > 0 ) THEN
90         CALL ctl_stop( 'trc_dta_init: unable to allocate n_trc_index' )   ;   RETURN
91      ENDIF
92      nb_trcdta      = 0
93      n_trc_index(:) = 0
94      DO jn = 1, ntrc
95         IF( ln_trc_ini(jn) ) THEN
96             nb_trcdta       = nb_trcdta + 1 
97             n_trc_index(jn) = nb_trcdta 
98         ENDIF
99      ENDDO
100      !
101      ntra = MAX( 1, nb_trcdta )   ! To avoid compilation error with bounds checking
102      IF(lwp) THEN
103         WRITE(numout,*) ' '
104         WRITE(numout,*) 'trc_dta_init : Passive tracers Initial Conditions '
105         WRITE(numout,*) '~~~~~~~~~~~~~~ '
106         WRITE(numout,*) ' number of passive tracers to be initialize by data :', ntra
107         WRITE(numout,*) ' '
108      ENDIF
109      !
110      REWIND( numnat_ref )              ! Namelist namtrc_dta in reference namelist : Passive tracer input data
111      READ  ( numnat_ref, namtrc_dta, IOSTAT = ios, ERR = 901)
112901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtrc_dta in reference namelist', lwp )
113
114      REWIND( numnat_cfg )              ! Namelist namtrc_dta in configuration namelist : Passive tracer input data
115      READ  ( numnat_cfg, namtrc_dta, IOSTAT = ios, ERR = 902 )
116902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtrc_dta in configuration namelist', lwp )
117      IF(lwm) WRITE ( numont, namtrc_dta )
118
119      IF( lwp ) THEN
120         DO jn = 1, ntrc
121            IF( ln_trc_ini(jn) )  THEN    ! open input file only if ln_trc_ini(jn) is true
122               clndta = TRIM( sn_trcdta(jn)%clvar )
123               if (jn > jptra) then
124                  clntrc='Dummy' ! By pass weird formats in ocean.output if ntrc > jptra
125               else
126                  clntrc = TRIM( ctrcnm   (jn)       )
127               endif
128               zfact  = rn_trfac(jn)
129               IF( clndta /=  clntrc ) THEN
130                  CALL ctl_warn( 'trc_dta_init: passive tracer data initialisation    ',   &
131                  &              'Input name of data file : '//TRIM(clndta)//   &
132                  &              ' differs from that of tracer : '//TRIM(clntrc)//' ')
133               ENDIF
134               WRITE(numout,'(a, i4,3a,e11.3)') ' Read IC file for tracer number :', &
135               &            jn, ', name : ', TRIM(clndta), ', Multiplicative Scaling factor : ', zfact
136            ENDIF
137         END DO
138      ENDIF
139      !
140      IF( nb_trcdta > 0 ) THEN       !  allocate only if the number of tracer to initialise is greater than zero
141         ALLOCATE( sf_trcdta(nb_trcdta), rf_trfac(nb_trcdta), STAT=ierr1 )
142         IF( ierr1 > 0 ) THEN
143            CALL ctl_stop( 'trc_dta_init: unable to allocate  sf_trcdta structure' )   ;   RETURN
144         ENDIF
145         !
146         DO jn = 1, ntrc
147            IF( ln_trc_ini(jn) ) THEN      ! update passive tracers arrays with input data read from file
148               jl = n_trc_index(jn)
149               slf_i(jl)    = sn_trcdta(jn)
150               rf_trfac(jl) = rn_trfac(jn)
151                                            ALLOCATE( sf_trcdta(jl)%fnow(jpi,jpj,jpk)   , STAT=ierr2 )
152               IF( sn_trcdta(jn)%ln_tint )  ALLOCATE( sf_trcdta(jl)%fdta(jpi,jpj,jpk,2) , STAT=ierr3 )
153               IF( ierr2 + ierr3 > 0 ) THEN
154                 CALL ctl_stop( 'trc_dta_init : unable to allocate passive tracer data arrays' )   ;   RETURN
155               ENDIF
156            ENDIF
157            !   
158         ENDDO
159         !                         ! fill sf_trcdta with slf_i and control print
160         CALL fld_fill( sf_trcdta, slf_i, cn_dir, 'trc_dta_init', 'Passive tracer data', 'namtrc' )
161         !
162      ENDIF
163      !
164      DEALLOCATE( slf_i )          ! deallocate local field structure
165      IF( nn_timing == 1 )  CALL timing_stop('trc_dta_init')
166      !
167      IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
168   END SUBROUTINE trc_dta_init
169
170
171   SUBROUTINE trc_dta( kt, sf_dta, ptrfac, ptrc)
172      !!----------------------------------------------------------------------
173      !!                   ***  ROUTINE trc_dta  ***
174      !!                   
175      !! ** Purpose :   provides passive tracer data at kt
176      !!
177      !! ** Method  : - call fldread routine
178      !!              - s- or mixed z-s coordinate: vertical interpolation on model mesh
179      !!              - ln_trcdmp=F: deallocates the data structure as they are not used
180      !!
181      !! ** Action  :   sf_dta   passive tracer data on medl mesh and interpolated at time-step kt
182      !!----------------------------------------------------------------------
183      INTEGER                     , INTENT(in   ) ::   kt     ! ocean time-step
184      TYPE(FLD), DIMENSION(1)     , INTENT(inout) ::   sf_dta     ! array of information on the field to read
185      REAL(wp)                    , INTENT(in   ) ::   ptrfac  ! multiplication factor
186      REAL(wp), DIMENSION(jpi,jpj,jpk), OPTIONAL  , INTENT(out  ) ::   ptrc
187      !
188      INTEGER ::   ji, jj, jk, jl, jkk, ik    ! dummy loop indices
189      REAL(wp)::   zl, zi
190      REAL(wp), DIMENSION(jpk) ::  ztp                ! 1D workspace
191      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrcdta   ! 3D  workspace
192      CHARACTER(len=100) :: clndta
193      INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0
194      INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
195      REAL(KIND=jprb)               :: zhook_handle
196
197      CHARACTER(LEN=*), PARAMETER :: RoutineName='TRC_DTA'
198
199      IF (lhook) CALL dr_hook(RoutineName,zhook_in,zhook_handle)
200
201      !!----------------------------------------------------------------------
202      !
203      IF( nn_timing == 1 )  CALL timing_start('trc_dta')
204      !
205      IF( nb_trcdta > 0 ) THEN
206         !
207         CALL wrk_alloc( jpi, jpj, jpk, ztrcdta )    ! Memory allocation
208         !
209         CALL fld_read( kt, 1, sf_dta )      !==   read data at kt time step   ==!
210         ztrcdta(:,:,:) = sf_dta(1)%fnow(:,:,:) * tmask(:,:,:)    ! Mask
211         !
212         IF( ln_sco ) THEN                   !==   s- or mixed s-zps-coordinate   ==!
213            !
214            IF( kt == nit000 .AND. lwp )THEN
215               WRITE(numout,*)
216               WRITE(numout,*) 'trc_dta: interpolates passive tracer data onto the s- or mixed s-z-coordinate mesh'
217            ENDIF
218            !
219            DO jj = 1, jpj                         ! vertical interpolation of T & S
220               DO ji = 1, jpi
221                  DO jk = 1, jpk                        ! determines the intepolated T-S profiles at each (i,j) points
222                     zl = fsdept_n(ji,jj,jk)
223                     IF(     zl < gdept_1d(1  ) ) THEN         ! above the first level of data
224                        ztp(jk) = ztrcdta(ji,jj,1)
225                     ELSEIF( zl > gdept_1d(jpk) ) THEN         ! below the last level of data
226                        ztp(jk) =  ztrcdta(ji,jj,jpkm1)
227                     ELSE                                      ! inbetween : vertical interpolation between jkk & jkk+1
228                        DO jkk = 1, jpkm1                                  ! when  gdept(jkk) < zl < gdept(jkk+1)
229                           IF( (zl-gdept_1d(jkk)) * (zl-gdept_1d(jkk+1)) <= 0._wp ) THEN
230                              zi = ( zl - gdept_1d(jkk) ) / (gdept_1d(jkk+1)-gdept_1d(jkk))
231                              ztp(jk) = ztrcdta(ji,jj,jkk) + ( ztrcdta(ji,jj,jkk+1) - &
232                                        ztrcdta(ji,jj,jkk) ) * zi 
233                           ENDIF
234                        END DO
235                     ENDIF
236                  END DO
237                  DO jk = 1, jpkm1
238                    ztrcdta(ji,jj,jk) = ztp(jk) * tmask(ji,jj,jk)     ! mask required for mixed zps-s-coord
239                  END DO
240                  ztrcdta(ji,jj,jpk) = 0._wp
241                END DO
242            END DO
243            !
244         ELSE                                !==   z- or zps- coordinate   ==!
245            !
246            IF( ln_zps ) THEN                      ! zps-coordinate (partial steps) interpolation at the last ocean level
247               DO jj = 1, jpj
248                  DO ji = 1, jpi
249                     ik = mbkt(ji,jj) 
250                     IF( ik > 1 ) THEN
251                        zl = ( gdept_1d(ik) - fsdept_n(ji,jj,ik) ) / ( gdept_1d(ik) - gdept_1d(ik-1) )
252                        ztrcdta(ji,jj,ik) = (1.-zl) * ztrcdta(ji,jj,ik) + zl * ztrcdta(ji,jj,ik-1)
253                     ENDIF
254                     ik = mikt(ji,jj)
255                     IF( ik > 1 ) THEN
256                        zl = ( fsdept_n(ji,jj,ik) - gdept_1d(ik) ) / ( gdept_1d(ik+1) - gdept_1d(ik) )
257                        ztrcdta(ji,jj,ik) = (1.-zl) * ztrcdta(ji,jj,ik) + zl * ztrcdta(ji,jj,ik+1)
258                     ENDIF
259                  END DO
260               END DO
261            ENDIF
262            !
263         ENDIF
264         !
265         ! Add multiplicative factor
266         ztrcdta(:,:,:) = ztrcdta(:,:,:) * ptrfac
267         !
268         ! Data structure for trc_ini (and BFMv5.1 coupling)
269         IF( .NOT. PRESENT(ptrc) ) sf_dta(1)%fnow(:,:,:) = ztrcdta(:,:,:)
270         !
271         ! Data structure for trc_dmp
272         IF( PRESENT(ptrc) )  ptrc(:,:,:) = ztrcdta(:,:,:)
273         !
274         IF( lwp .AND. kt == nit000 ) THEN
275               clndta = TRIM( sf_dta(1)%clvar ) 
276               WRITE(numout,*) ''//clndta//' data '
277               WRITE(numout,*)
278               WRITE(numout,*)'  level = 1'
279               CALL prihre( ztrcdta(:,:,1), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1., numout )
280               WRITE(numout,*)'  level = ', jpk/2
281               CALL prihre( ztrcdta(:,:,jpk/2), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1., numout )
282               WRITE(numout,*)'  level = ', jpkm1
283               CALL prihre( ztrcdta(:,:,jpkm1), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1., numout )
284               WRITE(numout,*)
285         ENDIF
286         !
287         CALL wrk_dealloc( jpi, jpj, jpk, ztrcdta )
288         !
289      ENDIF
290      !
291      IF( nn_timing == 1 )  CALL timing_stop('trc_dta')
292      !
293      IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
294   END SUBROUTINE trc_dta
295#else
296   !!----------------------------------------------------------------------
297   !!   Dummy module                              NO 3D passive tracer data
298   !!----------------------------------------------------------------------
299CONTAINS
300   SUBROUTINE trc_dta( kt, sf_dta, ptrfac, ptrc)        ! Empty routine
301   INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0
302   INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
303   REAL(KIND=jprb)               :: zhook_handle
304
305   CHARACTER(LEN=*), PARAMETER :: RoutineName='TRC_DTA'
306
307   IF (lhook) CALL dr_hook(RoutineName,zhook_in,zhook_handle)
308
309      WRITE(*,*) 'trc_dta: You should not have seen this print! error?', kt
310   IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
311   END SUBROUTINE trc_dta
312#endif
313   !!======================================================================
314END MODULE trcdta
Note: See TracBrowser for help on using the repository browser.