Changeset 14789 for NEMO/branches/2021/dev_r13747_HPC-11_mcastril_HPDAonline_DiagGPU/src/OFF/dtadyn.F90
- Timestamp:
- 2021-05-05T13:18:04+02:00 (3 years ago)
- Location:
- NEMO/branches/2021/dev_r13747_HPC-11_mcastril_HPDAonline_DiagGPU
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
src/OFF/dtadyn.F90 (modified) (16 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2021/dev_r13747_HPC-11_mcastril_HPDAonline_DiagGPU
- Property svn:externals
-
old new 3 3 ^/utils/build/mk@HEAD mk 4 4 ^/utils/tools@HEAD tools 5 ^/vendors/AGRIF/dev _r12970_AGRIF_CMEMSext/AGRIF5 ^/vendors/AGRIF/dev@HEAD ext/AGRIF 6 6 ^/vendors/FCM@HEAD ext/FCM 7 7 ^/vendors/IOIPSL@HEAD ext/IOIPSL 8 ^/vendors/PPR@HEAD ext/PPR 8 9 9 10 # SETTE 10 ^/utils/CI/sette@1 3559sette11 ^/utils/CI/sette@14244 sette
-
- Property svn:externals
-
NEMO/branches/2021/dev_r13747_HPC-11_mcastril_HPDAonline_DiagGPU/src/OFF/dtadyn.F90
r13497 r14789 23 23 USE c1d ! 1D configuration: lk_c1d 24 24 USE dom_oce ! ocean domain: variables 25 #if !defined key_qco26 USE dom vvl! variable volume25 #if defined key_qco 26 USE domqco ! variable volume 27 27 #else 28 USE dom qco28 USE domvvl 29 29 #endif 30 30 USE zdf_oce ! ocean vertical physics: variables … … 46 46 USE fldread ! read input fields 47 47 USE timing ! Timing 48 USE trc, ONLY : ln_rsttr, numrtr, numrtw,lrst_trc48 USE trc, ONLY : ln_rsttr, lrst_trc 49 49 50 50 IMPLICIT NONE … … 53 53 PUBLIC dta_dyn_init ! called by nemo_init 54 54 PUBLIC dta_dyn ! called by nemo_gcm 55 PUBLIC dta_dyn_sed_init ! called by nemo_init56 PUBLIC dta_dyn_sed ! called by nemo_gcm57 55 PUBLIC dta_dyn_atf ! called by nemo_gcm 58 56 #if ! defined key_qco 59 57 PUBLIC dta_dyn_sf_interp ! called by nemo_gcm 60 58 #endif 59 #if defined key_sed_off 60 PUBLIC dta_dyn_sed_init ! called by nemo_init 61 PUBLIC dta_dyn_sed ! called by nemo_gcm 62 #endif 61 63 62 64 CHARACTER(len=100) :: cn_dir !: Root directory for location of ssr files 63 65 LOGICAL :: ln_dynrnf !: read runoff data in file (T) or set to zero (F) 64 66 LOGICAL :: ln_dynrnf_depth !: read runoff data in file (T) or set to zero (F) 65 REAL(wp) :: fwbcorr66 67 67 68 … … 97 98 !! * Substitutions 98 99 # include "do_loop_substitute.h90" 100 # include "domzgr_substitute.h90" 101 99 102 !!---------------------------------------------------------------------- 100 103 !! NEMO/OFF 4.0 , NEMO Consortium (2018) … … 136 139 ts(:,:,:,jp_tem,Kmm) = sf_dyn(jf_tem)%fnow(:,:,:) * tmask(:,:,:) ! temperature 137 140 ts(:,:,:,jp_sal,Kmm) = sf_dyn(jf_sal)%fnow(:,:,:) * tmask(:,:,:) ! salinity 138 wndm(:,:) = sf_dyn(jf_wnd)%fnow(:,:,1) * tmask(:,:,1) ! wind speed - needed for gas exchange139 fmmflx(:,:) = sf_dyn(jf_fmf)%fnow(:,:,1) * tmask(:,:,1) ! downward salt flux (v3.5+)140 fr_i(:,:) = sf_dyn(jf_ice)%fnow(:,:,1) * tmask(:,:,1) ! Sea-ice fraction141 qsr (:,:) = sf_dyn(jf_qsr)%fnow(:,:,1) * tmask(:,:,1) ! solar radiation142 emp (:,:) = sf_dyn(jf_emp)%fnow(:,:,1) * tmask(:,:,1) ! E-P141 wndm(:,:) = sf_dyn(jf_wnd)%fnow(:,:,1) * tmask(:,:,1) ! wind speed - needed for gas exchange 142 fmmflx(:,:) = sf_dyn(jf_fmf)%fnow(:,:,1) * tmask(:,:,1) ! downward salt flux (v3.5+) 143 fr_i(:,:) = sf_dyn(jf_ice)%fnow(:,:,1) * tmask(:,:,1) ! Sea-ice fraction 144 qsr (:,:) = sf_dyn(jf_qsr)%fnow(:,:,1) * tmask(:,:,1) ! solar radiation 145 emp (:,:) = sf_dyn(jf_emp)%fnow(:,:,1) * tmask(:,:,1) ! E-P 143 146 IF( ln_dynrnf ) THEN 144 rnf (:,:) = sf_dyn(jf_rnf)%fnow(:,:,1) * tmask(:,:,1) ! E-P145 IF( ln_dynrnf_depth .AND. .NOT. ln_linssh ) CALL dta_dyn_ hrnf(Kmm)147 rnf (:,:) = sf_dyn(jf_rnf)%fnow(:,:,1) * tmask(:,:,1) ! E-P 148 IF( ln_dynrnf_depth .AND. .NOT. ln_linssh ) CALL dta_dyn_rnf( Kmm ) 146 149 ENDIF 147 150 ! 148 151 uu(:,:,:,Kmm) = sf_dyn(jf_uwd)%fnow(:,:,:) * umask(:,:,:) ! effective u-transport 149 152 vv(:,:,:,Kmm) = sf_dyn(jf_vwd)%fnow(:,:,:) * vmask(:,:,:) ! effective v-transport 150 ww(:,:,:) = sf_dyn(jf_wwd)%fnow(:,:,:) * tmask(:,:,:) ! effective v-transport153 ww(:,:,:) = sf_dyn(jf_wwd)%fnow(:,:,:) * tmask(:,:,:) ! effective v-transport 151 154 ! 152 155 IF( .NOT.ln_linssh ) THEN … … 154 157 zhdivtr(:,:,:) = sf_dyn(jf_div)%fnow(:,:,:) * tmask(:,:,:) ! effective u-transport 155 158 emp_b (:,:) = sf_dyn(jf_empb)%fnow(:,:,1) * tmask(:,:,1) ! E-P 156 zemp (:,:) = ( 0.5_wp * ( emp(:,:) + emp_b(:,:) ) + rnf(:,:) + fwbcorr) * tmask(:,:,1)159 zemp (:,:) = ( 0.5_wp * ( emp(:,:) + emp_b(:,:) ) + rnf(:,:) ) * tmask(:,:,1) 157 160 #if defined key_qco 158 161 CALL dta_dyn_ssh( kt, zhdivtr, ssh(:,:,Kbb), zemp, ssh(:,:,Kaa) ) … … 223 226 INTEGER :: ios ! Local integer output status for namelist read 224 227 INTEGER :: ji, jj, jk 225 REAL(wp) :: zcoef226 INTEGER :: nkrnf_max227 REAL(wp) :: hrnf_max228 228 !! 229 229 CHARACTER(len=100) :: cn_dir ! Root directory for location of core files … … 235 235 TYPE(FLD_N) :: sn_div ! informations about the fields to be read 236 236 !! 237 NAMELIST/namdta_dyn/cn_dir, ln_dynrnf, ln_dynrnf_depth, fwbcorr,&237 NAMELIST/namdta_dyn/cn_dir, ln_dynrnf, ln_dynrnf_depth, & 238 238 & sn_uwd, sn_vwd, sn_wwd, sn_emp, & 239 239 & sn_avt, sn_tem, sn_sal, sn_mld , sn_qsr , & … … 257 257 WRITE(numout,*) ' runoffs option enabled (T) or not (F) ln_dynrnf = ', ln_dynrnf 258 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 = ', fwbcorr260 259 WRITE(numout,*) 261 260 ENDIF … … 353 352 DO jk = 1, jpkm1 354 353 e3t(:,:,jk,Kmm) = e3t_0(:,:,jk) * ( 1._wp + ssh(:,:,Kmm) * r1_ht_0(:,:) * tmask(:,:,jk) ) 354 e3t(:,:,jk,Kbb) = e3t_0(:,:,jk) * ( 1._wp + ssh(:,:,Kbb) * r1_ht_0(:,:) * tmask(:,:,jk) ) 355 355 ENDDO 356 e3t(:,:,jpk,Kaa) = e3t_0(:,:,jpk) 357 358 ! Horizontal scale factor interpolations 359 ! -------------------------------------- 360 CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3u(:,:,:,Kmm), 'U' ) 361 CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3v(:,:,:,Kmm), 'V' ) 362 363 ! Vertical scale factor interpolations 364 ! ------------------------------------ 365 CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3w(:,:,:,Kmm), 'W' ) 366 !!gm this should be computed from ssh(Kbb) 367 e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm) 368 e3u(:,:,:,Kbb) = e3u(:,:,:,Kmm) 369 e3v(:,:,:,Kbb) = e3v(:,:,:,Kmm) 370 371 ! t- and w- points depth 372 ! ---------------------- 373 gdept(:,:,1,Kmm) = 0.5_wp * e3w(:,:,1,Kmm) 374 gdepw(:,:,1,Kmm) = 0.0_wp 375 376 DO_3D( 1, 1, 1, 1, 2, jpk ) 377 ! zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) ! 0 everywhere 378 ! tmask = wmask, ie everywhere expect at jk = mikt 379 ! 1 for jk = 380 ! mikt 381 zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 382 gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm) 383 gdept(ji,jj,jk,Kmm) = zcoef * ( gdepw(ji,jj,jk ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm)) & 384 & + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) + e3w(ji,jj,jk,Kmm)) 385 END_3D 386 387 gdept(:,:,:,Kbb) = gdept(:,:,:,Kmm) 388 gdepw(:,:,:,Kbb) = gdepw(:,:,:,Kmm) 389 ! 390 ENDIF 356 357 CALL dta_dyn_sf_interp( nit000, Kmm ) 358 CALL dta_dyn_sf_interp( nit000, Kbb ) 391 359 #endif 392 ! 360 ENDIF 361 ! 362 CALL dta_dyn_rnf_init( Kmm ) 363 ! 364 CALL dta_dyn( nit000, Kbb, Kmm, Kaa ) 365 ! 366 END SUBROUTINE dta_dyn_init 367 368 SUBROUTINE dta_dyn_atf( kt, Kbb, Kmm, Kaa ) 369 !!--------------------------------------------------------------------- 370 !! *** ROUTINE dta_dyn_swp *** 371 !! 372 !! ** Purpose : Asselin time filter of now SSH 373 !!--------------------------------------------------------------------- 374 INTEGER, INTENT(in) :: kt ! time step 375 INTEGER, INTENT(in) :: Kbb, Kmm, Kaa ! ocean time level indices 376 ! 377 !!--------------------------------------------------------------------- 378 379 IF( kt == nit000 ) THEN 380 IF(lwp) WRITE(numout,*) 381 IF(lwp) WRITE(numout,*) 'dta_dyn_atf : Asselin time filter of sea surface height' 382 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ ' 383 ENDIF 384 385 ssh(:,:,Kmm) = ssh(:,:,Kmm) + rn_atfp * ( ssh(:,:,Kbb) - 2 * ssh(:,:,Kmm) + ssh(:,:,Kaa)) 386 387 ! 388 END SUBROUTINE dta_dyn_atf 389 390 391 #if ! defined key_qco 392 393 SUBROUTINE dta_dyn_sf_interp( kt, Kmm ) 394 !!--------------------------------------------------------------------- 395 !! *** ROUTINE dta_dyn_sf_interp *** 396 !! 397 !! ** Purpose : Calculate scale factors at U/V/W points and depths 398 !! given the after e3t field 399 !!--------------------------------------------------------------------- 400 INTEGER, INTENT(in) :: kt ! time step 401 INTEGER, INTENT(in) :: Kmm ! ocean time level indices 402 ! 403 INTEGER :: ji, jj, jk 404 REAL(wp) :: zcoef 405 !!--------------------------------------------------------------------- 406 407 ! Horizontal scale factor interpolations 408 ! -------------------------------------- 409 CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3u(:,:,:,Kmm), 'U' ) 410 CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3v(:,:,:,Kmm), 'V' ) 411 412 ! Vertical scale factor interpolations 413 ! ------------------------------------ 414 CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3w (:,:,:,Kmm), 'W' ) 415 416 ! t- and w- points depth 417 ! ---------------------- 418 gdept(:,:,1,Kmm) = 0.5_wp * e3w(:,:,1,Kmm) 419 gdepw(:,:,1,Kmm) = 0.0_wp 420 ! 421 DO_3D( 1, 1, 1, 1, 2, jpk ) 422 zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 423 gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm) 424 gdept(ji,jj,jk,Kmm) = zcoef * ( gdepw(ji,jj,jk ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm)) & 425 & + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) + e3w(ji,jj,jk,Kmm)) 426 END_3D 427 ! 428 END SUBROUTINE dta_dyn_sf_interp 429 430 #endif 431 432 SUBROUTINE dta_dyn_ssh( kt, phdivtr, psshb, pemp, pssha, pe3ta ) 433 !!---------------------------------------------------------------------- 434 !! *** ROUTINE dta_dyn_wzv *** 435 !! 436 !! ** Purpose : compute the after ssh (ssh(:,:,Kaa)) and the now vertical velocity 437 !! 438 !! ** Method : Using the incompressibility hypothesis, 439 !! - the ssh increment is computed by integrating the horizontal divergence 440 !! and multiply by the time step. 441 !! 442 !! - compute the after scale factor : repartition of ssh INCREMENT proportionnaly 443 !! to the level thickness ( z-star case ) 444 !! 445 !! - the vertical velocity is computed by integrating the horizontal divergence 446 !! from the bottom to the surface minus the scale factor evolution. 447 !! The boundary conditions are w=0 at the bottom (no flux) 448 !! 449 !! ** action : ssh(:,:,Kaa) / e3t(:,:,k,Kaa) / ww 450 !! 451 !! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling. 452 !!---------------------------------------------------------------------- 453 INTEGER, INTENT(in ) :: kt ! time-step 454 REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT(in ) :: phdivtr ! horizontal divergence transport 455 REAL(wp), DIMENSION(jpi,jpj) , OPTIONAL, INTENT(in ) :: psshb ! now ssh 456 REAL(wp), DIMENSION(jpi,jpj) , OPTIONAL, INTENT(in ) :: pemp ! evaporation minus precipitation 457 REAL(wp), DIMENSION(jpi,jpj) , OPTIONAL, INTENT(inout) :: pssha ! after ssh 458 REAL(wp), DIMENSION(jpi,jpj,jpk), OPTIONAL, INTENT(out) :: pe3ta ! after vertical scale factor 459 ! 460 INTEGER :: jk 461 REAL(wp), DIMENSION(jpi,jpj) :: zhdiv 462 REAL(wp) :: z2dt 463 !!---------------------------------------------------------------------- 464 ! 465 z2dt = 2._wp * rn_Dt 466 ! 467 zhdiv(:,:) = 0._wp 468 DO jk = 1, jpkm1 469 zhdiv(:,:) = zhdiv(:,:) + phdivtr(:,:,jk) * tmask(:,:,jk) 470 END DO 471 ! ! Sea surface elevation time-stepping 472 pssha(:,:) = ( psshb(:,:) - z2dt * ( r1_rho0 * pemp(:,:) + zhdiv(:,:) ) ) * ssmask(:,:) 473 ! 474 IF( PRESENT( pe3ta ) ) THEN ! After acale factors at t-points ( z_star coordinate ) 475 DO jk = 1, jpkm1 476 pe3ta(:,:,jk) = e3t_0(:,:,jk) * ( 1._wp + pssha(:,:) * r1_ht_0(:,:) * tmask(:,:,jk) ) 477 END DO 478 ENDIF 479 ! 480 END SUBROUTINE dta_dyn_ssh 481 482 SUBROUTINE dta_dyn_rnf_init( Kmm ) 483 !!---------------------------------------------------------------------- 484 !! *** ROUTINE dta_dyn_rnf_init *** 485 !! 486 !! ** Purpose : Initialisation of the runoffs if (ln_rnf=T) 487 !! 488 !!---------------------------------------------------------------------- 489 INTEGER, INTENT(in) :: Kmm ! ocean time level indices 490 ! 491 INTEGER :: inum ! local integer 492 INTEGER :: ji, jj, jk 493 REAL(wp) :: zcoef 494 INTEGER :: nkrnf_max 495 REAL(wp) :: hrnf_max 496 393 497 IF( ln_dynrnf .AND. ln_dynrnf_depth ) THEN ! read depht over which runoffs are distributed 394 498 IF(lwp) WRITE(numout,*) … … 412 516 ENDIF 413 517 END_2D 518 ! 414 519 DO_2D( 1, 1, 1, 1 ) ! set the associated depth 415 520 h_rnf(ji,jj) = 0._wp … … 432 537 IF(lwp) WRITE(numout,*) ' ' 433 538 ! 434 CALL dta_dyn( nit000, Kbb, Kmm, Kaa ) 435 ! 436 END SUBROUTINE dta_dyn_init 437 438 439 SUBROUTINE dta_dyn_sed( kt, Kmm ) 440 !!---------------------------------------------------------------------- 441 !! *** ROUTINE dta_dyn *** 442 !! 443 !! ** Purpose : Prepares dynamics and physics fields from a NEMO run 444 !! for an off-line simulation of passive tracers 445 !! 446 !! ** Method : calculates the position of data 447 !! - computes slopes if needed 448 !! - interpolates data if needed 449 !!---------------------------------------------------------------------- 450 INTEGER, INTENT(in) :: kt ! ocean time-step index 451 INTEGER, INTENT(in) :: Kmm ! ocean time level index 452 ! 453 !!---------------------------------------------------------------------- 454 ! 455 IF( ln_timing ) CALL timing_start( 'dta_dyn_sed') 456 ! 457 nsecdyn = nsec_year + nsec1jan000 ! number of seconds between Jan. 1st 00h of nit000 year and the middle of time step 458 ! 459 IF( kt == nit000 ) THEN ; nprevrec = 0 460 ELSE ; nprevrec = sf_dyn(jf_tem)%nrec(2,sf_dyn(jf_tem)%naa) 461 ENDIF 462 CALL fld_read( kt, 1, sf_dyn ) != read data at kt time step ==! 463 ! 464 ts(:,:,:,jp_tem,Kmm) = sf_dyn(jf_tem)%fnow(:,:,:) * tmask(:,:,:) ! temperature 465 ts(:,:,:,jp_sal,Kmm) = sf_dyn(jf_sal)%fnow(:,:,:) * tmask(:,:,:) ! salinity 466 ! 467 CALL eos ( ts(:,:,:,:,Kmm), rhd, rhop, gdept_0(:,:,:) ) ! In any case, we need rhop 468 469 IF(sn_cfctl%l_prtctl) THEN ! print control 470 CALL prt_ctl(tab3d_1=ts(:,:,:,jp_tem,Kmm), clinfo1=' tn - : ', mask1=tmask, kdim=jpk ) 471 CALL prt_ctl(tab3d_1=ts(:,:,:,jp_sal,Kmm), clinfo1=' sn - : ', mask1=tmask, kdim=jpk ) 472 ENDIF 473 ! 474 IF( ln_timing ) CALL timing_stop( 'dta_dyn_sed') 475 ! 476 END SUBROUTINE dta_dyn_sed 477 478 479 SUBROUTINE dta_dyn_sed_init( Kmm ) 480 !!---------------------------------------------------------------------- 481 !! *** ROUTINE dta_dyn_init *** 482 !! 483 !! ** Purpose : Initialisation of the dynamical data 484 !! ** Method : - read the data namdta_dyn namelist 485 !!---------------------------------------------------------------------- 486 INTEGER, INTENT( in ) :: Kmm ! ocean time level index 487 ! 488 INTEGER :: ierr, ierr0, ierr1, ierr2, ierr3 ! return error code 489 INTEGER :: ifpr ! dummy loop indice 490 INTEGER :: jfld ! dummy loop arguments 491 INTEGER :: inum, idv, idimv ! local integer 492 INTEGER :: ios ! Local integer output status for namelist read 493 !! 494 CHARACTER(len=100) :: cn_dir ! Root directory for location of core files 495 TYPE(FLD_N), DIMENSION(2) :: slf_d ! array of namelist informations on the fields to read 496 TYPE(FLD_N) :: sn_tem , sn_sal ! " " 497 !! 498 NAMELIST/namdta_dyn/cn_dir, ln_dynrnf, ln_dynrnf_depth, fwbcorr, & 499 & sn_tem, sn_sal 500 !!---------------------------------------------------------------------- 501 ! 502 READ ( numnam_ref, namdta_dyn, IOSTAT = ios, ERR = 901) 503 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdta_dyn in reference namelist' ) 504 READ ( numnam_cfg, namdta_dyn, IOSTAT = ios, ERR = 902 ) 505 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namdta_dyn in configuration namelist' ) 506 IF(lwm) WRITE ( numond, namdta_dyn ) 507 ! ! store namelist information in an array 508 ! ! Control print 509 IF(lwp) THEN 510 WRITE(numout,*) 511 WRITE(numout,*) 'dta_dyn : offline dynamics ' 512 WRITE(numout,*) '~~~~~~~ ' 513 WRITE(numout,*) ' Namelist namdta_dyn' 514 WRITE(numout,*) ' runoffs option enabled (T) or not (F) ln_dynrnf = ', ln_dynrnf 515 WRITE(numout,*) ' runoffs is spread in vertical ln_dynrnf_depth = ', ln_dynrnf_depth 516 WRITE(numout,*) ' annual global mean of empmr for ssh correction fwbcorr = ', fwbcorr 517 WRITE(numout,*) 518 ENDIF 519 ! 520 jf_tem = 1 ; jf_sal = 2 ; jfld = jf_sal 521 ! 522 slf_d(jf_tem) = sn_tem ; slf_d(jf_sal) = sn_sal 523 ! 524 ALLOCATE( sf_dyn(jfld), STAT=ierr ) ! set sf structure 525 IF( ierr > 0 ) THEN 526 CALL ctl_stop( 'dta_dyn: unable to allocate sf structure' ) ; RETURN 527 ENDIF 528 ! ! fill sf with slf_i and control print 529 CALL fld_fill( sf_dyn, slf_d, cn_dir, 'dta_dyn_init', 'Data in file', 'namdta_dyn' ) 530 ! 531 ! Open file for each variable to get his number of dimension 532 DO ifpr = 1, jfld 533 CALL fld_def( sf_dyn(ifpr) ) 534 CALL iom_open( sf_dyn(ifpr)%clname, sf_dyn(ifpr)%num ) 535 idv = iom_varid( sf_dyn(ifpr)%num , slf_d(ifpr)%clvar ) ! id of the variable sdjf%clvar 536 idimv = iom_file ( sf_dyn(ifpr)%num )%ndims(idv) ! number of dimension for variable sdjf%clvar 537 CALL iom_close( sf_dyn(ifpr)%num ) ! close file if already open 538 ierr1=0 539 IF( idimv == 3 ) THEN ! 2D variable 540 ALLOCATE( sf_dyn(ifpr)%fnow(jpi,jpj,1) , STAT=ierr0 ) 541 IF( slf_d(ifpr)%ln_tint ) ALLOCATE( sf_dyn(ifpr)%fdta(jpi,jpj,1,2) , STAT=ierr1 ) 542 ELSE ! 3D variable 543 ALLOCATE( sf_dyn(ifpr)%fnow(jpi,jpj,jpk) , STAT=ierr0 ) 544 IF( slf_d(ifpr)%ln_tint ) ALLOCATE( sf_dyn(ifpr)%fdta(jpi,jpj,jpk,2), STAT=ierr1 ) 545 ENDIF 546 IF( ierr0 + ierr1 > 0 ) THEN 547 CALL ctl_stop( 'dta_dyn_init : unable to allocate sf_dyn array structure' ) ; RETURN 548 ENDIF 549 END DO 550 ! 551 CALL dta_dyn_sed( nit000, Kmm ) 552 ! 553 END SUBROUTINE dta_dyn_sed_init 554 555 556 SUBROUTINE dta_dyn_atf( kt, Kbb, Kmm, Kaa ) 557 !!--------------------------------------------------------------------- 558 !! *** ROUTINE dta_dyn_swp *** 559 !! 560 !! ** Purpose : Asselin time filter of now SSH 561 !!--------------------------------------------------------------------- 562 INTEGER, INTENT(in) :: kt ! time step 563 INTEGER, INTENT(in) :: Kbb, Kmm, Kaa ! ocean time level indices 564 ! 565 !!--------------------------------------------------------------------- 566 567 IF( kt == nit000 ) THEN 568 IF(lwp) WRITE(numout,*) 569 IF(lwp) WRITE(numout,*) 'dta_dyn_atf : Asselin time filter of sea surface height' 570 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ ' 571 ENDIF 572 573 ssh(:,:,Kmm) = ssh(:,:,Kmm) + rn_atfp * ( ssh(:,:,Kbb) - 2 * ssh(:,:,Kmm) + ssh(:,:,Kaa)) 574 575 !! Do we also need to time filter e3t?? 576 ! 577 END SUBROUTINE dta_dyn_atf 578 579 580 #if ! defined key_qco 581 SUBROUTINE dta_dyn_sf_interp( kt, Kmm ) 582 !!--------------------------------------------------------------------- 583 !! *** ROUTINE dta_dyn_sf_interp *** 584 !! 585 !! ** Purpose : Calculate scale factors at U/V/W points and depths 586 !! given the after e3t field 587 !!--------------------------------------------------------------------- 588 INTEGER, INTENT(in) :: kt ! time step 589 INTEGER, INTENT(in) :: Kmm ! ocean time level indices 590 ! 591 INTEGER :: ji, jj, jk 592 REAL(wp) :: zcoef 593 !!--------------------------------------------------------------------- 594 595 ! Horizontal scale factor interpolations 596 ! -------------------------------------- 597 CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3u(:,:,:,Kmm), 'U' ) 598 CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3v(:,:,:,Kmm), 'V' ) 599 600 ! Vertical scale factor interpolations 601 ! ------------------------------------ 602 CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3w (:,:,:,Kmm), 'W' ) 603 604 ! t- and w- points depth 605 ! ---------------------- 606 gdept(:,:,1,Kmm) = 0.5_wp * e3w(:,:,1,Kmm) 607 gdepw(:,:,1,Kmm) = 0.0_wp 608 ! 609 DO_3D( 1, 1, 1, 1, 2, jpk ) 610 zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 611 gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm) 612 gdept(ji,jj,jk,Kmm) = zcoef * ( gdepw(ji,jj,jk ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm)) & 613 & + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) + e3w(ji,jj,jk,Kmm)) 614 END_3D 615 ! 616 END SUBROUTINE dta_dyn_sf_interp 617 #endif 618 619 620 SUBROUTINE dta_dyn_ssh( kt, phdivtr, psshb, pemp, pssha, pe3ta ) 621 !!---------------------------------------------------------------------- 622 !! *** ROUTINE dta_dyn_wzv *** 623 !! 624 !! ** Purpose : compute the after ssh (ssh(:,:,Kaa)) and the now vertical velocity 625 !! 626 !! ** Method : Using the incompressibility hypothesis, 627 !! - the ssh increment is computed by integrating the horizontal divergence 628 !! and multiply by the time step. 629 !! 630 !! - compute the after scale factor : repartition of ssh INCREMENT proportionnaly 631 !! to the level thickness ( z-star case ) 632 !! 633 !! - the vertical velocity is computed by integrating the horizontal divergence 634 !! from the bottom to the surface minus the scale factor evolution. 635 !! The boundary conditions are w=0 at the bottom (no flux) 636 !! 637 !! ** action : ssh(:,:,Kaa) / e3t(:,:,k,Kaa) / ww 638 !! 639 !! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling. 640 !!---------------------------------------------------------------------- 641 INTEGER, INTENT(in ) :: kt ! time-step 642 REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT(in ) :: phdivtr ! horizontal divergence transport 643 REAL(wp), DIMENSION(jpi,jpj) , OPTIONAL, INTENT(in ) :: psshb ! now ssh 644 REAL(wp), DIMENSION(jpi,jpj) , OPTIONAL, INTENT(in ) :: pemp ! evaporation minus precipitation 645 REAL(wp), DIMENSION(jpi,jpj) , OPTIONAL, INTENT(inout) :: pssha ! after ssh 646 REAL(wp), DIMENSION(jpi,jpj,jpk), OPTIONAL, INTENT(out) :: pe3ta ! after vertical scale factor 647 ! 648 INTEGER :: jk 649 REAL(wp), DIMENSION(jpi,jpj) :: zhdiv 650 REAL(wp) :: z2dt 651 !!---------------------------------------------------------------------- 652 ! 653 z2dt = 2._wp * rn_Dt 654 ! 655 zhdiv(:,:) = 0._wp 656 DO jk = 1, jpkm1 657 zhdiv(:,:) = zhdiv(:,:) + phdivtr(:,:,jk) * tmask(:,:,jk) 658 END DO 659 ! ! Sea surface elevation time-stepping 660 pssha(:,:) = ( psshb(:,:) - z2dt * ( r1_rho0 * pemp(:,:) + zhdiv(:,:) ) ) * ssmask(:,:) 661 ! 662 IF( PRESENT( pe3ta ) ) THEN ! After acale factors at t-points ( z_star coordinate ) 663 DO jk = 1, jpkm1 664 pe3ta(:,:,jk) = e3t_0(:,:,jk) * ( 1._wp + pssha(:,:) * r1_ht_0(:,:) * tmask(:,:,jk) ) 665 END DO 666 ENDIF 667 ! 668 END SUBROUTINE dta_dyn_ssh 669 670 671 SUBROUTINE dta_dyn_hrnf( Kmm ) 672 !!---------------------------------------------------------------------- 673 !! *** ROUTINE sbc_rnf *** 539 END SUBROUTINE dta_dyn_rnf_init 540 541 SUBROUTINE dta_dyn_rnf( Kmm ) 542 !!---------------------------------------------------------------------- 543 !! *** ROUTINE dta_dyn_rnf *** 674 544 !! 675 545 !! ** Purpose : update the horizontal divergence with the runoff inflow 676 546 !! 677 !! ** Method :678 !! CAUTION : rnf is positive (inflow) decreasing the679 !! divergence and expressed in m/s680 !!681 !! ** Action : phdivn decreased by the runoff inflow682 547 !!---------------------------------------------------------------------- 683 548 !! … … 694 559 END_2D 695 560 ! 696 END SUBROUTINE dta_dyn_hrnf 697 698 561 END SUBROUTINE dta_dyn_rnf 699 562 700 563 SUBROUTINE dta_dyn_slp( kt, Kbb, Kmm ) … … 787 650 END SUBROUTINE dta_dyn_slp 788 651 789 790 652 SUBROUTINE compute_slopes( kt, pts, puslp, pvslp, pwslpi, pwslpj, Kbb, Kmm ) 791 653 !!--------------------------------------------------------------------- … … 795 657 !!--------------------------------------------------------------------- 796 658 INTEGER , INTENT(in ) :: kt ! time step 797 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in 659 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(inout) :: pts ! temperature/salinity 798 660 REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT(out) :: puslp ! zonal isopycnal slopes 799 661 REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT(out) :: pvslp ! meridional isopycnal slopes … … 832 694 END SUBROUTINE compute_slopes 833 695 696 #if defined key_sed_off 697 698 SUBROUTINE dta_dyn_sed( kt, Kmm ) 699 !!---------------------------------------------------------------------- 700 !! *** ROUTINE dta_dyn *** 701 !! 702 !! ** Purpose : Prepares dynamics and physics fields from a NEMO run 703 !! for an off-line simulation of passive tracers 704 !! 705 !! ** Method : calculates the position of data 706 !! - computes slopes if needed 707 !! - interpolates data if needed 708 !!---------------------------------------------------------------------- 709 INTEGER, INTENT(in) :: kt ! ocean time-step index 710 INTEGER, INTENT(in) :: Kmm ! ocean time level index 711 ! 712 !!---------------------------------------------------------------------- 713 ! 714 IF( ln_timing ) CALL timing_start( 'dta_dyn_sed') 715 ! 716 nsecdyn = nsec_year + nsec1jan000 ! number of seconds between Jan. 1st 00h of nit000 year and the middle of time step 717 ! 718 IF( kt == nit000 ) THEN ; nprevrec = 0 719 ELSE ; nprevrec = sf_dyn(jf_tem)%nrec(2,sf_dyn(jf_tem)%naa) 720 ENDIF 721 CALL fld_read( kt, 1, sf_dyn ) != read data at kt time step ==! 722 ! 723 ts(:,:,:,jp_tem,Kmm) = sf_dyn(jf_tem)%fnow(:,:,:) * tmask(:,:,:) ! temperature 724 ts(:,:,:,jp_sal,Kmm) = sf_dyn(jf_sal)%fnow(:,:,:) * tmask(:,:,:) ! salinity 725 ! 726 CALL eos ( ts(:,:,:,:,Kmm), rhd, rhop, gdept_0(:,:,:) ) ! In any case, we need rhop 727 728 IF(sn_cfctl%l_prtctl) THEN ! print control 729 CALL prt_ctl(tab3d_1=ts(:,:,:,jp_tem,Kmm), clinfo1=' tn - : ', mask1=tmask, kdim=jpk ) 730 CALL prt_ctl(tab3d_1=ts(:,:,:,jp_sal,Kmm), clinfo1=' sn - : ', mask1=tmask, kdim=jpk ) 731 ENDIF 732 ! 733 IF( ln_timing ) CALL timing_stop( 'dta_dyn_sed') 734 ! 735 END SUBROUTINE dta_dyn_sed 736 737 738 SUBROUTINE dta_dyn_sed_init( Kmm ) 739 !!---------------------------------------------------------------------- 740 !! *** ROUTINE dta_dyn_init *** 741 !! 742 !! ** Purpose : Initialisation of the dynamical data 743 !! ** Method : - read the data namdta_dyn namelist 744 !!---------------------------------------------------------------------- 745 INTEGER, INTENT( in ) :: Kmm ! ocean time level index 746 ! 747 INTEGER :: ierr, ierr0, ierr1, ierr2, ierr3 ! return error code 748 INTEGER :: ifpr ! dummy loop indice 749 INTEGER :: jfld ! dummy loop arguments 750 INTEGER :: inum, idv, idimv ! local integer 751 INTEGER :: ios ! Local integer output status for namelist read 752 !! 753 CHARACTER(len=100) :: cn_dir ! Root directory for location of core files 754 TYPE(FLD_N), DIMENSION(2) :: slf_d ! array of namelist informations on the fields to read 755 TYPE(FLD_N) :: sn_tem , sn_sal ! " " 756 !! 757 NAMELIST/namdta_dyn/cn_dir, ln_dynrnf, ln_dynrnf_depth, & 758 & sn_tem, sn_sal 759 !!---------------------------------------------------------------------- 760 ! 761 READ ( numnam_ref, namdta_dyn, IOSTAT = ios, ERR = 901) 762 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdta_dyn in reference namelist' ) 763 READ ( numnam_cfg, namdta_dyn, IOSTAT = ios, ERR = 902 ) 764 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namdta_dyn in configuration namelist' ) 765 IF(lwm) WRITE ( numond, namdta_dyn ) 766 ! ! store namelist information in an array 767 ! ! Control print 768 IF(lwp) THEN 769 WRITE(numout,*) 770 WRITE(numout,*) 'dta_dyn : offline dynamics ' 771 WRITE(numout,*) '~~~~~~~ ' 772 WRITE(numout,*) ' Namelist namdta_dyn' 773 WRITE(numout,*) ' runoffs option enabled (T) or not (F) ln_dynrnf = ', ln_dynrnf 774 WRITE(numout,*) ' runoffs is spread in vertical ln_dynrnf_depth = ', ln_dynrnf_depth 775 WRITE(numout,*) 776 ENDIF 777 ! 778 jf_tem = 1 ; jf_sal = 2 ; jfld = jf_sal 779 ! 780 slf_d(jf_tem) = sn_tem ; slf_d(jf_sal) = sn_sal 781 ! 782 ALLOCATE( sf_dyn(jfld), STAT=ierr ) ! set sf structure 783 IF( ierr > 0 ) THEN 784 CALL ctl_stop( 'dta_dyn: unable to allocate sf structure' ) ; RETURN 785 ENDIF 786 ! ! fill sf with slf_i and control print 787 CALL fld_fill( sf_dyn, slf_d, cn_dir, 'dta_dyn_init', 'Data in file', 'namdta_dyn' ) 788 ! 789 ! Open file for each variable to get his number of dimension 790 DO ifpr = 1, jfld 791 CALL fld_def( sf_dyn(ifpr) ) 792 CALL iom_open( sf_dyn(ifpr)%clname, sf_dyn(ifpr)%num ) 793 idv = iom_varid( sf_dyn(ifpr)%num , slf_d(ifpr)%clvar ) ! id of the variable sdjf%clvar 794 idimv = iom_file ( sf_dyn(ifpr)%num )%ndims(idv) ! number of dimension for variable sdjf%clvar 795 CALL iom_close( sf_dyn(ifpr)%num ) ! close file if already open 796 ierr1=0 797 IF( idimv == 3 ) THEN ! 2D variable 798 ALLOCATE( sf_dyn(ifpr)%fnow(jpi,jpj,1) , STAT=ierr0 ) 799 IF( slf_d(ifpr)%ln_tint ) ALLOCATE( sf_dyn(ifpr)%fdta(jpi,jpj,1,2) , STAT=ierr1 ) 800 ELSE ! 3D variable 801 ALLOCATE( sf_dyn(ifpr)%fnow(jpi,jpj,jpk) , STAT=ierr0 ) 802 IF( slf_d(ifpr)%ln_tint ) ALLOCATE( sf_dyn(ifpr)%fdta(jpi,jpj,jpk,2), STAT=ierr1 ) 803 ENDIF 804 IF( ierr0 + ierr1 > 0 ) THEN 805 CALL ctl_stop( 'dta_dyn_init : unable to allocate sf_dyn array structure' ) ; RETURN 806 ENDIF 807 END DO 808 ! 809 CALL dta_dyn_sed( nit000, Kmm ) 810 ! 811 END SUBROUTINE dta_dyn_sed_init 812 #endif 834 813 !!====================================================================== 835 814 END MODULE dtadyn
Note: See TracChangeset
for help on using the changeset viewer.
