- Timestamp:
- 2020-05-14T21:46:00+02:00 (4 years ago)
- Location:
- NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
src/OFF/dtadyn.F90 (modified) (39 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser
- Property svn:externals
-
old new 6 6 ^/vendors/FCM@HEAD ext/FCM 7 7 ^/vendors/IOIPSL@HEAD ext/IOIPSL 8 9 # SETTE 10 ^/utils/CI/sette@HEAD sette
-
- Property svn:externals
-
NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OFF/dtadyn.F90
r12178 r12928 13 13 !! 3.3 ! 2010-11 (C. Ethe) Full reorganization of the off-line: phasing with the on-line 14 14 !! 3.4 ! 2011-05 (C. Ethe) Use of fldread 15 !! 4.1 ! 2019-08 (A. Coward, D. Storkey) split dta_dyn_sf_swp -> dta_dyn_sf_atf and dta_dyn_sf_interp 15 16 !!---------------------------------------------------------------------- 16 17 … … 46 47 PRIVATE 47 48 48 PUBLIC dta_dyn_init ! called by opa.F90 49 PUBLIC dta_dyn ! called by step.F90 50 PUBLIC dta_dyn_sed_init ! called by opa.F90 51 PUBLIC dta_dyn_sed ! called by step.F90 52 PUBLIC dta_dyn_swp ! called by step.F90 49 PUBLIC dta_dyn_init ! called by nemo_init 50 PUBLIC dta_dyn ! called by nemo_gcm 51 PUBLIC dta_dyn_sed_init ! called by nemo_init 52 PUBLIC dta_dyn_sed ! called by nemo_gcm 53 PUBLIC dta_dyn_atf ! called by nemo_gcm 54 PUBLIC dta_dyn_sf_interp ! called by nemo_gcm 53 55 54 56 CHARACTER(len=100) :: cn_dir !: Root directory for location of ssr files … … 87 89 INTEGER, SAVE :: nprevrec, nsecdyn 88 90 91 !! * Substitutions 92 # include "do_loop_substitute.h90" 89 93 !!---------------------------------------------------------------------- 90 94 !! NEMO/OFF 4.0 , NEMO Consortium (2018) … … 94 98 CONTAINS 95 99 96 SUBROUTINE dta_dyn( kt )100 SUBROUTINE dta_dyn( kt, Kbb, Kmm, Kaa ) 97 101 !!---------------------------------------------------------------------- 98 102 !! *** ROUTINE dta_dyn *** … … 105 109 !! - interpolates data if needed 106 110 !!---------------------------------------------------------------------- 107 INTEGER, INTENT(in) :: kt ! ocean time-step index 111 INTEGER, INTENT(in) :: kt ! ocean time-step index 112 INTEGER, INTENT(in) :: Kbb, Kmm, Kaa ! ocean time level indices 108 113 ! 109 114 INTEGER :: ji, jj, jk … … 121 126 CALL fld_read( kt, 1, sf_dyn ) != read data at kt time step ==! 122 127 ! 123 IF( l_ldfslp .AND. .NOT.lk_c1d ) CALL dta_dyn_slp( kt ) ! Computation of slopes124 ! 125 ts n(:,:,:,jp_tem) = sf_dyn(jf_tem)%fnow(:,:,:) * tmask(:,:,:) ! temperature126 ts n(:,:,:,jp_sal) = sf_dyn(jf_sal)%fnow(:,:,:) * tmask(:,:,:) ! salinity128 IF( l_ldfslp .AND. .NOT.lk_c1d ) CALL dta_dyn_slp( kt, Kbb, Kmm ) ! Computation of slopes 129 ! 130 ts(:,:,:,jp_tem,Kmm) = sf_dyn(jf_tem)%fnow(:,:,:) * tmask(:,:,:) ! temperature 131 ts(:,:,:,jp_sal,Kmm) = sf_dyn(jf_sal)%fnow(:,:,:) * tmask(:,:,:) ! salinity 127 132 wndm(:,:) = sf_dyn(jf_wnd)%fnow(:,:,1) * tmask(:,:,1) ! wind speed - needed for gas exchange 128 133 fmmflx(:,:) = sf_dyn(jf_fmf)%fnow(:,:,1) * tmask(:,:,1) ! downward salt flux (v3.5+) … … 132 137 IF( ln_dynrnf ) THEN 133 138 rnf (:,:) = sf_dyn(jf_rnf)%fnow(:,:,1) * tmask(:,:,1) ! E-P 134 IF( ln_dynrnf_depth .AND. .NOT. ln_linssh ) CALL dta_dyn_hrnf 135 ENDIF 136 ! 137 u n(:,:,:) = sf_dyn(jf_uwd)%fnow(:,:,:) * umask(:,:,:) ! effective u-transport138 v n(:,:,:) = sf_dyn(jf_vwd)%fnow(:,:,:) * vmask(:,:,:) ! effective v-transport139 w n(:,:,:) = sf_dyn(jf_wwd)%fnow(:,:,:) * tmask(:,:,:) ! effective v-transport139 IF( ln_dynrnf_depth .AND. .NOT. ln_linssh ) CALL dta_dyn_hrnf(Kmm) 140 ENDIF 141 ! 142 uu(:,:,:,Kmm) = sf_dyn(jf_uwd)%fnow(:,:,:) * umask(:,:,:) ! effective u-transport 143 vv(:,:,:,Kmm) = sf_dyn(jf_vwd)%fnow(:,:,:) * vmask(:,:,:) ! effective v-transport 144 ww(:,:,:) = sf_dyn(jf_wwd)%fnow(:,:,:) * tmask(:,:,:) ! effective v-transport 140 145 ! 141 146 IF( .NOT.ln_linssh ) THEN … … 144 149 emp_b (:,:) = sf_dyn(jf_empb)%fnow(:,:,1) * tmask(:,:,1) ! E-P 145 150 zemp (:,:) = ( 0.5_wp * ( emp(:,:) + emp_b(:,:) ) + rnf(:,:) + fwbcorr ) * tmask(:,:,1) 146 CALL dta_dyn_ssh( kt, zhdivtr, ssh b, zemp, ssha, e3t_a(:,:,:) ) != ssh, vertical scale factor & vertical transport151 CALL dta_dyn_ssh( kt, zhdivtr, ssh(:,:,Kbb), zemp, ssh(:,:,Kaa), e3t(:,:,:,Kaa) ) != ssh, vertical scale factor & vertical transport 147 152 DEALLOCATE( zemp , zhdivtr ) 148 153 ! Write in the tracer restart file … … 152 157 IF(lwp) WRITE(numout,*) 'dta_dyn_ssh : ssh field written in tracer restart file at it= ', kt,' date= ', ndastp 153 158 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 154 CALL iom_rstput( kt, nitrst, numrtw, 'sshn', ssh a)155 CALL iom_rstput( kt, nitrst, numrtw, 'sshb', ssh n)159 CALL iom_rstput( kt, nitrst, numrtw, 'sshn', ssh(:,:,Kaa) ) 160 CALL iom_rstput( kt, nitrst, numrtw, 'sshb', ssh(:,:,Kmm) ) 156 161 ENDIF 157 162 ENDIF 158 163 ! 159 CALL eos ( ts n, rhd, rhop, gdept_0(:,:,:) ) ! In any case, we need rhop160 CALL eos_rab( ts n, rab_n) ! now local thermal/haline expension ratio at T-points161 CALL bn2 ( ts n, rab_n, rn2 )! before Brunt-Vaisala frequency need for zdfmxl162 163 rn2b(:,:,:) = rn2(:,:,:) ! need for zdfmxl164 CALL zdf_mxl( kt )! In any case, we need mxl164 CALL eos ( ts(:,:,:,:,Kmm), rhd, rhop, gdept_0(:,:,:) ) ! In any case, we need rhop 165 CALL eos_rab( ts(:,:,:,:,Kmm), rab_n, Kmm ) ! now local thermal/haline expension ratio at T-points 166 CALL bn2 ( ts(:,:,:,:,Kmm), rab_n, rn2, Kmm ) ! before Brunt-Vaisala frequency need for zdfmxl 167 168 rn2b(:,:,:) = rn2(:,:,:) ! needed for zdfmxl 169 CALL zdf_mxl( kt, Kmm ) ! In any case, we need mxl 165 170 ! 166 171 hmld(:,:) = sf_dyn(jf_mld)%fnow(:,:,1) * tmask(:,:,1) ! mixed layer depht … … 174 179 ! 175 180 ! 176 CALL eos( ts n, rhd, rhop, gdept_0(:,:,:) ) ! In any case, we need rhop177 ! 178 IF( ln_ctl) THEN! print control179 CALL prt_ctl(tab3d_1=ts n(:,:,:,jp_tem), clinfo1=' tn - : ', mask1=tmask, kdim=jpk )180 CALL prt_ctl(tab3d_1=ts n(:,:,:,jp_sal), clinfo1=' sn - : ', mask1=tmask, kdim=jpk )181 CALL prt_ctl(tab3d_1=u n , clinfo1=' un- : ', mask1=umask, kdim=jpk )182 CALL prt_ctl(tab3d_1=v n , clinfo1=' vn- : ', mask1=vmask, kdim=jpk )183 CALL prt_ctl(tab3d_1=w n , clinfo1=' wn- : ', mask1=tmask, kdim=jpk )181 CALL eos( ts(:,:,:,:,Kmm), rhd, rhop, gdept_0(:,:,:) ) ! In any case, we need rhop 182 ! 183 IF(sn_cfctl%l_prtctl) THEN ! print control 184 CALL prt_ctl(tab3d_1=ts(:,:,:,jp_tem,Kmm), clinfo1=' tn - : ', mask1=tmask, kdim=jpk ) 185 CALL prt_ctl(tab3d_1=ts(:,:,:,jp_sal,Kmm), clinfo1=' sn - : ', mask1=tmask, kdim=jpk ) 186 CALL prt_ctl(tab3d_1=uu(:,:,:,Kmm) , clinfo1=' uu(:,:,:,Kmm) - : ', mask1=umask, kdim=jpk ) 187 CALL prt_ctl(tab3d_1=vv(:,:,:,Kmm) , clinfo1=' vv(:,:,:,Kmm) - : ', mask1=vmask, kdim=jpk ) 188 CALL prt_ctl(tab3d_1=ww , clinfo1=' ww - : ', mask1=tmask, kdim=jpk ) 184 189 CALL prt_ctl(tab3d_1=avt , clinfo1=' kz - : ', mask1=tmask, kdim=jpk ) 185 190 CALL prt_ctl(tab3d_1=uslp , clinfo1=' slp - u : ', tab3d_2=vslp, clinfo2=' v : ', kdim=jpk) … … 192 197 193 198 194 SUBROUTINE dta_dyn_init 199 SUBROUTINE dta_dyn_init( Kbb, Kmm, Kaa ) 195 200 !!---------------------------------------------------------------------- 196 201 !! *** ROUTINE dta_dyn_init *** … … 199 204 !! ** Method : - read the data namdta_dyn namelist 200 205 !!---------------------------------------------------------------------- 206 INTEGER, INTENT(in) :: Kbb, Kmm, Kaa ! ocean time level indices 207 ! 201 208 INTEGER :: ierr, ierr0, ierr1, ierr2, ierr3 ! return error code 202 209 INTEGER :: ifpr ! dummy loop indice … … 225 232 !!---------------------------------------------------------------------- 226 233 ! 227 REWIND( numnam_ref ) ! Namelist namdta_dyn in reference namelist : Offline: init. of dynamical data228 234 READ ( numnam_ref, namdta_dyn, IOSTAT = ios, ERR = 901) 229 235 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdta_dyn in reference namelist' ) 230 REWIND( numnam_cfg ) ! Namelist namdta_dyn in configuration namelist : Offline: init. of dynamical data231 236 READ ( numnam_cfg, namdta_dyn, IOSTAT = ios, ERR = 902 ) 232 237 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namdta_dyn in configuration namelist' ) … … 281 286 ! Open file for each variable to get his number of dimension 282 287 DO ifpr = 1, jfld 283 CALL fld_clopn( sf_dyn(ifpr), nyear, nmonth, nday ) 288 CALL fld_def( sf_dyn(ifpr) ) 289 CALL iom_open( sf_dyn(ifpr)%clname, sf_dyn(ifpr)%num ) 284 290 idv = iom_varid( sf_dyn(ifpr)%num , slf_d(ifpr)%clvar ) ! id of the variable sdjf%clvar 285 291 idimv = iom_file ( sf_dyn(ifpr)%num )%ndims(idv) ! number of dimension for variable sdjf%clvar 286 IF( sf_dyn(ifpr)%num /= 0 ) CALL iom_close( sf_dyn(ifpr)%num )! close file if already open292 CALL iom_close( sf_dyn(ifpr)%num ) ! close file if already open 287 293 ierr1=0 288 294 IF( idimv == 3 ) THEN ! 2D variable … … 312 318 IF( .NOT. sf_dyn(jf_uwd)%ln_clim .AND. ln_rsttr .AND. & ! Restart: read in restart file 313 319 iom_varid( numrtr, 'sshn', ldstop = .FALSE. ) > 0 ) THEN 314 IF(lwp) WRITE(numout,*) ' ssh nforcing fields read in the restart file for initialisation'315 CALL iom_get( numrtr, jpdom_autoglo, 'sshn', ssh n(:,:) )316 CALL iom_get( numrtr, jpdom_autoglo, 'sshb', ssh b(:,:) )320 IF(lwp) WRITE(numout,*) ' ssh(:,:,Kmm) forcing fields read in the restart file for initialisation' 321 CALL iom_get( numrtr, jpdom_autoglo, 'sshn', ssh(:,:,Kmm) ) 322 CALL iom_get( numrtr, jpdom_autoglo, 'sshb', ssh(:,:,Kbb) ) 317 323 ELSE 318 IF(lwp) WRITE(numout,*) ' ssh nforcing fields read in the restart file for initialisation'324 IF(lwp) WRITE(numout,*) ' ssh(:,:,Kmm) forcing fields read in the restart file for initialisation' 319 325 CALL iom_open( 'restart', inum ) 320 CALL iom_get( inum, jpdom_autoglo, 'sshn', ssh n(:,:) )321 CALL iom_get( inum, jpdom_autoglo, 'sshb', ssh b(:,:) )326 CALL iom_get( inum, jpdom_autoglo, 'sshn', ssh(:,:,Kmm) ) 327 CALL iom_get( inum, jpdom_autoglo, 'sshb', ssh(:,:,Kbb) ) 322 328 CALL iom_close( inum ) ! close file 323 329 ENDIF 324 330 ! 325 331 DO jk = 1, jpkm1 326 e3t _n(:,:,jk) = e3t_0(:,:,jk) * ( 1._wp + sshn(:,:) * tmask(:,:,1) / ( ht_0(:,:) + 1.0 - tmask(:,:,1) ) )332 e3t(:,:,jk,Kmm) = e3t_0(:,:,jk) * ( 1._wp + ssh(:,:,Kmm) * tmask(:,:,1) / ( ht_0(:,:) + 1.0 - tmask(:,:,1) ) ) 327 333 ENDDO 328 e3t _a(:,:,jpk) = e3t_0(:,:,jpk)334 e3t(:,:,jpk,Kaa) = e3t_0(:,:,jpk) 329 335 330 336 ! Horizontal scale factor interpolations 331 337 ! -------------------------------------- 332 CALL dom_vvl_interpol( e3t _n(:,:,:), e3u_n(:,:,:), 'U' )333 CALL dom_vvl_interpol( e3t _n(:,:,:), e3v_n(:,:,:), 'V' )338 CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3u(:,:,:,Kmm), 'U' ) 339 CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3v(:,:,:,Kmm), 'V' ) 334 340 335 341 ! Vertical scale factor interpolations 336 342 ! ------------------------------------ 337 CALL dom_vvl_interpol( e3t _n(:,:,:), e3w_n(:,:,:), 'W' )343 CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3w(:,:,:,Kmm), 'W' ) 338 344 339 e3t _b(:,:,:) = e3t_n(:,:,:)340 e3u _b(:,:,:) = e3u_n(:,:,:)341 e3v _b(:,:,:) = e3v_n(:,:,:)345 e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm) 346 e3u(:,:,:,Kbb) = e3u(:,:,:,Kmm) 347 e3v(:,:,:,Kbb) = e3v(:,:,:,Kmm) 342 348 343 349 ! t- and w- points depth 344 350 ! ---------------------- 345 gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1) 346 gdepw_n(:,:,1) = 0.0_wp 347 348 DO jk = 2, jpk 349 DO jj = 1,jpj 350 DO ji = 1,jpi 351 ! zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) ! 0 everywhere 352 ! tmask = wmask, ie everywhere expect at jk = mikt 353 ! 1 for jk = 354 ! mikt 355 zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 356 gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1) 357 gdept_n(ji,jj,jk) = zcoef * ( gdepw_n(ji,jj,jk ) + 0.5 * e3w_n(ji,jj,jk)) & 358 & + (1-zcoef) * ( gdept_n(ji,jj,jk-1) + e3w_n(ji,jj,jk)) 359 END DO 360 END DO 361 END DO 362 363 gdept_b(:,:,:) = gdept_n(:,:,:) 364 gdepw_b(:,:,:) = gdepw_n(:,:,:) 351 gdept(:,:,1,Kmm) = 0.5_wp * e3w(:,:,1,Kmm) 352 gdepw(:,:,1,Kmm) = 0.0_wp 353 354 DO_3D_11_11( 2, jpk ) 355 ! zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) ! 0 everywhere 356 ! tmask = wmask, ie everywhere expect at jk = mikt 357 ! 1 for jk = 358 ! mikt 359 zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 360 gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm) 361 gdept(ji,jj,jk,Kmm) = zcoef * ( gdepw(ji,jj,jk ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm)) & 362 & + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) + e3w(ji,jj,jk,Kmm)) 363 END_3D 364 365 gdept(:,:,:,Kbb) = gdept(:,:,:,Kmm) 366 gdepw(:,:,:,Kbb) = gdepw(:,:,:,Kmm) 365 367 ! 366 368 ENDIF … … 374 376 ! 375 377 nk_rnf(:,:) = 0 ! set the number of level over which river runoffs are applied 376 DO jj = 1, jpj 377 DO ji = 1, jpi 378 IF( h_rnf(ji,jj) > 0._wp ) THEN 379 jk = 2 380 DO WHILE ( jk /= mbkt(ji,jj) .AND. gdept_0(ji,jj,jk) < h_rnf(ji,jj) ) ; jk = jk + 1 381 END DO 382 nk_rnf(ji,jj) = jk 383 ELSEIF( h_rnf(ji,jj) == -1._wp ) THEN ; nk_rnf(ji,jj) = 1 384 ELSEIF( h_rnf(ji,jj) == -999._wp ) THEN ; nk_rnf(ji,jj) = mbkt(ji,jj) 385 ELSE 386 CALL ctl_stop( 'sbc_rnf_init: runoff depth not positive, and not -999 or -1, rnf value in file fort.999' ) 387 WRITE(999,*) 'ji, jj, h_rnf(ji,jj) :', ji, jj, h_rnf(ji,jj) 388 ENDIF 378 DO_2D_11_11 379 IF( h_rnf(ji,jj) > 0._wp ) THEN 380 jk = 2 381 DO WHILE ( jk /= mbkt(ji,jj) .AND. gdept_0(ji,jj,jk) < h_rnf(ji,jj) ) ; jk = jk + 1 382 END DO 383 nk_rnf(ji,jj) = jk 384 ELSEIF( h_rnf(ji,jj) == -1._wp ) THEN ; nk_rnf(ji,jj) = 1 385 ELSEIF( h_rnf(ji,jj) == -999._wp ) THEN ; nk_rnf(ji,jj) = mbkt(ji,jj) 386 ELSE 387 CALL ctl_stop( 'sbc_rnf_init: runoff depth not positive, and not -999 or -1, rnf value in file fort.999' ) 388 WRITE(999,*) 'ji, jj, h_rnf(ji,jj) :', ji, jj, h_rnf(ji,jj) 389 ENDIF 390 END_2D 391 DO_2D_11_11 392 h_rnf(ji,jj) = 0._wp 393 DO jk = 1, nk_rnf(ji,jj) 394 h_rnf(ji,jj) = h_rnf(ji,jj) + e3t(ji,jj,jk,Kmm) 389 395 END DO 390 END DO 391 DO jj = 1, jpj ! set the associated depth 392 DO ji = 1, jpi 393 h_rnf(ji,jj) = 0._wp 394 DO jk = 1, nk_rnf(ji,jj) 395 h_rnf(ji,jj) = h_rnf(ji,jj) + e3t_n(ji,jj,jk) 396 END DO 397 END DO 398 END DO 396 END_2D 399 397 ELSE ! runoffs applied at the surface 400 398 nk_rnf(:,:) = 1 401 h_rnf (:,:) = e3t _n(:,:,1)399 h_rnf (:,:) = e3t(:,:,1,Kmm) 402 400 ENDIF 403 401 nkrnf_max = MAXVAL( nk_rnf(:,:) ) … … 411 409 IF(lwp) WRITE(numout,*) ' ' 412 410 ! 413 CALL dta_dyn( nit000 )411 CALL dta_dyn( nit000, Kbb, Kmm, Kaa ) 414 412 ! 415 413 END SUBROUTINE dta_dyn_init 416 414 417 SUBROUTINE dta_dyn_sed( kt )415 SUBROUTINE dta_dyn_sed( kt, Kmm ) 418 416 !!---------------------------------------------------------------------- 419 417 !! *** ROUTINE dta_dyn *** … … 427 425 !!---------------------------------------------------------------------- 428 426 INTEGER, INTENT(in) :: kt ! ocean time-step index 427 INTEGER, INTENT(in) :: Kmm ! ocean time level index 429 428 ! 430 429 !!---------------------------------------------------------------------- … … 439 438 CALL fld_read( kt, 1, sf_dyn ) != read data at kt time step ==! 440 439 ! 441 ts n(:,:,:,jp_tem) = sf_dyn(jf_tem)%fnow(:,:,:) * tmask(:,:,:) ! temperature442 ts n(:,:,:,jp_sal) = sf_dyn(jf_sal)%fnow(:,:,:) * tmask(:,:,:) ! salinity443 ! 444 CALL eos ( ts n, rhd, rhop, gdept_0(:,:,:) ) ! In any case, we need rhop445 446 IF( ln_ctl) THEN! print control447 CALL prt_ctl(tab3d_1=ts n(:,:,:,jp_tem), clinfo1=' tn - : ', mask1=tmask, kdim=jpk )448 CALL prt_ctl(tab3d_1=ts n(:,:,:,jp_sal), clinfo1=' sn - : ', mask1=tmask, kdim=jpk )440 ts(:,:,:,jp_tem,Kmm) = sf_dyn(jf_tem)%fnow(:,:,:) * tmask(:,:,:) ! temperature 441 ts(:,:,:,jp_sal,Kmm) = sf_dyn(jf_sal)%fnow(:,:,:) * tmask(:,:,:) ! salinity 442 ! 443 CALL eos ( ts(:,:,:,:,Kmm), rhd, rhop, gdept_0(:,:,:) ) ! In any case, we need rhop 444 445 IF(sn_cfctl%l_prtctl) THEN ! print control 446 CALL prt_ctl(tab3d_1=ts(:,:,:,jp_tem,Kmm), clinfo1=' tn - : ', mask1=tmask, kdim=jpk ) 447 CALL prt_ctl(tab3d_1=ts(:,:,:,jp_sal,Kmm), clinfo1=' sn - : ', mask1=tmask, kdim=jpk ) 449 448 ENDIF 450 449 ! … … 454 453 455 454 456 SUBROUTINE dta_dyn_sed_init 455 SUBROUTINE dta_dyn_sed_init( Kmm ) 457 456 !!---------------------------------------------------------------------- 458 457 !! *** ROUTINE dta_dyn_init *** … … 461 460 !! ** Method : - read the data namdta_dyn namelist 462 461 !!---------------------------------------------------------------------- 462 INTEGER, INTENT( in ) :: Kmm ! ocean time level index 463 ! 463 464 INTEGER :: ierr, ierr0, ierr1, ierr2, ierr3 ! return error code 464 465 INTEGER :: ifpr ! dummy loop indice … … 475 476 !!---------------------------------------------------------------------- 476 477 ! 477 REWIND( numnam_ref ) ! Namelist namdta_dyn in reference namelist : Offline: init. of dynamical data478 478 READ ( numnam_ref, namdta_dyn, IOSTAT = ios, ERR = 901) 479 479 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdta_dyn in reference namelist' ) 480 REWIND( numnam_cfg ) ! Namelist namdta_dyn in configuration namelist : Offline: init. of dynamical data481 480 READ ( numnam_cfg, namdta_dyn, IOSTAT = ios, ERR = 902 ) 482 481 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namdta_dyn in configuration namelist' ) … … 508 507 ! Open file for each variable to get his number of dimension 509 508 DO ifpr = 1, jfld 510 CALL fld_clopn( sf_dyn(ifpr), nyear, nmonth, nday ) 509 CALL fld_def( sf_dyn(ifpr) ) 510 CALL iom_open( sf_dyn(ifpr)%clname, sf_dyn(ifpr)%num ) 511 511 idv = iom_varid( sf_dyn(ifpr)%num , slf_d(ifpr)%clvar ) ! id of the variable sdjf%clvar 512 512 idimv = iom_file ( sf_dyn(ifpr)%num )%ndims(idv) ! number of dimension for variable sdjf%clvar 513 IF( sf_dyn(ifpr)%num /= 0 ) CALL iom_close( sf_dyn(ifpr)%num )! close file if already open513 CALL iom_close( sf_dyn(ifpr)%num ) ! close file if already open 514 514 ierr1=0 515 515 IF( idimv == 3 ) THEN ! 2D variable … … 525 525 END DO 526 526 ! 527 CALL dta_dyn_sed( nit000 )527 CALL dta_dyn_sed( nit000, Kmm ) 528 528 ! 529 529 END SUBROUTINE dta_dyn_sed_init 530 530 531 SUBROUTINE dta_dyn_ swp( kt)531 SUBROUTINE dta_dyn_atf( kt, Kbb, Kmm, Kaa ) 532 532 !!--------------------------------------------------------------------- 533 533 !! *** ROUTINE dta_dyn_swp *** 534 534 !! 535 !! ** Purpose : Swap and the data and compute the vertical scale factor 536 !! at U/V/W pointand the depht 537 !!--------------------------------------------------------------------- 538 INTEGER, INTENT(in) :: kt ! time step 535 !! ** Purpose : Asselin time filter of now SSH 536 !!--------------------------------------------------------------------- 537 INTEGER, INTENT(in) :: kt ! time step 538 INTEGER, INTENT(in) :: Kbb, Kmm, Kaa ! ocean time level indices 539 ! 540 !!--------------------------------------------------------------------- 541 542 IF( kt == nit000 ) THEN 543 IF(lwp) WRITE(numout,*) 544 IF(lwp) WRITE(numout,*) 'dta_dyn_atf : Asselin time filter of sea surface height' 545 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ ' 546 ENDIF 547 548 ssh(:,:,Kmm) = ssh(:,:,Kmm) + rn_atfp * ( ssh(:,:,Kbb) - 2 * ssh(:,:,Kmm) + ssh(:,:,Kaa)) 549 550 !! Do we also need to time filter e3t?? 551 ! 552 END SUBROUTINE dta_dyn_atf 553 554 SUBROUTINE dta_dyn_sf_interp( kt, Kmm ) 555 !!--------------------------------------------------------------------- 556 !! *** ROUTINE dta_dyn_sf_interp *** 557 !! 558 !! ** Purpose : Calculate scale factors at U/V/W points and depths 559 !! given the after e3t field 560 !!--------------------------------------------------------------------- 561 INTEGER, INTENT(in) :: kt ! time step 562 INTEGER, INTENT(in) :: Kmm ! ocean time level indices 539 563 ! 540 564 INTEGER :: ji, jj, jk … … 542 566 !!--------------------------------------------------------------------- 543 567 544 IF( kt == nit000 ) THEN545 IF(lwp) WRITE(numout,*)546 IF(lwp) WRITE(numout,*) 'ssh_swp : Asselin time filter and swap of sea surface height'547 IF(lwp) WRITE(numout,*) '~~~~~~~ '548 ENDIF549 550 sshb(:,:) = sshn(:,:) + atfp * ( sshb(:,:) - 2 * sshn(:,:) + ssha(:,:)) ! before <-- now filtered551 sshn(:,:) = ssha(:,:)552 553 e3t_n(:,:,:) = e3t_a(:,:,:)554 555 ! Reconstruction of all vertical scale factors at now and before time steps556 ! =============================================================================557 558 568 ! Horizontal scale factor interpolations 559 569 ! -------------------------------------- 560 CALL dom_vvl_interpol( e3t _n(:,:,:), e3u_n(:,:,:), 'U' )561 CALL dom_vvl_interpol( e3t _n(:,:,:), e3v_n(:,:,:), 'V' )570 CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3u(:,:,:,Kmm), 'U' ) 571 CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3v(:,:,:,Kmm), 'V' ) 562 572 563 573 ! Vertical scale factor interpolations 564 574 ! ------------------------------------ 565 CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n (:,:,:), 'W' ) 566 567 e3t_b(:,:,:) = e3t_n(:,:,:) 568 e3u_b(:,:,:) = e3u_n(:,:,:) 569 e3v_b(:,:,:) = e3v_n(:,:,:) 575 CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3w (:,:,:,Kmm), 'W' ) 570 576 571 577 ! t- and w- points depth 572 578 ! ---------------------- 573 gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1) 574 gdepw_n(:,:,1) = 0.0_wp 575 ! 576 DO jk = 2, jpk 577 DO jj = 1,jpj 578 DO ji = 1,jpi 579 zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 580 gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1) 581 gdept_n(ji,jj,jk) = zcoef * ( gdepw_n(ji,jj,jk ) + 0.5 * e3w_n(ji,jj,jk)) & 582 & + (1-zcoef) * ( gdept_n(ji,jj,jk-1) + e3w_n(ji,jj,jk)) 583 END DO 584 END DO 585 END DO 586 ! 587 gdept_b(:,:,:) = gdept_n(:,:,:) 588 gdepw_b(:,:,:) = gdepw_n(:,:,:) 589 ! 590 END SUBROUTINE dta_dyn_swp 591 579 gdept(:,:,1,Kmm) = 0.5_wp * e3w(:,:,1,Kmm) 580 gdepw(:,:,1,Kmm) = 0.0_wp 581 ! 582 DO_3D_11_11( 2, jpk ) 583 zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 584 gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm) 585 gdept(ji,jj,jk,Kmm) = zcoef * ( gdepw(ji,jj,jk ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm)) & 586 & + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) + e3w(ji,jj,jk,Kmm)) 587 END_3D 588 ! 589 END SUBROUTINE dta_dyn_sf_interp 592 590 593 591 SUBROUTINE dta_dyn_ssh( kt, phdivtr, psshb, pemp, pssha, pe3ta ) … … 595 593 !! *** ROUTINE dta_dyn_wzv *** 596 594 !! 597 !! ** Purpose : compute the after ssh (ssh a) and the now vertical velocity595 !! ** Purpose : compute the after ssh (ssh(:,:,Kaa)) and the now vertical velocity 598 596 !! 599 597 !! ** Method : Using the incompressibility hypothesis, … … 608 606 !! The boundary conditions are w=0 at the bottom (no flux) 609 607 !! 610 !! ** action : ssh a / e3t_a / wn608 !! ** action : ssh(:,:,Kaa) / e3t(:,:,:,Kaa) / ww 611 609 !! 612 610 !! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling. … … 624 622 !!---------------------------------------------------------------------- 625 623 ! 626 z2dt = 2._wp * r dt624 z2dt = 2._wp * rn_Dt 627 625 ! 628 626 zhdiv(:,:) = 0._wp … … 631 629 END DO 632 630 ! ! Sea surface elevation time-stepping 633 pssha(:,:) = ( psshb(:,:) - z2dt * ( r1_r au0 * pemp(:,:) + zhdiv(:,:) ) ) * ssmask(:,:)631 pssha(:,:) = ( psshb(:,:) - z2dt * ( r1_rho0 * pemp(:,:) + zhdiv(:,:) ) ) * ssmask(:,:) 634 632 ! ! 635 633 ! ! After acale factors at t-points ( z_star coordinate ) … … 641 639 642 640 643 SUBROUTINE dta_dyn_hrnf 641 SUBROUTINE dta_dyn_hrnf( Kmm ) 644 642 !!---------------------------------------------------------------------- 645 643 !! *** ROUTINE sbc_rnf *** … … 654 652 !!---------------------------------------------------------------------- 655 653 !! 656 INTEGER :: ji, jj, jk ! dummy loop indices657 ! !----------------------------------------------------------------------658 !659 DO jj = 1, jpj ! update the depth over which runoffs are distributed660 DO ji = 1, jpi661 h_rnf(ji,jj) = 0._wp662 DO jk = 1, nk_rnf(ji,jj) ! recalculates h_rnf to be the depth in metres663 h_rnf(ji,jj) = h_rnf(ji,jj) + e3t_n(ji,jj,jk) ! to the bottom of the relevant grid box664 END DO665 END DO666 END DO654 INTEGER, INTENT(in) :: Kmm ! ocean time level index 655 ! 656 INTEGER :: ji, jj, jk ! dummy loop indices 657 !!---------------------------------------------------------------------- 658 ! 659 DO_2D_11_11 660 h_rnf(ji,jj) = 0._wp 661 DO jk = 1, nk_rnf(ji,jj) ! recalculates h_rnf to be the depth in metres 662 h_rnf(ji,jj) = h_rnf(ji,jj) + e3t(ji,jj,jk,Kmm) ! to the bottom of the relevant grid box 663 END DO 664 END_2D 667 665 ! 668 666 END SUBROUTINE dta_dyn_hrnf … … 670 668 671 669 672 SUBROUTINE dta_dyn_slp( kt )670 SUBROUTINE dta_dyn_slp( kt, Kbb, Kmm ) 673 671 !!--------------------------------------------------------------------- 674 672 !! *** ROUTINE dta_dyn_slp *** … … 678 676 !!--------------------------------------------------------------------- 679 677 INTEGER, INTENT(in) :: kt ! time step 678 INTEGER, INTENT(in) :: Kbb, Kmm ! ocean time level indices 680 679 ! 681 680 INTEGER :: ji, jj ! dummy loop indices … … 693 692 zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fdta(:,:,:,1) * tmask(:,:,:) ! salinity 694 693 avt(:,:,:) = sf_dyn(jf_avt)%fdta(:,:,:,1) * tmask(:,:,:) ! vertical diffusive coef. 695 CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj )694 CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj, Kbb, Kmm ) 696 695 uslpdta (:,:,:,1) = zuslp (:,:,:) 697 696 vslpdta (:,:,:,1) = zvslp (:,:,:) … … 702 701 zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fdta(:,:,:,2) * tmask(:,:,:) ! salinity 703 702 avt(:,:,:) = sf_dyn(jf_avt)%fdta(:,:,:,2) * tmask(:,:,:) ! vertical diffusive coef. 704 CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj )703 CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj, Kbb, Kmm ) 705 704 uslpdta (:,:,:,2) = zuslp (:,:,:) 706 705 vslpdta (:,:,:,2) = zvslp (:,:,:) … … 721 720 zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fdta(:,:,:,2) * tmask(:,:,:) ! salinity 722 721 avt(:,:,:) = sf_dyn(jf_avt)%fdta(:,:,:,2) * tmask(:,:,:) ! vertical diffusive coef. 723 CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj )722 CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj, Kbb, Kmm ) 724 723 ! 725 724 uslpdta (:,:,:,2) = zuslp (:,:,:) … … 745 744 zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fnow(:,:,:) * tmask(:,:,:) ! salinity 746 745 avt(:,:,:) = sf_dyn(jf_avt)%fnow(:,:,:) * tmask(:,:,:) ! vertical diffusive coef. 747 CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj )746 CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj, Kbb, Kmm ) 748 747 ! 749 748 IF( l_ldfslp .AND. .NOT.lk_c1d ) THEN ! Computes slopes (here avt is used as workspace) … … 758 757 759 758 760 SUBROUTINE compute_slopes( kt, pts, puslp, pvslp, pwslpi, pwslpj )759 SUBROUTINE compute_slopes( kt, pts, puslp, pvslp, pwslpi, pwslpj, Kbb, Kmm ) 761 760 !!--------------------------------------------------------------------- 762 761 !! *** ROUTINE dta_dyn_slp *** … … 770 769 REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT(out) :: pwslpi ! zonal diapycnal slopes 771 770 REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT(out) :: pwslpj ! meridional diapycnal slopes 771 INTEGER , INTENT(in ) :: Kbb, Kmm ! ocean time level indices 772 772 !!--------------------------------------------------------------------- 773 773 ! 774 774 IF( l_ldfslp .AND. .NOT.lk_c1d ) THEN ! Computes slopes (here avt is used as workspace) 775 775 CALL eos ( pts, rhd, rhop, gdept_0(:,:,:) ) 776 CALL eos_rab( pts, rab_n ) ! now local thermal/haline expension ratio at T-points777 CALL bn2 ( pts, rab_n, rn2 ) ! now Brunt-Vaisala776 CALL eos_rab( pts, rab_n, Kmm ) ! now local thermal/haline expension ratio at T-points 777 CALL bn2 ( pts, rab_n, rn2, Kmm ) ! now Brunt-Vaisala 778 778 779 779 ! Partial steps: before Horizontal DErivative 780 780 IF( ln_zps .AND. .NOT. ln_isfcav) & 781 & CALL zps_hde ( kt, jpts, pts, gtsu, gtsv, & ! Partial steps: before horizontal gradient781 & CALL zps_hde ( kt, Kmm, jpts, pts, gtsu, gtsv, & ! Partial steps: before horizontal gradient 782 782 & rhd, gru , grv ) ! of t, s, rd at the last ocean level 783 783 IF( ln_zps .AND. ln_isfcav) & 784 & CALL zps_hde_isf( kt, jpts, pts, gtsu, gtsv, gtui, gtvi, & ! Partial steps for top cell (ISF)784 & CALL zps_hde_isf( kt, Kmm, jpts, pts, gtsu, gtsv, gtui, gtvi, & ! Partial steps for top cell (ISF) 785 785 & rhd, gru , grv , grui, grvi ) ! of t, s, rd at the first ocean level 786 786 787 rn2b(:,:,:) = rn2(:,:,:) ! need for zdfmxl788 CALL zdf_mxl( kt )! mixed layer depth789 CALL ldf_slp( kt, rhd, rn2 ) ! slopes787 rn2b(:,:,:) = rn2(:,:,:) ! needed for zdfmxl 788 CALL zdf_mxl( kt, Kmm ) ! mixed layer depth 789 CALL ldf_slp( kt, rhd, rn2, Kbb, Kmm ) ! slopes 790 790 puslp (:,:,:) = uslp (:,:,:) 791 791 pvslp (:,:,:) = vslp (:,:,:)
Note: See TracChangeset
for help on using the changeset viewer.
