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/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OFF_SRC – NEMO

source: branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OFF_SRC/dtadyn.F90 @ 5870

Last change on this file since 5870 was 5870, checked in by acc, 8 years ago

Branch 2015/dev_r5803_NOC_WAD. Merge in trunk changes from 5803 to 5869 in preparation for merge. Also tidied and reorganised some wetting and drying code. Renamed wadlmt.F90 to wetdry.F90. Wetting drying code changes restricted to domzgr.F90, domvvl.F90 nemogcm.F90 sshwzv.F90, dynspg_ts.F90, wetdry.F90 and dynhpg.F90. Code passes full SETTE tests with ln_wd=.false.. Still awaiting test case for checking with ln_wd=.false.
  • Property svn:keywords set to Id
File size: 28.1 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 ldftra          ! lateral diffusivity coefficients
29   USE trabbl          ! active tracer: bottom boundary layer
30   USE ldfslp          ! lateral diffusion: iso-neutral slopes
31   USE zdfmxl          ! vertical physics: mixed layer depth
32   USE eosbn2          ! equation of state - Brunt Vaisala frequency
33   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
34   USE zpshde          ! z-coord. with partial steps: horizontal derivatives
35   USE in_out_manager  ! I/O manager
36   USE iom             ! I/O library
37   USE lib_mpp         ! distributed memory computing library
38   USE prtctl          ! print control
39   USE fldread         ! read input fields
40   USE timing          ! Timing
41   USE wrk_nemo
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_dynrnf    !: read runoff data in file (T) or set to zero (F)
53
54   INTEGER  , PARAMETER ::   jpfld = 15     ! 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_rnf         ! index of river runoff
67   INTEGER  , SAVE      ::   jf_ubl         ! index of u-bbl coef
68   INTEGER  , SAVE      ::   jf_vbl         ! index of v-bbl coef
69   INTEGER  , SAVE      ::   jf_fmf         ! index of downward salt flux
70
71   TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_dyn  ! structure of input fields (file informations, fields read)
72   !                                               !
73   REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: wdta       ! vertical velocity at 2 time step
74   REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:  ) :: wnow       ! vertical velocity at 2 time step
75   REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: uslpdta    ! zonal isopycnal slopes
76   REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: vslpdta    ! meridional isopycnal slopes
77   REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: wslpidta   ! zonal diapycnal slopes
78   REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: wslpjdta   ! meridional diapycnal slopes
79   REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:)   :: uslpnow    ! zonal isopycnal slopes
80   REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:)   :: vslpnow    ! meridional isopycnal slopes
81   REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:)   :: wslpinow   ! zonal diapycnal slopes
82   REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:)   :: wslpjnow   ! meridional diapycnal slopes
83
84   INTEGER :: nrecprev_tem , nrecprev_uwd
85
86   !! * Substitutions
87#  include "domzgr_substitute.h90"
88#  include "vectopt_loop_substitute.h90"
89   !!----------------------------------------------------------------------
90   !! NEMO/OFF 3.3 , NEMO Consortium (2010)
91   !! $Id$
92   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
93   !!----------------------------------------------------------------------
94CONTAINS
95
96   SUBROUTINE dta_dyn( kt )
97      !!----------------------------------------------------------------------
98      !!                  ***  ROUTINE dta_dyn  ***
99      !!
100      !! ** Purpose :  Prepares dynamics and physics fields from a NEMO run
101      !!               for an off-line simulation of passive tracers
102      !!
103      !! ** Method : calculates the position of data
104      !!             - computes slopes if needed
105      !!             - interpolates data if needed
106      !!----------------------------------------------------------------------
107      USE oce, ONLY:  zts    => tsa
108      USE oce, ONLY:  zuslp  => ua   , zvslp  => va
109      USE oce, ONLY:  zwslpi => ua_sv , zwslpj => va_sv
110      USE oce, ONLY:  zu     => ub   , zv     => vb,  zw => rke
111      !
112      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
113      !
114!      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts)  :: zts
115!      REAL(wp), DIMENSION(jpi,jpj,jpk     )  :: zuslp, zvslp, zwslpi, zwslpj
116!      REAL(wp), DIMENSION(jpi,jpj,jpk     )  :: zu, zv, zw
117      !
118      !
119      INTEGER  ::   ji, jj     ! dummy loop indices
120      INTEGER  ::   isecsbc    ! number of seconds between Jan. 1st 00h of nit000 year and the middle of time step
121      REAL(wp) ::   ztinta     ! ratio applied to after  records when doing time interpolation
122      REAL(wp) ::   ztintb     ! ratio applied to before records when doing time interpolation
123      INTEGER  ::   iswap_tem, iswap_uwd    !
124      !!----------------------------------------------------------------------
125     
126      !
127      IF( nn_timing == 1 )  CALL timing_start( 'dta_dyn')
128      !
129      isecsbc = nsec_year + nsec1jan000 
130      !
131      IF( kt == nit000 ) THEN
132         nrecprev_tem = 0
133         nrecprev_uwd = 0
134         !
135         CALL fld_read( kt, 1, sf_dyn )      !==   read data at kt time step   ==!
136         !
137         IF( l_ldfslp .AND. .NOT.lk_c1d .AND. sf_dyn(jf_tem)%ln_tint ) THEN    ! Computes slopes (here avt is used as workspace)                       
138            zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fdta(:,:,:,1) * tmask(:,:,:)   ! temperature
139            zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fdta(:,:,:,1) * tmask(:,:,:)   ! salinity
140            avt(:,:,:)        = sf_dyn(jf_avt)%fdta(:,:,:,1) * tmask(:,:,:)   ! vertical diffusive coef.
141            CALL dta_dyn_slp( kt, zts, zuslp, zvslp, zwslpi, zwslpj )
142            uslpdta (:,:,:,1) = zuslp (:,:,:) 
143            vslpdta (:,:,:,1) = zvslp (:,:,:) 
144            wslpidta(:,:,:,1) = zwslpi(:,:,:) 
145            wslpjdta(:,:,:,1) = zwslpj(:,:,:) 
146         ENDIF
147         IF( ln_dynwzv .AND. sf_dyn(jf_uwd)%ln_tint )  THEN    ! compute vertical velocity from u/v
148            zu(:,:,:) = sf_dyn(jf_uwd)%fdta(:,:,:,1)
149            zv(:,:,:) = sf_dyn(jf_vwd)%fdta(:,:,:,1)
150            CALL dta_dyn_wzv( zu, zv, zw )
151            wdta(:,:,:,1) = zw(:,:,:) * tmask(:,:,:)
152         ENDIF
153      ELSE
154         nrecprev_tem = sf_dyn(jf_tem)%nrec_a(2)
155         nrecprev_uwd = sf_dyn(jf_uwd)%nrec_a(2)
156         !
157         CALL fld_read( kt, 1, sf_dyn )      !==   read data at kt time step   ==!
158         !
159      ENDIF
160      !
161      IF( l_ldfslp .AND. .NOT.lk_c1d ) THEN    ! Computes slopes (here avt is used as workspace)                       
162         iswap_tem = 0
163         IF(  kt /= nit000 .AND. ( sf_dyn(jf_tem)%nrec_a(2) - nrecprev_tem ) /= 0 )  iswap_tem = 1
164         IF( ( isecsbc > sf_dyn(jf_tem)%nrec_b(2) .AND. iswap_tem == 1 ) .OR. kt == nit000 )  THEN    ! read/update the after data
165            IF(lwp) WRITE(numout,*) ' Compute new slopes at kt = ', kt
166            IF( sf_dyn(jf_tem)%ln_tint ) THEN                 ! time interpolation of data
167               IF( kt /= nit000 ) THEN
168                  uslpdta (:,:,:,1) =  uslpdta (:,:,:,2)         ! swap the data
169                  vslpdta (:,:,:,1) =  vslpdta (:,:,:,2) 
170                  wslpidta(:,:,:,1) =  wslpidta(:,:,:,2) 
171                  wslpjdta(:,:,:,1) =  wslpjdta(:,:,:,2) 
172               ENDIF
173               !
174               zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fdta(:,:,:,2) * tmask(:,:,:)   ! temperature
175               zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fdta(:,:,:,2) * tmask(:,:,:)   ! salinity
176               avt(:,:,:)        = sf_dyn(jf_avt)%fdta(:,:,:,2) * tmask(:,:,:)   ! vertical diffusive coef.
177               CALL dta_dyn_slp( kt, zts, zuslp, zvslp, zwslpi, zwslpj )
178               !
179               uslpdta (:,:,:,2) = zuslp (:,:,:) 
180               vslpdta (:,:,:,2) = zvslp (:,:,:) 
181               wslpidta(:,:,:,2) = zwslpi(:,:,:) 
182               wslpjdta(:,:,:,2) = zwslpj(:,:,:) 
183            ELSE
184               zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fnow(:,:,:) * tmask(:,:,:)
185               zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fnow(:,:,:) * tmask(:,:,:)
186               avt(:,:,:)        = sf_dyn(jf_avt)%fnow(:,:,:) * tmask(:,:,:)
187               CALL dta_dyn_slp( kt, zts, zuslp, zvslp, zwslpi, zwslpj )
188               uslpnow (:,:,:)   = zuslp (:,:,:) 
189               vslpnow (:,:,:)   = zvslp (:,:,:) 
190               wslpinow(:,:,:)   = zwslpi(:,:,:) 
191               wslpjnow(:,:,:)   = zwslpj(:,:,:) 
192            ENDIF
193         ENDIF
194         IF( sf_dyn(jf_tem)%ln_tint )  THEN
195            ztinta =  REAL( isecsbc - sf_dyn(jf_tem)%nrec_b(2), wp )  &
196               &    / REAL( sf_dyn(jf_tem)%nrec_a(2) - sf_dyn(jf_tem)%nrec_b(2), wp )
197            ztintb =  1. - ztinta
198            uslp (:,:,:) = ztintb * uslpdta (:,:,:,1)  + ztinta * uslpdta (:,:,:,2) 
199            vslp (:,:,:) = ztintb * vslpdta (:,:,:,1)  + ztinta * vslpdta (:,:,:,2) 
200            wslpi(:,:,:) = ztintb * wslpidta(:,:,:,1)  + ztinta * wslpidta(:,:,:,2) 
201            wslpj(:,:,:) = ztintb * wslpjdta(:,:,:,1)  + ztinta * wslpjdta(:,:,:,2) 
202         ELSE
203            uslp (:,:,:) = uslpnow (:,:,:)
204            vslp (:,:,:) = vslpnow (:,:,:)
205            wslpi(:,:,:) = wslpinow(:,:,:)
206            wslpj(:,:,:) = wslpjnow(:,:,:)
207         ENDIF
208      ENDIF
209      !
210      IF( ln_dynwzv )  THEN    ! compute vertical velocity from u/v
211         iswap_uwd = 0
212         IF(  kt /= nit000 .AND. ( sf_dyn(jf_uwd)%nrec_a(2) - nrecprev_uwd ) /= 0 )  iswap_uwd = 1
213         IF( ( isecsbc > sf_dyn(jf_uwd)%nrec_b(2) .AND. iswap_uwd == 1 ) .OR. kt == nit000 )  THEN    ! read/update the after data
214            IF(lwp) WRITE(numout,*) ' Compute new vertical velocity at kt = ', kt
215            IF(lwp) WRITE(numout,*)
216            IF( sf_dyn(jf_uwd)%ln_tint ) THEN                 ! time interpolation of data
217               IF( kt /= nit000 )  THEN
218                  wdta(:,:,:,1) =  wdta(:,:,:,2)     ! swap the data for initialisation
219               ENDIF
220               zu(:,:,:) = sf_dyn(jf_uwd)%fdta(:,:,:,2)
221               zv(:,:,:) = sf_dyn(jf_vwd)%fdta(:,:,:,2)
222               CALL dta_dyn_wzv( zu, zv, zw )
223               wdta(:,:,:,2) = zw(:,:,:) * tmask(:,:,:)
224            ELSE
225               zu(:,:,:) = sf_dyn(jf_uwd)%fnow(:,:,:) 
226               zv(:,:,:) = sf_dyn(jf_vwd)%fnow(:,:,:)
227               CALL dta_dyn_wzv( zu, zv, zw )
228               wnow(:,:,:)  = zw(:,:,:) * tmask(:,:,:)
229            ENDIF
230         ENDIF
231         IF( sf_dyn(jf_uwd)%ln_tint )  THEN
232            ztinta =  REAL( isecsbc - sf_dyn(jf_uwd)%nrec_b(2), wp )  &
233               &    / REAL( sf_dyn(jf_uwd)%nrec_a(2) - sf_dyn(jf_uwd)%nrec_b(2), wp )
234            ztintb =  1. - ztinta
235            wn(:,:,:) = ztintb * wdta(:,:,:,1)  + ztinta * wdta(:,:,:,2) 
236         ELSE
237            wn(:,:,:) = wnow(:,:,:)
238         ENDIF
239      ENDIF
240      !
241      tsn(:,:,:,jp_tem) = sf_dyn(jf_tem)%fnow(:,:,:)  * tmask(:,:,:)    ! temperature
242      tsn(:,:,:,jp_sal) = sf_dyn(jf_sal)%fnow(:,:,:)  * tmask(:,:,:)    ! salinity
243      !
244      !
245      CALL eos    ( tsn, rhd, rhop, gdept_0(:,:,:) ) ! In any case, we need rhop
246      CALL eos_rab( tsn, rab_n )       ! now    local thermal/haline expension ratio at T-points
247      CALL bn2    ( tsn, rab_n, rn2 ) ! before Brunt-Vaisala frequency need for zdfmxl
248
249      rn2b(:,:,:) = rn2(:,:,:)         ! need for zdfmxl
250      CALL zdf_mxl( kt )                                                   ! In any case, we need mxl
251      !
252      avt(:,:,:)       = sf_dyn(jf_avt)%fnow(:,:,:)  * tmask(:,:,:)    ! vertical diffusive coefficient
253      un (:,:,:)       = sf_dyn(jf_uwd)%fnow(:,:,:)  * umask(:,:,:)    ! u-velocity
254      vn (:,:,:)       = sf_dyn(jf_vwd)%fnow(:,:,:)  * vmask(:,:,:)    ! v-velocity
255      IF( .NOT.ln_dynwzv ) &                                          ! w-velocity read in file
256         wn (:,:,:)    = sf_dyn(jf_wwd)%fnow(:,:,:) * tmask(:,:,:)   
257      hmld(:,:)        = sf_dyn(jf_mld)%fnow(:,:,1) * tmask(:,:,1)    ! mixed layer depht
258      wndm(:,:)        = sf_dyn(jf_wnd)%fnow(:,:,1) * tmask(:,:,1)    ! wind speed - needed for gas exchange
259      emp (:,:)        = sf_dyn(jf_emp)%fnow(:,:,1) * tmask(:,:,1)    ! E-P
260      fmmflx(:,:)      = sf_dyn(jf_fmf)%fnow(:,:,1) * tmask(:,:,1)    ! downward salt flux (v3.5+)
261      fr_i(:,:)        = sf_dyn(jf_ice)%fnow(:,:,1) * tmask(:,:,1)    ! Sea-ice fraction
262      qsr (:,:)        = sf_dyn(jf_qsr)%fnow(:,:,1) * tmask(:,:,1)    ! solar radiation
263      IF( ln_dynrnf ) &
264      rnf (:,:)        = sf_dyn(jf_rnf)%fnow(:,:,1) * tmask(:,:,1)    ! river runoffs
265
266      !                                               ! update eddy diffusivity coeff. and/or eiv coeff. at kt
267      IF( l_ldftra_time .OR. l_ldfeiv_time )   CALL ldf_tra( kt ) 
268      !                                                      ! bbl diffusive coef
269#if defined key_trabbl && ! defined key_c1d
270      IF( ln_dynbbl ) THEN                                        ! read in a file
271         ahu_bbl(:,:)  = sf_dyn(jf_ubl)%fnow(:,:,1) * umask(:,:,1)
272         ahv_bbl(:,:)  = sf_dyn(jf_vbl)%fnow(:,:,1) * vmask(:,:,1)
273      ELSE                                                        ! Compute bbl coefficients if needed
274         tsb(:,:,:,:) = tsn(:,:,:,:)
275         CALL bbl( kt, nit000, 'TRC')
276      END IF
277#endif
278      !
279      IF(ln_ctl) THEN                  ! print control
280         CALL prt_ctl(tab3d_1=tsn(:,:,:,jp_tem), clinfo1=' tn      - : ', mask1=tmask, ovlap=1, kdim=jpk   )
281         CALL prt_ctl(tab3d_1=tsn(:,:,:,jp_sal), clinfo1=' sn      - : ', mask1=tmask, ovlap=1, kdim=jpk   )
282         CALL prt_ctl(tab3d_1=un               , clinfo1=' un      - : ', mask1=umask, ovlap=1, kdim=jpk   )
283         CALL prt_ctl(tab3d_1=vn               , clinfo1=' vn      - : ', mask1=vmask, ovlap=1, kdim=jpk   )
284         CALL prt_ctl(tab3d_1=wn               , clinfo1=' wn      - : ', mask1=tmask, ovlap=1, kdim=jpk   )
285         CALL prt_ctl(tab3d_1=avt              , clinfo1=' kz      - : ', mask1=tmask, ovlap=1, kdim=jpk   )
286         CALL prt_ctl(tab2d_1=fr_i             , clinfo1=' fr_i    - : ', mask1=tmask, ovlap=1 )
287         CALL prt_ctl(tab2d_1=hmld             , clinfo1=' hmld    - : ', mask1=tmask, ovlap=1 )
288         CALL prt_ctl(tab2d_1=fmmflx           , clinfo1=' fmmflx  - : ', mask1=tmask, ovlap=1 )
289         CALL prt_ctl(tab2d_1=emp              , clinfo1=' emp     - : ', mask1=tmask, ovlap=1 )
290         CALL prt_ctl(tab2d_1=wndm             , clinfo1=' wspd    - : ', mask1=tmask, ovlap=1 )
291         CALL prt_ctl(tab2d_1=qsr              , clinfo1=' qsr     - : ', mask1=tmask, ovlap=1 )
292      ENDIF
293      !
294      IF( nn_timing == 1 )  CALL timing_stop( 'dta_dyn')
295      !
296   END SUBROUTINE dta_dyn
297
298
299   SUBROUTINE dta_dyn_init
300      !!----------------------------------------------------------------------
301      !!                  ***  ROUTINE dta_dyn_init  ***
302      !!
303      !! ** Purpose :   Initialisation of the dynamical data     
304      !! ** Method  : - read the data namdta_dyn namelist
305      !!
306      !! ** Action  : - read parameters
307      !!----------------------------------------------------------------------
308      INTEGER  :: ierr, ierr0, ierr1, ierr2, ierr3   ! return error code
309      INTEGER  :: ifpr                               ! dummy loop indice
310      INTEGER  :: jfld                               ! dummy loop arguments
311      INTEGER  :: inum, idv, idimv                   ! local integer
312      INTEGER  :: ios                                ! Local integer output status for namelist read
313      !!
314      CHARACTER(len=100)            ::  cn_dir   !   Root directory for location of core files
315      TYPE(FLD_N), DIMENSION(jpfld) ::  slf_d    ! array of namelist informations on the fields to read
316      TYPE(FLD_N) :: sn_tem, sn_sal, sn_mld, sn_emp, sn_ice, sn_qsr, sn_wnd, sn_rnf  ! informations about the fields to be read
317      TYPE(FLD_N) :: sn_uwd, sn_vwd, sn_wwd, sn_avt, sn_ubl, sn_vbl, sn_fmf          !   "                                 "
318      NAMELIST/namdta_dyn/cn_dir, ln_dynwzv, ln_dynbbl, ln_dynrnf,    &
319         &                sn_tem, sn_sal, sn_mld, sn_emp, sn_ice, sn_qsr, sn_wnd, sn_rnf,  &
320         &                sn_uwd, sn_vwd, sn_wwd, sn_avt, sn_ubl, sn_vbl, sn_fmf 
321      !!----------------------------------------------------------------------
322      !
323      REWIND( numnam_ref )              ! Namelist namdta_dyn in reference namelist : Offline: init. of dynamical data
324      READ  ( numnam_ref, namdta_dyn, IOSTAT = ios, ERR = 901)
325901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdta_dyn in reference namelist', lwp )
326
327      REWIND( numnam_cfg )              ! Namelist namdta_dyn in configuration namelist : Offline: init. of dynamical data
328      READ  ( numnam_cfg, namdta_dyn, IOSTAT = ios, ERR = 902 )
329902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdta_dyn in configuration namelist', lwp )
330      IF(lwm) WRITE ( numond, namdta_dyn )
331      !                                         ! store namelist information in an array
332      !                                         ! Control print
333      IF(lwp) THEN
334         WRITE(numout,*)
335         WRITE(numout,*) 'dta_dyn : offline dynamics '
336         WRITE(numout,*) '~~~~~~~ '
337         WRITE(numout,*) '   Namelist namdta_dyn'
338         WRITE(numout,*) '      vertical velocity read from file (T) or computed (F) ln_dynwzv  = ', ln_dynwzv
339         WRITE(numout,*) '      bbl coef read from file (T) or computed (F)          ln_dynbbl  = ', ln_dynbbl
340         WRITE(numout,*) '      river runoff option enabled (T) or not (F)           ln_dynrnf  = ', ln_dynrnf
341         WRITE(numout,*)
342      ENDIF
343      !
344      IF( ln_dynbbl .AND. ( .NOT.lk_trabbl .OR. lk_c1d ) ) THEN
345         CALL ctl_warn( 'dta_dyn_init: bbl option requires key_trabbl activated ; force ln_dynbbl to false' )
346         ln_dynbbl = .FALSE.
347      ENDIF
348
349      jf_tem = 1   ;   jf_sal = 2   ;  jf_mld = 3   ;  jf_emp = 4   ;   jf_fmf  = 5   ;  jf_ice = 6   ;   jf_qsr = 7
350      jf_wnd = 8   ;   jf_uwd = 9   ;  jf_vwd = 10  ;  jf_wwd = 11  ;   jf_avt  = 12  ;  jfld  = jf_avt
351      !
352      slf_d(jf_tem) = sn_tem   ;   slf_d(jf_sal)  = sn_sal   ;   slf_d(jf_mld) = sn_mld
353      slf_d(jf_emp) = sn_emp   ;   slf_d(jf_fmf ) = sn_fmf   ;   slf_d(jf_ice) = sn_ice 
354      slf_d(jf_qsr) = sn_qsr   ;   slf_d(jf_wnd)  = sn_wnd   ;   slf_d(jf_avt) = sn_avt 
355      slf_d(jf_uwd) = sn_uwd   ;   slf_d(jf_vwd)  = sn_vwd   ;   slf_d(jf_wwd) = sn_wwd
356
357      !
358      IF( ln_dynrnf ) THEN
359                jf_rnf = jfld + 1  ;  jfld  = jf_rnf
360         slf_d(jf_rnf) = sn_rnf
361      ELSE
362         rnf (:,:) = 0._wp
363      ENDIF
364
365      IF( ln_dynbbl ) THEN         ! eiv & bbl
366                 jf_ubl  = jfld + 1 ;        jf_vbl  = jfld + 2 ;  jfld = jf_vbl
367           slf_d(jf_ubl) = sn_ubl  ;   slf_d(jf_vbl) = sn_vbl
368      ENDIF
369
370
371      ALLOCATE( sf_dyn(jfld), STAT=ierr )         ! set sf structure
372      IF( ierr > 0 ) THEN
373         CALL ctl_stop( 'dta_dyn: unable to allocate sf structure' )   ;   RETURN
374      ENDIF
375      !                                         ! fill sf with slf_i and control print
376      CALL fld_fill( sf_dyn, slf_d, cn_dir, 'dta_dyn_init', 'Data in file', 'namdta_dyn' )
377      ! Open file for each variable to get his number of dimension
378      DO ifpr = 1, jfld
379         CALL fld_clopn( sf_dyn(ifpr), nyear, nmonth, nday )
380         idv   = iom_varid( sf_dyn(ifpr)%num , slf_d(ifpr)%clvar )        ! id of the variable sdjf%clvar
381         idimv = iom_file ( sf_dyn(ifpr)%num )%ndims(idv)                 ! number of dimension for variable sdjf%clvar
382         IF( sf_dyn(ifpr)%num /= 0 )   CALL iom_close( sf_dyn(ifpr)%num ) ! close file if already open
383         ierr1=0
384         IF( idimv == 3 ) THEN    ! 2D variable
385                                      ALLOCATE( sf_dyn(ifpr)%fnow(jpi,jpj,1)    , STAT=ierr0 )
386            IF( slf_d(ifpr)%ln_tint ) ALLOCATE( sf_dyn(ifpr)%fdta(jpi,jpj,1,2)  , STAT=ierr1 )
387         ELSE                     ! 3D variable
388                                      ALLOCATE( sf_dyn(ifpr)%fnow(jpi,jpj,jpk)  , STAT=ierr0 )
389            IF( slf_d(ifpr)%ln_tint ) ALLOCATE( sf_dyn(ifpr)%fdta(jpi,jpj,jpk,2), STAT=ierr1 )
390         ENDIF
391         IF( ierr0 + ierr1 > 0 ) THEN
392            CALL ctl_stop( 'dta_dyn_init : unable to allocate sf_dyn array structure' )   ;   RETURN
393         ENDIF
394      END DO
395      !
396      IF( l_ldfslp .AND. .NOT.lk_c1d ) THEN                  ! slopes
397         IF( sf_dyn(jf_tem)%ln_tint ) THEN      ! time interpolation
398            ALLOCATE( uslpdta (jpi,jpj,jpk,2), vslpdta (jpi,jpj,jpk,2),    &
399            &         wslpidta(jpi,jpj,jpk,2), wslpjdta(jpi,jpj,jpk,2), STAT=ierr2 )
400         ELSE
401            ALLOCATE( uslpnow (jpi,jpj,jpk)  , vslpnow (jpi,jpj,jpk)  ,    &
402            &         wslpinow(jpi,jpj,jpk)  , wslpjnow(jpi,jpj,jpk)  , STAT=ierr2 )
403         ENDIF
404         IF( ierr2 > 0 ) THEN
405            CALL ctl_stop( 'dta_dyn_init : unable to allocate slope arrays' )   ;   RETURN
406         ENDIF
407      ENDIF
408      IF( ln_dynwzv ) THEN                  ! slopes
409         IF( sf_dyn(jf_uwd)%ln_tint ) THEN      ! time interpolation
410            ALLOCATE( wdta(jpi,jpj,jpk,2), STAT=ierr3 )
411         ELSE
412            ALLOCATE( wnow(jpi,jpj,jpk)  , STAT=ierr3 )
413         ENDIF
414         IF( ierr3 > 0 ) THEN
415            CALL ctl_stop( 'dta_dyn_init : unable to allocate wdta arrays' )   ;   RETURN
416         ENDIF
417      ENDIF
418      !
419      CALL dta_dyn( nit000 )
420      !
421   END SUBROUTINE dta_dyn_init
422
423   SUBROUTINE dta_dyn_wzv( pu, pv, pw )
424      !!----------------------------------------------------------------------
425      !!                    ***  ROUTINE wzv  ***
426      !!
427      !! ** Purpose :   Compute the now vertical velocity after the array swap
428      !!
429      !! ** Method  : - compute the now divergence given by :
430      !!         * z-coordinate ONLY !!!!
431      !!         hdiv = 1/(e1t*e2t) [ di(e2u  u) + dj(e1v  v) ]
432      !!     - Using the incompressibility hypothesis, the vertical
433      !!      velocity is computed by integrating the horizontal divergence
434      !!      from the bottom to the surface.
435      !!        The boundary conditions are w=0 at the bottom (no flux).
436      !!----------------------------------------------------------------------
437      USE oce, ONLY:  zhdiv => hdivn
438      !
439      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) :: pu, pv    !:  horizontal velocities
440      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(  out) :: pw        !:  vertical velocity
441      !!
442      INTEGER  ::  ji, jj, jk
443      REAL(wp) ::  zu, zu1, zv, zv1, zet
444      !!----------------------------------------------------------------------
445      !
446      ! Computation of vertical velocity using horizontal divergence
447      zhdiv(:,:,:) = 0._wp
448      DO jk = 1, jpkm1
449         DO jj = 2, jpjm1
450            DO ji = fs_2, fs_jpim1   ! vector opt.
451               zu  = pu(ji  ,jj  ,jk) * umask(ji  ,jj  ,jk) * e2u(ji  ,jj  ) * fse3u(ji  ,jj  ,jk)
452               zu1 = pu(ji-1,jj  ,jk) * umask(ji-1,jj  ,jk) * e2u(ji-1,jj  ) * fse3u(ji-1,jj  ,jk)
453               zv  = pv(ji  ,jj  ,jk) * vmask(ji  ,jj  ,jk) * e1v(ji  ,jj  ) * fse3v(ji  ,jj  ,jk)
454               zv1 = pv(ji  ,jj-1,jk) * vmask(ji  ,jj-1,jk) * e1v(ji  ,jj-1) * fse3v(ji  ,jj-1,jk)
455               zet = 1. / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) )
456               zhdiv(ji,jj,jk) = ( zu - zu1 + zv - zv1 ) * zet 
457            END DO
458         END DO
459      END DO
460      !                              !  update the horizontal divergence with the runoff inflow
461      IF( ln_dynrnf )  zhdiv(:,:,1) = zhdiv(:,:,1) - rnf(:,:) * r1_rau0 / fse3t(:,:,1)
462      !
463      CALL lbc_lnk( zhdiv, 'T', 1. )      ! Lateral boundary conditions on zhdiv
464      ! computation of vertical velocity from the bottom
465      pw(:,:,jpk) = 0._wp
466      DO jk = jpkm1, 1, -1
467         pw(:,:,jk) = pw(:,:,jk+1) - fse3t(:,:,jk) * zhdiv(:,:,jk)
468      END DO
469      !
470   END SUBROUTINE dta_dyn_wzv
471
472   SUBROUTINE dta_dyn_slp( kt, pts, puslp, pvslp, pwslpi, pwslpj )
473      !!---------------------------------------------------------------------
474      !!                    ***  ROUTINE dta_dyn_slp  ***
475      !!
476      !! ** Purpose : Computation of slope
477      !!
478      !!---------------------------------------------------------------------
479      INTEGER ,                              INTENT(in ) :: kt       ! time step
480      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: pts      ! temperature/salinity
481      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(out) :: puslp    ! zonal isopycnal slopes
482      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(out) :: pvslp    ! meridional isopycnal slopes
483      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(out) :: pwslpi   ! zonal diapycnal slopes
484      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(out) :: pwslpj   ! meridional diapycnal slopes
485      !!---------------------------------------------------------------------
486      IF( l_ldfslp .AND. .NOT.lk_c1d ) THEN    ! Computes slopes (here avt is used as workspace)                       
487         CALL eos    ( pts, rhd, rhop, gdept_0(:,:,:) )
488         CALL eos_rab( pts, rab_n )       ! now local thermal/haline expension ratio at T-points
489         CALL bn2    ( pts, rab_n, rn2  ) ! now    Brunt-Vaisala
490
491         ! Partial steps: before Horizontal DErivative
492         IF( ln_zps  .AND. .NOT. ln_isfcav)                            &
493            &            CALL zps_hde    ( kt, jpts, pts, gtsu, gtsv,  &  ! Partial steps: before horizontal gradient
494            &                                        rhd, gru , grv    )  ! of t, s, rd at the last ocean level
495         IF( ln_zps .AND.        ln_isfcav)                            &
496            &            CALL zps_hde_isf( kt, jpts, pts, gtsu, gtsv,  &    ! Partial steps for top cell (ISF)
497            &                                        rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   &
498            &                                 gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi    ) ! of t, s, rd at the first ocean level
499
500         rn2b(:,:,:) = rn2(:,:,:)         ! need for zdfmxl
501         CALL zdf_mxl( kt )            ! mixed layer depth
502         CALL ldf_slp( kt, rhd, rn2 )  ! slopes
503         puslp (:,:,:) = uslp (:,:,:) 
504         pvslp (:,:,:) = vslp (:,:,:) 
505         pwslpi(:,:,:) = wslpi(:,:,:) 
506         pwslpj(:,:,:) = wslpj(:,:,:) 
507     ELSE
508         puslp (:,:,:) = 0.            ! to avoid warning when compiling
509         pvslp (:,:,:) = 0.
510         pwslpi(:,:,:) = 0.
511         pwslpj(:,:,:) = 0.
512     ENDIF
513      !
514   END SUBROUTINE dta_dyn_slp
515   !!======================================================================
516END MODULE dtadyn
Note: See TracBrowser for help on using the repository browser.