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

source: trunk/NEMOGCM/NEMO/OFF_SRC/dtadyn.F90 @ 4624

Last change on this file since 4624 was 4624, checked in by acc, 10 years ago

#1305. Fix slow start-up problems on some systems by introducing and using lwm logical to restrict output of merged namelists to the first (or only) processor. lwm is true only on the first processor regardless of ln_ctl. Small changes to all flavours of nemogcm.F90 are also required to write namctl and namcfg after the call to mynode which now opens output.namelist.dyn and writes nammpp.

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