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/2014/dev_r4704_NOC5_MPP_BDY_UPDATE/NEMOGCM/NEMO/OFF_SRC – NEMO

source: branches/2014/dev_r4704_NOC5_MPP_BDY_UPDATE/NEMOGCM/NEMO/OFF_SRC/dtadyn.F90 @ 6225

Last change on this file since 6225 was 6225, checked in by jamesharle, 8 years ago

Update MPP_BDY_UPDATE branch to be consistent with head of trunk

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