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/2012/dev_r3379_CMCC6_topbfm/NEMOGCM/NEMO/OFF_SRC – NEMO

source: branches/2012/dev_r3379_CMCC6_topbfm/NEMOGCM/NEMO/OFF_SRC/dtadyn.F90 @ 3561

Last change on this file since 3561 was 3561, checked in by vichi, 11 years ago

Corrected bug in directory of ocean input files

It was not possible to read ocean input files from an external directory.
Also updated the header of ocean.output with the list of institutes
  • Property svn:keywords set to Id
File size: 31.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  = .true.  !: vertical velocity read in a file (T) or computed from u/v (F)
51   LOGICAL            ::   ln_dynbbl  = .true.  !: bbl coef read in a file (T) or computed (F)
52   LOGICAL            ::   ln_degrad  = .false. !: degradation option enabled or not
53
54   INTEGER  , PARAMETER ::   jpfld = 19     ! maximum number of files 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
75   TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_dyn  ! structure of input fields (file informations, fields read)
76   !                                               !
77   REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: wdta       ! vertical velocity at 2 time step
78   REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:  ) :: wnow       ! vertical velocity at 2 time step
79   REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: uslpdta    ! zonal isopycnal slopes
80   REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: vslpdta    ! meridional isopycnal slopes
81   REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: wslpidta   ! zonal diapycnal slopes
82   REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: wslpjdta   ! meridional diapycnal slopes
83   REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:)   :: uslpnow    ! zonal isopycnal slopes
84   REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:)   :: vslpnow    ! meridional isopycnal slopes
85   REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:)   :: wslpinow   ! zonal diapycnal slopes
86   REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:)   :: wslpjnow   ! meridional diapycnal slopes
87
88   INTEGER :: nrecprev_tem , nrecprev_uwd
89
90   !! * Substitutions
91#  include "domzgr_substitute.h90"
92#  include "vectopt_loop_substitute.h90"
93   !!----------------------------------------------------------------------
94   !! NEMO/OFF 3.3 , NEMO Consortium (2010)
95   !! $Id$
96   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
97   !!----------------------------------------------------------------------
98CONTAINS
99
100   SUBROUTINE dta_dyn( kt )
101      !!----------------------------------------------------------------------
102      !!                  ***  ROUTINE dta_dyn  ***
103      !!
104      !! ** Purpose :  Prepares dynamics and physics fields from a NEMO run
105      !!               for an off-line simulation of passive tracers
106      !!
107      !! ** Method : calculates the position of data
108      !!             - computes slopes if needed
109      !!             - interpolates data if needed
110      !!----------------------------------------------------------------------
111      !
112      USE oce, ONLY:  zts    => tsa 
113      USE oce, ONLY:  zuslp  => ua   , zvslp  => va
114      USE oce, ONLY:  zwslpi => rotb , zwslpj => rotn
115      USE oce, ONLY:  zu     => ub   , zv     => vb,  zw => hdivb
116      !
117      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
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( lk_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( lk_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            write(numout,*)
166            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            write(numout,*)
216            write(numout,*) ' Compute new vertical velocity at kt = ', kt
217            write(numout,*)
218            IF( sf_dyn(jf_uwd)%ln_tint ) THEN                 ! time interpolation of data
219               IF( kt /= nit000 )  THEN
220                  wdta(:,:,:,1) =  wdta(:,:,:,2)     ! swap the data for initialisation
221               ENDIF
222               zu(:,:,:) = sf_dyn(jf_uwd)%fdta(:,:,:,2)
223               zv(:,:,:) = sf_dyn(jf_vwd)%fdta(:,:,:,2)
224               CALL dta_dyn_wzv( zu, zv, zw )
225               wdta(:,:,:,2) = zw(:,:,:) * tmask(:,:,:)
226            ELSE
227               zu(:,:,:) = sf_dyn(jf_uwd)%fnow(:,:,:) 
228               zv(:,:,:) = sf_dyn(jf_vwd)%fnow(:,:,:)
229               CALL dta_dyn_wzv( zu, zv, zw )
230               wnow(:,:,:)  = zw(:,:,:) * tmask(:,:,:)
231            ENDIF
232         ENDIF
233         IF( sf_dyn(jf_uwd)%ln_tint )  THEN
234            ztinta =  REAL( isecsbc - sf_dyn(jf_uwd)%nrec_b(2), wp )  &
235               &    / REAL( sf_dyn(jf_uwd)%nrec_a(2) - sf_dyn(jf_uwd)%nrec_b(2), wp )
236            ztintb =  1. - ztinta
237            wn(:,:,:) = ztintb * wdta(:,:,:,1)  + ztinta * wdta(:,:,:,2) 
238         ELSE
239            wn(:,:,:) = wnow(:,:,:)
240         ENDIF
241      ENDIF
242      !
243      tsn(:,:,:,jp_tem) = sf_dyn(jf_tem)%fnow(:,:,:) * tmask(:,:,:)    ! temperature
244      tsn(:,:,:,jp_sal) = sf_dyn(jf_sal)%fnow(:,:,:) * tmask(:,:,:)    ! salinity
245      !
246      CALL eos    ( tsn, rhd, rhop )                                       ! In any case, we need rhop
247      CALL zdf_mxl( kt )                                                   ! In any case, we need mxl
248      !
249      avt(:,:,:)       = sf_dyn(jf_avt)%fnow(:,:,:) * tmask(:,:,:)    ! vertical diffusive coefficient
250      un (:,:,:)       = sf_dyn(jf_uwd)%fnow(:,:,:) * umask(:,:,:)    ! u-velocity
251      vn (:,:,:)       = sf_dyn(jf_vwd)%fnow(:,:,:) * vmask(:,:,:)    ! v-velocity
252      IF( .NOT.ln_dynwzv ) &                                           ! w-velocity read in file
253         wn (:,:,:)    = sf_dyn(jf_wwd)%fnow(:,:,:) * tmask(:,:,:)   
254      hmld(:,:)        = sf_dyn(jf_mld)%fnow(:,:,1) * tmask(:,:,1)    ! mixed layer depht
255      wndm(:,:)        = sf_dyn(jf_wnd)%fnow(:,:,1) * tmask(:,:,1)    ! wind speed - needed for gas exchange
256      emp (:,:)        = sf_dyn(jf_emp)%fnow(:,:,1) * tmask(:,:,1)    ! E-P
257      emps(:,:)        = emp(:,:) 
258      fr_i(:,:)        = sf_dyn(jf_ice)%fnow(:,:,1) * tmask(:,:,1)     ! Sea-ice fraction
259      qsr (:,:)        = sf_dyn(jf_qsr)%fnow(:,:,1) * tmask(:,:,1)    ! solar radiation
260
261      !                                                      ! bbl diffusive coef
262#if defined key_trabbl && ! defined key_c1d
263      IF( ln_dynbbl ) THEN                                        ! read in a file
264         ahu_bbl(:,:)  = sf_dyn(jf_ubl)%fnow(:,:,1) * umask(:,:,1)
265         ahv_bbl(:,:)  = sf_dyn(jf_vbl)%fnow(:,:,1) * vmask(:,:,1)
266      ELSE                                                        ! Compute bbl coefficients if needed
267         tsb(:,:,:,:) = tsn(:,:,:,:)
268         CALL bbl( kt, nit000, 'TRC')
269      END IF
270#endif
271#if ( ! defined key_degrad && defined key_traldf_c2d && defined key_traldf_eiv ) && ! defined key_c1d
272      aeiw(:,:)        = sf_dyn(jf_eiw)%fnow(:,:,1) * tmask(:,:,1)    ! w-eiv
273      !                                                           ! Computes the horizontal values from the vertical value
274      DO jj = 2, jpjm1
275         DO ji = fs_2, fs_jpim1   ! vector opt.
276            aeiu(ji,jj) = .5 * ( aeiw(ji,jj) + aeiw(ji+1,jj  ) )  ! Average the diffusive coefficient at u- v- points
277            aeiv(ji,jj) = .5 * ( aeiw(ji,jj) + aeiw(ji  ,jj+1) )  ! at u- v- points
278         END DO
279      END DO
280      CALL lbc_lnk( aeiu, 'U', 1. )   ;   CALL lbc_lnk( aeiv, 'V', 1. )    ! lateral boundary condition
281#endif
282     
283#if defined key_degrad && ! defined key_c1d
284      !                                          ! degrad option : diffusive and eiv coef are 3D
285      ahtu(:,:,:) = sf_dyn(jf_ahu)%fnow(:,:,:) * umask(:,:,:)
286      ahtv(:,:,:) = sf_dyn(jf_ahv)%fnow(:,:,:) * vmask(:,:,:)
287      ahtw(:,:,:) = sf_dyn(jf_ahw)%fnow(:,:,:) * tmask(:,:,:)
288#  if defined key_traldf_eiv 
289      aeiu(:,:,:) = sf_dyn(jf_eiu)%fnow(:,:,:) * umask(:,:,:)
290      aeiv(:,:,:) = sf_dyn(jf_eiv)%fnow(:,:,:) * vmask(:,:,:)
291      aeiw(:,:,:) = sf_dyn(jf_eiw)%fnow(:,:,:) * tmask(:,:,:)
292#  endif
293#endif
294      !
295      IF(ln_ctl) THEN                  ! print control
296         CALL prt_ctl(tab3d_1=tsn(:,:,:,jp_tem), clinfo1=' tn      - : ', mask1=tmask, ovlap=1, kdim=jpk   )
297         CALL prt_ctl(tab3d_1=tsn(:,:,:,jp_sal), clinfo1=' sn      - : ', mask1=tmask, ovlap=1, kdim=jpk   )
298         CALL prt_ctl(tab3d_1=un               , clinfo1=' un      - : ', mask1=umask, ovlap=1, kdim=jpk   )
299         CALL prt_ctl(tab3d_1=vn               , clinfo1=' vn      - : ', mask1=vmask, ovlap=1, kdim=jpk   )
300         CALL prt_ctl(tab3d_1=wn               , clinfo1=' wn      - : ', mask1=tmask, ovlap=1, kdim=jpk   )
301         CALL prt_ctl(tab3d_1=avt              , clinfo1=' kz      - : ', mask1=tmask, ovlap=1, kdim=jpk   )
302         CALL prt_ctl(tab2d_1=fr_i             , clinfo1=' fr_i    - : ', mask1=tmask, ovlap=1 )
303         CALL prt_ctl(tab2d_1=hmld             , clinfo1=' hmld    - : ', mask1=tmask, ovlap=1 )
304         CALL prt_ctl(tab2d_1=emps             , clinfo1=' emps    - : ', 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      !!
328      CHARACTER(len=100)            ::  cn_dir   !   Root directory for location of core files
329      TYPE(FLD_N), DIMENSION(jpfld) ::  slf_d    ! array of namelist informations on the fields to read
330      TYPE(FLD_N) :: sn_tem, sn_sal, sn_mld, sn_emp, sn_ice, sn_qsr, sn_wnd  ! informations about the fields to be read
331      TYPE(FLD_N) :: sn_uwd, sn_vwd, sn_wwd, sn_avt, sn_ubl, sn_vbl          !   "                                 "
332      TYPE(FLD_N) :: sn_ahu, sn_ahv, sn_ahw, sn_eiu, sn_eiv, sn_eiw          !   "                                 "
333      !
334      NAMELIST/namdta_dyn/cn_dir, ln_dynwzv, ln_dynbbl, ln_degrad,    &
335         &                sn_tem, sn_sal, sn_mld, sn_emp, sn_ice, sn_qsr, sn_wnd,  &
336         &                sn_uwd, sn_vwd, sn_wwd, sn_avt, sn_ubl, sn_vbl,          &
337         &                sn_ahu, sn_ahv, sn_ahw, sn_eiu, sn_eiv, sn_eiw
338
339      !!----------------------------------------------------------------------
340      !                                   ! ============
341      !                                   !   Namelist
342      !                                   ! ============
343      ! (NB: frequency positive => hours, negative => months)
344      !                !   file      ! frequency !  variable  ! time intep !  clim  ! 'yearly' or ! weights  ! rotation   !
345      !                !   name      !  (hours)  !   name     !   (T/F)    !  (T/F) !  'monthly'  ! filename ! pairs      !
346      sn_tem  = FLD_N( 'dyna_grid_T' ,    120    , 'votemper' ,  .true.    , .true. ,   'yearly'  , ''       , ''         )
347      sn_sal  = FLD_N( 'dyna_grid_T' ,    120    , 'vosaline' ,  .true.    , .true. ,   'yearly'  , ''       , ''         )
348      sn_mld  = FLD_N( 'dyna_grid_T' ,    120    , 'somixght' ,  .true.    , .true. ,   'yearly'  , ''       , ''         )
349      sn_emp  = FLD_N( 'dyna_grid_T' ,    120    , 'sowaflcd' ,  .true.    , .true. ,   'yearly'  , ''       , ''         )
350      sn_ice  = FLD_N( 'dyna_grid_T' ,    120    , 'soicecov' ,  .true.    , .true. ,   'yearly'  , ''       , ''         )
351      sn_qsr  = FLD_N( 'dyna_grid_T' ,    120    , 'soshfldo' ,  .true.    , .true. ,   'yearly'  , ''       , ''         )
352      sn_wnd  = FLD_N( 'dyna_grid_T' ,    120    , 'sowindsp' ,  .true.    , .true. ,   'yearly'  , ''       , ''         )
353      sn_uwd  = FLD_N( 'dyna_grid_U' ,    120    , 'vozocrtx' ,  .true.    , .true. ,   'yearly'  , ''       , ''         )
354      sn_vwd  = FLD_N( 'dyna_grid_V' ,    120    , 'vomecrty' ,  .true.    , .true. ,   'yearly'  , ''       , ''         )
355      sn_wwd  = FLD_N( 'dyna_grid_W' ,    120    , 'vovecrtz' ,  .true.    , .true. ,   'yearly'  , ''       , ''         )
356      sn_avt  = FLD_N( 'dyna_grid_W' ,    120    , 'votkeavt' ,  .true.    , .true. ,   'yearly'  , ''       , ''         )
357      sn_ubl  = FLD_N( 'dyna_grid_U' ,    120    , 'sobblcox' ,  .true.    , .true. ,   'yearly'  , ''       , ''         )
358      sn_vbl  = FLD_N( 'dyna_grid_V' ,    120    , 'sobblcoy' ,  .true.    , .true. ,   'yearly'  , ''       , ''         )
359      sn_ahu  = FLD_N( 'dyna_grid_U' ,    120    , 'vozoahtu' ,  .true.    , .true. ,   'yearly'  , ''       , ''         )
360      sn_ahv  = FLD_N( 'dyna_grid_V' ,    120    , 'vomeahtv' ,  .true.    , .true. ,   'yearly'  , ''       , ''         )
361      sn_ahw  = FLD_N( 'dyna_grid_W' ,    120    , 'voveahtz' ,  .true.    , .true. ,   'yearly'  , ''       , ''         )
362      sn_eiu  = FLD_N( 'dyna_grid_U' ,    120    , 'vozoaeiu' ,  .true.    , .true. ,   'yearly'  , ''       , ''         )
363      sn_eiv  = FLD_N( 'dyna_grid_V' ,    120    , 'vomeaeiv' ,  .true.    , .true. ,   'yearly'  , ''       , ''         )
364      sn_eiw  = FLD_N( 'dyna_grid_W' ,    120    , 'voveaeiw' ,  .true.    , .true. ,   'yearly'  , ''       , ''         )
365      !
366      REWIND( numnam )                          ! read in namlist namdta_dyn
367      READ  ( numnam, namdta_dyn )
368      !                                         ! store namelist information in an array
369      !                                         ! Control print
370      IF(lwp) THEN
371         WRITE(numout,*)
372         WRITE(numout,*) 'dta_dyn : offline dynamics '
373         WRITE(numout,*) '~~~~~~~ '
374         WRITE(numout,*) '   Namelist namdta_dyn'
375         WRITE(numout,*) '      vertical velocity read from file (T) or computed (F) ln_dynwzv  = ', ln_dynwzv
376         WRITE(numout,*) '      bbl coef read from file (T) or computed (F)          ln_dynbbl  = ', ln_dynbbl
377         WRITE(numout,*) '      degradation option enabled (T) or not (F)            ln_degrad  = ', ln_degrad
378         WRITE(numout,*)
379      ENDIF
380      !
381      IF( ln_degrad .AND. .NOT.lk_degrad ) THEN
382         CALL ctl_warn( 'dta_dyn_init: degradation option requires key_degrad activated ; force ln_degrad to false' )
383         ln_degrad = .FALSE.
384      ENDIF
385      IF( ln_dynbbl .AND. ( .NOT.lk_trabbl .OR. lk_c1d ) ) THEN
386         CALL ctl_warn( 'dta_dyn_init: bbl option requires key_trabbl activated ; force ln_dynbbl to false' )
387         ln_dynbbl = .FALSE.
388      ENDIF
389
390      jf_tem = 1   ;   jf_sal = 2   ;  jf_mld = 3   ;  jf_emp = 4   ;   jf_ice = 5   ;   jf_qsr = 6 
391      jf_wnd = 7   ;   jf_uwd = 8   ;  jf_vwd = 9   ;  jf_wwd = 10  ;   jf_avt = 11  ;   jfld  = 11
392      !
393      slf_d(jf_tem) = sn_tem   ;   slf_d(jf_sal) = sn_sal   ;   slf_d(jf_mld) = sn_mld
394      slf_d(jf_emp) = sn_emp   ;   slf_d(jf_ice) = sn_ice   ;   slf_d(jf_qsr) = sn_qsr
395      slf_d(jf_wnd) = sn_wnd   ;   slf_d(jf_uwd) = sn_uwd   ;   slf_d(jf_vwd) = sn_vwd
396      slf_d(jf_wwd) = sn_wwd   ;   slf_d(jf_avt) = sn_avt 
397      !
398      IF( .NOT.ln_degrad ) THEN     ! no degrad option
399         IF( lk_traldf_eiv .AND. ln_dynbbl ) THEN        ! eiv & bbl
400                 jf_ubl  = 12      ;         jf_vbl  = 13      ;         jf_eiw  = 14   ;   jfld = 14
401           slf_d(jf_ubl) = sn_ubl  ;   slf_d(jf_vbl) = sn_vbl  ;   slf_d(jf_eiw) = sn_eiw
402         ENDIF
403         IF( .NOT.lk_traldf_eiv .AND. ln_dynbbl ) THEN   ! no eiv & bbl
404                 jf_ubl  = 12      ;         jf_vbl  = 13      ;   jfld = 13
405           slf_d(jf_ubl) = sn_ubl  ;   slf_d(jf_vbl) = sn_vbl
406         ENDIF
407         IF( lk_traldf_eiv .AND. .NOT.ln_dynbbl ) THEN   ! eiv & no bbl
408           jf_eiw = 12   ;   jfld = 12   ;   slf_d(jf_eiw) = sn_eiw
409         ENDIF
410      ELSE
411              jf_ahu  = 12      ;         jf_ahv  = 13      ;         jf_ahw  = 14   ;   jfld = 14
412        slf_d(jf_ahu) = sn_ahu  ;   slf_d(jf_ahv) = sn_ahv  ;   slf_d(jf_ahw) = sn_ahw
413        IF( lk_traldf_eiv .AND. ln_dynbbl ) THEN         ! eiv & bbl
414                 jf_ubl  = 15      ;         jf_vbl  = 16     
415           slf_d(jf_ubl) = sn_ubl  ;   slf_d(jf_vbl) = sn_vbl 
416                 jf_eiu  = 17      ;         jf_eiv  = 18      ;          jf_eiw  = 19   ;   jfld = 19
417           slf_d(jf_eiu) = sn_eiu  ;   slf_d(jf_eiv) = sn_eiv  ;    slf_d(jf_eiw) = sn_eiw
418        ENDIF
419        IF( .NOT.lk_traldf_eiv .AND. ln_dynbbl ) THEN    ! no eiv & bbl
420                 jf_ubl  = 15      ;         jf_vbl  = 16      ;   jfld = 16
421           slf_d(jf_ubl) = sn_ubl  ;   slf_d(jf_vbl) = sn_vbl
422        ENDIF
423        IF( lk_traldf_eiv .AND. .NOT.ln_dynbbl ) THEN    ! eiv & no bbl
424                 jf_eiu  = 15      ;         jf_eiv  = 16      ;         jf_eiw  = 17   ;   jfld = 17
425           slf_d(jf_eiu) = sn_eiu  ;   slf_d(jf_eiv) = sn_eiv  ;   slf_d(jf_eiw) = sn_eiw
426        ENDIF
427      ENDIF
428 
429      ALLOCATE( sf_dyn(jfld), STAT=ierr )         ! set sf structure
430      IF( ierr > 0 ) THEN
431         CALL ctl_stop( 'dta_dyn: unable to allocate sf structure' )   ;   RETURN
432      ENDIF
433      ! Open file for each variable to get his number of dimension
434      DO ifpr = 1, jfld
435         CALL iom_open( TRIM( cn_dir )//TRIM( slf_d(ifpr)%clname ), inum )
436         idv   = iom_varid( inum , slf_d(ifpr)%clvar )  ! id of the variable sdjf%clvar
437         idimv = iom_file ( inum )%ndims(idv)             ! number of dimension for variable sdjf%clvar
438         IF( inum /= 0 )   CALL iom_close( inum )       ! close file if already open
439         IF( idimv == 3 ) THEN    ! 2D variable
440                                      ALLOCATE( sf_dyn(ifpr)%fnow(jpi,jpj,1)    , STAT=ierr0 )
441            IF( slf_d(ifpr)%ln_tint ) ALLOCATE( sf_dyn(ifpr)%fdta(jpi,jpj,1,2)  , STAT=ierr1 )
442         ELSE                     ! 3D variable
443                                      ALLOCATE( sf_dyn(ifpr)%fnow(jpi,jpj,jpk)  , STAT=ierr0 )
444            IF( slf_d(ifpr)%ln_tint ) ALLOCATE( sf_dyn(ifpr)%fdta(jpi,jpj,jpk,2), STAT=ierr1 )
445         ENDIF
446         IF( ierr0 + ierr1 > 0 ) THEN
447            CALL ctl_stop( 'dta_dyn_init : unable to allocate sf_dyn array structure' )   ;   RETURN
448         ENDIF
449      END DO
450      !                                         ! fill sf with slf_i and control print
451      CALL fld_fill( sf_dyn, slf_d, cn_dir, 'dta_dyn_init', 'Data in file', 'namdta_dyn' )
452      !
453      IF( lk_ldfslp .AND. .NOT.lk_c1d ) THEN                  ! slopes
454         IF( sf_dyn(jf_tem)%ln_tint ) THEN      ! time interpolation
455            ALLOCATE( uslpdta (jpi,jpj,jpk,2), vslpdta (jpi,jpj,jpk,2),    &
456            &         wslpidta(jpi,jpj,jpk,2), wslpjdta(jpi,jpj,jpk,2), STAT=ierr2 )
457         ELSE
458            ALLOCATE( uslpnow (jpi,jpj,jpk)  , vslpnow (jpi,jpj,jpk)  ,    &
459            &         wslpinow(jpi,jpj,jpk)  , wslpjnow(jpi,jpj,jpk)  , STAT=ierr2 )
460         ENDIF
461         IF( ierr2 > 0 ) THEN
462            CALL ctl_stop( 'dta_dyn_init : unable to allocate slope arrays' )   ;   RETURN
463         ENDIF
464      ENDIF
465      IF( ln_dynwzv ) THEN                  ! slopes
466         IF( sf_dyn(jf_uwd)%ln_tint ) THEN      ! time interpolation
467            ALLOCATE( wdta(jpi,jpj,jpk,2), STAT=ierr3 )
468         ELSE
469            ALLOCATE( wnow(jpi,jpj,jpk)  , STAT=ierr3 )
470         ENDIF
471         IF( ierr3 > 0 ) THEN
472            CALL ctl_stop( 'dta_dyn_init : unable to allocate wdta arrays' )   ;   RETURN
473         ENDIF
474      ENDIF
475      !
476      CALL dta_dyn( nit000 )
477      !
478   END SUBROUTINE dta_dyn_init
479
480   SUBROUTINE dta_dyn_wzv( pu, pv, pw )
481      !!----------------------------------------------------------------------
482      !!                    ***  ROUTINE wzv  ***
483      !!
484      !! ** Purpose :   Compute the now vertical velocity after the array swap
485      !!
486      !! ** Method  : - compute the now divergence given by :
487      !!         * z-coordinate ONLY !!!!
488      !!         hdiv = 1/(e1t*e2t) [ di(e2u  u) + dj(e1v  v) ]
489      !!     - Using the incompressibility hypothesis, the vertical
490      !!      velocity is computed by integrating the horizontal divergence
491      !!      from the bottom to the surface.
492      !!        The boundary conditions are w=0 at the bottom (no flux).
493      !!----------------------------------------------------------------------
494      USE oce, ONLY:  zhdiv => hdivn
495      !
496      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) :: pu, pv    !:  horizontal velocities
497      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(  out) :: pw        !:  vertical velocity
498      !!
499      INTEGER  ::  ji, jj, jk
500      REAL(wp) ::  zu, zu1, zv, zv1, zet
501      !!----------------------------------------------------------------------
502      !
503      ! Computation of vertical velocity using horizontal divergence
504      zhdiv(:,:,:) = 0._wp
505      DO jk = 1, jpkm1
506         DO jj = 2, jpjm1
507            DO ji = fs_2, fs_jpim1   ! vector opt.
508               zu  = pu(ji  ,jj  ,jk) * umask(ji  ,jj  ,jk) * e2u(ji  ,jj  ) * fse3u(ji  ,jj  ,jk)
509               zu1 = pu(ji-1,jj  ,jk) * umask(ji-1,jj  ,jk) * e2u(ji-1,jj  ) * fse3u(ji-1,jj  ,jk)
510               zv  = pv(ji  ,jj  ,jk) * vmask(ji  ,jj  ,jk) * e1v(ji  ,jj  ) * fse3v(ji  ,jj  ,jk)
511               zv1 = pv(ji  ,jj-1,jk) * vmask(ji  ,jj-1,jk) * e1v(ji  ,jj-1) * fse3v(ji  ,jj-1,jk)
512               zet = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
513               zhdiv(ji,jj,jk) = ( zu - zu1 + zv - zv1 ) * zet 
514            END DO
515         END DO
516      END DO
517      CALL lbc_lnk( zhdiv, 'T', 1. )      ! Lateral boundary conditions on zhdiv
518      !
519      ! computation of vertical velocity from the bottom
520      pw(:,:,jpk) = 0._wp
521      DO jk = jpkm1, 1, -1
522         pw(:,:,jk) = pw(:,:,jk+1) - fse3t(:,:,jk) * zhdiv(:,:,jk)
523      END DO
524      !
525   END SUBROUTINE dta_dyn_wzv
526
527   SUBROUTINE dta_dyn_slp( kt, pts, puslp, pvslp, pwslpi, pwslpj )
528      !!---------------------------------------------------------------------
529      !!                    ***  ROUTINE dta_dyn_slp  ***
530      !!
531      !! ** Purpose : Computation of slope
532      !!
533      !!---------------------------------------------------------------------
534      INTEGER ,                              INTENT(in ) :: kt       ! time step
535      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: pts      ! temperature/salinity
536      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(out) :: puslp    ! zonal isopycnal slopes
537      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(out) :: pvslp    ! meridional isopycnal slopes
538      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(out) :: pwslpi   ! zonal diapycnal slopes
539      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(out) :: pwslpj   ! meridional diapycnal slopes
540      !!---------------------------------------------------------------------
541#if defined key_ldfslp && ! defined key_c1d
542      CALL eos( pts, rhd, rhop )   ! Time-filtered in situ density
543      CALL bn2( pts, rn2 )         ! before Brunt-Vaisala frequency
544      IF( ln_zps )   &
545         &  CALL zps_hde( kt, jpts, pts, gtsu, gtsv, rhd, gru, grv )  ! Partial steps: before Horizontal DErivative
546         !                                                            ! of t, s, rd at the bottom ocean level
547      CALL zdf_mxl( kt )            ! mixed layer depth
548      CALL ldf_slp( kt, rhd, rn2 )  ! slopes
549      puslp (:,:,:) = uslp (:,:,:) 
550      pvslp (:,:,:) = vslp (:,:,:) 
551      pwslpi(:,:,:) = wslpi(:,:,:) 
552      pwslpj(:,:,:) = wslpj(:,:,:) 
553#else
554      puslp (:,:,:) = 0.            ! to avoid warning when compiling
555      pvslp (:,:,:) = 0.
556      pwslpi(:,:,:) = 0.
557      pwslpj(:,:,:) = 0.
558#endif
559      !
560   END SUBROUTINE dta_dyn_slp
561   !!======================================================================
562END MODULE dtadyn
Note: See TracBrowser for help on using the repository browser.