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.
dtatsd.F90 in branches/UKMO/r8395_India_uncoupled/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

source: branches/UKMO/r8395_India_uncoupled/NEMOGCM/NEMO/OPA_SRC/DOM/dtatsd.F90 @ 10812

Last change on this file since 10812 was 10812, checked in by jcastill, 5 years ago

Add again the changes that were reverted before: they are needed after all
File size: 16.2 KB
Line 
1MODULE dtatsd
2   !!======================================================================
3   !!                     ***  MODULE  dtatsd  ***
4   !! Ocean data  :  read ocean Temperature & Salinity Data from gridded data
5   !!======================================================================
6   !! History :  OPA  ! 1991-03  ()  Original code
7   !!             -   ! 1992-07  (M. Imbard)
8   !!            8.0  ! 1999-10  (M.A. Foujols, M. Imbard)  NetCDF FORMAT
9   !!   NEMO     1.0  ! 2002-06  (G. Madec)  F90: Free form and module
10   !!            3.3  ! 2010-10  (C. Bricaud, S. Masson)  use of fldread
11   !!            3.4  ! 2010-11  (G. Madec, C. Ethe) Merge of dtatem and dtasal + suppression of CPP keys
12   !!----------------------------------------------------------------------
13
14   !!----------------------------------------------------------------------
15   !!   dta_tsd      : read and time interpolated ocean Temperature & Salinity Data
16   !!----------------------------------------------------------------------
17   USE oce             ! ocean dynamics and tracers
18   USE dom_oce         ! ocean space and time domain
19   USE fldread         ! read input fields
20   USE in_out_manager  ! I/O manager
21   USE phycst          ! physical constants
22   USE lib_mpp         ! MPP library
23   USE wrk_nemo        ! Memory allocation
24   USE timing          ! Timing
25   USE iom
26
27   IMPLICIT NONE
28   PRIVATE
29
30   PUBLIC   dta_tsd_init   ! called by opa.F90
31   PUBLIC   dta_tsd        ! called by istate.F90 and tradmp.90
32
33   LOGICAL , PUBLIC ::   ln_tsd_init      !: T & S data flag
34   LOGICAL , PUBLIC ::   ln_tsd_interp    !: vertical interpolation flag
35   LOGICAL , PUBLIC ::   ln_tsd_tradmp    !: internal damping toward input data flag
36
37   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_tsd   ! structure of input SST (file informations, fields read)
38   INTEGER                                 ::   jpk_init , inum_dta
39   INTEGER                                 ::   id ,linum   ! local integers
40   INTEGER                                 ::   zdim(4)
41
42   !!----------------------------------------------------------------------
43   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
44   !! $Id$
45   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
46   !!----------------------------------------------------------------------
47CONTAINS
48
49   SUBROUTINE dta_tsd_init( ld_tradmp )
50      !!----------------------------------------------------------------------
51      !!                   ***  ROUTINE dta_tsd_init  ***
52      !!                   
53      !! ** Purpose :   initialisation of T & S input data
54      !!
55      !! ** Method  : - Read namtsd namelist
56      !!              - allocates T & S data structure
57      !!----------------------------------------------------------------------
58      LOGICAL, INTENT(in), OPTIONAL ::   ld_tradmp   ! force the initialization when tradp is used
59      !
60      INTEGER ::   ios, ierr0, ierr1, ierr2, ierr3, ierr4, ierr5   ! local integers
61      !!
62      CHARACTER(len=100)            ::   cn_dir          ! Root directory for location of ssr files
63      TYPE(FLD_N), DIMENSION(jpts+2)::   slf_i           ! array of namelist informations on the fields to read
64      TYPE(FLD_N)                   ::   sn_tem, sn_sal, sn_dep, sn_msk
65     
66      !!
67      NAMELIST/namtsd/   ln_tsd_init, ln_tsd_interp, ln_tsd_tradmp, cn_dir, sn_tem, sn_sal, sn_dep, sn_msk
68      !!----------------------------------------------------------------------
69      !
70      IF( nn_timing == 1 )  CALL timing_start('dta_tsd_init')
71      !
72      !  Initialisation
73      ierr0 = 0  ;  ierr1 = 0  ;  ierr2 = 0  ;  ierr3 = 0  ; ierr4 = 0  ;  ierr5 = 0 
74      !
75      REWIND( numnam_ref )              ! Namelist namtsd in reference namelist :
76      READ  ( numnam_ref, namtsd, IOSTAT = ios, ERR = 901)
77901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtsd in reference namelist', lwp )
78
79      REWIND( numnam_cfg )              ! Namelist namtsd in configuration namelist : Parameters of the run
80      READ  ( numnam_cfg, namtsd, IOSTAT = ios, ERR = 902 )
81902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtsd in configuration namelist', lwp )
82      IF(lwm) WRITE ( numond, namtsd )
83
84      IF( PRESENT( ld_tradmp ) )   ln_tsd_tradmp = .TRUE.     ! forces the initialization when tradmp is used
85     
86      IF(lwp) THEN                  ! control print
87         WRITE(numout,*)
88         WRITE(numout,*) 'dta_tsd_init : Temperature & Salinity data '
89         WRITE(numout,*) '~~~~~~~~~~~~ '
90         WRITE(numout,*) '   Namelist namtsd'
91         WRITE(numout,*) '      Initialisation of ocean T & S with T &S input data   ln_tsd_init   = ', ln_tsd_init
92         WRITE(numout,*) '      iInterpolation of initial conditions in the vertical ln_tsd_interp = ', ln_tsd_interp
93         WRITE(numout,*) '      damping of ocean T & S toward T &S input data        ln_tsd_tradmp = ', ln_tsd_tradmp
94         WRITE(numout,*)
95         IF( .NOT.ln_tsd_init .AND. .NOT.ln_tsd_tradmp ) THEN
96            WRITE(numout,*)
97            WRITE(numout,*) '   T & S data not used'
98         ENDIF
99      ENDIF
100      !
101      IF( ln_rstart .AND. ln_tsd_init ) THEN
102         CALL ctl_warn( 'dta_tsd_init: ocean restart and T & S data intialisation, ',   &
103            &           'we keep the restart T & S values and set ln_tsd_init to FALSE' )
104         ln_tsd_init = .FALSE.
105      ENDIF
106      IF( ln_tsd_interp .AND. ln_tsd_tradmp ) THEN
107            CALL ctl_stop( 'dta_tsd_init: Tracer damping and vertical interpolation not yet configured' )   ;   RETURN
108      ENDIF
109      IF( ln_tsd_interp .AND. LEN(TRIM(sn_msk%wname)) > 0 ) THEN
110            CALL ctl_stop( 'dta_tsd_init: Using vertical interpolation and weights files not recommended' )   ;   RETURN
111      ENDIF
112      !
113      !                             ! allocate the arrays (if necessary)
114      IF(  ln_tsd_init .OR. ln_tsd_tradmp  ) THEN
115         !
116         IF( ln_tsd_interp ) THEN
117           ALLOCATE( sf_tsd(jpts+2), STAT=ierr0 ) ! to carry the addtional depth information
118         ELSE
119           ALLOCATE( sf_tsd(jpts  ), STAT=ierr0 ) 
120         ENDIF
121         IF( ierr0 > 0 ) THEN
122            CALL ctl_stop( 'dta_tsd_init: unable to allocate sf_tsd structure' )   ;   RETURN
123         ENDIF
124         !
125         IF( ln_tsd_interp ) THEN
126            CALL iom_open ( trim(cn_dir) // trim(sn_dep%clname), inum_dta ) 
127            id = iom_varid( inum_dta, sn_dep%clvar, zdim )
128            jpk_init = zdim(3)
129            IF(lwp) WRITE(numout,*) 'Dimension of veritcal coordinate in ICs: ', jpk_init
130            CALL iom_close( inum_dta )   ! Close the input file
131            !
132                                 ALLOCATE( sf_tsd(jp_tem)%fnow(jpi,jpj,jpk_init  ) , STAT=ierr0 )
133            IF( sn_tem%ln_tint ) ALLOCATE( sf_tsd(jp_tem)%fdta(jpi,jpj,jpk_init,2) , STAT=ierr1 )
134                                 ALLOCATE( sf_tsd(jp_sal)%fnow(jpi,jpj,jpk_init  ) , STAT=ierr2 )
135            IF( sn_sal%ln_tint ) ALLOCATE( sf_tsd(jp_sal)%fdta(jpi,jpj,jpk_init,2) , STAT=ierr3 ) 
136                                 ALLOCATE( sf_tsd(jp_dep)%fnow(jpi,jpj,jpk_init  ) , STAT=ierr4 )
137                                 ALLOCATE( sf_tsd(jp_msk)%fnow(jpi,jpj,jpk_init  ) , STAT=ierr5 )
138         ELSE
139                                 ALLOCATE( sf_tsd(jp_tem)%fnow(jpi,jpj,jpk)   , STAT=ierr0 )
140            IF( sn_tem%ln_tint ) ALLOCATE( sf_tsd(jp_tem)%fdta(jpi,jpj,jpk,2) , STAT=ierr1 )
141                                 ALLOCATE( sf_tsd(jp_sal)%fnow(jpi,jpj,jpk)   , STAT=ierr2 )
142            IF( sn_sal%ln_tint ) ALLOCATE( sf_tsd(jp_sal)%fdta(jpi,jpj,jpk,2) , STAT=ierr3 ) 
143         ENDIF ! ln_tsd_interp
144
145         !
146         IF( ierr0 + ierr1 + ierr2 + ierr3 + ierr4 + ierr5 > 0 ) THEN
147            CALL ctl_stop( 'dta_tsd : unable to allocate T & S data arrays' )   ;   RETURN
148         ENDIF
149         !                         ! fill sf_tsd with sn_tem & sn_sal and control print
150         slf_i(jp_tem) = sn_tem   ;   slf_i(jp_sal) = sn_sal
151         IF( ln_tsd_interp ) slf_i(jp_dep) = sn_dep   ;   slf_i(jp_msk) = sn_msk
152         CALL fld_fill( sf_tsd, slf_i, cn_dir, 'dta_tsd', 'Temperature & Salinity data', 'namtsd', no_print )
153         !
154      ENDIF
155      !
156      IF( nn_timing == 1 )  CALL timing_stop('dta_tsd_init')
157      !
158   END SUBROUTINE dta_tsd_init
159
160
161   SUBROUTINE dta_tsd( kt, ptsd )
162      !!----------------------------------------------------------------------
163      !!                   ***  ROUTINE dta_tsd  ***
164      !!                   
165      !! ** Purpose :   provides T and S data at kt
166      !!
167      !! ** Method  : - call fldread routine
168      !!              - ORCA_R2: add some hand made alteration to read data 
169      !!              - 'key_orca_lev10' interpolates on 10 times more levels
170      !!              - s- or mixed z-s coordinate: vertical interpolation on model mesh
171      !!              - ln_tsd_tradmp=F: deallocates the T-S data structure
172      !!                as T-S data are no are used
173      !!
174      !! ** Action  :   ptsd   T-S data on medl mesh and interpolated at time-step kt
175      !!----------------------------------------------------------------------
176      INTEGER                              , INTENT(in   ) ::   kt     ! ocean time-step
177      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(  out) ::   ptsd   ! T & S data
178      !
179      INTEGER ::   ji, jj, jk, jl, jk_init   ! dummy loop indicies
180      INTEGER ::   ik, il0, il1, ii0, ii1, ij0, ij1        ! local integers
181      REAL(wp)::   zl, zi
182      !!----------------------------------------------------------------------
183      !
184      IF( nn_timing == 1 )  CALL timing_start('dta_tsd')
185      !
186      CALL fld_read( kt, 1, sf_tsd )      !==   read T & S data at kt time step   ==!
187      !
188      !
189!!gm  This should be removed from the code   ===>>>>  T & S files has to be changed
190      !
191      !                                   !==   ORCA_R2 configuration and T & S damping   ==!
192      IF( cn_cfg == "orca" .AND. nn_cfg == 2 .AND. ln_tsd_tradmp ) THEN    ! some hand made alterations
193         !
194         ij0 = 101   ;   ij1 = 109                       ! Reduced T & S in the Alboran Sea
195         ii0 = 141   ;   ii1 = 155
196         DO jj = mj0(ij0), mj1(ij1)
197            DO ji = mi0(ii0), mi1(ii1)
198               sf_tsd(jp_tem)%fnow(ji,jj,13:13) = sf_tsd(jp_tem)%fnow(ji,jj,13:13) - 0.20_wp
199               sf_tsd(jp_tem)%fnow(ji,jj,14:15) = sf_tsd(jp_tem)%fnow(ji,jj,14:15) - 0.35_wp
200               sf_tsd(jp_tem)%fnow(ji,jj,16:25) = sf_tsd(jp_tem)%fnow(ji,jj,16:25) - 0.40_wp
201               !
202               sf_tsd(jp_sal)%fnow(ji,jj,13:13) = sf_tsd(jp_sal)%fnow(ji,jj,13:13) - 0.15_wp
203               sf_tsd(jp_sal)%fnow(ji,jj,14:15) = sf_tsd(jp_sal)%fnow(ji,jj,14:15) - 0.25_wp
204               sf_tsd(jp_sal)%fnow(ji,jj,16:17) = sf_tsd(jp_sal)%fnow(ji,jj,16:17) - 0.30_wp
205               sf_tsd(jp_sal)%fnow(ji,jj,18:25) = sf_tsd(jp_sal)%fnow(ji,jj,18:25) - 0.35_wp
206            END DO
207         END DO
208         ij0 =  87   ;   ij1 =  96                          ! Reduced temperature in Red Sea
209         ii0 = 148   ;   ii1 = 160
210         sf_tsd(jp_tem)%fnow( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ,  4:10 ) = 7.0_wp
211         sf_tsd(jp_tem)%fnow( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 11:13 ) = 6.5_wp
212         sf_tsd(jp_tem)%fnow( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 14:20 ) = 6.0_wp
213      ENDIF
214!!gm end
215      !
216      IF( kt == nit000 .AND. lwp )THEN
217         WRITE(numout,*)
218         WRITE(numout,*) 'dta_tsd: interpolates T & S data onto current mesh'
219      ENDIF
220      !
221      IF( ln_tsd_interp ) THEN                 ! probably should use pointers in the following to make more readable
222      !
223         DO jk = 1, jpk                        ! determines the intepolated T-S profiles at each (i,j) points
224            DO jj= 1, jpj
225               DO ji= 1, jpi
226                  zl = gdept_0(ji,jj,jk)
227                  IF( zl < sf_tsd(jp_dep)%fnow(ji,jj,1) ) THEN                     ! above the first level of data
228                     ptsd(ji,jj,jk,jp_tem) = sf_tsd(jp_tem)%fnow(ji,jj,1) 
229                     ptsd(ji,jj,jk,jp_sal) = sf_tsd(jp_sal)%fnow(ji,jj,1)
230                  ELSEIF( zl > sf_tsd(jp_dep)%fnow(ji,jj,jpk_init) ) THEN          ! below the last level of data
231                     ptsd(ji,jj,jk,jp_tem) = sf_tsd(jp_tem)%fnow(ji,jj,jpk_init)
232                     ptsd(ji,jj,jk,jp_sal) = sf_tsd(jp_sal)%fnow(ji,jj,jpk_init)
233                  ELSE                                                             ! inbetween : vertical interpolation between jk_init & jk_init+1
234                     DO jk_init = 1, jpk_init-1                                    ! when  gdept(jk_init) < zl < gdept(jk_init+1)
235                        IF( sf_tsd(jp_msk)%fnow(ji,jj,jk_init+1) == 0 ) THEN       ! if there is no data fill down
236                           sf_tsd(jp_tem)%fnow(ji,jj,jk_init+1) = sf_tsd(jp_tem)%fnow(ji,jj,jk_init)
237                           sf_tsd(jp_sal)%fnow(ji,jj,jk_init+1) = sf_tsd(jp_sal)%fnow(ji,jj,jk_init)
238                        ENDIF
239                        IF( (zl-sf_tsd(jp_dep)%fnow(ji,jj,jk_init)) * (zl-sf_tsd(jp_dep)%fnow(ji,jj,jk_init+1)) <= 0._wp ) THEN
240                           zi = ( zl - sf_tsd(jp_dep)%fnow(ji,jj,jk_init) ) / &
241                        &       (sf_tsd(jp_dep)%fnow(ji,jj,jk_init+1)-sf_tsd(jp_dep)%fnow(ji,jj,jk_init))
242                           ptsd(ji,jj,jk,jp_tem) = sf_tsd(jp_tem)%fnow(ji,jj,jk_init) + &
243                        &                          (sf_tsd(jp_tem)%fnow(ji,jj,jk_init+1)-sf_tsd(jp_tem)%fnow(ji,jj,jk_init)) * zi
244                           ptsd(ji,jj,jk,jp_sal) = sf_tsd(jp_sal)%fnow(ji,jj,jk_init) + &
245                        &                          (sf_tsd(jp_sal)%fnow(ji,jj,jk_init+1)-sf_tsd(jp_sal)%fnow(ji,jj,jk_init)) * zi
246                        ENDIF
247                     END DO
248                  ENDIF
249               ENDDO
250            ENDDO
251         END DO
252         !
253         ptsd(:,:,:,jp_tem) = ptsd(:,:,:,jp_tem) *tmask(:,:,:)
254         ptsd(:,:,:,jp_sal) = ptsd(:,:,:,jp_sal) *tmask(:,:,:)
255      ELSE                                !==   z- or zps- coordinate   ==!
256         !                             
257         ptsd(:,:,:,jp_tem) = sf_tsd(jp_tem)%fnow(:,:,:)  * tmask(:,:,:)  ! Mask
258         ptsd(:,:,:,jp_sal) = sf_tsd(jp_sal)%fnow(:,:,:)  * tmask(:,:,:)
259         !
260         IF( ln_zps ) THEN                      ! zps-coordinate (partial steps) interpolation at the last ocean level
261            DO jj = 1, jpj
262               DO ji = 1, jpi
263                  ik = mbkt(ji,jj) 
264                  IF( ik > 1 ) THEN
265                     zl = ( gdept_1d(ik) - gdept_0(ji,jj,ik) ) / ( gdept_1d(ik) - gdept_1d(ik-1) )
266                     ptsd(ji,jj,ik,jp_tem) = (1.-zl) * ptsd(ji,jj,ik,jp_tem) + zl * ptsd(ji,jj,ik-1,jp_tem)
267                     ptsd(ji,jj,ik,jp_sal) = (1.-zl) * ptsd(ji,jj,ik,jp_sal) + zl * ptsd(ji,jj,ik-1,jp_sal)
268                  ENDIF
269                  ik = mikt(ji,jj)
270                  IF( ik > 1 ) THEN
271                     zl = ( gdept_0(ji,jj,ik) - gdept_1d(ik) ) / ( gdept_1d(ik+1) - gdept_1d(ik) ) 
272                     ptsd(ji,jj,ik,jp_tem) = (1.-zl) * ptsd(ji,jj,ik,jp_tem) + zl * ptsd(ji,jj,ik+1,jp_tem)
273                     ptsd(ji,jj,ik,jp_sal) = (1.-zl) * ptsd(ji,jj,ik,jp_sal) + zl * ptsd(ji,jj,ik+1,jp_sal)
274                  END IF
275               END DO
276            END DO
277         ENDIF
278         !
279      ENDIF
280      !
281      IF( .NOT.ln_tsd_tradmp ) THEN                   !==   deallocate T & S structure   ==!
282         !                                              (data used only for initialisation)
283         IF(lwp) WRITE(numout,*) 'dta_tsd: deallocte T & S arrays as they are only use to initialize the run'
284                                        DEALLOCATE( sf_tsd(jp_tem)%fnow )     ! T arrays in the structure
285         IF( sf_tsd(jp_tem)%ln_tint )   DEALLOCATE( sf_tsd(jp_tem)%fdta )
286                                        DEALLOCATE( sf_tsd(jp_sal)%fnow )     ! S arrays in the structure
287         IF( sf_tsd(jp_sal)%ln_tint )   DEALLOCATE( sf_tsd(jp_sal)%fdta )
288         IF( ln_tsd_interp )            DEALLOCATE( sf_tsd(jp_dep)%fnow )     ! T arrays in the structure
289         IF( ln_tsd_interp )            DEALLOCATE( sf_tsd(jp_msk)%fnow )     ! T arrays in the structure
290                                        DEALLOCATE( sf_tsd              )     ! the structure itself
291      ENDIF
292      !
293      IF( nn_timing == 1 )  CALL timing_stop('dta_tsd')
294      !
295   END SUBROUTINE dta_tsd
296
297   !!======================================================================
298END MODULE dtatsd
Note: See TracBrowser for help on using the repository browser.