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

source: branches/UKMO/dev_r10171_test_crs_AMM7/NEMOGCM/NEMO/OFF_SRC/dtadyn.F90 @ 10207

Last change on this file since 10207 was 10207, checked in by cmao, 6 years ago

remove svn keyword

File size: 35.7 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 domvvl          ! variable volume
25   USE zdf_oce         ! ocean vertical physics: variables
26   USE sbc_oce         ! surface module: variables
27   USE trc_oce         ! share ocean/biogeo variables
28   USE phycst          ! physical constants
29   USE trabbl          ! active tracer: bottom boundary layer
30   USE ldfslp          ! lateral diffusion: iso-neutral slopes
31   USE sbcrnf          ! river runoffs
32   USE ldfeiv          ! eddy induced velocity coef.
33   USE ldftra_oce      ! ocean tracer   lateral physics
34   USE zdfmxl          ! vertical physics: mixed layer depth
35   USE eosbn2          ! equation of state - Brunt Vaisala frequency
36   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
37   USE zpshde          ! z-coord. with partial steps: horizontal derivatives
38   USE in_out_manager  ! I/O manager
39   USE iom             ! I/O library
40   USE lib_mpp         ! distributed memory computing library
41   USE prtctl          ! print control
42   USE fldread         ! read input fields
43   USE wrk_nemo        ! Memory allocation
44   USE timing          ! Timing
45   USE trc, ONLY : ln_rsttr, numrtr, numrtw, lrst_trc
46
47   IMPLICIT NONE
48   PRIVATE
49
50   PUBLIC   dta_dyn_init   ! called by opa.F90
51   PUBLIC   dta_dyn        ! called by step.F90
52   PUBLIC   dta_dyn_swp   ! called by step.F90
53
54   CHARACTER(len=100) ::   cn_dir           !: Root directory for location of ssr files
55   LOGICAL            ::   ln_ssh_ini       !: initial ssh from dyn file (T) or not (F) - ssh is then read from passive tracer restart
56   LOGICAL            ::   ln_dynrnf        !: read runoff data in file (T) or set to zero (F)
57   LOGICAL            ::   ln_dynrnf_depth  !: read runoff data in file (T) or set to zero (F)
58   REAL(wp)           ::   fwbcorr
59
60
61   INTEGER  , PARAMETER ::   jpfld = 20     ! maximum number of fields to read
62   INTEGER  , SAVE      ::   jf_tem         ! index of temperature
63   INTEGER  , SAVE      ::   jf_sal         ! index of salinity
64   INTEGER  , SAVE      ::   jf_uwd         ! index of u-transport
65   INTEGER  , SAVE      ::   jf_vwd         ! index of v-transport
66   INTEGER  , SAVE      ::   jf_wwd         ! index of v-transport
67   INTEGER  , SAVE      ::   jf_avt         ! index of Kz
68   INTEGER  , SAVE      ::   jf_mld         ! index of mixed layer deptht
69   INTEGER  , SAVE      ::   jf_emp         ! index of water flux
70   INTEGER  , SAVE      ::   jf_empb        ! index of water flux
71   INTEGER  , SAVE      ::   jf_qsr         ! index of solar radiation
72   INTEGER  , SAVE      ::   jf_wnd         ! index of wind speed
73   INTEGER  , SAVE      ::   jf_ice         ! index of sea ice cover
74   INTEGER  , SAVE      ::   jf_rnf         ! index of river runoff
75   INTEGER  , SAVE      ::   jf_fmf         ! index of downward salt flux
76   INTEGER  , SAVE      ::   jf_ubl         ! index of u-bbl coef
77   INTEGER  , SAVE      ::   jf_vbl         ! index of v-bbl coef
78   INTEGER  , SAVE      ::   jf_div         ! index of e3t
79
80
81   TYPE(FLD), ALLOCATABLE, SAVE, DIMENSION(:) :: sf_dyn  ! structure of input fields (file informations, fields read)
82   !                                               !
83   REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: uslpdta    ! zonal isopycnal slopes
84   REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: vslpdta    ! meridional isopycnal slopes
85   REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: wslpidta   ! zonal diapycnal slopes
86   REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: wslpjdta   ! meridional diapycnal slopes
87
88   INTEGER, SAVE  :: nprevrec, nsecdyn
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:  zhdivtr => ua
113      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
114      INTEGER             ::   ji, jj, jk
115      REAL(wp), POINTER, DIMENSION(:,:)   :: zemp
116      !
117      !!----------------------------------------------------------------------
118     
119      !
120      IF( nn_timing == 1 )  CALL timing_start( 'dta_dyn')
121      !
122      !
123      nsecdyn = nsec_year + nsec1jan000   ! number of seconds between Jan. 1st 00h of nit000 year and the middle of time step
124      !
125      IF( kt == nit000 ) THEN    ;    nprevrec = 0
126      ELSE                       ;    nprevrec = sf_dyn(jf_tem)%nrec_a(2)
127      ENDIF
128      !
129      CALL fld_read( kt, 1, sf_dyn )      !=  read data at kt time step   ==!
130      !
131      IF( lk_ldfslp .AND. .NOT.lk_c1d )   CALL  dta_dyn_slp( kt )    ! Computation of slopes
132      !
133      tsn(:,:,:,jp_tem) = sf_dyn(jf_tem)%fnow(:,:,:)  * tmask(:,:,:)    ! temperature
134      tsn(:,:,:,jp_sal) = sf_dyn(jf_sal)%fnow(:,:,:)  * tmask(:,:,:)    ! salinity
135      wndm(:,:)         = sf_dyn(jf_wnd)%fnow(:,:,1)  * tmask(:,:,1)    ! wind speed - needed for gas exchange
136      fmmflx(:,:)       = sf_dyn(jf_fmf)%fnow(:,:,1)  * tmask(:,:,1)    ! downward salt flux (v3.5+)
137      fr_i(:,:)         = sf_dyn(jf_ice)%fnow(:,:,1)  * tmask(:,:,1)    ! Sea-ice fraction
138      qsr (:,:)         = sf_dyn(jf_qsr)%fnow(:,:,1)  * tmask(:,:,1)    ! solar radiation
139      emp (:,:)         = sf_dyn(jf_emp)%fnow(:,:,1)  * tmask(:,:,1)    ! E-P
140      IF( ln_dynrnf ) THEN
141         rnf (:,:)      = sf_dyn(jf_rnf)%fnow(:,:,1) * tmask(:,:,1)    ! E-P
142         IF( ln_dynrnf_depth .AND. lk_vvl )    CALL  dta_dyn_hrnf
143      ENDIF
144      !
145      un(:,:,:)        = sf_dyn(jf_uwd)%fnow(:,:,:) * umask(:,:,:)    ! effective u-transport
146      vn(:,:,:)        = sf_dyn(jf_vwd)%fnow(:,:,:) * vmask(:,:,:)    ! effective v-transport
147      wn(:,:,:)        = sf_dyn(jf_wwd)%fnow(:,:,:) * tmask(:,:,:)    ! effective v-transport
148      !
149      IF( lk_vvl ) THEN
150         CALL wrk_alloc(jpi, jpj, zemp )
151         zhdivtr(:,:,:)    = sf_dyn(jf_div)%fnow(:,:,:)  * tmask(:,:,:)    ! effective u-transport
152         emp_b (:,:)       = sf_dyn(jf_empb)%fnow(:,:,1) * tmask(:,:,1)    ! E-P
153         zemp(:,:) = 0.5_wp * ( emp(:,:) + emp_b(:,:) ) + rnf(:,:) + fwbcorr * tmask(:,:,1)
154         CALL dta_dyn_ssh( kt, zhdivtr, sshb, zemp, ssha, fse3t_a(:,:,:) )  !=  ssh, vertical scale factor & vertical transport
155         CALL wrk_dealloc(jpi, jpj, zemp )
156         !                                           Write in the tracer restart file
157         !                                          *******************************
158         IF( lrst_trc ) THEN
159            IF(lwp) WRITE(numout,*)
160            IF(lwp) WRITE(numout,*) 'dta_dyn_ssh : ssh field written in tracer restart file ',   &
161               &                    'at it= ', kt,' date= ', ndastp
162            IF(lwp) WRITE(numout,*) '~~~~'
163            CALL iom_rstput( kt, nitrst, numrtw, 'sshn', ssha )
164            CALL iom_rstput( kt, nitrst, numrtw, 'sshb', sshn )
165         ENDIF
166      ENDIF
167      !
168      CALL eos    ( tsn, rhd, rhop, gdept_0(:,:,:) ) ! In any case, we need rhop
169      CALL eos_rab( tsn, rab_n )       ! now    local thermal/haline expension ratio at T-points
170      CALL bn2    ( tsn, rab_n, rn2 ) ! before Brunt-Vaisala frequency need for zdfmxl
171
172      rn2b(:,:,:) = rn2(:,:,:)         ! need for zdfmxl
173      CALL zdf_mxl( kt )                                                   ! In any case, we need mxl
174      !
175      hmld(:,:)         = sf_dyn(jf_mld)%fnow(:,:,1) * tmask(:,:,1)    ! mixed layer depht
176      avt(:,:,:)        = sf_dyn(jf_avt)%fnow(:,:,:) * tmask(:,:,:)    ! vertical diffusive coefficient
177      !
178#if defined key_trabbl && ! defined key_c1d
179      ahu_bbl(:,:)      = sf_dyn(jf_ubl)%fnow(:,:,1) * umask(:,:,1)    ! bbl diffusive coef
180      ahv_bbl(:,:)      = sf_dyn(jf_vbl)%fnow(:,:,1) * vmask(:,:,1)
181#endif
182      !
183      !
184      CALL eos( tsn, rhd, rhop, gdept_0(:,:,:) ) ! In any case, we need rhop
185      !
186      IF(ln_ctl) THEN                  ! print control
187         CALL prt_ctl(tab3d_1=tsn(:,:,:,jp_tem), clinfo1=' tn      - : ', mask1=tmask, ovlap=1, kdim=jpk   )
188         CALL prt_ctl(tab3d_1=tsn(:,:,:,jp_sal), clinfo1=' sn      - : ', mask1=tmask, ovlap=1, kdim=jpk   )
189         CALL prt_ctl(tab3d_1=un               , clinfo1=' un      - : ', mask1=umask, ovlap=1, kdim=jpk   )
190         CALL prt_ctl(tab3d_1=vn               , clinfo1=' vn      - : ', mask1=vmask, ovlap=1, kdim=jpk   )
191         CALL prt_ctl(tab3d_1=wn               , clinfo1=' wn      - : ', mask1=tmask, ovlap=1, kdim=jpk   )
192         CALL prt_ctl(tab3d_1=avt              , clinfo1=' kz      - : ', mask1=tmask, ovlap=1, kdim=jpk   )
193!         CALL prt_ctl(tab2d_1=fr_i             , clinfo1=' fr_i    - : ', mask1=tmask, ovlap=1 )
194!         CALL prt_ctl(tab2d_1=hmld             , clinfo1=' hmld    - : ', mask1=tmask, ovlap=1 )
195!         CALL prt_ctl(tab2d_1=fmmflx           , clinfo1=' fmmflx  - : ', mask1=tmask, ovlap=1 )
196!         CALL prt_ctl(tab2d_1=emp              , clinfo1=' emp     - : ', mask1=tmask, ovlap=1 )
197!         CALL prt_ctl(tab2d_1=wndm             , clinfo1=' wspd    - : ', mask1=tmask, ovlap=1 )
198!         CALL prt_ctl(tab2d_1=qsr              , clinfo1=' qsr     - : ', mask1=tmask, ovlap=1 )
199      ENDIF
200      !
201      IF( nn_timing == 1 )  CALL timing_stop( 'dta_dyn')
202      !
203   END SUBROUTINE dta_dyn
204
205
206   SUBROUTINE dta_dyn_init
207      !!----------------------------------------------------------------------
208      !!                  ***  ROUTINE dta_dyn_init  ***
209      !!
210      !! ** Purpose :   Initialisation of the dynamical data     
211      !! ** Method  : - read the data namdta_dyn namelist
212      !!
213      !! ** Action  : - read parameters
214      !!----------------------------------------------------------------------
215      INTEGER  :: ierr, ierr0, ierr1, ierr2, ierr3   ! return error code
216      INTEGER  :: ifpr                               ! dummy loop indice
217      INTEGER  :: jfld                               ! dummy loop arguments
218      INTEGER  :: inum, idv, idimv                   ! local integer
219      INTEGER  :: ios                                ! Local integer output status for namelist read
220      INTEGER  :: ji, jj, jk
221      REAL(wp) :: zcoef
222      INTEGER  :: nkrnf_max
223      REAL(wp) :: hrnf_max
224      !!
225      CHARACTER(len=100)            ::  cn_dir        !   Root directory for location of core files
226      TYPE(FLD_N), DIMENSION(jpfld) ::  slf_d         ! array of namelist informations on the fields to read
227      TYPE(FLD_N) :: sn_uwd, sn_vwd, sn_wwd, sn_empb, sn_emp  ! informations about the fields to be read
228      TYPE(FLD_N) :: sn_tem , sn_sal , sn_avt   !   "                 "
229      TYPE(FLD_N) :: sn_mld, sn_qsr, sn_wnd , sn_ice , sn_fmf   !   "               "
230      TYPE(FLD_N) :: sn_ubl, sn_vbl, sn_rnf    !   "              "
231      TYPE(FLD_N) :: sn_div  ! informations about the fields to be read
232      !!----------------------------------------------------------------------
233
234      NAMELIST/namdta_dyn/cn_dir, ln_dynrnf, ln_dynrnf_depth,  ln_ssh_ini, fwbcorr, &
235         &                sn_uwd, sn_vwd, sn_wwd, sn_emp,    &
236         &                sn_avt, sn_tem, sn_sal, sn_mld , sn_qsr ,   &
237         &                sn_wnd, sn_ice, sn_fmf,                    &
238         &                sn_ubl, sn_vbl, sn_rnf,                   &
239         &                sn_empb, sn_div 
240      !
241      REWIND( numnam_ref )              ! Namelist namdta_dyn in reference namelist : Offline: init. of dynamical data
242      READ  ( numnam_ref, namdta_dyn, IOSTAT = ios, ERR = 901)
243901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdta_dyn in reference namelist', lwp )
244
245      REWIND( numnam_cfg )              ! Namelist namdta_dyn in configuration namelist : Offline: init. of dynamical data
246      READ  ( numnam_cfg, namdta_dyn, IOSTAT = ios, ERR = 902 )
247902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdta_dyn in configuration namelist', lwp )
248      IF(lwm) WRITE ( numond, namdta_dyn )
249      !                                         ! store namelist information in an array
250      !                                         ! Control print
251      IF(lwp) THEN
252         WRITE(numout,*)
253         WRITE(numout,*) 'dta_dyn : offline dynamics '
254         WRITE(numout,*) '~~~~~~~ '
255         WRITE(numout,*) '   Namelist namdta_dyn'
256         WRITE(numout,*) '      ssh initialised from dyn file (T) or not (F)     ln_ssh_ini       = ', ln_ssh_ini
257         WRITE(numout,*) '      runoffs option enabled (T) or not (F)            ln_dynrnf        = ', ln_dynrnf
258         WRITE(numout,*) '      runoffs is spread in vertical                    ln_dynrnf_depth  = ', ln_dynrnf_depth
259         WRITE(numout,*) '      annual global mean of empmr for ssh correction   fwbcorr          = ', fwbcorr
260         WRITE(numout,*)
261      ENDIF
262      !
263
264      jf_uwd  = 1     ;   jf_vwd  = 2    ;   jf_wwd = 3    ;   jf_emp = 4    ;   jf_avt = 5
265      jf_tem  = 6     ;   jf_sal  = 7    ;   jf_mld = 8    ;   jf_qsr = 9
266      jf_wnd  = 10    ;   jf_ice  = 11   ;   jf_fmf = 12   ;   jfld   = jf_fmf
267
268      !
269      slf_d(jf_uwd)  = sn_uwd    ;   slf_d(jf_vwd)  = sn_vwd   ;   slf_d(jf_wwd) = sn_wwd
270      slf_d(jf_emp)  = sn_emp    ;   slf_d(jf_avt)  = sn_avt
271      slf_d(jf_tem)  = sn_tem    ;   slf_d(jf_sal)  = sn_sal   ;   slf_d(jf_mld) = sn_mld
272      slf_d(jf_qsr)  = sn_qsr    ;   slf_d(jf_wnd)  = sn_wnd   ;   slf_d(jf_ice) = sn_ice
273      slf_d(jf_fmf)  = sn_fmf
274
275
276      !
277      IF( lk_vvl ) THEN
278                 jf_div  = jfld + 1    ;         jf_empb  = jfld + 2      ;      jfld = jf_empb
279           slf_d(jf_div) = sn_div      ;   slf_d(jf_empb) = sn_empb
280      ENDIF
281      !
282      IF( lk_trabbl ) THEN
283                 jf_ubl  = jfld + 1    ;         jf_vbl  = jfld + 2       ;      jfld = jf_vbl
284           slf_d(jf_ubl) = sn_ubl      ;   slf_d(jf_vbl) = sn_vbl
285      ENDIF
286      !
287      IF( ln_dynrnf ) THEN
288                jf_rnf  = jfld + 1     ;     jfld  = jf_rnf
289          slf_d(jf_rnf) = sn_rnf   
290      ELSE
291         rnf(:,:) = 0._wp
292      ENDIF
293
294 
295      ALLOCATE( sf_dyn(jfld), STAT=ierr )         ! set sf structure
296      IF( ierr > 0 )  THEN
297         CALL ctl_stop( 'dta_dyn: unable to allocate sf structure' )   ;   RETURN
298      ENDIF
299      !                                         ! fill sf with slf_i and control print
300      CALL fld_fill( sf_dyn, slf_d, cn_dir, 'dta_dyn_init', 'Data in file', 'namdta_dyn' )
301      !
302      ! Open file for each variable to get his number of dimension
303      DO ifpr = 1, jfld
304         CALL fld_clopn( sf_dyn(ifpr), nyear, nmonth, nday )
305         idv   = iom_varid( sf_dyn(ifpr)%num , slf_d(ifpr)%clvar )        ! id of the variable sdjf%clvar
306         idimv = iom_file ( sf_dyn(ifpr)%num )%ndims(idv)                 ! number of dimension for variable sdjf%clvar
307         IF( sf_dyn(ifpr)%num /= 0 )   CALL iom_close( sf_dyn(ifpr)%num ) ! close file if already open
308         ierr1=0
309         IF( idimv == 3 ) THEN    ! 2D variable
310                                      ALLOCATE( sf_dyn(ifpr)%fnow(jpi,jpj,1)    , STAT=ierr0 )
311            IF( slf_d(ifpr)%ln_tint ) ALLOCATE( sf_dyn(ifpr)%fdta(jpi,jpj,1,2)  , STAT=ierr1 )
312         ELSE                     ! 3D variable
313                                      ALLOCATE( sf_dyn(ifpr)%fnow(jpi,jpj,jpk)  , STAT=ierr0 )
314            IF( slf_d(ifpr)%ln_tint ) ALLOCATE( sf_dyn(ifpr)%fdta(jpi,jpj,jpk,2), STAT=ierr1 )
315         ENDIF
316         IF( ierr0 + ierr1 > 0 ) THEN
317            CALL ctl_stop( 'dta_dyn_init : unable to allocate sf_dyn array structure' )   ;   RETURN
318         ENDIF
319      END DO
320      !
321      IF( lk_ldfslp .AND. .NOT.lk_c1d ) THEN                  ! slopes
322         IF( sf_dyn(jf_tem)%ln_tint ) THEN      ! time interpolation
323            ALLOCATE( uslpdta (jpi,jpj,jpk,2), vslpdta (jpi,jpj,jpk,2),    &
324            &         wslpidta(jpi,jpj,jpk,2), wslpjdta(jpi,jpj,jpk,2), STAT=ierr2 )
325            !
326            IF( ierr2 > 0 )  THEN
327               CALL ctl_stop( 'dta_dyn_init : unable to allocate slope arrays' )   ;   RETURN
328            ENDIF
329         ENDIF
330      ENDIF
331      !
332      IF( lk_vvl ) THEN
333        IF( ln_ssh_ini ) THEN                     ! Restart: read in restart file
334           IF(lwp) WRITE(numout,*) ' sshn forcing fields read in the dynamics restart file for initialisation'
335           CALL iom_open( 'restart', inum )
336           CALL iom_get( inum, jpdom_autoglo, 'sshn', sshn(:,:)   )
337           CALL iom_get( inum, jpdom_autoglo, 'sshb', sshb(:,:)   )
338           CALL iom_close( inum )                                        ! close file
339        ELSE
340           IF(lwp) WRITE(numout,*) ' sshn forcing fields read in passive tracers restart file for initialisation'
341           CALL iom_get( numrtr, jpdom_autoglo, 'sshn', sshn(:,:)   )
342           CALL iom_get( numrtr, jpdom_autoglo, 'sshb', sshb(:,:)   )
343        ENDIF
344        !
345        DO jk = 1, jpkm1
346           fse3t_n(:,:,jk) = e3t_0(:,:,jk) * ( 1._wp + sshn(:,:) * tmask(:,:,1) / ( ht_0(:,:) + 1.0 - tmask(:,:,1) ) )
347        ENDDO
348        fse3t_a(:,:,jpk) = e3t_0(:,:,jpk)
349
350        ! Horizontal scale factor interpolations
351        ! --------------------------------------
352        CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3u_n(:,:,:), 'U' )
353        CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3v_n(:,:,:), 'V' )
354
355        ! Vertical scale factor interpolations
356        ! ------------------------------------
357        CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3w_n(:,:,:), 'W' )
358 
359        fse3t_b(:,:,:)  = fse3t_n(:,:,:)
360        fse3u_b(:,:,:)  = fse3u_n(:,:,:)
361        fse3v_b(:,:,:)  = fse3v_n(:,:,:)
362
363        ! t- and w- points depth
364        ! ----------------------
365        fsdept_n(:,:,1) = 0.5_wp * fse3w_n(:,:,1)
366        fsdepw_n(:,:,1) = 0.0_wp
367
368        DO jk = 2, jpk
369           DO jj = 1,jpj
370              DO ji = 1,jpi
371                !    zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))   ! 0 everywhere
372                !    tmask = wmask, ie everywhere expect at jk = mikt
373                                                                   ! 1 for jk =
374                                                                   ! mikt
375                 zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
376                 fsdepw_n(ji,jj,jk) = fsdepw_n(ji,jj,jk-1) + fse3t_n(ji,jj,jk-1)
377                 fsdept_n(ji,jj,jk) =      zcoef  * ( fsdepw_n(ji,jj,jk  ) + 0.5 * fse3w_n(ji,jj,jk))  &
378                     &                + (1-zcoef) * ( fsdept_n(ji,jj,jk-1) + fse3w_n(ji,jj,jk))
379              END DO
380           END DO
381        END DO
382
383        fsdept_b(:,:,:) = fsdept_n(:,:,:)
384        fsdepw_b(:,:,:) = fsdepw_n(:,:,:)
385        !
386      ENDIF
387      !
388      IF( ln_dynrnf .AND. ln_dynrnf_depth ) THEN       ! read depht over which runoffs are distributed
389         IF(lwp) WRITE(numout,*) 
390         IF(lwp) WRITE(numout,*) ' read in the file depht over which runoffs are distributed'
391         CALL iom_open ( "runoffs", inum )                           ! open file
392         CALL iom_get  ( inum, jpdom_data, 'rodepth', h_rnf )   ! read the river mouth array
393         CALL iom_close( inum )                                        ! close file
394         !
395         nk_rnf(:,:) = 0                               ! set the number of level over which river runoffs are applied
396         DO jj = 1, jpj
397            DO ji = 1, jpi
398               IF( h_rnf(ji,jj) > 0._wp ) THEN
399                  jk = 2
400                  DO WHILE ( jk /= mbkt(ji,jj) .AND. gdept_0(ji,jj,jk) < h_rnf(ji,jj) ) ;  jk = jk + 1
401                  END DO
402                  nk_rnf(ji,jj) = jk
403               ELSEIF( h_rnf(ji,jj) == -1._wp   ) THEN   ;  nk_rnf(ji,jj) = 1
404               ELSEIF( h_rnf(ji,jj) == -999._wp ) THEN   ;  nk_rnf(ji,jj) = mbkt(ji,jj)
405               ELSE
406                  CALL ctl_stop( 'sbc_rnf_init: runoff depth not positive, and not -999 or -1, rnf value in file fort.999'  )
407                  WRITE(999,*) 'ji, jj, h_rnf(ji,jj) :', ji, jj, h_rnf(ji,jj)
408               ENDIF
409            END DO
410         END DO
411         DO jj = 1, jpj                                ! set the associated depth
412            DO ji = 1, jpi
413               h_rnf(ji,jj) = 0._wp
414               DO jk = 1, nk_rnf(ji,jj)
415                  h_rnf(ji,jj) = h_rnf(ji,jj) + fse3t(ji,jj,jk)
416               END DO
417            END DO
418         END DO
419      ELSE                                       ! runoffs applied at the surface
420         nk_rnf(:,:) = 1
421         h_rnf (:,:) = fse3t(:,:,1)
422      ENDIF
423      nkrnf_max = MAXVAL( nk_rnf(:,:) )
424      hrnf_max = MAXVAL( h_rnf(:,:) )
425      IF( lk_mpp )  THEN
426         CALL mpp_max( nkrnf_max )                 ! max over the  global domain
427         CALL mpp_max( hrnf_max )                 ! max over the  global domain
428      ENDIF
429      IF(lwp) WRITE(numout,*) ' '
430      IF(lwp) WRITE(numout,*) ' max depht of runoff : ', hrnf_max,'    max level  : ', nkrnf_max
431      IF(lwp) WRITE(numout,*) ' '
432      !
433      CALL dta_dyn( nit000 )
434      !
435   END SUBROUTINE dta_dyn_init
436
437   SUBROUTINE dta_dyn_swp( kt )
438     !!---------------------------------------------------------------------
439      !!                    ***  ROUTINE dta_dyn_swp  ***
440      !!
441      !! ** Purpose : Swap and the data and compute the vertical scale factor at U/V/W point
442      !!              and the depht
443      !!
444      !!---------------------------------------------------------------------
445      INTEGER, INTENT(in) :: kt       ! time step
446      INTEGER             :: ji, jj, jk
447      REAL(wp)            :: zcoef
448      !
449      !!---------------------------------------------------------------------
450
451      IF( kt == nit000 ) THEN
452         IF(lwp) WRITE(numout,*)
453         IF(lwp) WRITE(numout,*) 'ssh_swp : Asselin time filter and swap of sea surface height'
454         IF(lwp) WRITE(numout,*) '~~~~~~~ '
455      ENDIF
456
457      sshb(:,:) = sshn(:,:) + atfp * ( sshb(:,:) - 2 * sshn(:,:) + ssha(:,:))  ! before <-- now filtered
458      sshn(:,:) = ssha(:,:)
459
460      fse3t_n(:,:,:) = fse3t_a(:,:,:)
461
462      ! Reconstruction of all vertical scale factors at now and before time steps
463      ! =============================================================================
464
465      ! Horizontal scale factor interpolations
466      ! --------------------------------------
467      CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3u_n(:,:,:), 'U' )
468      CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3v_n(:,:,:), 'V' )
469
470      ! Vertical scale factor interpolations
471      ! ------------------------------------
472      CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3w_n (:,:,:), 'W' )
473
474      fse3t_b(:,:,:)  = fse3t_n(:,:,:)
475      fse3u_b(:,:,:)  = fse3u_n(:,:,:)
476      fse3v_b(:,:,:)  = fse3v_n(:,:,:)
477
478      ! t- and w- points depth
479      ! ----------------------
480      fsdept_n(:,:,1) = 0.5_wp * fse3w_n(:,:,1)
481      fsdepw_n(:,:,1) = 0.0_wp
482
483      DO jk = 2, jpk
484         DO jj = 1,jpj
485            DO ji = 1,jpi
486                 zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
487                 fsdepw_n(ji,jj,jk) = fsdepw_n(ji,jj,jk-1) + fse3t_n(ji,jj,jk-1)
488                 fsdept_n(ji,jj,jk) =      zcoef  * ( fsdepw_n(ji,jj,jk  ) + 0.5 * fse3w_n(ji,jj,jk))  &
489                     &                + (1-zcoef) * ( fsdept_n(ji,jj,jk-1) + fse3w_n(ji,jj,jk))
490              END DO
491           END DO
492        END DO
493
494      fsdept_b(:,:,:) = fsdept_n(:,:,:)
495      fsdepw_b(:,:,:) = fsdepw_n(:,:,:)
496
497      !
498   END SUBROUTINE dta_dyn_swp
499
500   SUBROUTINE dta_dyn_ssh( kt, phdivtr, psshb,  pemp, pssha, pe3ta )
501      !!----------------------------------------------------------------------
502      !!                ***  ROUTINE dta_dyn_wzv  ***
503      !!                   
504      !! ** Purpose :   compute the after ssh (ssha) and the now vertical velocity
505      !!
506      !! ** Method  : Using the incompressibility hypothesis,
507      !!        - the ssh increment is computed by integrating the horizontal divergence
508      !!          and multiply by the time step.
509      !!
510      !!        - compute the after scale factor : repartition of ssh INCREMENT proportionnaly
511      !!                                           to the level thickness ( z-star case )
512      !!
513      !!        - the vertical velocity is computed by integrating the horizontal divergence 
514      !!          from the bottom to the surface minus the scale factor evolution.
515      !!          The boundary conditions are w=0 at the bottom (no flux)
516      !!
517      !! ** action  :   ssha / e3t_a / wn
518      !!
519      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
520      !!----------------------------------------------------------------------
521      !! * Arguments
522      INTEGER,                                   INTENT(in )    :: kt        !  time-step
523      REAL(wp), DIMENSION(jpi,jpj,jpk)          , INTENT(in )   :: phdivtr   ! horizontal divergence transport
524      REAL(wp), DIMENSION(jpi,jpj)    , OPTIONAL, INTENT(in )   :: psshb     ! now ssh
525      REAL(wp), DIMENSION(jpi,jpj)    , OPTIONAL, INTENT(in )   :: pemp      ! evaporation minus precipitation
526      REAL(wp), DIMENSION(jpi,jpj)    , OPTIONAL, INTENT(inout) :: pssha     ! after ssh
527      REAL(wp), DIMENSION(jpi,jpj,jpk), OPTIONAL, INTENT(out)   :: pe3ta     ! after vertical scale factor
528      !! * Local declarations
529      INTEGER                       :: jk
530      REAL(wp), DIMENSION(jpi,jpj)  :: zhdiv 
531      REAL(wp)  :: z2dt 
532      !!----------------------------------------------------------------------
533     
534      !
535      z2dt = 2._wp * rdt
536      !
537      zhdiv(:,:) = 0._wp
538      DO jk = 1, jpkm1
539         zhdiv(:,:) = zhdiv(:,:) +  phdivtr(:,:,jk) * tmask(:,:,jk)
540      END DO
541      !                                                ! Sea surface  elevation time-stepping
542      pssha(:,:) = ( psshb(:,:) - z2dt * ( r1_rau0 * pemp(:,:)  + zhdiv(:,:) ) ) * ssmask(:,:)
543      !                                                 !
544      !                                                 ! After acale factors at t-points ( z_star coordinate )
545      DO jk = 1, jpkm1
546        pe3ta(:,:,jk) = e3t_0(:,:,jk) * ( 1._wp + pssha(:,:) * tmask(:,:,1) / ( ht_0(:,:) + 1.0 - tmask(:,:,1) ) )
547      END DO
548      !
549   END SUBROUTINE dta_dyn_ssh
550
551
552   SUBROUTINE dta_dyn_hrnf
553      !!----------------------------------------------------------------------
554      !!                  ***  ROUTINE sbc_rnf  ***
555      !!
556      !! ** Purpose :   update the horizontal divergence with the runoff inflow
557      !!
558      !! ** Method  :
559      !!                CAUTION : rnf is positive (inflow) decreasing the
560      !!                          divergence and expressed in m/s
561      !!
562      !! ** Action  :   phdivn   decreased by the runoff inflow
563      !!----------------------------------------------------------------------
564      !!
565      INTEGER  ::   ji, jj, jk   ! dummy loop indices
566      !!----------------------------------------------------------------------
567      !
568      DO jj = 1, jpj                   ! update the depth over which runoffs are distributed
569         DO ji = 1, jpi
570            h_rnf(ji,jj) = 0._wp
571            DO jk = 1, nk_rnf(ji,jj)                           ! recalculates h_rnf to be the depth in metres
572                h_rnf(ji,jj) = h_rnf(ji,jj) + fse3t(ji,jj,jk)   ! to the bottom of the relevant grid box
573            END DO
574        END DO
575      END DO
576      !
577   END SUBROUTINE dta_dyn_hrnf
578
579
580
581   SUBROUTINE dta_dyn_slp( kt )
582      !!---------------------------------------------------------------------
583      !!                    ***  ROUTINE dta_dyn_slp  ***
584      !!
585      !! ** Purpose : Computation of slope
586      !!
587      !!---------------------------------------------------------------------
588      USE oce, ONLY:  zts => tsa 
589      !
590      INTEGER,  INTENT(in) :: kt       ! time step
591      !
592      INTEGER  ::   ji, jj     ! dummy loop indices
593      REAL(wp) ::   ztinta     ! ratio applied to after  records when doing time interpolation
594      REAL(wp) ::   ztintb     ! ratio applied to before records when doing time interpolation
595      INTEGER  ::   iswap
596      REAL(wp), POINTER, DIMENSION(:,:,:) :: zuslp, zvslp, zwslpi, zwslpj
597      !!---------------------------------------------------------------------
598      !
599      CALL wrk_alloc(jpi, jpj, jpk, zuslp, zvslp, zwslpi, zwslpj )
600      !
601      IF( sf_dyn(jf_tem)%ln_tint ) THEN    ! Computes slopes (here avt is used as workspace)                       
602         IF( kt == nit000 ) THEN
603            zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fdta(:,:,:,1) * tmask(:,:,:)   ! temperature
604            zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fdta(:,:,:,1) * tmask(:,:,:)   ! salinity
605            avt(:,:,:)        = sf_dyn(jf_avt)%fdta(:,:,:,1) * tmask(:,:,:)   ! vertical diffusive coef.
606            CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj )
607            uslpdta (:,:,:,1) = zuslp (:,:,:) 
608            vslpdta (:,:,:,1) = zvslp (:,:,:) 
609            wslpidta(:,:,:,1) = zwslpi(:,:,:) 
610            wslpjdta(:,:,:,1) = zwslpj(:,:,:) 
611            !
612            zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fdta(:,:,:,2) * tmask(:,:,:)   ! temperature
613            zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fdta(:,:,:,2) * tmask(:,:,:)   ! salinity
614            avt(:,:,:)        = sf_dyn(jf_avt)%fdta(:,:,:,2) * tmask(:,:,:)   ! vertical diffusive coef.
615            CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj )
616            uslpdta (:,:,:,2) = zuslp (:,:,:) 
617            vslpdta (:,:,:,2) = zvslp (:,:,:) 
618            wslpidta(:,:,:,2) = zwslpi(:,:,:) 
619            wslpjdta(:,:,:,2) = zwslpj(:,:,:) 
620         ELSE
621           !
622           iswap = 0
623           IF( sf_dyn(jf_tem)%nrec_a(2) - nprevrec /= 0 )  iswap = 1
624           IF( nsecdyn > sf_dyn(jf_tem)%nrec_b(2) .AND. iswap == 1 )  THEN    ! read/update the after data
625              IF(lwp) WRITE(numout,*) ' Compute new slopes at kt = ', kt
626              uslpdta (:,:,:,1) =  uslpdta (:,:,:,2)         ! swap the data
627              vslpdta (:,:,:,1) =  vslpdta (:,:,:,2) 
628              wslpidta(:,:,:,1) =  wslpidta(:,:,:,2) 
629              wslpjdta(:,:,:,1) =  wslpjdta(:,:,:,2) 
630              !
631              zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fdta(:,:,:,2) * tmask(:,:,:)   ! temperature
632              zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fdta(:,:,:,2) * tmask(:,:,:)   ! salinity
633              avt(:,:,:)        = sf_dyn(jf_avt)%fdta(:,:,:,2) * tmask(:,:,:)   ! vertical diffusive coef.
634              CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj )
635              !
636              uslpdta (:,:,:,2) = zuslp (:,:,:) 
637              vslpdta (:,:,:,2) = zvslp (:,:,:) 
638              wslpidta(:,:,:,2) = zwslpi(:,:,:) 
639              wslpjdta(:,:,:,2) = zwslpj(:,:,:) 
640            ENDIF
641         ENDIF
642      ENDIF
643      !
644      IF( sf_dyn(jf_tem)%ln_tint )  THEN
645         ztinta =  REAL( nsecdyn - sf_dyn(jf_tem)%nrec_b(2), wp )  &
646            &    / REAL( sf_dyn(jf_tem)%nrec_a(2) - sf_dyn(jf_tem)%nrec_b(2), wp )
647         ztintb =  1. - ztinta
648#if defined key_ldfslp && ! defined key_c1d
649         uslp (:,:,:) = ztintb * uslpdta (:,:,:,1)  + ztinta * uslpdta (:,:,:,2) 
650         vslp (:,:,:) = ztintb * vslpdta (:,:,:,1)  + ztinta * vslpdta (:,:,:,2) 
651         wslpi(:,:,:) = ztintb * wslpidta(:,:,:,1)  + ztinta * wslpidta(:,:,:,2) 
652         wslpj(:,:,:) = ztintb * wslpjdta(:,:,:,1)  + ztinta * wslpjdta(:,:,:,2) 
653#endif
654      ELSE
655         zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fnow(:,:,:) * tmask(:,:,:)   ! temperature
656         zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fnow(:,:,:) * tmask(:,:,:)   ! salinity
657         avt(:,:,:)        = sf_dyn(jf_avt)%fnow(:,:,:) * tmask(:,:,:)   ! vertical diffusive coef.
658         CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj )
659         !
660#if defined key_ldfslp && ! defined key_c1d
661         uslp (:,:,:) = zuslp (:,:,:)
662         vslp (:,:,:) = zvslp (:,:,:)
663         wslpi(:,:,:) = zwslpi(:,:,:)
664         wslpj(:,:,:) = zwslpj(:,:,:)
665#endif
666      ENDIF
667      !
668      CALL wrk_dealloc(jpi, jpj, jpk, zuslp, zvslp, zwslpi, zwslpj )
669      !
670   END SUBROUTINE dta_dyn_slp
671
672   SUBROUTINE compute_slopes( kt, pts, puslp, pvslp, pwslpi, pwslpj )
673      !!---------------------------------------------------------------------
674      !!                    ***  ROUTINE dta_dyn_slp  ***
675      !!
676      !! ** Purpose : Computation of slope
677      !!
678      !!---------------------------------------------------------------------
679      INTEGER ,                              INTENT(in ) :: kt       ! time step
680      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: pts      ! temperature/salinity
681      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(out) :: puslp    ! zonal isopycnal slopes
682      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(out) :: pvslp    ! meridional isopycnal slopes
683      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(out) :: pwslpi   ! zonal diapycnal slopes
684      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(out) :: pwslpj   ! meridional diapycnal slopes
685      !!---------------------------------------------------------------------
686#if defined key_ldfslp && ! defined key_c1d
687      CALL eos    ( pts, rhd, rhop, gdept_0(:,:,:) )
688      CALL eos_rab( pts, rab_n )       ! now local thermal/haline expension ratio at T-points
689      CALL bn2    ( pts, rab_n, rn2  ) ! now    Brunt-Vaisala
690
691      ! Partial steps: before Horizontal DErivative
692      IF( ln_zps  .AND. .NOT. ln_isfcav)                            &
693         &            CALL zps_hde    ( kt, jpts, pts, gtsu, gtsv,  &  ! Partial steps: before horizontal gradient
694         &                                        rhd, gru , grv    )  ! of t, s, rd at the last ocean level
695      IF( ln_zps .AND.        ln_isfcav)                            &
696         &            CALL zps_hde_isf( kt, jpts, pts, gtsu, gtsv,  &    ! Partial steps for top cell (ISF)
697         &                                        rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   &
698         &                                 gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi    ) ! of t, s, rd at the first ocean level
699
700      rn2b(:,:,:) = rn2(:,:,:)         ! need for zdfmxl
701      CALL zdf_mxl( kt )            ! mixed layer depth
702      CALL ldf_slp( kt, rhd, rn2 )  ! slopes
703      puslp (:,:,:) = uslp (:,:,:) 
704      pvslp (:,:,:) = vslp (:,:,:) 
705      pwslpi(:,:,:) = wslpi(:,:,:) 
706      pwslpj(:,:,:) = wslpj(:,:,:) 
707#else
708      puslp (:,:,:) = 0.            ! to avoid warning when compiling
709      pvslp (:,:,:) = 0.
710      pwslpi(:,:,:) = 0.
711      pwslpj(:,:,:) = 0.
712#endif
713      !
714   END SUBROUTINE compute_slopes
715
716   !!======================================================================
717END MODULE dtadyn
Note: See TracBrowser for help on using the repository browser.