Changeset 14255 for NEMO/trunk/src
- Timestamp:
- 2021-01-04T12:35:00+01:00 (3 years ago)
- Location:
- NEMO/trunk/src
- Files:
-
- 4 edited
-
OCE/DOM/domain.F90 (modified) (1 diff)
-
OFF/dtadyn.F90 (modified) (11 diffs)
-
OFF/nemogcm.F90 (modified) (2 diffs)
-
TOP/trcwri.F90 (modified) (3 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk/src/OCE/DOM/domain.F90
r14239 r14255 182 182 ! != ssh initialization 183 183 ! 184 IF( l_ offline .OR. l_SAS ) THEN !* No ocean dynamics calculation : set to 0184 IF( l_SAS ) THEN !* No ocean dynamics calculation : set to 0 185 185 ssh(:,:,:) = 0._wp 186 186 #if defined key_agrif -
NEMO/trunk/src/OFF/dtadyn.F90
r14053 r14255 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 … … 138 139 ts(:,:,:,jp_tem,Kmm) = sf_dyn(jf_tem)%fnow(:,:,:) * tmask(:,:,:) ! temperature 139 140 ts(:,:,:,jp_sal,Kmm) = sf_dyn(jf_sal)%fnow(:,:,:) * tmask(:,:,:) ! salinity 140 wndm(:,:) = sf_dyn(jf_wnd)%fnow(:,:,1) * tmask(:,:,1) ! wind speed - needed for gas exchange141 fmmflx(:,:) = sf_dyn(jf_fmf)%fnow(:,:,1) * tmask(:,:,1) ! downward salt flux (v3.5+)142 fr_i(:,:) = sf_dyn(jf_ice)%fnow(:,:,1) * tmask(:,:,1) ! Sea-ice fraction143 qsr (:,:) = sf_dyn(jf_qsr)%fnow(:,:,1) * tmask(:,:,1) ! solar radiation144 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 145 146 IF( ln_dynrnf ) THEN 146 rnf (:,:) = sf_dyn(jf_rnf)%fnow(:,:,1) * tmask(:,:,1) ! E-P147 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 ) 148 149 ENDIF 149 150 ! 150 151 uu(:,:,:,Kmm) = sf_dyn(jf_uwd)%fnow(:,:,:) * umask(:,:,:) ! effective u-transport 151 152 vv(:,:,:,Kmm) = sf_dyn(jf_vwd)%fnow(:,:,:) * vmask(:,:,:) ! effective v-transport 152 ww(:,:,:) = sf_dyn(jf_wwd)%fnow(:,:,:) * tmask(:,:,:) ! effective v-transport153 ww(:,:,:) = sf_dyn(jf_wwd)%fnow(:,:,:) * tmask(:,:,:) ! effective v-transport 153 154 ! 154 155 IF( .NOT.ln_linssh ) THEN … … 156 157 zhdivtr(:,:,:) = sf_dyn(jf_div)%fnow(:,:,:) * tmask(:,:,:) ! effective u-transport 157 158 emp_b (:,:) = sf_dyn(jf_empb)%fnow(:,:,1) * tmask(:,:,1) ! E-P 158 zemp (:,:) = ( 0.5_wp * ( emp(:,:) + emp_b(:,:) ) + rnf(:,:) + fwbcorr) * tmask(:,:,1)159 zemp (:,:) = ( 0.5_wp * ( emp(:,:) + emp_b(:,:) ) + rnf(:,:) ) * tmask(:,:,1) 159 160 #if defined key_qco 160 161 CALL dta_dyn_ssh( kt, zhdivtr, ssh(:,:,Kbb), zemp, ssh(:,:,Kaa) ) … … 225 226 INTEGER :: ios ! Local integer output status for namelist read 226 227 INTEGER :: ji, jj, jk 227 REAL(wp) :: zcoef228 INTEGER :: nkrnf_max229 REAL(wp) :: hrnf_max230 228 !! 231 229 CHARACTER(len=100) :: cn_dir ! Root directory for location of core files … … 237 235 TYPE(FLD_N) :: sn_div ! informations about the fields to be read 238 236 !! 239 NAMELIST/namdta_dyn/cn_dir, ln_dynrnf, ln_dynrnf_depth, fwbcorr,&237 NAMELIST/namdta_dyn/cn_dir, ln_dynrnf, ln_dynrnf_depth, & 240 238 & sn_uwd, sn_vwd, sn_wwd, sn_emp, & 241 239 & sn_avt, sn_tem, sn_sal, sn_mld , sn_qsr , & … … 259 257 WRITE(numout,*) ' runoffs option enabled (T) or not (F) ln_dynrnf = ', ln_dynrnf 260 258 WRITE(numout,*) ' runoffs is spread in vertical ln_dynrnf_depth = ', ln_dynrnf_depth 261 WRITE(numout,*) ' annual global mean of empmr for ssh correction fwbcorr = ', fwbcorr262 259 WRITE(numout,*) 263 260 ENDIF … … 355 352 DO jk = 1, jpkm1 356 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) ) 357 355 ENDDO 358 e3t(:,:,jpk,Kaa) = e3t_0(:,:,jpk) 359 360 ! Horizontal scale factor interpolations 361 ! -------------------------------------- 362 CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3u(:,:,:,Kmm), 'U' ) 363 CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3v(:,:,:,Kmm), 'V' ) 364 365 ! Vertical scale factor interpolations 366 ! ------------------------------------ 367 CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3w(:,:,:,Kmm), 'W' ) 368 !!gm this should be computed from ssh(Kbb) 369 e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm) 370 e3u(:,:,:,Kbb) = e3u(:,:,:,Kmm) 371 e3v(:,:,:,Kbb) = e3v(:,:,:,Kmm) 372 373 ! t- and w- points depth 374 ! ---------------------- 375 gdept(:,:,1,Kmm) = 0.5_wp * e3w(:,:,1,Kmm) 376 gdepw(:,:,1,Kmm) = 0.0_wp 377 378 DO_3D( 1, 1, 1, 1, 2, jpk ) 379 ! zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) ! 0 everywhere 380 ! tmask = wmask, ie everywhere expect at jk = mikt 381 ! 1 for jk = 382 ! mikt 383 zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 384 gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm) 385 gdept(ji,jj,jk,Kmm) = zcoef * ( gdepw(ji,jj,jk ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm)) & 386 & + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) + e3w(ji,jj,jk,Kmm)) 387 END_3D 388 389 gdept(:,:,:,Kbb) = gdept(:,:,:,Kmm) 390 gdepw(:,:,:,Kbb) = gdepw(:,:,:,Kmm) 391 ! 356 357 CALL dta_dyn_sf_interp( kt, Kmm ) 358 CALL dta_dyn_sf_interp( kt, Kbb ) 392 359 #endif 393 360 ENDIF 394 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 395 497 IF( ln_dynrnf .AND. ln_dynrnf_depth ) THEN ! read depht over which runoffs are distributed 396 498 IF(lwp) WRITE(numout,*) … … 435 537 IF(lwp) WRITE(numout,*) ' ' 436 538 ! 437 CALL dta_dyn( nit000, Kbb, Kmm, Kaa ) 438 ! 439 END SUBROUTINE dta_dyn_init 440 441 442 SUBROUTINE dta_dyn_sed( kt, Kmm ) 443 !!---------------------------------------------------------------------- 444 !! *** ROUTINE dta_dyn *** 445 !! 446 !! ** Purpose : Prepares dynamics and physics fields from a NEMO run 447 !! for an off-line simulation of passive tracers 448 !! 449 !! ** Method : calculates the position of data 450 !! - computes slopes if needed 451 !! - interpolates data if needed 452 !!---------------------------------------------------------------------- 453 INTEGER, INTENT(in) :: kt ! ocean time-step index 454 INTEGER, INTENT(in) :: Kmm ! ocean time level index 455 ! 456 !!---------------------------------------------------------------------- 457 ! 458 IF( ln_timing ) CALL timing_start( 'dta_dyn_sed') 459 ! 460 nsecdyn = nsec_year + nsec1jan000 ! number of seconds between Jan. 1st 00h of nit000 year and the middle of time step 461 ! 462 IF( kt == nit000 ) THEN ; nprevrec = 0 463 ELSE ; nprevrec = sf_dyn(jf_tem)%nrec(2,sf_dyn(jf_tem)%naa) 464 ENDIF 465 CALL fld_read( kt, 1, sf_dyn ) != read data at kt time step ==! 466 ! 467 ts(:,:,:,jp_tem,Kmm) = sf_dyn(jf_tem)%fnow(:,:,:) * tmask(:,:,:) ! temperature 468 ts(:,:,:,jp_sal,Kmm) = sf_dyn(jf_sal)%fnow(:,:,:) * tmask(:,:,:) ! salinity 469 ! 470 CALL eos ( ts(:,:,:,:,Kmm), rhd, rhop, gdept_0(:,:,:) ) ! In any case, we need rhop 471 472 IF(sn_cfctl%l_prtctl) THEN ! print control 473 CALL prt_ctl(tab3d_1=ts(:,:,:,jp_tem,Kmm), clinfo1=' tn - : ', mask1=tmask, kdim=jpk ) 474 CALL prt_ctl(tab3d_1=ts(:,:,:,jp_sal,Kmm), clinfo1=' sn - : ', mask1=tmask, kdim=jpk ) 475 ENDIF 476 ! 477 IF( ln_timing ) CALL timing_stop( 'dta_dyn_sed') 478 ! 479 END SUBROUTINE dta_dyn_sed 480 481 482 SUBROUTINE dta_dyn_sed_init( Kmm ) 483 !!---------------------------------------------------------------------- 484 !! *** ROUTINE dta_dyn_init *** 485 !! 486 !! ** Purpose : Initialisation of the dynamical data 487 !! ** Method : - read the data namdta_dyn namelist 488 !!---------------------------------------------------------------------- 489 INTEGER, INTENT( in ) :: Kmm ! ocean time level index 490 ! 491 INTEGER :: ierr, ierr0, ierr1, ierr2, ierr3 ! return error code 492 INTEGER :: ifpr ! dummy loop indice 493 INTEGER :: jfld ! dummy loop arguments 494 INTEGER :: inum, idv, idimv ! local integer 495 INTEGER :: ios ! Local integer output status for namelist read 496 !! 497 CHARACTER(len=100) :: cn_dir ! Root directory for location of core files 498 TYPE(FLD_N), DIMENSION(2) :: slf_d ! array of namelist informations on the fields to read 499 TYPE(FLD_N) :: sn_tem , sn_sal ! " " 500 !! 501 NAMELIST/namdta_dyn/cn_dir, ln_dynrnf, ln_dynrnf_depth, fwbcorr, & 502 & sn_tem, sn_sal 503 !!---------------------------------------------------------------------- 504 ! 505 READ ( numnam_ref, namdta_dyn, IOSTAT = ios, ERR = 901) 506 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdta_dyn in reference namelist' ) 507 READ ( numnam_cfg, namdta_dyn, IOSTAT = ios, ERR = 902 ) 508 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namdta_dyn in configuration namelist' ) 509 IF(lwm) WRITE ( numond, namdta_dyn ) 510 ! ! store namelist information in an array 511 ! ! Control print 512 IF(lwp) THEN 513 WRITE(numout,*) 514 WRITE(numout,*) 'dta_dyn : offline dynamics ' 515 WRITE(numout,*) '~~~~~~~ ' 516 WRITE(numout,*) ' Namelist namdta_dyn' 517 WRITE(numout,*) ' runoffs option enabled (T) or not (F) ln_dynrnf = ', ln_dynrnf 518 WRITE(numout,*) ' runoffs is spread in vertical ln_dynrnf_depth = ', ln_dynrnf_depth 519 WRITE(numout,*) ' annual global mean of empmr for ssh correction fwbcorr = ', fwbcorr 520 WRITE(numout,*) 521 ENDIF 522 ! 523 jf_tem = 1 ; jf_sal = 2 ; jfld = jf_sal 524 ! 525 slf_d(jf_tem) = sn_tem ; slf_d(jf_sal) = sn_sal 526 ! 527 ALLOCATE( sf_dyn(jfld), STAT=ierr ) ! set sf structure 528 IF( ierr > 0 ) THEN 529 CALL ctl_stop( 'dta_dyn: unable to allocate sf structure' ) ; RETURN 530 ENDIF 531 ! ! fill sf with slf_i and control print 532 CALL fld_fill( sf_dyn, slf_d, cn_dir, 'dta_dyn_init', 'Data in file', 'namdta_dyn' ) 533 ! 534 ! Open file for each variable to get his number of dimension 535 DO ifpr = 1, jfld 536 CALL fld_def( sf_dyn(ifpr) ) 537 CALL iom_open( sf_dyn(ifpr)%clname, sf_dyn(ifpr)%num ) 538 idv = iom_varid( sf_dyn(ifpr)%num , slf_d(ifpr)%clvar ) ! id of the variable sdjf%clvar 539 idimv = iom_file ( sf_dyn(ifpr)%num )%ndims(idv) ! number of dimension for variable sdjf%clvar 540 CALL iom_close( sf_dyn(ifpr)%num ) ! close file if already open 541 ierr1=0 542 IF( idimv == 3 ) THEN ! 2D variable 543 ALLOCATE( sf_dyn(ifpr)%fnow(jpi,jpj,1) , STAT=ierr0 ) 544 IF( slf_d(ifpr)%ln_tint ) ALLOCATE( sf_dyn(ifpr)%fdta(jpi,jpj,1,2) , STAT=ierr1 ) 545 ELSE ! 3D variable 546 ALLOCATE( sf_dyn(ifpr)%fnow(jpi,jpj,jpk) , STAT=ierr0 ) 547 IF( slf_d(ifpr)%ln_tint ) ALLOCATE( sf_dyn(ifpr)%fdta(jpi,jpj,jpk,2), STAT=ierr1 ) 548 ENDIF 549 IF( ierr0 + ierr1 > 0 ) THEN 550 CALL ctl_stop( 'dta_dyn_init : unable to allocate sf_dyn array structure' ) ; RETURN 551 ENDIF 552 END DO 553 ! 554 CALL dta_dyn_sed( nit000, Kmm ) 555 ! 556 END SUBROUTINE dta_dyn_sed_init 557 558 559 SUBROUTINE dta_dyn_atf( kt, Kbb, Kmm, Kaa ) 560 !!--------------------------------------------------------------------- 561 !! *** ROUTINE dta_dyn_swp *** 562 !! 563 !! ** Purpose : Asselin time filter of now SSH 564 !!--------------------------------------------------------------------- 565 INTEGER, INTENT(in) :: kt ! time step 566 INTEGER, INTENT(in) :: Kbb, Kmm, Kaa ! ocean time level indices 567 ! 568 !!--------------------------------------------------------------------- 569 570 IF( kt == nit000 ) THEN 571 IF(lwp) WRITE(numout,*) 572 IF(lwp) WRITE(numout,*) 'dta_dyn_atf : Asselin time filter of sea surface height' 573 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ ' 574 ENDIF 575 576 ssh(:,:,Kmm) = ssh(:,:,Kmm) + rn_atfp * ( ssh(:,:,Kbb) - 2 * ssh(:,:,Kmm) + ssh(:,:,Kaa)) 577 578 !! Do we also need to time filter e3t?? 579 ! 580 END SUBROUTINE dta_dyn_atf 581 582 583 #if ! defined key_qco 584 SUBROUTINE dta_dyn_sf_interp( kt, Kmm ) 585 !!--------------------------------------------------------------------- 586 !! *** ROUTINE dta_dyn_sf_interp *** 587 !! 588 !! ** Purpose : Calculate scale factors at U/V/W points and depths 589 !! given the after e3t field 590 !!--------------------------------------------------------------------- 591 INTEGER, INTENT(in) :: kt ! time step 592 INTEGER, INTENT(in) :: Kmm ! ocean time level indices 593 ! 594 INTEGER :: ji, jj, jk 595 REAL(wp) :: zcoef 596 !!--------------------------------------------------------------------- 597 598 ! Horizontal scale factor interpolations 599 ! -------------------------------------- 600 CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3u(:,:,:,Kmm), 'U' ) 601 CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3v(:,:,:,Kmm), 'V' ) 602 603 ! Vertical scale factor interpolations 604 ! ------------------------------------ 605 CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3w (:,:,:,Kmm), 'W' ) 606 607 ! t- and w- points depth 608 ! ---------------------- 609 gdept(:,:,1,Kmm) = 0.5_wp * e3w(:,:,1,Kmm) 610 gdepw(:,:,1,Kmm) = 0.0_wp 611 ! 612 DO_3D( 1, 1, 1, 1, 2, jpk ) 613 zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 614 gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm) 615 gdept(ji,jj,jk,Kmm) = zcoef * ( gdepw(ji,jj,jk ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm)) & 616 & + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) + e3w(ji,jj,jk,Kmm)) 617 END_3D 618 ! 619 END SUBROUTINE dta_dyn_sf_interp 620 #endif 621 622 623 SUBROUTINE dta_dyn_ssh( kt, phdivtr, psshb, pemp, pssha, pe3ta ) 624 !!---------------------------------------------------------------------- 625 !! *** ROUTINE dta_dyn_wzv *** 626 !! 627 !! ** Purpose : compute the after ssh (ssh(:,:,Kaa)) and the now vertical velocity 628 !! 629 !! ** Method : Using the incompressibility hypothesis, 630 !! - the ssh increment is computed by integrating the horizontal divergence 631 !! and multiply by the time step. 632 !! 633 !! - compute the after scale factor : repartition of ssh INCREMENT proportionnaly 634 !! to the level thickness ( z-star case ) 635 !! 636 !! - the vertical velocity is computed by integrating the horizontal divergence 637 !! from the bottom to the surface minus the scale factor evolution. 638 !! The boundary conditions are w=0 at the bottom (no flux) 639 !! 640 !! ** action : ssh(:,:,Kaa) / e3t(:,:,k,Kaa) / ww 641 !! 642 !! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling. 643 !!---------------------------------------------------------------------- 644 INTEGER, INTENT(in ) :: kt ! time-step 645 REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT(in ) :: phdivtr ! horizontal divergence transport 646 REAL(wp), DIMENSION(jpi,jpj) , OPTIONAL, INTENT(in ) :: psshb ! now ssh 647 REAL(wp), DIMENSION(jpi,jpj) , OPTIONAL, INTENT(in ) :: pemp ! evaporation minus precipitation 648 REAL(wp), DIMENSION(jpi,jpj) , OPTIONAL, INTENT(inout) :: pssha ! after ssh 649 REAL(wp), DIMENSION(jpi,jpj,jpk), OPTIONAL, INTENT(out) :: pe3ta ! after vertical scale factor 650 ! 651 INTEGER :: jk 652 REAL(wp), DIMENSION(jpi,jpj) :: zhdiv 653 REAL(wp) :: z2dt 654 !!---------------------------------------------------------------------- 655 ! 656 z2dt = 2._wp * rn_Dt 657 ! 658 zhdiv(:,:) = 0._wp 659 DO jk = 1, jpkm1 660 zhdiv(:,:) = zhdiv(:,:) + phdivtr(:,:,jk) * tmask(:,:,jk) 661 END DO 662 ! ! Sea surface elevation time-stepping 663 pssha(:,:) = ( psshb(:,:) - z2dt * ( r1_rho0 * pemp(:,:) + zhdiv(:,:) ) ) * ssmask(:,:) 664 ! 665 IF( PRESENT( pe3ta ) ) THEN ! After acale factors at t-points ( z_star coordinate ) 666 DO jk = 1, jpkm1 667 pe3ta(:,:,jk) = e3t_0(:,:,jk) * ( 1._wp + pssha(:,:) * r1_ht_0(:,:) * tmask(:,:,jk) ) 668 END DO 669 ENDIF 670 ! 671 END SUBROUTINE dta_dyn_ssh 672 673 674 SUBROUTINE dta_dyn_hrnf( Kmm ) 675 !!---------------------------------------------------------------------- 676 !! *** ROUTINE sbc_rnf *** 539 END SUBROUTINE dta_dyn_rnf_init 540 541 SUBROUTINE dta_dyn_rnf( Kmm ) 542 !!---------------------------------------------------------------------- 543 !! *** ROUTINE dta_dyn_rnf *** 677 544 !! 678 545 !! ** Purpose : update the horizontal divergence with the runoff inflow 679 546 !! 680 !! ** Method :681 !! CAUTION : rnf is positive (inflow) decreasing the682 !! divergence and expressed in m/s683 !!684 !! ** Action : phdivn decreased by the runoff inflow685 547 !!---------------------------------------------------------------------- 686 548 !! … … 697 559 END_2D 698 560 ! 699 END SUBROUTINE dta_dyn_hrnf 700 701 561 END SUBROUTINE dta_dyn_rnf 702 562 703 563 SUBROUTINE dta_dyn_slp( kt, Kbb, Kmm ) … … 790 650 END SUBROUTINE dta_dyn_slp 791 651 792 793 652 SUBROUTINE compute_slopes( kt, pts, puslp, pvslp, pwslpi, pwslpj, Kbb, Kmm ) 794 653 !!--------------------------------------------------------------------- … … 835 694 END SUBROUTINE compute_slopes 836 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 837 813 !!====================================================================== 838 814 END MODULE dtadyn -
NEMO/trunk/src/OFF/nemogcm.F90
r14239 r14255 138 138 IF( istp /= nit000 ) CALL day ( istp ) ! Calendar (day was already called at nit000 in day_init) 139 139 CALL iom_setkt ( istp - nit000 + 1, cxios_context ) ! say to iom that we are at time step kstp 140 #if defined key_sed_off 141 CALL dta_dyn_sed( istp, Nnn ) ! Interpolation of the dynamical fields 142 #else 140 #if ! defined key_sed_off 143 141 CALL dta_dyn ( istp, Nbb, Nnn, Naa ) ! Interpolation of the dynamical fields 144 #endif145 #if ! defined key_sed_off146 142 IF( .NOT.ln_linssh ) THEN 147 143 CALL dta_dyn_atf( istp, Nbb, Nnn, Naa ) ! time filter of sea surface height and vertical scale factors … … 151 147 ENDIF 152 148 CALL trc_stp ( istp, Nbb, Nnn, Nrhs, Naa ) ! time-stepping 153 # if defined key_qco154 !r3t(:,:,Nnn) = r3t_f(:,:) ! update ssh to h0 ratio155 !r3u(:,:,Nnn) = r3u_f(:,:)156 !r3v(:,:,Nnn) = r3v_f(:,:)157 # endif158 #endif159 149 ! Swap time levels 160 150 Nrhs = Nbb 161 Nbb = Nnn162 Nnn = Naa163 Naa = Nrhs151 Nbb = Nnn 152 Nnn = Naa 153 Naa = Nrhs 164 154 ! 165 155 #if ! defined key_qco 166 # if ! defined key_sed_off167 156 IF( .NOT.ln_linssh ) CALL dta_dyn_sf_interp( istp, Nnn ) ! calculate now grid parameters 168 # endif 169 #endif 157 #endif 158 159 #else 160 CALL dta_dyn_sed( istp, Nnn ) ! Interpolation of the dynamical fields 161 162 #endif 170 163 CALL stp_ctl ( istp ) ! Time loop: control and print 171 164 istp = istp + 1 -
NEMO/trunk/src/TOP/trcwri.F90
r14239 r14255 28 28 PUBLIC trc_wri 29 29 30 !! * Substitutions 31 # include "do_loop_substitute.h90" 32 # include "domzgr_substitute.h90" 33 30 34 CONTAINS 31 35 … … 39 43 INTEGER, INTENT( in ) :: Kmm ! time level indices 40 44 ! 41 INTEGER :: j n45 INTEGER :: jk, jn 42 46 CHARACTER (len=20) :: cltra 43 47 CHARACTER (len=40) :: clhstnam 44 48 INTEGER :: inum = 11 ! temporary logical unit 49 REAL(wp), DIMENSION(jpi,jpj,jpk) :: z3d ! 3D workspace 45 50 !!--------------------------------------------------------------------- 46 51 ! … … 53 58 WRITE(inum,*) clhstnam 54 59 CLOSE(inum) 55 ENDIF60 ENDIF 56 61 57 ! Output of initial vertical scale factor 58 CALL iom_put( "e3t_0", e3t_0(:,:,:) ) 59 CALL iom_put( "e3u_0", e3u_0(:,:,:) ) 60 CALL iom_put( "e3v_0", e3v_0(:,:,:) ) 61 ! 62 #if ! defined key_qco 63 CALL iom_put( "e3t" , e3t(:,:,:,Kmm) ) 64 CALL iom_put( "e3u" , e3u(:,:,:,Kmm) ) 65 CALL iom_put( "e3v" , e3v(:,:,:,Kmm) ) 66 #endif 67 ! 62 ! Output of initial vertical scale factor 63 CALL iom_put( "e3t_0", e3t_0(:,:,:) ) 64 CALL iom_put( "e3u_0", e3u_0(:,:,:) ) 65 CALL iom_put( "e3v_0", e3v_0(:,:,:) ) 66 ! 67 IF( .NOT.ln_linssh ) CALL iom_put( "ssh" , ssh(:,:,Kmm) ) ! sea surface height 68 ! 69 IF ( iom_use("e3t") ) THEN ! time-varying e3t 70 DO jk = 1, jpk 71 z3d(:,:,jk) = e3t(:,:,jk,Kmm) 72 END DO 73 CALL iom_put( "e3t", z3d(:,:,:) ) 74 ENDIF 75 IF ( iom_use("e3u") ) THEN ! time-varying e3u 76 DO jk = 1, jpk 77 z3d(:,:,jk) = e3u(:,:,jk,Kmm) 78 END DO 79 CALL iom_put( "e3u", z3d(:,:,:) ) 80 ENDIF 81 IF ( iom_use("e3v") ) THEN ! time-varying e3v 82 DO jk = 1, jpk 83 z3d(:,:,jk) = e3v(:,:,jk,Kmm) 84 END DO 85 CALL iom_put( "e3v", z3d(:,:,:) ) 86 ENDIF 87 ! 68 88 ENDIF 89 ! 69 90 ! write the tracer concentrations in the file 70 91 ! ---------------------------------------
Note: See TracChangeset
for help on using the changeset viewer.
