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.
dtadyn.F90 in branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/OFF_SRC – NEMO

source: branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/OFF_SRC/dtadyn.F90 @ 6812

Last change on this file since 6812 was 4148, checked in by cetlod, 11 years ago

merge in trunk changes between r3853 and r3940 and commit the changes, see ticket #1169
  • Property svn:keywords set to Id
File size: 29.8 KB
Line 
1MODULE dtadyn
2   !!======================================================================
3   !!                       ***  MODULE  dtadyn  ***
4   !! Off-line : interpolation of the physical fields
5   !!======================================================================
6   !! History :   OPA  ! 1992-01 (M. Imbard) Original code
7   !!             8.0  ! 1998-04 (L.Bopp MA Foujols) slopes for isopyc.
8   !!              -   ! 1998-05 (L. Bopp) read output of coupled run
9   !!             8.2  ! 2001-01 (M. Levy et M. Benjelloul) add netcdf FORMAT
10   !!   NEMO      1.0  ! 2005-03 (O. Aumont and A. El Moussaoui) F90
11   !!              -   ! 2005-12 (C. Ethe) Adapted for DEGINT
12   !!             3.0  ! 2007-06 (C. Ethe) use of iom module
13   !!             3.3  ! 2010-11 (C. Ethe) Full reorganization of the off-line: phasing with the on-line
14   !!             3.4  ! 2011-05 (C. Ethe) Use of fldread
15   !!----------------------------------------------------------------------
16
17   !!----------------------------------------------------------------------
18   !!   dta_dyn_init : initialization, namelist read, and SAVEs control
19   !!   dta_dyn      : Interpolation of the fields
20   !!----------------------------------------------------------------------
21   USE oce             ! ocean dynamics and tracers variables
22   USE c1d             ! 1D configuration: lk_c1d
23   USE dom_oce         ! ocean domain: variables
24   USE zdf_oce         ! ocean vertical physics: variables
25   USE sbc_oce         ! surface module: variables
26   USE trc_oce         ! share ocean/biogeo variables
27   USE phycst          ! physical constants
28   USE trabbl          ! active tracer: bottom boundary layer
29   USE ldfslp          ! lateral diffusion: iso-neutral slopes
30   USE ldfeiv          ! eddy induced velocity coef.
31   USE ldftra_oce      ! ocean tracer   lateral physics
32   USE zdfmxl          ! vertical physics: mixed layer depth
33   USE eosbn2          ! equation of state - Brunt Vaisala frequency
34   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
35   USE zpshde          ! z-coord. with partial steps: horizontal derivatives
36   USE in_out_manager  ! I/O manager
37   USE iom             ! I/O library
38   USE lib_mpp         ! distributed memory computing library
39   USE prtctl          ! print control
40   USE fldread         ! read input fields
41   USE timing          ! Timing
42
43   IMPLICIT NONE
44   PRIVATE
45
46   PUBLIC   dta_dyn_init   ! called by opa.F90
47   PUBLIC   dta_dyn        ! called by step.F90
48
49   CHARACTER(len=100) ::   cn_dir       !: Root directory for location of ssr files
50   LOGICAL            ::   ln_dynwzv    !: vertical velocity read in a file (T) or computed from u/v (F)
51   LOGICAL            ::   ln_dynbbl    !: bbl coef read in a file (T) or computed (F)
52   LOGICAL            ::   ln_degrad    !: degradation option enabled or not
53
54   INTEGER  , PARAMETER ::   jpfld = 20     ! maximum number of fields to read
55   INTEGER  , SAVE      ::   jf_tem         ! index of temperature
56   INTEGER  , SAVE      ::   jf_sal         ! index of salinity
57   INTEGER  , SAVE      ::   jf_uwd         ! index of u-wind
58   INTEGER  , SAVE      ::   jf_vwd         ! index of v-wind
59   INTEGER  , SAVE      ::   jf_wwd         ! index of w-wind
60   INTEGER  , SAVE      ::   jf_avt         ! index of Kz
61   INTEGER  , SAVE      ::   jf_mld         ! index of mixed layer deptht
62   INTEGER  , SAVE      ::   jf_emp         ! index of water flux
63   INTEGER  , SAVE      ::   jf_qsr         ! index of solar radiation
64   INTEGER  , SAVE      ::   jf_wnd         ! index of wind speed
65   INTEGER  , SAVE      ::   jf_ice         ! index of sea ice cover
66   INTEGER  , SAVE      ::   jf_ubl         ! index of u-bbl coef
67   INTEGER  , SAVE      ::   jf_vbl         ! index of v-bbl coef
68   INTEGER  , SAVE      ::   jf_ahu         ! index of u-diffusivity coef
69   INTEGER  , SAVE      ::   jf_ahv         ! index of v-diffusivity coef
70   INTEGER  , SAVE      ::   jf_ahw         ! index of w-diffusivity coef
71   INTEGER  , SAVE      ::   jf_eiu         ! index of u-eiv
72   INTEGER  , SAVE      ::   jf_eiv         ! index of v-eiv
73   INTEGER  , SAVE      ::   jf_eiw         ! index of w-eiv
74   INTEGER  , SAVE      ::   jf_fmf         ! index of downward salt flux
75
76   TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_dyn  ! structure of input fields (file informations, fields read)
77   !                                               !
78   REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: wdta       ! vertical velocity at 2 time step
79   REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:  ) :: wnow       ! vertical velocity at 2 time step
80   REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: uslpdta    ! zonal isopycnal slopes
81   REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: vslpdta    ! meridional isopycnal slopes
82   REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: wslpidta   ! zonal diapycnal slopes
83   REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: wslpjdta   ! meridional diapycnal slopes
84   REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:)   :: uslpnow    ! zonal isopycnal slopes
85   REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:)   :: vslpnow    ! meridional isopycnal slopes
86   REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:)   :: wslpinow   ! zonal diapycnal slopes
87   REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:)   :: wslpjnow   ! meridional diapycnal slopes
88
89   INTEGER :: nrecprev_tem , nrecprev_uwd
90
91   !! * Substitutions
92#  include "domzgr_substitute.h90"
93#  include "vectopt_loop_substitute.h90"
94   !!----------------------------------------------------------------------
95   !! NEMO/OFF 3.3 , NEMO Consortium (2010)
96   !! $Id$
97   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
98   !!----------------------------------------------------------------------
99CONTAINS
100
101   SUBROUTINE dta_dyn( kt )
102      !!----------------------------------------------------------------------
103      !!                  ***  ROUTINE dta_dyn  ***
104      !!
105      !! ** Purpose :  Prepares dynamics and physics fields from a NEMO run
106      !!               for an off-line simulation of passive tracers
107      !!
108      !! ** Method : calculates the position of data
109      !!             - computes slopes if needed
110      !!             - interpolates data if needed
111      !!----------------------------------------------------------------------
112      !
113      USE oce, ONLY:  zts    => tsa 
114      USE oce, ONLY:  zuslp  => ua   , zvslp  => va
115      USE oce, ONLY:  zwslpi => rotb , zwslpj => rotn
116      USE oce, ONLY:  zu     => ub   , zv     => vb,  zw => hdivb
117      !
118      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
119      !
120      INTEGER  ::   ji, jj     ! dummy loop indices
121      INTEGER  ::   isecsbc    ! number of seconds between Jan. 1st 00h of nit000 year and the middle of time step
122      REAL(wp) ::   ztinta     ! ratio applied to after  records when doing time interpolation
123      REAL(wp) ::   ztintb     ! ratio applied to before records when doing time interpolation
124      INTEGER  ::   iswap_tem, iswap_uwd    !
125      !!----------------------------------------------------------------------
126     
127      !
128      IF( nn_timing == 1 )  CALL timing_start( 'dta_dyn')
129      !
130      isecsbc = nsec_year + nsec1jan000 
131      !
132      IF( kt == nit000 ) THEN
133         nrecprev_tem = 0
134         nrecprev_uwd = 0
135         !
136         CALL fld_read( kt, 1, sf_dyn )      !==   read data at kt time step   ==!
137         !
138         IF( lk_ldfslp .AND. .NOT.lk_c1d .AND. sf_dyn(jf_tem)%ln_tint ) THEN    ! Computes slopes (here avt is used as workspace)                       
139            zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fdta(:,:,:,1) * tmask(:,:,:)   ! temperature
140            zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fdta(:,:,:,1) * tmask(:,:,:)   ! salinity
141            avt(:,:,:)        = sf_dyn(jf_avt)%fdta(:,:,:,1) * tmask(:,:,:)   ! vertical diffusive coef.
142            CALL dta_dyn_slp( kt, zts, zuslp, zvslp, zwslpi, zwslpj )
143            uslpdta (:,:,:,1) = zuslp (:,:,:) 
144            vslpdta (:,:,:,1) = zvslp (:,:,:) 
145            wslpidta(:,:,:,1) = zwslpi(:,:,:) 
146            wslpjdta(:,:,:,1) = zwslpj(:,:,:) 
147         ENDIF
148         IF( ln_dynwzv .AND. sf_dyn(jf_uwd)%ln_tint )  THEN    ! compute vertical velocity from u/v
149            zu(:,:,:) = sf_dyn(jf_uwd)%fdta(:,:,:,1)
150            zv(:,:,:) = sf_dyn(jf_vwd)%fdta(:,:,:,1)
151            CALL dta_dyn_wzv( zu, zv, zw )
152            wdta(:,:,:,1) = zw(:,:,:) * tmask(:,:,:)
153         ENDIF
154      ELSE
155         nrecprev_tem = sf_dyn(jf_tem)%nrec_a(2)
156         nrecprev_uwd = sf_dyn(jf_uwd)%nrec_a(2)
157         !
158         CALL fld_read( kt, 1, sf_dyn )      !==   read data at kt time step   ==!
159         !
160      ENDIF
161      !
162      IF( lk_ldfslp .AND. .NOT.lk_c1d ) THEN    ! Computes slopes (here avt is used as workspace)                       
163         iswap_tem = 0
164         IF(  kt /= nit000 .AND. ( sf_dyn(jf_tem)%nrec_a(2) - nrecprev_tem ) /= 0 )  iswap_tem = 1
165         IF( ( isecsbc > sf_dyn(jf_tem)%nrec_b(2) .AND. iswap_tem == 1 ) .OR. kt == nit000 )  THEN    ! read/update the after data
166            IF(lwp) WRITE(numout,*) ' Compute new slopes at kt = ', kt
167            IF( sf_dyn(jf_tem)%ln_tint ) THEN                 ! time interpolation of data
168               IF( kt /= nit000 ) THEN
169                  uslpdta (:,:,:,1) =  uslpdta (:,:,:,2)         ! swap the data
170                  vslpdta (:,:,:,1) =  vslpdta (:,:,:,2) 
171                  wslpidta(:,:,:,1) =  wslpidta(:,:,:,2) 
172                  wslpjdta(:,:,:,1) =  wslpjdta(:,:,:,2) 
173               ENDIF
174               !
175               zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fdta(:,:,:,2) * tmask(:,:,:)   ! temperature
176               zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fdta(:,:,:,2) * tmask(:,:,:)   ! salinity
177               avt(:,:,:)        = sf_dyn(jf_avt)%fdta(:,:,:,2) * tmask(:,:,:)   ! vertical diffusive coef.
178               CALL dta_dyn_slp( kt, zts, zuslp, zvslp, zwslpi, zwslpj )
179               !
180               uslpdta (:,:,:,2) = zuslp (:,:,:) 
181               vslpdta (:,:,:,2) = zvslp (:,:,:) 
182               wslpidta(:,:,:,2) = zwslpi(:,:,:) 
183               wslpjdta(:,:,:,2) = zwslpj(:,:,:) 
184            ELSE
185               zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fnow(:,:,:) * tmask(:,:,:)
186               zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fnow(:,:,:) * tmask(:,:,:)
187               avt(:,:,:)        = sf_dyn(jf_avt)%fnow(:,:,:) * tmask(:,:,:)
188               CALL dta_dyn_slp( kt, zts, zuslp, zvslp, zwslpi, zwslpj )
189               uslpnow (:,:,:)   = zuslp (:,:,:) 
190               vslpnow (:,:,:)   = zvslp (:,:,:) 
191               wslpinow(:,:,:)   = zwslpi(:,:,:) 
192               wslpjnow(:,:,:)   = zwslpj(:,:,:) 
193            ENDIF
194         ENDIF
195         IF( sf_dyn(jf_tem)%ln_tint )  THEN
196            ztinta =  REAL( isecsbc - sf_dyn(jf_tem)%nrec_b(2), wp )  &
197               &    / REAL( sf_dyn(jf_tem)%nrec_a(2) - sf_dyn(jf_tem)%nrec_b(2), wp )
198            ztintb =  1. - ztinta
199            uslp (:,:,:) = ztintb * uslpdta (:,:,:,1)  + ztinta * uslpdta (:,:,:,2) 
200            vslp (:,:,:) = ztintb * vslpdta (:,:,:,1)  + ztinta * vslpdta (:,:,:,2) 
201            wslpi(:,:,:) = ztintb * wslpidta(:,:,:,1)  + ztinta * wslpidta(:,:,:,2) 
202            wslpj(:,:,:) = ztintb * wslpjdta(:,:,:,1)  + ztinta * wslpjdta(:,:,:,2) 
203         ELSE
204            uslp (:,:,:) = uslpnow (:,:,:)
205            vslp (:,:,:) = vslpnow (:,:,:)
206            wslpi(:,:,:) = wslpinow(:,:,:)
207            wslpj(:,:,:) = wslpjnow(:,:,:)
208         ENDIF
209      ENDIF
210      !
211      IF( ln_dynwzv )  THEN    ! compute vertical velocity from u/v
212         iswap_uwd = 0
213         IF(  kt /= nit000 .AND. ( sf_dyn(jf_uwd)%nrec_a(2) - nrecprev_uwd ) /= 0 )  iswap_uwd = 1
214         IF( ( isecsbc > sf_dyn(jf_uwd)%nrec_b(2) .AND. iswap_uwd == 1 ) .OR. kt == nit000 )  THEN    ! read/update the after data
215            IF(lwp) WRITE(numout,*) ' Compute new vertical velocity at kt = ', kt
216            IF(lwp) WRITE(numout,*)
217            IF( sf_dyn(jf_uwd)%ln_tint ) THEN                 ! time interpolation of data
218               IF( kt /= nit000 )  THEN
219                  wdta(:,:,:,1) =  wdta(:,:,:,2)     ! swap the data for initialisation
220               ENDIF
221               zu(:,:,:) = sf_dyn(jf_uwd)%fdta(:,:,:,2)
222               zv(:,:,:) = sf_dyn(jf_vwd)%fdta(:,:,:,2)
223               CALL dta_dyn_wzv( zu, zv, zw )
224               wdta(:,:,:,2) = zw(:,:,:) * tmask(:,:,:)
225            ELSE
226               zu(:,:,:) = sf_dyn(jf_uwd)%fnow(:,:,:) 
227               zv(:,:,:) = sf_dyn(jf_vwd)%fnow(:,:,:)
228               CALL dta_dyn_wzv( zu, zv, zw )
229               wnow(:,:,:)  = zw(:,:,:) * tmask(:,:,:)
230            ENDIF
231         ENDIF
232         IF( sf_dyn(jf_uwd)%ln_tint )  THEN
233            ztinta =  REAL( isecsbc - sf_dyn(jf_uwd)%nrec_b(2), wp )  &
234               &    / REAL( sf_dyn(jf_uwd)%nrec_a(2) - sf_dyn(jf_uwd)%nrec_b(2), wp )
235            ztintb =  1. - ztinta
236            wn(:,:,:) = ztintb * wdta(:,:,:,1)  + ztinta * wdta(:,:,:,2) 
237         ELSE
238            wn(:,:,:) = wnow(:,:,:)
239         ENDIF
240      ENDIF
241      !
242      tsn(:,:,:,jp_tem) = sf_dyn(jf_tem)%fnow(:,:,:)  * tmask(:,:,:)    ! temperature
243      tsn(:,:,:,jp_sal) = sf_dyn(jf_sal)%fnow(:,:,:)  * tmask(:,:,:)    ! salinity
244      !
245      CALL eos    ( tsn, rhd, rhop )                                       ! In any case, we need rhop
246      CALL zdf_mxl( kt )                                                   ! In any case, we need mxl
247      !
248      avt(:,:,:)       = sf_dyn(jf_avt)%fnow(:,:,:)  * tmask(:,:,:)    ! vertical diffusive coefficient
249      un (:,:,:)       = sf_dyn(jf_uwd)%fnow(:,:,:)  * umask(:,:,:)    ! u-velocity
250      vn (:,:,:)       = sf_dyn(jf_vwd)%fnow(:,:,:)  * vmask(:,:,:)    ! v-velocity
251      IF( .NOT.ln_dynwzv ) &                                          ! w-velocity read in file
252         wn (:,:,:)    = sf_dyn(jf_wwd)%fnow(:,:,:) * tmask(:,:,:)   
253      hmld(:,:)        = sf_dyn(jf_mld)%fnow(:,:,1) * tmask(:,:,1)    ! mixed layer depht
254      wndm(:,:)        = sf_dyn(jf_wnd)%fnow(:,:,1) * tmask(:,:,1)    ! wind speed - needed for gas exchange
255      emp (:,:)        = sf_dyn(jf_emp)%fnow(:,:,1) * tmask(:,:,1)    ! E-P
256      fmmflx(:,:)      = sf_dyn(jf_fmf)%fnow(:,:,1) * tmask(:,:,1)    ! downward salt flux (v3.5+)
257      fr_i(:,:)        = sf_dyn(jf_ice)%fnow(:,:,1) * tmask(:,:,1)    ! Sea-ice fraction
258      qsr (:,:)        = sf_dyn(jf_qsr)%fnow(:,:,1) * tmask(:,:,1)    ! solar radiation
259
260      !                                                      ! bbl diffusive coef
261#if defined key_trabbl && ! defined key_c1d
262      IF( ln_dynbbl ) THEN                                        ! read in a file
263         ahu_bbl(:,:)  = sf_dyn(jf_ubl)%fnow(:,:,1) * umask(:,:,1)
264         ahv_bbl(:,:)  = sf_dyn(jf_vbl)%fnow(:,:,1) * vmask(:,:,1)
265      ELSE                                                        ! Compute bbl coefficients if needed
266         tsb(:,:,:,:) = tsn(:,:,:,:)
267         CALL bbl( kt, nit000, 'TRC')
268      END IF
269#endif
270#if ( ! defined key_degrad && defined key_traldf_c2d && defined key_traldf_eiv ) && ! defined key_c1d
271      aeiw(:,:)        = sf_dyn(jf_eiw)%fnow(:,:,1) * tmask(:,:,1)    ! w-eiv
272      !                                                           ! Computes the horizontal values from the vertical value
273      DO jj = 2, jpjm1
274         DO ji = fs_2, fs_jpim1   ! vector opt.
275            aeiu(ji,jj) = .5 * ( aeiw(ji,jj) + aeiw(ji+1,jj  ) )  ! Average the diffusive coefficient at u- v- points
276            aeiv(ji,jj) = .5 * ( aeiw(ji,jj) + aeiw(ji  ,jj+1) )  ! at u- v- points
277         END DO
278      END DO
279      CALL lbc_lnk( aeiu, 'U', 1. )   ;   CALL lbc_lnk( aeiv, 'V', 1. )    ! lateral boundary condition
280#endif
281     
282#if defined key_degrad && ! defined key_c1d
283      !                                          ! degrad option : diffusive and eiv coef are 3D
284      ahtu(:,:,:) = sf_dyn(jf_ahu)%fnow(:,:,:) * umask(:,:,:)
285      ahtv(:,:,:) = sf_dyn(jf_ahv)%fnow(:,:,:) * vmask(:,:,:)
286      ahtw(:,:,:) = sf_dyn(jf_ahw)%fnow(:,:,:) * tmask(:,:,:)
287#  if defined key_traldf_eiv 
288      aeiu(:,:,:) = sf_dyn(jf_eiu)%fnow(:,:,:) * umask(:,:,:)
289      aeiv(:,:,:) = sf_dyn(jf_eiv)%fnow(:,:,:) * vmask(:,:,:)
290      aeiw(:,:,:) = sf_dyn(jf_eiw)%fnow(:,:,:) * tmask(:,:,:)
291#  endif
292#endif
293      !
294      IF(ln_ctl) THEN                  ! print control
295         CALL prt_ctl(tab3d_1=tsn(:,:,:,jp_tem), clinfo1=' tn      - : ', mask1=tmask, ovlap=1, kdim=jpk   )
296         CALL prt_ctl(tab3d_1=tsn(:,:,:,jp_sal), clinfo1=' sn      - : ', mask1=tmask, ovlap=1, kdim=jpk   )
297         CALL prt_ctl(tab3d_1=un               , clinfo1=' un      - : ', mask1=umask, ovlap=1, kdim=jpk   )
298         CALL prt_ctl(tab3d_1=vn               , clinfo1=' vn      - : ', mask1=vmask, ovlap=1, kdim=jpk   )
299         CALL prt_ctl(tab3d_1=wn               , clinfo1=' wn      - : ', mask1=tmask, ovlap=1, kdim=jpk   )
300         CALL prt_ctl(tab3d_1=avt              , clinfo1=' kz      - : ', mask1=tmask, ovlap=1, kdim=jpk   )
301         CALL prt_ctl(tab2d_1=fr_i             , clinfo1=' fr_i    - : ', mask1=tmask, ovlap=1 )
302         CALL prt_ctl(tab2d_1=hmld             , clinfo1=' hmld    - : ', mask1=tmask, ovlap=1 )
303         CALL prt_ctl(tab2d_1=fmmflx           , clinfo1=' fmmflx  - : ', mask1=tmask, ovlap=1 )
304         CALL prt_ctl(tab2d_1=emp              , clinfo1=' emp     - : ', mask1=tmask, ovlap=1 )
305         CALL prt_ctl(tab2d_1=wndm             , clinfo1=' wspd    - : ', mask1=tmask, ovlap=1 )
306         CALL prt_ctl(tab2d_1=qsr              , clinfo1=' qsr     - : ', mask1=tmask, ovlap=1 )
307      ENDIF
308      !
309      IF( nn_timing == 1 )  CALL timing_stop( 'dta_dyn')
310      !
311   END SUBROUTINE dta_dyn
312
313
314   SUBROUTINE dta_dyn_init
315      !!----------------------------------------------------------------------
316      !!                  ***  ROUTINE dta_dyn_init  ***
317      !!
318      !! ** Purpose :   Initialisation of the dynamical data     
319      !! ** Method  : - read the data namdta_dyn namelist
320      !!
321      !! ** Action  : - read parameters
322      !!----------------------------------------------------------------------
323      INTEGER  :: ierr, ierr0, ierr1, ierr2, ierr3   ! return error code
324      INTEGER  :: ifpr                               ! dummy loop indice
325      INTEGER  :: jfld                               ! dummy loop arguments
326      INTEGER  :: inum, idv, idimv                   ! local integer
327      INTEGER  :: ios                                ! Local integer output status for namelist read
328      !!
329      CHARACTER(len=100)            ::  cn_dir   !   Root directory for location of core files
330      TYPE(FLD_N), DIMENSION(jpfld) ::  slf_d    ! array of namelist informations on the fields to read
331      TYPE(FLD_N) :: sn_tem, sn_sal, sn_mld, sn_emp, sn_ice, sn_qsr, sn_wnd  ! informations about the fields to be read
332      TYPE(FLD_N) :: sn_uwd, sn_vwd, sn_wwd, sn_avt, sn_ubl, sn_vbl          !   "                                 "
333      TYPE(FLD_N) :: sn_ahu, sn_ahv, sn_ahw, sn_eiu, sn_eiv, sn_eiw, sn_fmf  !   "                                 "
334      !!----------------------------------------------------------------------
335      !
336      NAMELIST/namdta_dyn/cn_dir, ln_dynwzv, ln_dynbbl, ln_degrad,    &
337         &                sn_tem, sn_sal, sn_mld, sn_emp, sn_ice, sn_qsr, sn_wnd,  &
338         &                sn_uwd, sn_vwd, sn_wwd, sn_avt, sn_ubl, sn_vbl,          &
339         &                sn_ahu, sn_ahv, sn_ahw, sn_eiu, sn_eiv, sn_eiw, sn_fmf
340      !
341      REWIND( numnam_ref )              ! Namelist namdta_dyn in reference namelist : Offline: init. of dynamical data
342      READ  ( numnam_ref, namdta_dyn, IOSTAT = ios, ERR = 901)
343901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdta_dyn in reference namelist', lwp )
344
345      REWIND( numnam_cfg )              ! Namelist namdta_dyn in configuration namelist : Offline: init. of dynamical data
346      READ  ( numnam_cfg, namdta_dyn, IOSTAT = ios, ERR = 902 )
347902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdta_dyn in configuration namelist', lwp )
348      WRITE ( numond, namdta_dyn )
349
350      !                                         ! store namelist information in an array
351      !                                         ! Control print
352      IF(lwp) THEN
353         WRITE(numout,*)
354         WRITE(numout,*) 'dta_dyn : offline dynamics '
355         WRITE(numout,*) '~~~~~~~ '
356         WRITE(numout,*) '   Namelist namdta_dyn'
357         WRITE(numout,*) '      vertical velocity read from file (T) or computed (F) ln_dynwzv  = ', ln_dynwzv
358         WRITE(numout,*) '      bbl coef read from file (T) or computed (F)          ln_dynbbl  = ', ln_dynbbl
359         WRITE(numout,*) '      degradation option enabled (T) or not (F)            ln_degrad  = ', ln_degrad
360         WRITE(numout,*)
361      ENDIF
362      !
363      IF( ln_degrad .AND. .NOT.lk_degrad ) THEN
364         CALL ctl_warn( 'dta_dyn_init: degradation option requires key_degrad activated ; force ln_degrad to false' )
365         ln_degrad = .FALSE.
366      ENDIF
367      IF( ln_dynbbl .AND. ( .NOT.lk_trabbl .OR. lk_c1d ) ) THEN
368         CALL ctl_warn( 'dta_dyn_init: bbl option requires key_trabbl activated ; force ln_dynbbl to false' )
369         ln_dynbbl = .FALSE.
370      ENDIF
371
372      jf_tem = 1   ;   jf_sal = 2   ;  jf_mld = 3   ;  jf_emp = 4   ;   jf_fmf  = 5   ;  jf_ice = 6   ;   jf_qsr = 7
373      jf_wnd = 8   ;   jf_uwd = 9   ;  jf_vwd = 10  ;  jf_wwd = 11  ;   jf_avt  = 12  ;  jfld  = 12
374      !
375      slf_d(jf_tem) = sn_tem   ;   slf_d(jf_sal)  = sn_sal   ;   slf_d(jf_mld) = sn_mld
376      slf_d(jf_emp) = sn_emp   ;   slf_d(jf_fmf ) = sn_fmf   ;   slf_d(jf_ice) = sn_ice 
377      slf_d(jf_qsr) = sn_qsr   ;   slf_d(jf_wnd)  = sn_wnd   ;   slf_d(jf_avt) = sn_avt 
378      slf_d(jf_uwd) = sn_uwd   ;   slf_d(jf_vwd)  = sn_vwd   ;   slf_d(jf_wwd) = sn_wwd
379      !
380      IF( .NOT.ln_degrad ) THEN     ! no degrad option
381         IF( lk_traldf_eiv .AND. ln_dynbbl ) THEN        ! eiv & bbl
382                 jf_ubl  = 13      ;         jf_vbl  = 14      ;         jf_eiw  = 15   ;   jfld = 15
383           slf_d(jf_ubl) = sn_ubl  ;   slf_d(jf_vbl) = sn_vbl  ;   slf_d(jf_eiw) = sn_eiw
384         ENDIF
385         IF( .NOT.lk_traldf_eiv .AND. ln_dynbbl ) THEN   ! no eiv & bbl
386                 jf_ubl  = 13      ;         jf_vbl  = 14      ;   jfld = 14
387           slf_d(jf_ubl) = sn_ubl  ;   slf_d(jf_vbl) = sn_vbl
388         ENDIF
389         IF( lk_traldf_eiv .AND. .NOT.ln_dynbbl ) THEN   ! eiv & no bbl
390           jf_eiw = 13   ;   jfld = 13   ;   slf_d(jf_eiw) = sn_eiw
391         ENDIF
392      ELSE
393              jf_ahu  = 13      ;         jf_ahv  = 14      ;         jf_ahw  = 15   ;   jfld = 15
394        slf_d(jf_ahu) = sn_ahu  ;   slf_d(jf_ahv) = sn_ahv  ;   slf_d(jf_ahw) = sn_ahw
395        IF( lk_traldf_eiv .AND. ln_dynbbl ) THEN         ! eiv & bbl
396                 jf_ubl  = 16      ;         jf_vbl  = 17     
397           slf_d(jf_ubl) = sn_ubl  ;   slf_d(jf_vbl) = sn_vbl 
398                 jf_eiu  = 18      ;         jf_eiv  = 19      ;          jf_eiw  = 20   ;   jfld = 20
399           slf_d(jf_eiu) = sn_eiu  ;   slf_d(jf_eiv) = sn_eiv  ;    slf_d(jf_eiw) = sn_eiw
400        ENDIF
401        IF( .NOT.lk_traldf_eiv .AND. ln_dynbbl ) THEN    ! no eiv & bbl
402                 jf_ubl  = 16      ;         jf_vbl  = 17      ;   jfld = 17
403           slf_d(jf_ubl) = sn_ubl  ;   slf_d(jf_vbl) = sn_vbl
404        ENDIF
405        IF( lk_traldf_eiv .AND. .NOT.ln_dynbbl ) THEN    ! eiv & no bbl
406                 jf_eiu  = 16      ;         jf_eiv  = 17      ;         jf_eiw  = 18   ;   jfld = 18
407           slf_d(jf_eiu) = sn_eiu  ;   slf_d(jf_eiv) = sn_eiv  ;   slf_d(jf_eiw) = sn_eiw
408        ENDIF
409      ENDIF
410 
411      ALLOCATE( sf_dyn(jfld), STAT=ierr )         ! set sf structure
412      IF( ierr > 0 ) THEN
413         CALL ctl_stop( 'dta_dyn: unable to allocate sf structure' )   ;   RETURN
414      ENDIF
415      ! Open file for each variable to get his number of dimension
416      DO ifpr = 1, jfld
417         CALL iom_open( TRIM( cn_dir )//TRIM( slf_d(ifpr)%clname ), inum )
418         idv   = iom_varid( inum , slf_d(ifpr)%clvar )  ! id of the variable sdjf%clvar
419         idimv = iom_file ( inum )%ndims(idv)             ! number of dimension for variable sdjf%clvar
420         IF( inum /= 0 )   CALL iom_close( inum )       ! close file if already open
421         IF( idimv == 3 ) THEN    ! 2D variable
422                                      ALLOCATE( sf_dyn(ifpr)%fnow(jpi,jpj,1)    , STAT=ierr0 )
423            IF( slf_d(ifpr)%ln_tint ) ALLOCATE( sf_dyn(ifpr)%fdta(jpi,jpj,1,2)  , STAT=ierr1 )
424         ELSE                     ! 3D variable
425                                      ALLOCATE( sf_dyn(ifpr)%fnow(jpi,jpj,jpk)  , STAT=ierr0 )
426            IF( slf_d(ifpr)%ln_tint ) ALLOCATE( sf_dyn(ifpr)%fdta(jpi,jpj,jpk,2), STAT=ierr1 )
427         ENDIF
428         IF( ierr0 + ierr1 > 0 ) THEN
429            CALL ctl_stop( 'dta_dyn_init : unable to allocate sf_dyn array structure' )   ;   RETURN
430         ENDIF
431      END DO
432      !                                         ! fill sf with slf_i and control print
433      CALL fld_fill( sf_dyn, slf_d, cn_dir, 'dta_dyn_init', 'Data in file', 'namdta_dyn' )
434      !
435      IF( lk_ldfslp .AND. .NOT.lk_c1d ) THEN                  ! slopes
436         IF( sf_dyn(jf_tem)%ln_tint ) THEN      ! time interpolation
437            ALLOCATE( uslpdta (jpi,jpj,jpk,2), vslpdta (jpi,jpj,jpk,2),    &
438            &         wslpidta(jpi,jpj,jpk,2), wslpjdta(jpi,jpj,jpk,2), STAT=ierr2 )
439         ELSE
440            ALLOCATE( uslpnow (jpi,jpj,jpk)  , vslpnow (jpi,jpj,jpk)  ,    &
441            &         wslpinow(jpi,jpj,jpk)  , wslpjnow(jpi,jpj,jpk)  , STAT=ierr2 )
442         ENDIF
443         IF( ierr2 > 0 ) THEN
444            CALL ctl_stop( 'dta_dyn_init : unable to allocate slope arrays' )   ;   RETURN
445         ENDIF
446      ENDIF
447      IF( ln_dynwzv ) THEN                  ! slopes
448         IF( sf_dyn(jf_uwd)%ln_tint ) THEN      ! time interpolation
449            ALLOCATE( wdta(jpi,jpj,jpk,2), STAT=ierr3 )
450         ELSE
451            ALLOCATE( wnow(jpi,jpj,jpk)  , STAT=ierr3 )
452         ENDIF
453         IF( ierr3 > 0 ) THEN
454            CALL ctl_stop( 'dta_dyn_init : unable to allocate wdta arrays' )   ;   RETURN
455         ENDIF
456      ENDIF
457      !
458      CALL dta_dyn( nit000 )
459      !
460   END SUBROUTINE dta_dyn_init
461
462   SUBROUTINE dta_dyn_wzv( pu, pv, pw )
463      !!----------------------------------------------------------------------
464      !!                    ***  ROUTINE wzv  ***
465      !!
466      !! ** Purpose :   Compute the now vertical velocity after the array swap
467      !!
468      !! ** Method  : - compute the now divergence given by :
469      !!         * z-coordinate ONLY !!!!
470      !!         hdiv = 1/(e1t*e2t) [ di(e2u  u) + dj(e1v  v) ]
471      !!     - Using the incompressibility hypothesis, the vertical
472      !!      velocity is computed by integrating the horizontal divergence
473      !!      from the bottom to the surface.
474      !!        The boundary conditions are w=0 at the bottom (no flux).
475      !!----------------------------------------------------------------------
476      USE oce, ONLY:  zhdiv => hdivn
477      !
478      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) :: pu, pv    !:  horizontal velocities
479      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(  out) :: pw        !:  vertical velocity
480      !!
481      INTEGER  ::  ji, jj, jk
482      REAL(wp) ::  zu, zu1, zv, zv1, zet
483      !!----------------------------------------------------------------------
484      !
485      ! Computation of vertical velocity using horizontal divergence
486      zhdiv(:,:,:) = 0._wp
487      DO jk = 1, jpkm1
488         DO jj = 2, jpjm1
489            DO ji = fs_2, fs_jpim1   ! vector opt.
490               zu  = pu(ji  ,jj  ,jk) * umask(ji  ,jj  ,jk) * e2u(ji  ,jj  ) * fse3u(ji  ,jj  ,jk)
491               zu1 = pu(ji-1,jj  ,jk) * umask(ji-1,jj  ,jk) * e2u(ji-1,jj  ) * fse3u(ji-1,jj  ,jk)
492               zv  = pv(ji  ,jj  ,jk) * vmask(ji  ,jj  ,jk) * e1v(ji  ,jj  ) * fse3v(ji  ,jj  ,jk)
493               zv1 = pv(ji  ,jj-1,jk) * vmask(ji  ,jj-1,jk) * e1v(ji  ,jj-1) * fse3v(ji  ,jj-1,jk)
494               zet = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
495               zhdiv(ji,jj,jk) = ( zu - zu1 + zv - zv1 ) * zet 
496            END DO
497         END DO
498      END DO
499      CALL lbc_lnk( zhdiv, 'T', 1. )      ! Lateral boundary conditions on zhdiv
500      !
501      ! computation of vertical velocity from the bottom
502      pw(:,:,jpk) = 0._wp
503      DO jk = jpkm1, 1, -1
504         pw(:,:,jk) = pw(:,:,jk+1) - fse3t(:,:,jk) * zhdiv(:,:,jk)
505      END DO
506      !
507   END SUBROUTINE dta_dyn_wzv
508
509   SUBROUTINE dta_dyn_slp( kt, pts, puslp, pvslp, pwslpi, pwslpj )
510      !!---------------------------------------------------------------------
511      !!                    ***  ROUTINE dta_dyn_slp  ***
512      !!
513      !! ** Purpose : Computation of slope
514      !!
515      !!---------------------------------------------------------------------
516      INTEGER ,                              INTENT(in ) :: kt       ! time step
517      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: pts      ! temperature/salinity
518      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(out) :: puslp    ! zonal isopycnal slopes
519      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(out) :: pvslp    ! meridional isopycnal slopes
520      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(out) :: pwslpi   ! zonal diapycnal slopes
521      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(out) :: pwslpj   ! meridional diapycnal slopes
522      !!---------------------------------------------------------------------
523#if defined key_ldfslp && ! defined key_c1d
524      CALL eos( pts, rhd, rhop )   ! Time-filtered in situ density
525      CALL bn2( pts, rn2 )         ! before Brunt-Vaisala frequency
526      IF( ln_zps )   &
527         &  CALL zps_hde( kt, jpts, pts, gtsu, gtsv, rhd, gru, grv )  ! Partial steps: before Horizontal DErivative
528         !                                                            ! of t, s, rd at the bottom ocean level
529      CALL zdf_mxl( kt )            ! mixed layer depth
530      CALL ldf_slp( kt, rhd, rn2 )  ! slopes
531      puslp (:,:,:) = uslp (:,:,:) 
532      pvslp (:,:,:) = vslp (:,:,:) 
533      pwslpi(:,:,:) = wslpi(:,:,:) 
534      pwslpj(:,:,:) = wslpj(:,:,:) 
535#else
536      puslp (:,:,:) = 0.            ! to avoid warning when compiling
537      pvslp (:,:,:) = 0.
538      pwslpi(:,:,:) = 0.
539      pwslpj(:,:,:) = 0.
540#endif
541      !
542   END SUBROUTINE dta_dyn_slp
543   !!======================================================================
544END MODULE dtadyn
Note: See TracBrowser for help on using the repository browser.