Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • nemo/nemo
  • sparonuz/nemo
  • hatfield/nemo
  • extdevs/nemo
4 results
Show changes
Hide whitespace changes
Inline Side-by-side
Showing
with 932 additions and 330 deletions
src/OCE/SBC/sbcice_if.F90
View file @ 3efbcb17
...@@ -85,8 +85,8 @@ CONTAINS ...@@ -85,8 +85,8 @@ CONTAINS
ALLOCATE( sf_ice(1), STAT=ierror ) ALLOCATE( sf_ice(1), STAT=ierror )
IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'sbc_ice_if: unable to allocate sf_ice structure' ) IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'sbc_ice_if: unable to allocate sf_ice structure' )
ALLOCATE( sf_ice(1)%fnow(jpi,jpj,1) ) ALLOCATE( sf_ice(1)%fnow(A2D(0),1) )
IF( sn_ice%ln_tint ) ALLOCATE( sf_ice(1)%fdta(jpi,jpj,1,2) ) IF( sn_ice%ln_tint ) ALLOCATE( sf_ice(1)%fdta(A2D(0),1,2) )
! fill sf_ice with sn_ice and control print ! fill sf_ice with sn_ice and control print
CALL fld_fill( sf_ice, (/ sn_ice /), cn_dir, 'sbc_ice_if', 'ice-if sea-ice model', 'namsbc_iif' ) CALL fld_fill( sf_ice, (/ sn_ice /), cn_dir, 'sbc_ice_if', 'ice-if sea-ice model', 'namsbc_iif' )
...@@ -108,8 +108,12 @@ CONTAINS ...@@ -108,8 +108,12 @@ CONTAINS
IF( ln_cpl ) a_i(:,:,1) = fr_i(:,:) IF( ln_cpl ) a_i(:,:,1) = fr_i(:,:)
! Flux and ice fraction computation ! Flux and ice fraction computation
DO_2D( 1, 1, 1, 1 ) DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
! zt_fzp = fr_i(ji,jj) ! freezing point temperature
ts(ji,jj,1,jp_tem,Kmm) = MAX( ts(ji,jj,1,jp_tem,Kmm), zt_fzp ) ! avoid over-freezing point temperature
END_2D
DO_2D( 0, 0, 0, 0 )
zt_fzp = fr_i(ji,jj) ! freezing point temperature zt_fzp = fr_i(ji,jj) ! freezing point temperature
zfr_obs = sf_ice(1)%fnow(ji,jj,1) ! observed ice cover zfr_obs = sf_ice(1)%fnow(ji,jj,1) ! observed ice cover
! ! ocean ice fraction (0/1) from the freezing point temperature ! ! ocean ice fraction (0/1) from the freezing point temperature
...@@ -117,8 +121,6 @@ CONTAINS ...@@ -117,8 +121,6 @@ CONTAINS
ELSE ; fr_i(ji,jj) = 0.e0 ELSE ; fr_i(ji,jj) = 0.e0
ENDIF ENDIF
ts(ji,jj,1,jp_tem,Kmm) = MAX( ts(ji,jj,1,jp_tem,Kmm), zt_fzp ) ! avoid over-freezing point temperature
qsr(ji,jj) = ( 1. - zfr_obs ) * qsr(ji,jj) ! solar heat flux : zero below observed ice cover qsr(ji,jj) = ( 1. - zfr_obs ) * qsr(ji,jj) ! solar heat flux : zero below observed ice cover
! ! non solar heat flux : add a damping term ! ! non solar heat flux : add a damping term
...@@ -127,7 +129,7 @@ CONTAINS ...@@ -127,7 +129,7 @@ CONTAINS
zqri = ztrp * ( ts(ji,jj,1,jp_tem,Kbb) - ( zt_fzp - 1.) ) zqri = ztrp * ( ts(ji,jj,1,jp_tem,Kbb) - ( zt_fzp - 1.) )
zqrj = ztrp * MIN( 0., ts(ji,jj,1,jp_tem,Kbb) - zt_fzp ) zqrj = ztrp * MIN( 0., ts(ji,jj,1,jp_tem,Kbb) - zt_fzp )
zqrp = ( zfr_obs * ( (1. - fr_i(ji,jj) ) * zqri & zqrp = ( zfr_obs * ( (1. - fr_i(ji,jj) ) * zqri &
& + fr_i(ji,jj) * zqrj ) ) * tmask(ji,jj,1) & + fr_i(ji,jj) * zqrj ) ) * smask0(ji,jj)
! ! non-solar heat flux ! ! non-solar heat flux
! # qns unchanged if no climatological ice (zfr_obs=0) ! # qns unchanged if no climatological ice (zfr_obs=0)
...@@ -136,7 +138,7 @@ CONTAINS ...@@ -136,7 +138,7 @@ CONTAINS
! (-2=arctic, -4=antarctic) ! (-2=arctic, -4=antarctic)
zqi = -3. + SIGN( 1._wp, ff_f(ji,jj) ) zqi = -3. + SIGN( 1._wp, ff_f(ji,jj) )
qns(ji,jj) = ( ( 1.- zfr_obs ) * qns(ji,jj) & qns(ji,jj) = ( ( 1.- zfr_obs ) * qns(ji,jj) &
& + zfr_obs * fr_i(ji,jj) * zqi ) * tmask(ji,jj,1) & & + zfr_obs * fr_i(ji,jj) * zqi ) * smask0(ji,jj) &
& + zqrp & + zqrp
END_2D END_2D
! !
......
src/OCE/SBC/sbcmod.F90
View file @ 3efbcb17
...@@ -51,6 +51,7 @@ MODULE sbcmod ...@@ -51,6 +51,7 @@ MODULE sbcmod
USE sbcfwb ! surface boundary condition: freshwater budget USE sbcfwb ! surface boundary condition: freshwater budget
USE icbstp ! Icebergs USE icbstp ! Icebergs
USE icb_oce , ONLY : ln_passive_mode ! iceberg interaction mode USE icb_oce , ONLY : ln_passive_mode ! iceberg interaction mode
USE isf_oce , ONLY : ln_isf, l_isfoasis, fwfisf_oasis
USE traqsr ! active tracers: light penetration USE traqsr ! active tracers: light penetration
USE sbcwave ! Wave module USE sbcwave ! Wave module
USE bdy_oce , ONLY: ln_bdy USE bdy_oce , ONLY: ln_bdy
...@@ -158,8 +159,8 @@ CONTAINS ...@@ -158,8 +159,8 @@ CONTAINS
! !
IF( .NOT.ln_usr ) THEN ! the model calendar needs some specificities (except in user defined case) IF( .NOT.ln_usr ) THEN ! the model calendar needs some specificities (except in user defined case)
IF( MOD( rday , rn_Dt ) /= 0. ) CALL ctl_stop( 'the time step must devide the number of second of in a day' ) IF( MOD( rday , rn_Dt ) /= 0. ) CALL ctl_stop( 'the time step must devide the number of second of in a day' )
IF( MOD( rday , 2. ) /= 0. ) CALL ctl_stop( 'the number of second of in a day must be an even number' ) IF( MOD( rday , 2. ) /= 0. ) CALL ctl_stop( 'the number of second of in a day must be an even number' )
IF( MOD( rn_Dt , 2. ) /= 0. ) CALL ctl_stop( 'the time step (in second) must be an even number' ) IF( MOD( rn_Dt, 2. ) /= 0. ) CALL ctl_stop( 'the time step (in second) must be an even number' )
ENDIF ENDIF
! !** check option consistency ! !** check option consistency
! !
...@@ -374,7 +375,7 @@ CONTAINS ...@@ -374,7 +375,7 @@ CONTAINS
!! !!
!! ** Action : - set the ocean surface boundary condition at before and now !! ** Action : - set the ocean surface boundary condition at before and now
!! time step, i.e. !! time step, i.e.
!! utau_b, vtau_b, qns_b, qsr_b, emp_n, sfx_b, qrp_b, erp_b !! utau_b, vtau_b, qns_b, qsr_b, emp_b, sfx_b
!! utau , vtau , qns , qsr , emp , sfx , qrp , erp !! utau , vtau , qns , qsr , emp , sfx , qrp , erp
!! - updte the ice fraction : fr_i !! - updte the ice fraction : fr_i
!!---------------------------------------------------------------------- !!----------------------------------------------------------------------
...@@ -382,12 +383,10 @@ CONTAINS ...@@ -382,12 +383,10 @@ CONTAINS
INTEGER, INTENT(in) :: Kbb, Kmm ! ocean time level indices INTEGER, INTENT(in) :: Kbb, Kmm ! ocean time level indices
INTEGER :: jj, ji ! dummy loop argument INTEGER :: jj, ji ! dummy loop argument
! !
LOGICAL :: ll_sas, ll_opa ! local logical LOGICAL :: ll_sas, ll_opa ! local logical
! !
REAL(wp) :: zthscl ! wd tanh scale REAL(wp) :: zthscl ! wd tanh scale
REAL(wp), DIMENSION(jpi,jpj) :: zwdht, zwght ! wd dep over wd limit, wgt REAL(wp) :: zwdht, zwght ! wd dep over wd limit, wgt
REAL(wp), DIMENSION(jpi,jpj) :: z2d ! temporary array used for iom_put
!!--------------------------------------------------------------------- !!---------------------------------------------------------------------
! !
IF( ln_timing ) CALL timing_start('sbc') IF( ln_timing ) CALL timing_start('sbc')
...@@ -395,20 +394,22 @@ CONTAINS ...@@ -395,20 +394,22 @@ CONTAINS
! ! ---------------------------------------- ! ! ! ---------------------------------------- !
IF( kt /= nit000 ) THEN ! Swap of forcing fields ! IF( kt /= nit000 ) THEN ! Swap of forcing fields !
! ! ---------------------------------------- ! ! ! ---------------------------------------- !
utau_b(:,:) = utau(:,:) ! Swap the ocean forcing fields utau_b(:,:) = utauU(:,:) ! Swap the ocean forcing fields
vtau_b(:,:) = vtau(:,:) ! (except at nit000 where before fields vtau_b(:,:) = vtauV(:,:) ! (except at nit000 where before fields
qns_b (:,:) = qns (:,:) ! are set at the end of the routine) qns_b (:,:) = qns (:,:) ! are set at the end of the routine)
emp_b (:,:) = emp (:,:) emp_b (:,:) = emp (:,:)
sfx_b (:,:) = sfx (:,:) sfx_b (:,:) = sfx (:,:)
IF( ln_rnf ) THEN IF( ln_rnf ) THEN
rnf_b (:,: ) = rnf (:,: ) rnf_b (:,: ) = rnf (:,: )
rnf_tsc_b(:,:,:) = rnf_tsc(:,:,:) rnf_tsc_b(:,:,:) = rnf_tsc(:,:,:)
ENDIF ENDIF
! !
ENDIF ENDIF
! ! ---------------------------------------- ! ! ! ---------------------------------------- !
! ! forcing field computation ! ! ! forcing field computation !
! ! ---------------------------------------- ! ! ! ---------------------------------------- !
! most of the following routines update fields only in the interior
! with the exception of sbcssm, sbcrnf and sbcwave modules
! !
ll_sas = nn_components == jp_iam_sas ! component flags ll_sas = nn_components == jp_iam_sas ! component flags
ll_opa = nn_components == jp_iam_oce ll_opa = nn_components == jp_iam_oce
...@@ -417,87 +418,95 @@ CONTAINS ...@@ -417,87 +418,95 @@ CONTAINS
! !
! !== sbc formulation ==! ! !== sbc formulation ==!
! !
!
SELECT CASE( nsbc ) ! Compute ocean surface boundary condition SELECT CASE( nsbc ) ! Compute ocean surface boundary condition
! ! (i.e. utau,vtau, qns, qsr, emp, sfx) ! ! (i.e. utau,vtau, qns, qsr, emp, sfx)
CASE( jp_usr ) ; CALL usrdef_sbc_oce( kt, Kbb ) ! user defined formulation CASE( jp_usr ) ; CALL usrdef_sbc_oce( kt, Kbb ) ! user defined formulation
CASE( jp_flx ) ; CALL sbc_flx ( kt ) ! flux formulation CASE( jp_flx ) ; CALL sbc_flx ( kt ) ! flux formulation
CASE( jp_blk ) CASE( jp_blk )
IF( ll_sas ) CALL sbc_cpl_rcv ( kt, nn_fsbc, nn_ice, Kbb, Kmm ) ! OCE-SAS coupling: SAS receiving fields from OCE IF( ll_sas ) CALL sbc_cpl_rcv ( kt, nn_fsbc, nn_ice, Kbb, Kmm ) ! OCE-SAS coupling: SAS receiving fields from OCE
!!!!!!!!!!! ATTENTION:ln_wave is not only used for oasis coupling !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!! ATTENTION:ln_wave is not only used for oasis coupling !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
IF( ln_wave ) THEN IF( ln_wave ) THEN
IF ( lk_oasis ) CALL sbc_cpl_rcv ( kt, nn_fsbc, nn_ice, Kbb, Kmm ) ! OCE-wave coupling IF ( lk_oasis ) CALL sbc_cpl_rcv ( kt, nn_fsbc, nn_ice, Kbb, Kmm ) ! OCE-wave coupling
CALL sbc_wave ( kt, Kmm ) CALL sbc_wave ( kt, Kmm )
ENDIF ENDIF
CALL sbc_blk ( kt ) ! bulk formulation for the ocean CALL sbc_blk ( kt ) ! bulk formulation for the ocean
! !
CASE( jp_abl ) CASE( jp_abl )
IF( ll_sas ) CALL sbc_cpl_rcv ( kt, nn_fsbc, nn_ice, Kbb, Kmm ) ! OCE-SAS coupling: SAS receiving fields from OCE IF( ll_sas ) CALL sbc_cpl_rcv ( kt, nn_fsbc, nn_ice, Kbb, Kmm ) ! OCE-SAS coupling: SAS receiving fields from OCE
CALL sbc_abl ( kt ) ! ABL formulation for the ocean CALL sbc_abl ( kt ) ! ABL formulation for the ocean
! !
CASE( jp_purecpl ) ; CALL sbc_cpl_rcv ( kt, nn_fsbc, nn_ice, Kbb, Kmm ) ! pure coupled formulation CASE( jp_purecpl ) ; CALL sbc_cpl_rcv ( kt, nn_fsbc, nn_ice, Kbb, Kmm ) ! pure coupled formulation
CASE( jp_none ) CASE( jp_none )
IF( ll_opa ) CALL sbc_cpl_rcv ( kt, nn_fsbc, nn_ice, Kbb, Kmm ) ! OCE-SAS coupling: OCE receiving fields from SAS IF( ll_opa ) CALL sbc_cpl_rcv ( kt, nn_fsbc, nn_ice, Kbb, Kmm ) ! OCE-SAS coupling: OCE receiving fields from SAS
END SELECT END SELECT
! !
IF( ln_mixcpl ) CALL sbc_cpl_rcv ( kt, nn_fsbc, nn_ice, Kbb, Kmm ) ! forced-coupled mixed formulation after forcing IF( ln_mixcpl ) CALL sbc_cpl_rcv ( kt, nn_fsbc, nn_ice, Kbb, Kmm ) ! forced-coupled mixed formulation after forcing
! !
IF( ln_wave .AND. ln_tauoc ) THEN ! Wave stress reduction IF( ln_wave .AND. ln_tauoc ) THEN ! Wave stress reduction
DO_2D( 0, 0, 0, 0)
utau(ji,jj) = utau(ji,jj) * ( tauoc_wave(ji,jj) + tauoc_wave(ji+1,jj) ) * 0.5_wp
vtau(ji,jj) = vtau(ji,jj) * ( tauoc_wave(ji,jj) + tauoc_wave(ji,jj+1) ) * 0.5_wp
END_2D
!
CALL lbc_lnk( 'sbcwave', utau, 'U', -1. )
CALL lbc_lnk( 'sbcwave', vtau, 'V', -1. )
! !
taum(:,:) = taum(:,:)*tauoc_wave(:,:) DO_2D( 0, 0, 0, 0 )
utau(ji,jj) = utau(ji,jj) * tauoc_wave(ji,jj)
vtau(ji,jj) = vtau(ji,jj) * tauoc_wave(ji,jj)
taum(ji,jj) = taum(ji,jj) * tauoc_wave(ji,jj)
END_2D
! !
IF( kt == nit000 ) CALL ctl_warn( 'sbc: You are subtracting the wave stress to the ocean.', & IF( kt == nit000 ) CALL ctl_warn( 'sbc: You are subtracting the wave stress to the ocean.', &
& 'If not requested select ln_tauoc=.false.' ) & 'If not requested select ln_tauoc=.false.' )
! !
ELSEIF( ln_wave .AND. ln_taw ) THEN ! Wave stress reduction ELSEIF( ln_wave .AND. ln_taw ) THEN ! Wave stress reduction
utau(:,:) = utau(:,:) - tawx(:,:) + twox(:,:) DO_2D( 0, 0, 0, 0 )
vtau(:,:) = vtau(:,:) - tawy(:,:) + twoy(:,:) utau(ji,jj) = utau(ji,jj) - tawx(ji,jj) + twox(ji,jj)
CALL lbc_lnk( 'sbcwave', utau, 'U', -1. ) vtau(ji,jj) = vtau(ji,jj) - tawy(ji,jj) + twoy(ji,jj)
CALL lbc_lnk( 'sbcwave', vtau, 'V', -1. ) taum(ji,jj) = SQRT( utau(ji,jj)*utau(ji,jj) + vtau(ji,jj)*vtau(ji,jj) )
!
DO_2D( 0, 0, 0, 0)
taum(ji,jj) = sqrt((.5*(utau(ji-1,jj)+utau(ji,jj)))**2 + (.5*(vtau(ji,jj-1)+vtau(ji,jj)))**2)
END_2D END_2D
! !
IF( kt == nit000 ) CALL ctl_warn( 'sbc: You are subtracting the wave stress to the ocean.', & IF( kt == nit000 ) CALL ctl_warn( 'sbc: You are subtracting the wave stress to the ocean.', &
& 'If not requested select ln_taw=.false.' ) & 'If not requested select ln_taw=.false.' )
! !
ENDIF ENDIF
CALL lbc_lnk( 'sbcmod', taum(:,:), 'T', 1. ) !
!clem: these calls are needed for sbccpl only => only for SAS I think?
IF( ll_sas .OR. ll_opa ) CALL lbc_lnk( 'sbcmod', sst_m, 'T', 1.0_wp, sss_m, 'T', 1.0_wp, ssh_m, 'T', 1.0_wp, &
& frq_m, 'T', 1.0_wp, e3t_m, 'T', 1.0_wp, fr_i , 'T', 1.0_wp )
!clem : these calls are needed for sbccpl => it needs an IF statement but it's complicated
IF( ln_isf .AND. l_isfoasis ) CALL lbc_lnk( 'sbcmod', fwfisf_oasis, 'T', 1.0_wp )
IF( ln_rnf .AND. l_rnfcpl ) CALL lbc_lnk( 'sbcmod', rnf, 'T', 1.0_wp, fwficb , 'T', 1.0_wp )
! !
IF( ln_icebergs ) THEN ! save pure stresses (with no ice-ocean stress) for use by icebergs IF( ln_icebergs ) THEN ! save pure stresses (with no ice-ocean stress) for use by icebergs
utau_icb(:,:) = utau(:,:) ; vtau_icb(:,:) = vtau(:,:) ! Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines
! and the use of MAX(tmask(i,j),tmask(i+1,j) is to mask tau over ice shelves
CALL lbc_lnk( 'sbcmod', utau, 'T', -1.0_wp, vtau, 'T', -1.0_wp )
DO_2D( 0, 0, 0, 0 )
utau_icb(ji,jj) = 0.5_wp * ( utau(ji,jj) + utau(ji+1,jj) ) * &
& ( 2. - umask(ji,jj,1) ) * MAX( tmask(ji,jj,1), tmask(ji+1,jj,1) )
vtau_icb(ji,jj) = 0.5_wp * ( vtau(ji,jj) + vtau(ji,jj+1) ) * &
& ( 2. - vmask(ji,jj,1) ) * MAX( tmask(ji,jj,1), tmask(ji,jj+1,1) )
END_2D
CALL lbc_lnk( 'sbcmod', utau_icb, 'U', -1.0_wp, vtau_icb, 'V', -1.0_wp )
ENDIF ENDIF
! !
! !== Misc. Options ==! ! !== Misc. Options ==!
! !
SELECT CASE( nn_ice ) ! Update heat and freshwater fluxes over sea-ice areas SELECT CASE( nn_ice ) ! Update heat and freshwater fluxes over sea-ice areas
CASE( 1 ) ; CALL sbc_ice_if ( kt, Kbb, Kmm ) ! Ice-cover climatology ("Ice-if" model) CASE( 1 ) ; CALL sbc_ice_if ( kt, Kbb, Kmm ) ! Ice-cover climatology ("Ice-if" model)
#if defined key_si3 #if defined key_si3
CASE( 2 ) ; CALL ice_stp ( kt, Kbb, Kmm, nsbc ) ! SI3 ice model CASE( 2 ) ; CALL ice_stp ( kt, Kbb, Kmm, nsbc ) ! SI3 ice model
#endif #endif
CASE( 3 ) ; CALL sbc_ice_cice ( kt, nsbc ) ! CICE ice model CASE( 3 ) ; CALL sbc_ice_cice ( kt, nsbc ) ! CICE ice model
END SELECT END SELECT
!==> clem: from here on, the following fields are ok on the halos: snwice_mass, snwice_mass_b, snwice_fmass
IF( ln_icebergs ) CALL icb_stp( kt, Kmm ) ! compute icebergs ! but not utau, vtau, emp (must be done later on)
IF( ln_icebergs ) CALL icb_stp( kt, Kmm ) ! compute icebergs
! Icebergs do not melt over the haloes. ! Icebergs do not melt over the haloes.
! So emp values over the haloes are no more consistent with the inner domain values. ! So emp values over the haloes are no more consistent with the inner domain values.
! A lbc_lnk is therefore needed to ensure reproducibility and restartability. ! A lbc_lnk is therefore needed to ensure reproducibility and restartability.
! see ticket #2113 for discussion about this lbc_lnk. ! see ticket #2113 for discussion about this lbc_lnk.
! The lbc_lnk is also needed for SI3 with nn_hls > 1 as emp is not yet defined for these points in iceupdate.F90 !!$ IF( ln_icebergs .AND. .NOT. ln_passive_mode ) CALL lbc_lnk( 'sbcmod', emp, 'T', 1.0_wp )
IF( (ln_icebergs .AND. .NOT. ln_passive_mode) .OR. (nn_ice == 2 .AND. nn_hls == 2) ) THEN !clem: not needed anymore since lbc is done afterwards
CALL lbc_lnk( 'sbcmod', emp, 'T', 1.0_wp )
ENDIF IF( ln_rnf ) CALL sbc_rnf( kt ) ! add runoffs to fresh water fluxes
IF( ln_rnf ) CALL sbc_rnf( kt ) ! add runoffs to fresh water fluxes
IF( ln_ssr ) CALL sbc_ssr( kt ) ! add SST/SSS damping term IF( ln_ssr ) CALL sbc_ssr( kt ) ! add SST/SSS damping term
...@@ -507,33 +516,49 @@ CONTAINS ...@@ -507,33 +516,49 @@ CONTAINS
! Should not be run if ln_diurnal_only ! Should not be run if ln_diurnal_only
IF( l_sbc_clo ) CALL sbc_clo( kt ) IF( l_sbc_clo ) CALL sbc_clo( kt )
!!$!RBbug do not understand why see ticket 667
!!$!clem: it looks like it is necessary for the north fold (in certain circumstances). Don't know why.
!!$ CALL lbc_lnk( 'sbcmod', emp, 'T', 1.0_wp )
IF( ll_wd ) THEN ! If near WAD point limit the flux for now IF( ll_wd ) THEN ! If near WAD point limit the flux for now
zthscl = atanh(rn_wd_sbcfra) ! taper frac default is .999 zthscl = atanh(rn_wd_sbcfra) ! taper frac default is .999
zwdht(:,:) = ssh(:,:,Kmm) + ht_0(:,:) - rn_wdmin1 ! do this calc of water DO_2D( 0, 0, 0, 0 )
! depth above wd limit once zwdht = ssh(ji,jj,Kmm) + ht_0(ji,jj) - rn_wdmin1 ! do this calc of water depth above wd limit once
WHERE( zwdht(:,:) <= 0.0 ) zwght = TANH(zthscl*zwdht)
taum(:,:) = 0.0 IF( zwdht <= 0.0 ) THEN
utau(:,:) = 0.0 taum(ji,jj) = 0.0
vtau(:,:) = 0.0 utau(ji,jj) = 0.0
qns (:,:) = 0.0 vtau(ji,jj) = 0.0
qsr (:,:) = 0.0 qns (ji,jj) = 0.0
emp (:,:) = min(emp(:,:),0.0) !can allow puddles to grow but not shrink qsr (ji,jj) = 0.0
sfx (:,:) = 0.0 emp (ji,jj) = MIN(emp(ji,jj),0.0) !can allow puddles to grow but not shrink
END WHERE sfx (ji,jj) = 0.0
zwght(:,:) = tanh(zthscl*zwdht(:,:)) ELSEIF( zwdht > 0.0 .AND. zwdht < rn_wd_sbcdep ) THEN ! 5 m hard limit here is arbitrary
WHERE( zwdht(:,:) > 0.0 .and. zwdht(:,:) < rn_wd_sbcdep ) ! 5 m hard limit here is arbitrary qsr (ji,jj) = qsr(ji,jj) * zwght
qsr (:,:) = qsr(:,:) * zwght(:,:) qns (ji,jj) = qns(ji,jj) * zwght
qns (:,:) = qns(:,:) * zwght(:,:) taum (ji,jj) = taum(ji,jj) * zwght
taum (:,:) = taum(:,:) * zwght(:,:) utau (ji,jj) = utau(ji,jj) * zwght
utau (:,:) = utau(:,:) * zwght(:,:) vtau (ji,jj) = vtau(ji,jj) * zwght
vtau (:,:) = vtau(:,:) * zwght(:,:) sfx (ji,jj) = sfx(ji,jj) * zwght
sfx (:,:) = sfx(:,:) * zwght(:,:) emp (ji,jj) = emp(ji,jj) * zwght
emp (:,:) = emp(:,:) * zwght(:,:) ENDIF
END WHERE END_2D
ENDIF
! clem: these should be the only fields that are needed over the entire domain
! (in addition to snwice_mass)
IF( ln_rnf ) THEN
CALL lbc_lnk( 'sbcmod', utau, 'T', -1.0_wp, vtau , 'T', -1.0_wp, emp, 'T', 1.0_wp, &
& rnf , 'T', 1.0_wp, fwficb, 'T', 1.0_wp ) ! fwficb is used on the halos in pisces (only)
ELSE
CALL lbc_lnk( 'sbcmod', utau, 'T', -1.0_wp, vtau , 'T', -1.0_wp, emp, 'T', 1.0_wp )
ENDIF ENDIF
! --- calculate utau and vtau on U,V-points --- !
! Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines
! and the use of MAX(tmask(i,j),tmask(i+1,j) is to mask tau over ice shelves
DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
utauU (ji,jj) = 0.5_wp * ( utau(ji,jj) + utau(ji+1,jj) ) * &
& ( 2. - umask(ji,jj,1) ) * MAX( tmask(ji,jj,1), tmask(ji+1,jj,1) )
vtauV (ji,jj) = 0.5_wp * ( vtau(ji,jj) + vtau(ji,jj+1) ) * &
& ( 2. - vmask(ji,jj,1) ) * MAX( tmask(ji,jj,1), tmask(ji,jj+1,1) )
END_2D
IF( nn_hls == 1 ) CALL lbc_lnk( 'sbcmod', utauU, 'U', -1.0_wp, vtauV, 'V', -1.0_wp )
! !
IF( kt == nit000 ) THEN ! set the forcing field at nit000 - 1 ! IF( kt == nit000 ) THEN ! set the forcing field at nit000 - 1 !
! ! ---------------------------------------- ! ! ! ---------------------------------------- !
...@@ -543,27 +568,28 @@ CONTAINS ...@@ -543,27 +568,28 @@ CONTAINS
IF( ln_rstart .AND. .NOT.l_1st_euler ) THEN !* MLF: Restart: read in restart file IF( ln_rstart .AND. .NOT.l_1st_euler ) THEN !* MLF: Restart: read in restart file
#endif #endif
IF(lwp) WRITE(numout,*) ' nit000-1 surface forcing fields read in the restart file' IF(lwp) WRITE(numout,*) ' nit000-1 surface forcing fields read in the restart file'
CALL iom_get( numror, jpdom_auto, 'utau_b', utau_b ) ! i-stress CALL iom_get( numror, jpdom_auto, 'utau_b', utau_b, cd_type = 'U', psgn = -1._wp ) ! i-stress
CALL iom_get( numror, jpdom_auto, 'vtau_b', vtau_b ) ! j-stress CALL iom_get( numror, jpdom_auto, 'vtau_b', vtau_b, cd_type = 'V', psgn = -1._wp ) ! j-stress
CALL iom_get( numror, jpdom_auto, 'qns_b', qns_b ) ! non solar heat flux CALL iom_get( numror, jpdom_auto, 'qns_b', qns_b, cd_type = 'T', psgn = 1._wp ) ! non solar heat flux
CALL iom_get( numror, jpdom_auto, 'emp_b', emp_b ) ! freshwater flux CALL iom_get( numror, jpdom_auto, 'emp_b', emp_b, cd_type = 'T', psgn = 1._wp ) ! freshwater flux
! NB: The 3D heat content due to qsr forcing (qsr_hc_b) is treated in traqsr ! NB: The 3D heat content due to qsr forcing (qsr_hc_b) is treated in traqsr
! To ensure restart capability with 3.3x/3.4 restart files !! to be removed in v3.6 ! To ensure restart capability with 3.3x/3.4 restart files !! to be removed in v3.6
IF( iom_varid( numror, 'sfx_b', ldstop = .FALSE. ) > 0 ) THEN IF( iom_varid( numror, 'sfx_b', ldstop = .FALSE. ) > 0 ) THEN
CALL iom_get( numror, jpdom_auto, 'sfx_b', sfx_b ) ! before salt flux (T-point) CALL iom_get( numror, jpdom_auto, 'sfx_b', sfx_b, cd_type = 'T', psgn = 1._wp ) ! before salt flux (T-point)
ELSE ELSE
sfx_b (:,:) = sfx(:,:) sfx_b (:,:) = sfx(:,:)
ENDIF ENDIF
ELSE !* no restart: set from nit000 values ELSE !* no restart: set from nit000 values
IF(lwp) WRITE(numout,*) ' nit000-1 surface forcing fields set to nit000' IF(lwp) WRITE(numout,*) ' nit000-1 surface forcing fields set to nit000'
utau_b(:,:) = utau(:,:) utau_b(:,:) = utauU(:,:)
vtau_b(:,:) = vtau(:,:) vtau_b(:,:) = vtauV(:,:)
qns_b (:,:) = qns (:,:) qns_b (:,:) = qns (:,:)
emp_b (:,:) = emp (:,:) emp_b (:,:) = emp (:,:)
sfx_b (:,:) = sfx (:,:) sfx_b (:,:) = sfx (:,:)
ENDIF ENDIF
ENDIF ENDIF
! !
!
#if defined key_RK3 #if defined key_RK3
! ! ---------------------------------------- ! ! ! ---------------------------------------- !
IF( lrst_oce .AND. lk_SWE ) THEN ! RK3: Write in the ocean restart file ! IF( lrst_oce .AND. lk_SWE ) THEN ! RK3: Write in the ocean restart file !
...@@ -578,41 +604,26 @@ CONTAINS ...@@ -578,41 +604,26 @@ CONTAINS
IF(lwp) WRITE(numout,*) 'sbc : ocean surface forcing fields written in ocean restart file ', & IF(lwp) WRITE(numout,*) 'sbc : ocean surface forcing fields written in ocean restart file ', &
& 'at it= ', kt,' date= ', ndastp & 'at it= ', kt,' date= ', ndastp
IF(lwp) WRITE(numout,*) '~~~~' IF(lwp) WRITE(numout,*) '~~~~'
CALL iom_rstput( kt, nitrst, numrow, 'utau_b' , utau ) CALL iom_rstput( kt, nitrst, numrow, 'utau_b' , utauU )
CALL iom_rstput( kt, nitrst, numrow, 'vtau_b' , vtau ) CALL iom_rstput( kt, nitrst, numrow, 'vtau_b' , vtauV )
CALL iom_rstput( kt, nitrst, numrow, 'qns_b' , qns ) CALL iom_rstput( kt, nitrst, numrow, 'qns_b' , qns )
! The 3D heat content due to qsr forcing is treated in traqsr ! The 3D heat content due to qsr forcing is treated in traqsr
! CALL iom_rstput( kt, nitrst, numrow, 'qsr_b' , qsr ) ! CALL iom_rstput( kt, nitrst, numrow, 'qsr_b' , qsr )
CALL iom_rstput( kt, nitrst, numrow, 'emp_b' , emp ) CALL iom_rstput( kt, nitrst, numrow, 'emp_b' , emp )
CALL iom_rstput( kt, nitrst, numrow, 'sfx_b' , sfx ) CALL iom_rstput( kt, nitrst, numrow, 'sfx_b' , sfx )
ENDIF ENDIF
! ! ---------------------------------------- ! ! ! ---------------------------------------- !
! ! Outputs and control print ! ! ! Outputs and control print !
! ! ---------------------------------------- ! ! ! ---------------------------------------- !
IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN
IF( iom_use("empmr") ) THEN CALL iom_put( "empmr" , emp(A2D(0))-rnf(A2D(0)) ) ! upward water flux
DO_2D( 0, 0, 0, 0 ) CALL iom_put( "empbmr" , emp_b(A2D(0))-rnf(A2D(0)) ) ! before upward water flux (for ssh in offline )
z2d(ji,jj) = emp(ji,jj) - rnf(ji,jj) CALL iom_put( "saltflx", sfx ) ! downward salt flux
END_2D
CALL iom_put( "empmr" , z2d ) ! upward water flux
ENDIF
IF( iom_use("empbmr") ) THEN
DO_2D( 0, 0, 0, 0 )
z2d(ji,jj) = emp_b(ji,jj) - rnf(ji,jj)
END_2D
CALL iom_put( "empbmr" , z2d ) ! before upward water flux ( needed to recalculate the time evolution of ssh in offline )
ENDIF
CALL iom_put( "saltflx", sfx ) ! downward salt flux (includes virtual salt flux beneath ice in linear free surface case)
CALL iom_put( "fmmflx" , fmmflx ) ! Freezing-melting water flux CALL iom_put( "fmmflx" , fmmflx ) ! Freezing-melting water flux
IF( iom_use("qt") ) THEN CALL iom_put( "qt" , qns+qsr ) ! total heat flux
DO_2D( 0, 0, 0, 0 )
z2d(ji,jj) = qns(ji,jj) + qsr(ji,jj)
END_2D
CALL iom_put( "qt" , z2d ) ! total heat flux
ENDIF
CALL iom_put( "qns" , qns ) ! solar heat flux CALL iom_put( "qns" , qns ) ! solar heat flux
CALL iom_put( "qsr" , qsr ) ! solar heat flux CALL iom_put( "qsr" , qsr ) ! solar heat flux
IF( nn_ice > 0 .OR. ll_opa ) CALL iom_put( "ice_cover", fr_i ) ! ice fraction IF( nn_ice > 0 .OR. ll_opa ) CALL iom_put( "ice_cover", fr_i(:,:) ) ! ice fraction
CALL iom_put( "taum" , taum ) ! wind stress module CALL iom_put( "taum" , taum ) ! wind stress module
CALL iom_put( "wspd" , wndm ) ! wind speed module over free ocean or leads in presence of sea-ice CALL iom_put( "wspd" , wndm ) ! wind speed module over free ocean or leads in presence of sea-ice
CALL iom_put( "qrp" , qrp ) ! heat flux damping CALL iom_put( "qrp" , qrp ) ! heat flux damping
...@@ -622,14 +633,14 @@ CONTAINS ...@@ -622,14 +633,14 @@ CONTAINS
IF(sn_cfctl%l_prtctl) THEN ! print mean trends (used for debugging) IF(sn_cfctl%l_prtctl) THEN ! print mean trends (used for debugging)
CALL prt_ctl(tab2d_1=fr_i , clinfo1=' fr_i - : ', mask1=tmask ) CALL prt_ctl(tab2d_1=fr_i , clinfo1=' fr_i - : ', mask1=tmask )
CALL prt_ctl(tab2d_1=(emp-rnf) , clinfo1=' emp-rnf - : ', mask1=tmask ) CALL prt_ctl(tab2d_1=(emp-rnf) , clinfo1=' emp-rnf - : ', mask1=tmask )
CALL prt_ctl(tab2d_1=(sfx-rnf) , clinfo1=' sfx-rnf - : ', mask1=tmask ) CALL prt_ctl(tab2d_1=(sfx-rnf(A2D(0))) , clinfo1=' sfx-rnf - : ', mask1=tmask )
CALL prt_ctl(tab2d_1=qns , clinfo1=' qns - : ', mask1=tmask ) CALL prt_ctl(tab2d_1=qns , clinfo1=' qns - : ', mask1=tmask )
CALL prt_ctl(tab2d_1=qsr , clinfo1=' qsr - : ', mask1=tmask ) CALL prt_ctl(tab2d_1=qsr , clinfo1=' qsr - : ', mask1=tmask )
CALL prt_ctl(tab3d_1=tmask , clinfo1=' tmask - : ', mask1=tmask, kdim=jpk ) CALL prt_ctl(tab3d_1=tmask , clinfo1=' tmask - : ', mask1=tmask, kdim=jpk )
CALL prt_ctl(tab3d_1=ts(:,:,:,jp_tem,Kmm), clinfo1=' sst - : ', mask1=tmask, kdim=1 ) CALL prt_ctl(tab2d_1=sst_m , clinfo1=' sst - : ', mask1=tmask )
CALL prt_ctl(tab3d_1=ts(:,:,:,jp_sal,Kmm), clinfo1=' sss - : ', mask1=tmask, kdim=1 ) CALL prt_ctl(tab2d_1=sss_m , clinfo1=' sss - : ', mask1=tmask )
CALL prt_ctl(tab2d_1=utau , clinfo1=' utau - : ', mask1=umask, & CALL prt_ctl(tab2d_1=utau , clinfo1=' utau - : ', mask1=tmask, &
& tab2d_2=vtau , clinfo2=' vtau - : ', mask2=vmask ) & tab2d_2=vtau , clinfo2=' vtau - : ', mask2=tmask )
ENDIF ENDIF
IF( kt == nitend ) CALL sbc_final ! Close down surface module if necessary IF( kt == nitend ) CALL sbc_final ! Close down surface module if necessary
......
src/OCE/SBC/sbcrnf.F90
View file @ 3efbcb17
...@@ -83,9 +83,9 @@ CONTAINS ...@@ -83,9 +83,9 @@ CONTAINS
!!---------------------------------------------------------------------- !!----------------------------------------------------------------------
!! *** ROUTINE sbc_rnf_alloc *** !! *** ROUTINE sbc_rnf_alloc ***
!!---------------------------------------------------------------------- !!----------------------------------------------------------------------
ALLOCATE( rnfmsk(jpi,jpj) , rnfmsk_z(jpk) , & ALLOCATE( rnfmsk(jpi,jpj) , rnfmsk_z(jpk) , &
& h_rnf (jpi,jpj) , nk_rnf (jpi,jpj) , & & h_rnf (jpi,jpj) , nk_rnf (jpi,jpj) , &
& rnf_tsc_b(jpi,jpj,jpts) , rnf_tsc (jpi,jpj,jpts) , STAT=sbc_rnf_alloc ) & rnf_tsc_b(A2D(0),jpts) , rnf_tsc (A2D(0),jpts) , STAT=sbc_rnf_alloc )
! !
CALL mpp_sum ( 'sbcrnf', sbc_rnf_alloc ) CALL mpp_sum ( 'sbcrnf', sbc_rnf_alloc )
IF( sbc_rnf_alloc > 0 ) CALL ctl_warn('sbc_rnf_alloc: allocation of arrays failed') IF( sbc_rnf_alloc > 0 ) CALL ctl_warn('sbc_rnf_alloc: allocation of arrays failed')
...@@ -109,8 +109,6 @@ CONTAINS ...@@ -109,8 +109,6 @@ CONTAINS
INTEGER :: ji, jj ! dummy loop indices INTEGER :: ji, jj ! dummy loop indices
INTEGER :: z_err = 0 ! dummy integer for error handling INTEGER :: z_err = 0 ! dummy integer for error handling
!!---------------------------------------------------------------------- !!----------------------------------------------------------------------
REAL(wp), DIMENSION(jpi,jpj) :: ztfrz ! freezing point used for temperature correction
!
! !
! !-------------------! ! !-------------------!
! ! Update runoff ! ! ! Update runoff !
...@@ -118,8 +116,8 @@ CONTAINS ...@@ -118,8 +116,8 @@ CONTAINS
! !
! !
IF( .NOT. l_rnfcpl ) THEN IF( .NOT. l_rnfcpl ) THEN
CALL fld_read ( kt, nn_fsbc, sf_rnf ) ! Read Runoffs data and provide it at kt ( runoffs + iceberg ) CALL fld_read ( kt, nn_fsbc, sf_rnf ) ! Read Runoffs data and provide it at kt ( runoffs + iceberg )
IF( ln_rnf_icb ) CALL fld_read ( kt, nn_fsbc, sf_i_rnf ) ! idem for iceberg flux if required IF( ln_rnf_icb ) CALL fld_read ( kt, nn_fsbc, sf_i_rnf ) ! idem for iceberg flux if required
ENDIF ENDIF
IF( ln_rnf_tem ) CALL fld_read ( kt, nn_fsbc, sf_t_rnf ) ! idem for runoffs temperature if required IF( ln_rnf_tem ) CALL fld_read ( kt, nn_fsbc, sf_t_rnf ) ! idem for runoffs temperature if required
IF( ln_rnf_sal ) CALL fld_read ( kt, nn_fsbc, sf_s_rnf ) ! idem for runoffs salinity if required IF( ln_rnf_sal ) CALL fld_read ( kt, nn_fsbc, sf_s_rnf ) ! idem for runoffs salinity if required
...@@ -127,33 +125,32 @@ CONTAINS ...@@ -127,33 +125,32 @@ CONTAINS
IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN
! !
IF( .NOT. l_rnfcpl ) THEN IF( .NOT. l_rnfcpl ) THEN
rnf(:,:) = rn_rfact * ( sf_rnf(1)%fnow(:,:,1) ) * tmask(:,:,1) ! updated runoff value at time step kt rnf(A2D(0)) = rn_rfact * ( sf_rnf(1)%fnow(:,:,1) ) * smask0(:,:) ! updated runoff value at time step kt
IF( ln_rnf_icb ) THEN IF( ln_rnf_icb ) THEN
fwficb(:,:) = rn_rfact * ( sf_i_rnf(1)%fnow(:,:,1) ) * tmask(:,:,1) ! updated runoff value at time step kt fwficb(A2D(0)) = rn_rfact * ( sf_i_rnf(1)%fnow(:,:,1) ) * smask0(:,:) ! updated runoff value at time step kt
rnf(:,:) = rnf(:,:) + fwficb(:,:) rnf(A2D(0)) = rnf(A2D(0)) + fwficb(A2D(0))
qns(:,:) = qns(:,:) - fwficb(:,:) * rLfus qns(:,:) = qns(:,:) - fwficb(A2D(0)) * rLfus
!!qns_tot(:,:) = qns_tot(:,:) - fwficb(:,:) * rLfus !!qns_tot(:,:) = qns_tot(:,:) - fwficb(:,:) * rLfus
!!qns_oce(:,:) = qns_oce(:,:) - fwficb(:,:) * rLfus !!qns_oce(:,:) = qns_oce(:,:) - fwficb(:,:) * rLfus
CALL iom_put( 'iceberg_cea' , fwficb(:,:) ) ! output iceberg flux CALL iom_put( 'iceberg_cea' , fwficb(A2D(0)) ) ! output iceberg flux
CALL iom_put( 'hflx_icb_cea' , -fwficb(:,:) * rLfus ) ! output Heat Flux into Sea Water due to Iceberg Thermodynamics --> CALL iom_put( 'hflx_icb_cea' , -fwficb(A2D(0)) * rLfus ) ! output Heat Flux into Sea Water due to Iceberg Thermodynamics -->
ENDIF ENDIF
ENDIF ENDIF
! !
! ! set temperature & salinity content of runoffs ! ! set temperature & salinity content of runoffs
IF( ln_rnf_tem ) THEN ! use runoffs temperature data IF( ln_rnf_tem ) THEN ! use runoffs temperature data
rnf_tsc(:,:,jp_tem) = ( sf_t_rnf(1)%fnow(:,:,1) ) * rnf(:,:) * r1_rho0 rnf_tsc(:,:,jp_tem) = ( sf_t_rnf(1)%fnow(:,:,1) ) * rnf(A2D(0)) * r1_rho0
CALL eos_fzp( sss_m(:,:), ztfrz(:,:) )
WHERE( sf_t_rnf(1)%fnow(:,:,1) == -999._wp ) ! if missing data value use SST as runoffs temperature WHERE( sf_t_rnf(1)%fnow(:,:,1) == -999._wp ) ! if missing data value use SST as runoffs temperature
rnf_tsc(:,:,jp_tem) = sst_m(:,:) * rnf(:,:) * r1_rho0 rnf_tsc(:,:,jp_tem) = sst_m(A2D(0)) * rnf(A2D(0)) * r1_rho0
END WHERE END WHERE
ELSE ! use SST as runoffs temperature ELSE ! use SST as runoffs temperature
!CEOD River is fresh water so must at least be 0 unless we consider ice !CEOD River is fresh water so must at least be 0 unless we consider ice
rnf_tsc(:,:,jp_tem) = MAX( sst_m(:,:), 0.0_wp ) * rnf(:,:) * r1_rho0 rnf_tsc(:,:,jp_tem) = MAX( sst_m(A2D(0)), 0.0_wp ) * rnf(A2D(0)) * r1_rho0
ENDIF ENDIF
! ! use runoffs salinity data ! ! use runoffs salinity data
IF( ln_rnf_sal ) rnf_tsc(:,:,jp_sal) = ( sf_s_rnf(1)%fnow(:,:,1) ) * rnf(:,:) * r1_rho0 IF( ln_rnf_sal ) rnf_tsc(:,:,jp_sal) = ( sf_s_rnf(1)%fnow(:,:,1) ) * rnf(A2D(0)) * r1_rho0
! ! else use S=0 for runoffs (done one for all in the init) ! ! else use S=0 for runoffs (done one for all in the init)
CALL iom_put( 'runoffs' , rnf(:,:) ) ! output runoff mass flux CALL iom_put( 'runoffs' , rnf(A2D(0)) ) ! output runoff mass flux
IF( iom_use('hflx_rnf_cea') ) CALL iom_put( 'hflx_rnf_cea', rnf_tsc(:,:,jp_tem) * rho0 * rcp ) ! output runoff sensible heat (W/m2) IF( iom_use('hflx_rnf_cea') ) CALL iom_put( 'hflx_rnf_cea', rnf_tsc(:,:,jp_tem) * rho0 * rcp ) ! output runoff sensible heat (W/m2)
IF( iom_use('sflx_rnf_cea') ) CALL iom_put( 'sflx_rnf_cea', rnf_tsc(:,:,jp_sal) * rho0 ) ! output runoff salt flux (g/m2/s) IF( iom_use('sflx_rnf_cea') ) CALL iom_put( 'sflx_rnf_cea', rnf_tsc(:,:,jp_sal) * rho0 ) ! output runoff salt flux (g/m2/s)
ENDIF ENDIF
...@@ -168,13 +165,15 @@ CONTAINS ...@@ -168,13 +165,15 @@ CONTAINS
CALL iom_get( numror, jpdom_auto, 'rnf_sc_b', rnf_tsc_b(:,:,jp_sal) ) ! before salinity content of runoff CALL iom_get( numror, jpdom_auto, 'rnf_sc_b', rnf_tsc_b(:,:,jp_sal) ) ! before salinity content of runoff
ELSE !* no restart: set from nit000 values ELSE !* no restart: set from nit000 values
IF(lwp) WRITE(numout,*) ' nit000-1 runoff forcing fields set to nit000' IF(lwp) WRITE(numout,*) ' nit000-1 runoff forcing fields set to nit000'
rnf_b (:,: ) = rnf (:,: ) CALL lbc_lnk( 'sbcrnf', rnf, 'T', 1.0_wp )
rnf_tsc_b(:,:,:) = rnf_tsc(:,:,:) rnf_b (:,:) = rnf (:,:)
rnf_tsc_b(:,:,:) = rnf_tsc(:,:,:)
ENDIF ENDIF
ENDIF ENDIF
! ! ---------------------------------------- ! ! ! ---------------------------------------- !
IF( lrst_oce ) THEN ! Write in the ocean restart file ! IF( lrst_oce ) THEN ! Write in the ocean restart file !
! ! ---------------------------------------- ! ! ! ---------------------------------------- !
!
IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) 'sbcrnf : runoff forcing fields written in ocean restart file ', & IF(lwp) WRITE(numout,*) 'sbcrnf : runoff forcing fields written in ocean restart file ', &
& 'at it= ', kt,' date= ', ndastp & 'at it= ', kt,' date= ', ndastp
...@@ -323,8 +322,8 @@ CONTAINS ...@@ -323,8 +322,8 @@ CONTAINS
IF( ierror > 0 ) THEN IF( ierror > 0 ) THEN
CALL ctl_stop( 'sbc_rnf_init: unable to allocate sf_rnf structure' ) ; RETURN CALL ctl_stop( 'sbc_rnf_init: unable to allocate sf_rnf structure' ) ; RETURN
ENDIF ENDIF
ALLOCATE( sf_rnf(1)%fnow(jpi,jpj,1) ) ALLOCATE( sf_rnf(1)%fnow(A2D(0),1) )
IF( sn_rnf%ln_tint ) ALLOCATE( sf_rnf(1)%fdta(jpi,jpj,1,2) ) IF( sn_rnf%ln_tint ) ALLOCATE( sf_rnf(1)%fdta(A2D(0),1,2) )
CALL fld_fill( sf_rnf, (/ sn_rnf /), cn_dir, 'sbc_rnf_init', 'read runoffs data', 'namsbc_rnf', no_print ) CALL fld_fill( sf_rnf, (/ sn_rnf /), cn_dir, 'sbc_rnf_init', 'read runoffs data', 'namsbc_rnf', no_print )
! !
IF( ln_rnf_icb ) THEN ! Create (if required) sf_i_rnf structure IF( ln_rnf_icb ) THEN ! Create (if required) sf_i_rnf structure
...@@ -334,11 +333,13 @@ CONTAINS ...@@ -334,11 +333,13 @@ CONTAINS
IF( ierror > 0 ) THEN IF( ierror > 0 ) THEN
CALL ctl_stop( 'sbc_rnf_init: unable to allocate sf_i_rnf structure' ) ; RETURN CALL ctl_stop( 'sbc_rnf_init: unable to allocate sf_i_rnf structure' ) ; RETURN
ENDIF ENDIF
ALLOCATE( sf_i_rnf(1)%fnow(jpi,jpj,1) ) ALLOCATE( sf_i_rnf(1)%fnow(A2D(0),1) )
IF( sn_i_rnf%ln_tint ) ALLOCATE( sf_i_rnf(1)%fdta(jpi,jpj,1,2) ) IF( sn_i_rnf%ln_tint ) ALLOCATE( sf_i_rnf(1)%fdta(A2D(0),1,2) )
CALL fld_fill (sf_i_rnf, (/ sn_i_rnf /), cn_dir, 'sbc_rnf_init', 'read iceberg flux data', 'namsbc_rnf' ) CALL fld_fill (sf_i_rnf, (/ sn_i_rnf /), cn_dir, 'sbc_rnf_init', 'read iceberg flux data', 'namsbc_rnf' )
ELSE ELSE
fwficb(:,:) = 0._wp DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
fwficb(ji,jj) = 0._wp
END_2D
ENDIF ENDIF
ENDIF ENDIF
...@@ -350,8 +351,8 @@ CONTAINS ...@@ -350,8 +351,8 @@ CONTAINS
IF( ierror > 0 ) THEN IF( ierror > 0 ) THEN
CALL ctl_stop( 'sbc_rnf_init: unable to allocate sf_t_rnf structure' ) ; RETURN CALL ctl_stop( 'sbc_rnf_init: unable to allocate sf_t_rnf structure' ) ; RETURN
ENDIF ENDIF
ALLOCATE( sf_t_rnf(1)%fnow(jpi,jpj,1) ) ALLOCATE( sf_t_rnf(1)%fnow(A2D(0),1) )
IF( sn_t_rnf%ln_tint ) ALLOCATE( sf_t_rnf(1)%fdta(jpi,jpj,1,2) ) IF( sn_t_rnf%ln_tint ) ALLOCATE( sf_t_rnf(1)%fdta(A2D(0),1,2) )
CALL fld_fill (sf_t_rnf, (/ sn_t_rnf /), cn_dir, 'sbc_rnf_init', 'read runoff temperature data', 'namsbc_rnf', no_print ) CALL fld_fill (sf_t_rnf, (/ sn_t_rnf /), cn_dir, 'sbc_rnf_init', 'read runoff temperature data', 'namsbc_rnf', no_print )
ENDIF ENDIF
! !
...@@ -362,8 +363,8 @@ CONTAINS ...@@ -362,8 +363,8 @@ CONTAINS
IF( ierror > 0 ) THEN IF( ierror > 0 ) THEN
CALL ctl_stop( 'sbc_rnf_init: unable to allocate sf_s_rnf structure' ) ; RETURN CALL ctl_stop( 'sbc_rnf_init: unable to allocate sf_s_rnf structure' ) ; RETURN
ENDIF ENDIF
ALLOCATE( sf_s_rnf(1)%fnow(jpi,jpj,1) ) ALLOCATE( sf_s_rnf(1)%fnow(A2D(0),1) )
IF( sn_s_rnf%ln_tint ) ALLOCATE( sf_s_rnf(1)%fdta(jpi,jpj,1,2) ) IF( sn_s_rnf%ln_tint ) ALLOCATE( sf_s_rnf(1)%fdta(A2D(0),1,2) )
CALL fld_fill (sf_s_rnf, (/ sn_s_rnf /), cn_dir, 'sbc_rnf_init', 'read runoff salinity data', 'namsbc_rnf', no_print ) CALL fld_fill (sf_s_rnf, (/ sn_s_rnf /), cn_dir, 'sbc_rnf_init', 'read runoff salinity data', 'namsbc_rnf', no_print )
ENDIF ENDIF
! !
...@@ -378,8 +379,7 @@ CONTAINS ...@@ -378,8 +379,7 @@ CONTAINS
CALL iom_get ( inum, jpdom_global, sn_dep_rnf%clvar, h_rnf, kfill = jpfillcopy ) ! read the river mouth. no 0 on halos! CALL iom_get ( inum, jpdom_global, sn_dep_rnf%clvar, h_rnf, kfill = jpfillcopy ) ! read the river mouth. no 0 on halos!
CALL iom_close( inum ) ! close file CALL iom_close( inum ) ! close file
! !
nk_rnf(:,:) = 0 ! set the number of level over which river runoffs are applied DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) ! set the number of level over which river runoffs are applied
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
IF( h_rnf(ji,jj) > 0._wp ) THEN IF( h_rnf(ji,jj) > 0._wp ) THEN
jk = 2 jk = 2
DO WHILE ( jk < mbkt(ji,jj) .AND. gdept_0(ji,jj,jk) < h_rnf(ji,jj) ) ; jk = jk + 1 DO WHILE ( jk < mbkt(ji,jj) .AND. gdept_0(ji,jj,jk) < h_rnf(ji,jj) ) ; jk = jk + 1
...@@ -392,7 +392,8 @@ CONTAINS ...@@ -392,7 +392,8 @@ CONTAINS
WRITE(999,*) 'ji, jj, h_rnf(ji,jj) :', ji, jj, h_rnf(ji,jj) WRITE(999,*) 'ji, jj, h_rnf(ji,jj) :', ji, jj, h_rnf(ji,jj)
ENDIF ENDIF
END_2D END_2D
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) ! set the associated depth !
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) ! set the associated depth
h_rnf(ji,jj) = 0._wp h_rnf(ji,jj) = 0._wp
DO jk = 1, nk_rnf(ji,jj) DO jk = 1, nk_rnf(ji,jj)
h_rnf(ji,jj) = h_rnf(ji,jj) + e3t(ji,jj,jk,Kmm) h_rnf(ji,jj) = h_rnf(ji,jj) + e3t(ji,jj,jk,Kmm)
...@@ -407,7 +408,7 @@ CONTAINS ...@@ -407,7 +408,7 @@ CONTAINS
IF(lwp) WRITE(numout,*) ' depth over which runoffs is spread rn_dep_max = ', rn_dep_max IF(lwp) WRITE(numout,*) ' depth over which runoffs is spread rn_dep_max = ', rn_dep_max
IF(lwp) WRITE(numout,*) ' create (=1) a runoff depth file or not (=0) nn_rnf_depth_file = ', nn_rnf_depth_file IF(lwp) WRITE(numout,*) ' create (=1) a runoff depth file or not (=0) nn_rnf_depth_file = ', nn_rnf_depth_file
CALL iom_open( TRIM( sn_rnf%clname ), inum ) ! open runoff file CALL iom_open( TRIM( sn_rnf%clname ), inum ) ! open runoff file
nbrec = iom_getszuld( inum ) nbrec = iom_getszuld( inum )
zrnfcl(:,:,1) = 0._wp ! init the max to 0. in 1 zrnfcl(:,:,1) = 0._wp ! init the max to 0. in 1
DO jm = 1, nbrec DO jm = 1, nbrec
...@@ -416,21 +417,20 @@ CONTAINS ...@@ -416,21 +417,20 @@ CONTAINS
END DO END DO
CALL iom_close( inum ) CALL iom_close( inum )
! !
h_rnf(:,:) = 1. zacoef = rn_dep_max / rn_rnf_max ! coef of linear relation between runoff and its depth (150m for max of runoff)
!
zacoef = rn_dep_max / rn_rnf_max ! coef of linear relation between runoff and its depth (150m for max of runoff)
! !
WHERE( zrnfcl(:,:,1) > 0._wp ) h_rnf(:,:) = zacoef * zrnfcl(:,:,1) ! compute depth for all runoffs WHERE( zrnfcl(:,:,1) > 0._wp ) ; h_rnf(:,:) = zacoef * zrnfcl(:,:,1) ! compute depth for all runoffs
ELSEWHERE ; h_rnf(:,:) = 1._wp
ENDWHERE
! !
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) ! take in account min depth of ocean rn_hmin DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) ! take in account min depth of ocean rn_hmin
IF( zrnfcl(ji,jj,1) > 0._wp ) THEN IF( zrnfcl(ji,jj,1) > 0._wp ) THEN
jk = mbkt(ji,jj) jk = mbkt(ji,jj)
h_rnf(ji,jj) = MIN( h_rnf(ji,jj), gdept_0(ji,jj,jk ) ) h_rnf(ji,jj) = MIN( h_rnf(ji,jj), gdept_0(ji,jj,jk ) )
ENDIF ENDIF
END_2D END_2D
! !
nk_rnf(:,:) = 0 ! number of levels on which runoffs are distributed DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) ! number of levels on which runoffs are distributed
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
IF( zrnfcl(ji,jj,1) > 0._wp ) THEN IF( zrnfcl(ji,jj,1) > 0._wp ) THEN
jk = 2 jk = 2
DO WHILE ( jk < mbkt(ji,jj) .AND. gdept_0(ji,jj,jk) < h_rnf(ji,jj) ) ; jk = jk + 1 DO WHILE ( jk < mbkt(ji,jj) .AND. gdept_0(ji,jj,jk) < h_rnf(ji,jj) ) ; jk = jk + 1
...@@ -441,27 +441,33 @@ CONTAINS ...@@ -441,27 +441,33 @@ CONTAINS
ENDIF ENDIF
END_2D END_2D
! !
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) ! set the associated depth DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) ! set the associated depth
h_rnf(ji,jj) = 0._wp h_rnf(ji,jj) = 0._wp
DO jk = 1, nk_rnf(ji,jj) DO jk = 1, nk_rnf(ji,jj)
h_rnf(ji,jj) = h_rnf(ji,jj) + e3t(ji,jj,jk,Kmm) h_rnf(ji,jj) = h_rnf(ji,jj) + e3t(ji,jj,jk,Kmm)
END DO END DO
END_2D END_2D
! !
IF( nn_rnf_depth_file == 1 ) THEN ! save output nb levels for runoff IF( nn_rnf_depth_file == 1 ) THEN ! save output nb levels for runoff
IF(lwp) WRITE(numout,*) ' ==>>> create runoff depht file' IF(lwp) WRITE(numout,*) ' ==>>> create runoff depht file'
CALL iom_open ( TRIM( sn_dep_rnf%clname ), inum, ldwrt = .TRUE. ) CALL iom_open ( TRIM( sn_dep_rnf%clname ), inum, ldwrt = .TRUE. )
CALL iom_rstput( 0, 0, inum, 'rodepth', h_rnf ) CALL iom_rstput( 0, 0, inum, 'rodepth', h_rnf )
CALL iom_close ( inum ) CALL iom_close ( inum )
ENDIF ENDIF
ELSE ! runoffs applied at the surface ELSE ! runoffs applied at the surface
nk_rnf(:,:) = 1 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
h_rnf (:,:) = e3t(:,:,1,Kmm) nk_rnf(ji,jj) = 1
h_rnf (ji,jj) = e3t(ji,jj,1,Kmm)
END_2D
ENDIF ENDIF
! !
rnf(:,:) = 0._wp ! runoff initialisation DO_2D( 0, 0, 0, 0 )
rnf_tsc(:,:,:) = 0._wp ! runoffs temperature & salinty contents initilisation rnf(ji,jj) = 0._wp ! runoff initialisation
! END_2D
DO_3D( 0, 0, 0, 0, 1, jpts )
rnf_tsc(ji,jj,jk) = 0._wp ! runoffs temperature & salinty contents initilisation
END_3D
!
! ! ======================== ! ! ========================
! ! River mouth vicinity ! ! River mouth vicinity
! ! ======================== ! ! ========================
...@@ -493,11 +499,13 @@ CONTAINS ...@@ -493,11 +499,13 @@ CONTAINS
ELSE ! No treatment at river mouths ELSE ! No treatment at river mouths
IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) ' ==>>> No specific treatment at river mouths' IF(lwp) WRITE(numout,*) ' ==>>> No specific treatment at river mouths'
rnfmsk (:,:) = 0._wp DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
rnfmsk(ji,jj) = 0._wp
END_2D
rnfmsk_z(:) = 0._wp rnfmsk_z(:) = 0._wp
nkrnf = 0 nkrnf = 0
ENDIF ENDIF
! !
END SUBROUTINE sbc_rnf_init END SUBROUTINE sbc_rnf_init
......
src/OCE/SBC/sbcssr.F90
View file @ 3efbcb17
...@@ -97,8 +97,8 @@ CONTAINS ...@@ -97,8 +97,8 @@ CONTAINS
erp(:,:) = 0._wp erp(:,:) = 0._wp
! !
IF( nn_sstr == 1 ) THEN !* Temperature restoring term IF( nn_sstr == 1 ) THEN !* Temperature restoring term
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) DO_2D( 0, 0, 0, 0 )
zqrp = rn_dqdt * ( sst_m(ji,jj) - sf_sst(1)%fnow(ji,jj,1) ) * tmask(ji,jj,1) zqrp = rn_dqdt * ( sst_m(ji,jj) - sf_sst(1)%fnow(ji,jj,1) ) * smask0(ji,jj)
qns(ji,jj) = qns(ji,jj) + zqrp qns(ji,jj) = qns(ji,jj) + zqrp
qrp(ji,jj) = zqrp qrp(ji,jj) = zqrp
END_2D END_2D
...@@ -107,7 +107,7 @@ CONTAINS ...@@ -107,7 +107,7 @@ CONTAINS
IF( nn_sssr /= 0 .AND. nn_sssr_ice /= 1 ) THEN IF( nn_sssr /= 0 .AND. nn_sssr_ice /= 1 ) THEN
! use fraction of ice ( fr_i ) to adjust relaxation under ice if nn_sssr_ice .ne. 1 ! use fraction of ice ( fr_i ) to adjust relaxation under ice if nn_sssr_ice .ne. 1
! n.b. coefice is initialised and fixed to 1._wp if nn_sssr_ice = 1 ! n.b. coefice is initialised and fixed to 1._wp if nn_sssr_ice = 1
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) DO_2D( 0, 0, 0, 0 )
SELECT CASE ( nn_sssr_ice ) SELECT CASE ( nn_sssr_ice )
CASE ( 0 ) ; coefice(ji,jj) = 1._wp - fr_i(ji,jj) ! no/reduced damping under ice CASE ( 0 ) ; coefice(ji,jj) = 1._wp - fr_i(ji,jj) ! no/reduced damping under ice
CASE DEFAULT ; coefice(ji,jj) = 1._wp + ( nn_sssr_ice - 1 ) * fr_i(ji,jj) ! reinforced damping (x nn_sssr_ice) under ice ) CASE DEFAULT ; coefice(ji,jj) = 1._wp + ( nn_sssr_ice - 1 ) * fr_i(ji,jj) ! reinforced damping (x nn_sssr_ice) under ice )
...@@ -117,10 +117,10 @@ CONTAINS ...@@ -117,10 +117,10 @@ CONTAINS
! !
IF( nn_sssr == 1 ) THEN !* Salinity damping term (salt flux only (sfx)) IF( nn_sssr == 1 ) THEN !* Salinity damping term (salt flux only (sfx))
zsrp = rn_deds / rday ! from [mm/day] to [kg/m2/s] zsrp = rn_deds / rday ! from [mm/day] to [kg/m2/s]
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) DO_2D( 0, 0, 0, 0 )
zerp = zsrp * ( 1. - 2.*rnfmsk(ji,jj) ) & ! No damping in vicinity of river mouths zerp = zsrp * ( 1. - 2.*rnfmsk(ji,jj) ) & ! No damping in vicinity of river mouths
& * coefice(ji,jj) & ! Optional control of damping under sea-ice & * coefice(ji,jj) & ! Optional control of damping under sea-ice
& * ( sss_m(ji,jj) - sf_sss(1)%fnow(ji,jj,1) ) * tmask(ji,jj,1) & * ( sss_m(ji,jj) - sf_sss(1)%fnow(ji,jj,1) ) * smask0(ji,jj)
sfx(ji,jj) = sfx(ji,jj) + zerp ! salt flux sfx(ji,jj) = sfx(ji,jj) + zerp ! salt flux
erp(ji,jj) = zerp / MAX( sss_m(ji,jj), 1.e-20 ) ! converted into an equivalent volume flux (diagnostic only) erp(ji,jj) = zerp / MAX( sss_m(ji,jj), 1.e-20 ) ! converted into an equivalent volume flux (diagnostic only)
END_2D END_2D
...@@ -128,21 +128,21 @@ CONTAINS ...@@ -128,21 +128,21 @@ CONTAINS
ELSEIF( nn_sssr == 2 ) THEN !* Salinity damping term (volume flux (emp) and associated heat flux (qns) ELSEIF( nn_sssr == 2 ) THEN !* Salinity damping term (volume flux (emp) and associated heat flux (qns)
zsrp = rn_deds / rday ! from [mm/day] to [kg/m2/s] zsrp = rn_deds / rday ! from [mm/day] to [kg/m2/s]
zerp_bnd = rn_sssr_bnd / rday ! - - zerp_bnd = rn_sssr_bnd / rday ! - -
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) DO_2D( 0, 0, 0, 0 )
zerp = zsrp * ( 1. - 2.*rnfmsk(ji,jj) ) & ! No damping in vicinity of river mouths zerp = zsrp * ( 1. - 2.*rnfmsk(ji,jj) ) & ! No damping in vicinity of river mouths
& * coefice(ji,jj) & ! Optional control of damping under sea-ice & * coefice(ji,jj) & ! Optional control of damping under sea-ice
& * ( sss_m(ji,jj) - sf_sss(1)%fnow(ji,jj,1) ) & & * ( sss_m(ji,jj) - sf_sss(1)%fnow(ji,jj,1) ) &
& / MAX( sss_m(ji,jj), 1.e-20 ) * tmask(ji,jj,1) & / MAX( sss_m(ji,jj), 1.e-20 ) * smask0(ji,jj)
IF( ln_sssr_bnd ) zerp = SIGN( 1.0_wp, zerp ) * MIN( zerp_bnd, ABS(zerp) ) IF( ln_sssr_bnd ) zerp = SIGN( 1.0_wp, zerp ) * MIN( zerp_bnd, ABS(zerp) )
emp(ji,jj) = emp (ji,jj) + zerp
qns(ji,jj) = qns(ji,jj) - zerp * rcp * sst_m(ji,jj)
erp(ji,jj) = zerp erp(ji,jj) = zerp
qrp(ji,jj) = qrp(ji,jj) - zerp * rcp * sst_m(ji,jj) emp(ji,jj) = emp(ji,jj) + erp(ji,jj)
qns(ji,jj) = qns(ji,jj) - erp(ji,jj) * rcp * sst_m(ji,jj)
qrp(ji,jj) = qrp(ji,jj) - erp(ji,jj) * rcp * sst_m(ji,jj)
END_2D END_2D
ENDIF ENDIF
! outputs ! outputs
CALL iom_put( 'hflx_ssr_cea', qrp(:,:) ) CALL iom_put( 'hflx_ssr_cea', qrp(:,:) )
IF( nn_sssr == 1 ) CALL iom_put( 'sflx_ssr_cea', erp(:,:) * sss_m(:,:) ) IF( nn_sssr == 1 ) CALL iom_put( 'sflx_ssr_cea', erp(:,:) * sss_m(A2D(0)) )
IF( nn_sssr == 2 ) CALL iom_put( 'vflx_ssr_cea', -erp(:,:) ) IF( nn_sssr == 2 ) CALL iom_put( 'vflx_ssr_cea', -erp(:,:) )
! !
ENDIF ENDIF
...@@ -207,12 +207,12 @@ CONTAINS ...@@ -207,12 +207,12 @@ CONTAINS
! !
ALLOCATE( sf_sst(1), STAT=ierror ) ALLOCATE( sf_sst(1), STAT=ierror )
IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'sbc_ssr: unable to allocate sf_sst structure' ) IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'sbc_ssr: unable to allocate sf_sst structure' )
ALLOCATE( sf_sst(1)%fnow(jpi,jpj,1), STAT=ierror ) ALLOCATE( sf_sst(1)%fnow(A2D(0),1), STAT=ierror )
IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'sbc_ssr: unable to allocate sf_sst now array' ) IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'sbc_ssr: unable to allocate sf_sst now array' )
! !
! fill sf_sst with sn_sst and control print ! fill sf_sst with sn_sst and control print
CALL fld_fill( sf_sst, (/ sn_sst /), cn_dir, 'sbc_ssr', 'SST restoring term toward SST data', 'namsbc_ssr', no_print ) CALL fld_fill( sf_sst, (/ sn_sst /), cn_dir, 'sbc_ssr', 'SST restoring term toward SST data', 'namsbc_ssr', no_print )
IF( sf_sst(1)%ln_tint ) ALLOCATE( sf_sst(1)%fdta(jpi,jpj,1,2), STAT=ierror ) IF( sf_sst(1)%ln_tint ) ALLOCATE( sf_sst(1)%fdta(A2D(0),1,2), STAT=ierror )
IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'sbc_ssr: unable to allocate sf_sst data array' ) IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'sbc_ssr: unable to allocate sf_sst data array' )
! !
ENDIF ENDIF
...@@ -221,12 +221,12 @@ CONTAINS ...@@ -221,12 +221,12 @@ CONTAINS
! !
ALLOCATE( sf_sss(1), STAT=ierror ) ALLOCATE( sf_sss(1), STAT=ierror )
IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'sbc_ssr: unable to allocate sf_sss structure' ) IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'sbc_ssr: unable to allocate sf_sss structure' )
ALLOCATE( sf_sss(1)%fnow(jpi,jpj,1), STAT=ierror ) ALLOCATE( sf_sss(1)%fnow(A2D(0),1), STAT=ierror )
IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'sbc_ssr: unable to allocate sf_sss now array' ) IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'sbc_ssr: unable to allocate sf_sss now array' )
! !
! fill sf_sss with sn_sss and control print ! fill sf_sss with sn_sss and control print
CALL fld_fill( sf_sss, (/ sn_sss /), cn_dir, 'sbc_ssr', 'SSS restoring term toward SSS data', 'namsbc_ssr', no_print ) CALL fld_fill( sf_sss, (/ sn_sss /), cn_dir, 'sbc_ssr', 'SSS restoring term toward SSS data', 'namsbc_ssr', no_print )
IF( sf_sss(1)%ln_tint ) ALLOCATE( sf_sss(1)%fdta(jpi,jpj,1,2), STAT=ierror ) IF( sf_sss(1)%ln_tint ) ALLOCATE( sf_sss(1)%fdta(A2D(0),1,2), STAT=ierror )
IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'sbc_ssr: unable to allocate sf_sss data array' ) IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'sbc_ssr: unable to allocate sf_sss data array' )
! !
ENDIF ENDIF
...@@ -244,7 +244,7 @@ CONTAINS ...@@ -244,7 +244,7 @@ CONTAINS
!!---------------------------------------------------------------------- !!----------------------------------------------------------------------
sbc_ssr_alloc = 0 ! set to zero if no array to be allocated sbc_ssr_alloc = 0 ! set to zero if no array to be allocated
IF( .NOT. ALLOCATED( erp ) ) THEN IF( .NOT. ALLOCATED( erp ) ) THEN
ALLOCATE( qrp(jpi,jpj), erp(jpi,jpj), coefice(jpi,jpj), STAT= sbc_ssr_alloc ) ALLOCATE( qrp(A2D(0)), erp(A2D(0)), coefice(A2D(0)), STAT= sbc_ssr_alloc )
! !
IF( lk_mpp ) CALL mpp_sum ( 'sbcssr', sbc_ssr_alloc ) IF( lk_mpp ) CALL mpp_sum ( 'sbcssr', sbc_ssr_alloc )
IF( sbc_ssr_alloc /= 0 ) CALL ctl_warn('sbc_ssr_alloc: failed to allocate arrays.') IF( sbc_ssr_alloc /= 0 ) CALL ctl_warn('sbc_ssr_alloc: failed to allocate arrays.')
......
src/OCE/SBC/sbcwave.F90
View file @ 3efbcb17
...@@ -115,14 +115,13 @@ CONTAINS ...@@ -115,14 +115,13 @@ CONTAINS
INTEGER :: jj, ji, jk ! dummy loop argument INTEGER :: jj, ji, jk ! dummy loop argument
INTEGER :: ik ! local integer INTEGER :: ik ! local integer
REAL(wp) :: ztransp, zfac, ztemp, zsp0, zsqrt, zbreiv16_w REAL(wp) :: ztransp, zfac, ztemp, zsp0, zsqrt, zbreiv16_w
REAL(wp) :: zdep_u, zdep_v, zkh_u, zkh_v, zda_u, zda_v, sdtrp REAL(wp) :: zdep_u, zdep_v, zkh_u, zkh_v, zda_u, zda_v, sdtrp, zInt_w0, zInt_w1
REAL(wp), DIMENSION(:,:) , ALLOCATABLE :: zk_t, zk_u, zk_v, zu0_sd, zv0_sd ! 2D workspace REAL(wp), DIMENSION(:,:) , ALLOCATABLE :: zk_t, zk_u, zk_v, zu0_sd, zv0_sd ! 2D workspace
REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ze3divh, zInt_w ! 3D workspace REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ze3divh ! 3D workspace
!!--------------------------------------------------------------------- !!---------------------------------------------------------------------
! !
ALLOCATE( ze3divh(jpi,jpj,jpkm1) ) ! jpkm1 -> avoid lbc_lnk on jpk that is not defined ALLOCATE( ze3divh(jpi,jpj,jpkm1) ) ! jpkm1 -> avoid lbc_lnk on jpk that is not defined
ALLOCATE( zInt_w(jpi,jpj,jpk) ) ALLOCATE( zk_t(A2D(1)), zk_u(A2D(0)), zk_v(A2D(0)), zu0_sd(A2D(1)), zv0_sd(A2D(1)) )
ALLOCATE( zk_t(jpi,jpj), zk_u(jpi,jpj), zk_v(jpi,jpj), zu0_sd(jpi,jpj), zv0_sd(jpi,jpj) )
zk_t (:,:) = 0._wp zk_t (:,:) = 0._wp
zk_u (:,:) = 0._wp zk_u (:,:) = 0._wp
zk_v (:,:) = 0._wp zk_v (:,:) = 0._wp
...@@ -138,7 +137,7 @@ CONTAINS ...@@ -138,7 +137,7 @@ CONTAINS
! sdtrp is the norm of Stokes transport ! sdtrp is the norm of Stokes transport
! !
zfac = 0.166666666667_wp zfac = 0.166666666667_wp
DO_2D( 1, 1, 1, 1 ) ! In the deep-water limit we have ke = ||ust0||/( 6 * ||transport|| ) DO_2D( 0, 1, 0, 1 ) ! In the deep-water limit we have ke = ||ust0||/( 6 * ||transport|| )
zsp0 = SQRT( ut0sd(ji,jj)*ut0sd(ji,jj) + vt0sd(ji,jj)*vt0sd(ji,jj) ) !<-- norm of Surface Stokes drift zsp0 = SQRT( ut0sd(ji,jj)*ut0sd(ji,jj) + vt0sd(ji,jj)*vt0sd(ji,jj) ) !<-- norm of Surface Stokes drift
tsd2d(ji,jj) = zsp0 tsd2d(ji,jj) = zsp0
IF( cpl_tusd .AND. cpl_tvsd ) THEN !stokes transport is provided in coupled mode IF( cpl_tusd .AND. cpl_tvsd ) THEN !stokes transport is provided in coupled mode
...@@ -150,35 +149,40 @@ CONTAINS ...@@ -150,35 +149,40 @@ CONTAINS
ENDIF ENDIF
zk_t (ji,jj) = zfac * zsp0 / MAX ( sdtrp, 0.0000001_wp ) !<-- ke = ||ust0||/( 6 * ||transport|| ) zk_t (ji,jj) = zfac * zsp0 / MAX ( sdtrp, 0.0000001_wp ) !<-- ke = ||ust0||/( 6 * ||transport|| )
END_2D END_2D
!# define zInt_w ze3divh !
DO_3D( 1, 1, 1, 1, 1, jpk ) ! Compute the primitive of Breivik 2016 function at W-points
zfac = - 2._wp * zk_t (ji,jj) * gdepw(ji,jj,jk,Kmm) !<-- zfac should be negative definite
ztemp = EXP ( zfac )
zsqrt = SQRT( -zfac )
zbreiv16_w = ztemp - SQRT(rpi)*zsqrt*ERFC(zsqrt) !Eq. 16 Breivik 2016
zInt_w(ji,jj,jk) = ztemp - 4._wp * zk_t (ji,jj) * gdepw(ji,jj,jk,Kmm) * zbreiv16_w
END_3D
!
DO jk = 1, jpkm1 DO jk = 1, jpkm1
zfac = 0.166666666667_wp zfac = 0.166666666667_wp
DO_2D( 1, 1, 1, 1 ) !++ Compute the FV Breivik 2016 function at T-points DO_2D( 0, 1, 0, 1 ) !++ Compute the FV Breivik 2016 function at T-points
! zInt at jk
zfac = - 2._wp * zk_t (ji,jj) * gdepw(ji,jj,jk,Kmm) !<-- zfac should be negative definite
ztemp = EXP ( zfac )
zsqrt = SQRT( -zfac )
zbreiv16_w = ztemp - SQRT(rpi)*zsqrt*ERFC(zsqrt) !Eq. 16 Breivik 2016
zInt_w0 = ztemp - 4._wp * zk_t (ji,jj) * gdepw(ji,jj,jk,Kmm) * zbreiv16_w
! zInt at jk+1
zfac = - 2._wp * zk_t (ji,jj) * gdepw(ji,jj,jk+1,Kmm) !<-- zfac should be negative definite
ztemp = EXP ( zfac )
zsqrt = SQRT( -zfac )
zbreiv16_w = ztemp - SQRT(rpi)*zsqrt*ERFC(zsqrt) !Eq. 16 Breivik 2016
zInt_w1 = ztemp - 4._wp * zk_t (ji,jj) * gdepw(ji,jj,jk+1,Kmm) * zbreiv16_w
!
!
zsp0 = zfac / MAX(zk_t (ji,jj),0.0000001_wp) zsp0 = zfac / MAX(zk_t (ji,jj),0.0000001_wp)
ztemp = zInt_w(ji,jj,jk) - zInt_w(ji,jj,jk+1) ztemp = zInt_w0 - zInt_w1
zu0_sd(ji,jj) = ut0sd(ji,jj) * zsp0 * ztemp * tmask(ji,jj,jk) zu0_sd(ji,jj) = ut0sd(ji,jj) * zsp0 * ztemp * tmask(ji,jj,jk)
zv0_sd(ji,jj) = vt0sd(ji,jj) * zsp0 * ztemp * tmask(ji,jj,jk) zv0_sd(ji,jj) = vt0sd(ji,jj) * zsp0 * ztemp * tmask(ji,jj,jk)
END_2D END_2D
DO_2D( 1, 0, 1, 0 ) ! ++ Interpolate at U/V points DO_2D( 0, 0, 0, 0 ) ! ++ Interpolate at U/V points
zfac = 1.0_wp / e3u(ji ,jj,jk,Kmm) zfac = 1.0_wp / e3u(ji ,jj,jk,Kmm)
usd(ji,jj,jk) = 0.5_wp * zfac * ( zu0_sd(ji,jj)+zu0_sd(ji+1,jj) ) * umask(ji,jj,jk) usd(ji,jj,jk) = 0.5_wp * zfac * ( zu0_sd(ji,jj)+zu0_sd(ji+1,jj) ) * umask(ji,jj,jk)
zfac = 1.0_wp / e3v(ji ,jj,jk,Kmm) zfac = 1.0_wp / e3v(ji ,jj,jk,Kmm)
vsd(ji,jj,jk) = 0.5_wp * zfac * ( zv0_sd(ji,jj)+zv0_sd(ji,jj+1) ) * vmask(ji,jj,jk) vsd(ji,jj,jk) = 0.5_wp * zfac * ( zv0_sd(ji,jj)+zv0_sd(ji,jj+1) ) * vmask(ji,jj,jk)
END_2D END_2D
ENDDO ENDDO
!# undef zInt_w !
!
ELSE ELSE
zfac = 2.0_wp * rpi / 16.0_wp zfac = 2.0_wp * rpi / 16.0_wp
DO_2D( 1, 1, 1, 1 ) DO_2D( 0, 1, 0, 1 )
! Stokes drift velocity estimated from Hs and Tmean ! Stokes drift velocity estimated from Hs and Tmean
ztransp = zfac * hsw(ji,jj)*hsw(ji,jj) / MAX( wmp(ji,jj), 0.0000001_wp ) ztransp = zfac * hsw(ji,jj)*hsw(ji,jj) / MAX( wmp(ji,jj), 0.0000001_wp )
! Stokes surface speed ! Stokes surface speed
...@@ -186,7 +190,7 @@ CONTAINS ...@@ -186,7 +190,7 @@ CONTAINS
! Wavenumber scale ! Wavenumber scale
zk_t(ji,jj) = ABS( tsd2d(ji,jj) ) / MAX( ABS( 5.97_wp*ztransp ), 0.0000001_wp ) zk_t(ji,jj) = ABS( tsd2d(ji,jj) ) / MAX( ABS( 5.97_wp*ztransp ), 0.0000001_wp )
END_2D END_2D
DO_2D( 1, 0, 1, 0 ) ! exp. wave number & Stokes drift velocity at u- & v-points DO_2D( 0, 0, 0, 0 ) ! exp. wave number & Stokes drift velocity at u- & v-points
zk_u(ji,jj) = 0.5_wp * ( zk_t(ji,jj) + zk_t(ji+1,jj) ) zk_u(ji,jj) = 0.5_wp * ( zk_t(ji,jj) + zk_t(ji+1,jj) )
zk_v(ji,jj) = 0.5_wp * ( zk_t(ji,jj) + zk_t(ji,jj+1) ) zk_v(ji,jj) = 0.5_wp * ( zk_t(ji,jj) + zk_t(ji,jj+1) )
! !
...@@ -248,7 +252,7 @@ CONTAINS ...@@ -248,7 +252,7 @@ CONTAINS
CALL iom_put( "vstokes", vsd ) CALL iom_put( "vstokes", vsd )
CALL iom_put( "wstokes", wsd ) CALL iom_put( "wstokes", wsd )
! ! ! !
DEALLOCATE( ze3divh, zInt_w ) DEALLOCATE( ze3divh )
DEALLOCATE( zk_t, zk_u, zk_v, zu0_sd, zv0_sd ) DEALLOCATE( zk_t, zk_u, zk_v, zu0_sd, zv0_sd )
! !
END SUBROUTINE sbc_stokes END SUBROUTINE sbc_stokes
...@@ -274,12 +278,12 @@ CONTAINS ...@@ -274,12 +278,12 @@ CONTAINS
! !
IF( ln_cdgw .AND. .NOT. cpl_wdrag ) THEN !== Neutral drag coefficient ==! IF( ln_cdgw .AND. .NOT. cpl_wdrag ) THEN !== Neutral drag coefficient ==!
CALL fld_read( kt, nn_fsbc, sf_cd ) ! read from external forcing CALL fld_read( kt, nn_fsbc, sf_cd ) ! read from external forcing
cdn_wave(:,:) = sf_cd(1)%fnow(:,:,1) * tmask(:,:,1) cdn_wave(:,:) = sf_cd(1)%fnow(:,:,1) * smask0(:,:)
ENDIF ENDIF
IF( ln_tauoc .AND. .NOT. cpl_wstrf ) THEN !== Wave induced stress ==! IF( ln_tauoc .AND. .NOT. cpl_wstrf ) THEN !== Wave induced stress ==!
CALL fld_read( kt, nn_fsbc, sf_tauoc ) ! read stress reduction factor due to wave from external forcing CALL fld_read( kt, nn_fsbc, sf_tauoc ) ! read stress reduction factor due to wave from external forcing
tauoc_wave(:,:) = sf_tauoc(1)%fnow(:,:,1) * tmask(:,:,1) tauoc_wave(:,:) = sf_tauoc(1)%fnow(:,:,1) * smask0(:,:)
ELSEIF ( ln_taw .AND. cpl_taw ) THEN ELSEIF ( ln_taw .AND. cpl_taw ) THEN
IF (kt < 1) THEN ! The first fields gave by OASIS have very high erroneous values .... IF (kt < 1) THEN ! The first fields gave by OASIS have very high erroneous values ....
twox(:,:)=0._wp twox(:,:)=0._wp
...@@ -315,7 +319,7 @@ CONTAINS ...@@ -315,7 +319,7 @@ CONTAINS
! coupling routines ! coupling routines
IF( ln_zdfswm .AND. .NOT. cpl_wnum ) THEN !==wavenumber==! IF( ln_zdfswm .AND. .NOT. cpl_wnum ) THEN !==wavenumber==!
CALL fld_read( kt, nn_fsbc, sf_wn ) ! read wave parameters from external forcing CALL fld_read( kt, nn_fsbc, sf_wn ) ! read wave parameters from external forcing
wnum(:,:) = sf_wn(1)%fnow(:,:,1) * tmask(:,:,1) wnum(:,:) = sf_wn(1)%fnow(:,:,1) * smask0(:,:)
ENDIF ENDIF
! !
...@@ -391,10 +395,10 @@ CONTAINS ...@@ -391,10 +395,10 @@ CONTAINS
! !== Allocate wave arrays ==! ! !== Allocate wave arrays ==!
ALLOCATE( ut0sd (jpi,jpj) , vt0sd (jpi,jpj) ) ALLOCATE( ut0sd (jpi,jpj) , vt0sd (jpi,jpj) )
ALLOCATE( hsw (jpi,jpj) , wmp (jpi,jpj) ) ALLOCATE( hsw (jpi,jpj) , wmp (jpi,jpj) )
ALLOCATE( wnum (jpi,jpj) )
ALLOCATE( tsd2d (jpi,jpj) , div_sd(jpi,jpj) , bhd_wave(jpi,jpj) ) ALLOCATE( tsd2d (jpi,jpj) , div_sd(jpi,jpj) , bhd_wave(jpi,jpj) )
ALLOCATE( usd (jpi,jpj,jpk), vsd (jpi,jpj,jpk), wsd (jpi,jpj,jpk) ) ALLOCATE( usd (jpi,jpj,jpk), vsd (jpi,jpj,jpk), wsd (jpi,jpj,jpk) )
ALLOCATE( tusd (jpi,jpj) , tvsd (jpi,jpj) , ZMX (jpi,jpj,jpk) ) ALLOCATE( tusd (jpi,jpj) , tvsd (jpi,jpj) )
ALLOCATE( wnum (A2D(0)) , ZMX (A2D(0),jpk) )
usd (:,:,:) = 0._wp usd (:,:,:) = 0._wp
vsd (:,:,:) = 0._wp vsd (:,:,:) = 0._wp
wsd (:,:,:) = 0._wp wsd (:,:,:) = 0._wp
...@@ -422,30 +426,30 @@ CONTAINS ...@@ -422,30 +426,30 @@ CONTAINS
ALLOCATE( sf_cd(1), STAT=ierror ) !* allocate and fill sf_wave with sn_cdg ALLOCATE( sf_cd(1), STAT=ierror ) !* allocate and fill sf_wave with sn_cdg
IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_wave structure' ) IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_wave structure' )
! !
ALLOCATE( sf_cd(1)%fnow(jpi,jpj,1) ) ALLOCATE( sf_cd(1)%fnow(A2D(0),1) )
IF( sn_cdg%ln_tint ) ALLOCATE( sf_cd(1)%fdta(jpi,jpj,1,2) ) IF( sn_cdg%ln_tint ) ALLOCATE( sf_cd(1)%fdta(A2D(0),1,2) )
CALL fld_fill( sf_cd, (/ sn_cdg /), cn_dir, 'sbc_wave_init', 'Wave module ', 'namsbc_wave' ) CALL fld_fill( sf_cd, (/ sn_cdg /), cn_dir, 'sbc_wave_init', 'Wave module ', 'namsbc_wave' )
ENDIF ENDIF
ALLOCATE( cdn_wave(jpi,jpj) ) ALLOCATE( cdn_wave(A2D(0)) )
cdn_wave(:,:) = 0._wp cdn_wave(:,:) = 0._wp
ENDIF ENDIF
IF( ln_charn ) THEN ! wave drag IF( ln_charn ) THEN ! wave drag
IF( .NOT. cpl_charn ) THEN IF( .NOT. cpl_charn ) THEN
CALL ctl_stop( 'STOP', 'Charnock based wind stress can be used in coupled mode only' ) CALL ctl_stop( 'STOP', 'Charnock based wind stress can be used in coupled mode only' )
ENDIF ENDIF
ALLOCATE( charn(jpi,jpj) ) ALLOCATE( charn(A2D(0)) )
charn(:,:) = 0._wp charn(:,:) = 0._wp
ENDIF ENDIF
IF( ln_taw ) THEN ! wind stress IF( ln_taw ) THEN ! wind stress
IF( .NOT. cpl_taw ) THEN IF( .NOT. cpl_taw ) THEN
CALL ctl_stop( 'STOP', 'wind stress from wave model can be used in coupled mode only, use ln_cdgw instead' ) CALL ctl_stop( 'STOP', 'wind stress from wave model can be used in coupled mode only, use ln_cdgw instead' )
ENDIF ENDIF
ALLOCATE( tawx(jpi,jpj) ) ALLOCATE( tawx(A2D(0)) )
ALLOCATE( tawy(jpi,jpj) ) ALLOCATE( tawy(A2D(0)) )
ALLOCATE( twox(jpi,jpj) ) ALLOCATE( twox(A2D(0)) )
ALLOCATE( twoy(jpi,jpj) ) ALLOCATE( twoy(A2D(0)) )
ALLOCATE( tauoc_wavex(jpi,jpj) ) ALLOCATE( tauoc_wavex(A2D(0)) )
ALLOCATE( tauoc_wavey(jpi,jpj) ) ALLOCATE( tauoc_wavey(A2D(0)) )
tawx(:,:) = 0._wp tawx(:,:) = 0._wp
tawy(:,:) = 0._wp tawy(:,:) = 0._wp
twox(:,:) = 0._wp twox(:,:) = 0._wp
...@@ -458,7 +462,7 @@ CONTAINS ...@@ -458,7 +462,7 @@ CONTAINS
IF( .NOT. cpl_phioc ) THEN IF( .NOT. cpl_phioc ) THEN
CALL ctl_stop( 'STOP', 'phioc can be used in coupled mode only' ) CALL ctl_stop( 'STOP', 'phioc can be used in coupled mode only' )
ENDIF ENDIF
ALLOCATE( phioc(jpi,jpj) ) ALLOCATE( phioc(A2D(0)) )
phioc(:,:) = 0._wp phioc(:,:) = 0._wp
ENDIF ENDIF
...@@ -467,11 +471,11 @@ CONTAINS ...@@ -467,11 +471,11 @@ CONTAINS
ALLOCATE( sf_tauoc(1), STAT=ierror ) !* allocate and fill sf_wave with sn_tauoc ALLOCATE( sf_tauoc(1), STAT=ierror ) !* allocate and fill sf_wave with sn_tauoc
IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_tauoc structure' ) IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_tauoc structure' )
! !
ALLOCATE( sf_tauoc(1)%fnow(jpi,jpj,1) ) ALLOCATE( sf_tauoc(1)%fnow(A2D(0),1) )
IF( sn_tauoc%ln_tint ) ALLOCATE( sf_tauoc(1)%fdta(jpi,jpj,1,2) ) IF( sn_tauoc%ln_tint ) ALLOCATE( sf_tauoc(1)%fdta(A2D(0),1,2) )
CALL fld_fill( sf_tauoc, (/ sn_tauoc /), cn_dir, 'sbc_wave_init', 'Wave module', 'namsbc_wave' ) CALL fld_fill( sf_tauoc, (/ sn_tauoc /), cn_dir, 'sbc_wave_init', 'Wave module', 'namsbc_wave' )
ENDIF ENDIF
ALLOCATE( tauoc_wave(jpi,jpj) ) ALLOCATE( tauoc_wave(A2D(0)) )
tauoc_wave(:,:) = 0._wp tauoc_wave(:,:) = 0._wp
ENDIF ENDIF
...@@ -518,8 +522,8 @@ CONTAINS ...@@ -518,8 +522,8 @@ CONTAINS
IF( .NOT. cpl_wnum ) THEN IF( .NOT. cpl_wnum ) THEN
ALLOCATE( sf_wn(1), STAT=ierror ) !* allocate and fill sf_wave with sn_wnum ALLOCATE( sf_wn(1), STAT=ierror ) !* allocate and fill sf_wave with sn_wnum
IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_wn structure' ) IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_wn structure' )
ALLOCATE( sf_wn(1)%fnow(jpi,jpj,1) ) ALLOCATE( sf_wn(1)%fnow(A2D(0),1) )
IF( sn_wnum%ln_tint ) ALLOCATE( sf_wn(1)%fdta(jpi,jpj,1,2) ) IF( sn_wnum%ln_tint ) ALLOCATE( sf_wn(1)%fdta(A2D(0),1,2) )
CALL fld_fill( sf_wn, (/ sn_wnum /), cn_dir, 'sbc_wave', 'Wave module', 'namsbc_wave' ) CALL fld_fill( sf_wn, (/ sn_wnum /), cn_dir, 'sbc_wave', 'Wave module', 'namsbc_wave' )
ENDIF ENDIF
! !
......
src/OCE/TRA/eosbn2.F90
View file @ 3efbcb17
...@@ -83,9 +83,11 @@ MODULE eosbn2 ...@@ -83,9 +83,11 @@ MODULE eosbn2
INTEGER , PARAMETER :: np_teos10 = -1 ! parameter for using TEOS10 INTEGER , PARAMETER :: np_teos10 = -1 ! parameter for using TEOS10
INTEGER , PARAMETER :: np_eos80 = 0 ! parameter for using EOS80 INTEGER , PARAMETER :: np_eos80 = 0 ! parameter for using EOS80
INTEGER , PARAMETER :: np_seos = 1 ! parameter for using Simplified Equation of state INTEGER , PARAMETER :: np_seos = 1 ! parameter for using Simplified Equation of state
! !!! simplified eos coefficients (default value: Vallis 2006) ! !!! simplified eos coefficients (default value: Vallis 2006)
REAL(wp), PUBLIC :: rn_T0 = 10._wp ! reference temperature
REAL(wp), PUBLIC :: rn_S0 = 35._wp ! reference salinity
REAL(wp), PUBLIC :: rn_a0 = 1.6550e-1_wp ! thermal expansion coeff. REAL(wp), PUBLIC :: rn_a0 = 1.6550e-1_wp ! thermal expansion coeff.
REAL(wp), PUBLIC :: rn_b0 = 7.6554e-1_wp ! saline expansion coeff. REAL(wp), PUBLIC :: rn_b0 = 7.6554e-1_wp ! saline expansion coeff.
REAL(wp) :: rn_lambda1 = 5.9520e-2_wp ! cabbeling coeff. in T^2 REAL(wp) :: rn_lambda1 = 5.9520e-2_wp ! cabbeling coeff. in T^2
...@@ -272,8 +274,8 @@ CONTAINS ...@@ -272,8 +274,8 @@ CONTAINS
CASE( np_seos ) !== simplified EOS ==! CASE( np_seos ) !== simplified EOS ==!
! !
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 ) DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )
zt = pts (ji,jj,jk,jp_tem,Knn) - 10._wp zt = pts (ji,jj,jk,jp_tem,Knn) - rn_T0
zs = pts (ji,jj,jk,jp_sal,Knn) - 35._wp zs = pts (ji,jj,jk,jp_sal,Knn) - rn_S0
zh = gdept(ji,jj,jk,Knn) zh = gdept(ji,jj,jk,Knn)
ztm = tmask(ji,jj,jk) ztm = tmask(ji,jj,jk)
! !
...@@ -391,8 +393,8 @@ CONTAINS ...@@ -391,8 +393,8 @@ CONTAINS
CASE( np_seos ) !== simplified EOS ==! CASE( np_seos ) !== simplified EOS ==!
! !
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 ) DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )
zt = pts (ji,jj,jk,jp_tem) - 10._wp zt = pts (ji,jj,jk,jp_tem) - rn_T0
zs = pts (ji,jj,jk,jp_sal) - 35._wp zs = pts (ji,jj,jk,jp_sal) - rn_S0
zh = pdep (ji,jj,jk) zh = pdep (ji,jj,jk)
ztm = tmask(ji,jj,jk) ztm = tmask(ji,jj,jk)
! !
...@@ -556,8 +558,8 @@ CONTAINS ...@@ -556,8 +558,8 @@ CONTAINS
CASE( np_seos ) !== simplified EOS ==! CASE( np_seos ) !== simplified EOS ==!
! !
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 ) DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )
zt = pts (ji,jj,jk,jp_tem) - 10._wp zt = pts (ji,jj,jk,jp_tem) - rn_T0
zs = pts (ji,jj,jk,jp_sal) - 35._wp zs = pts (ji,jj,jk,jp_sal) - rn_S0
zh = pdep (ji,jj,jk) zh = pdep (ji,jj,jk)
ztm = tmask(ji,jj,jk) ztm = tmask(ji,jj,jk)
! ! potential density referenced at the surface ! ! potential density referenced at the surface
...@@ -658,8 +660,8 @@ CONTAINS ...@@ -658,8 +660,8 @@ CONTAINS
! !
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
! !
zt = pts (ji,jj,jp_tem) - 10._wp zt = pts (ji,jj,jp_tem) - rn_T0
zs = pts (ji,jj,jp_sal) - 35._wp zs = pts (ji,jj,jp_sal) - rn_S0
zh = pdep (ji,jj) ! depth at the partial step level zh = pdep (ji,jj) ! depth at the partial step level
! !
zn = - rn_a0 * ( 1._wp + 0.5_wp*rn_lambda1*zt + rn_mu1*zh ) * zt & zn = - rn_a0 * ( 1._wp + 0.5_wp*rn_lambda1*zt + rn_mu1*zh ) * zt &
...@@ -742,8 +744,8 @@ CONTAINS ...@@ -742,8 +744,8 @@ CONTAINS
CASE( np_seos ) !== simplified EOS ==! CASE( np_seos ) !== simplified EOS ==!
! !
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
zt = pts (ji,jj,jp_tem) - 10._wp zt = pts (ji,jj,jp_tem) - rn_T0
zs = pts (ji,jj,jp_sal) - 35._wp zs = pts (ji,jj,jp_sal) - rn_S0
ztm = tmask(ji,jj,1) ztm = tmask(ji,jj,1)
! ! potential density referenced at the surface ! ! potential density referenced at the surface
zn = - rn_a0 * ( 1._wp + 0.5_wp*rn_lambda1*zt ) * zt & zn = - rn_a0 * ( 1._wp + 0.5_wp*rn_lambda1*zt ) * zt &
...@@ -853,9 +855,9 @@ CONTAINS ...@@ -853,9 +855,9 @@ CONTAINS
CASE( np_seos ) !== simplified EOS ==! CASE( np_seos ) !== simplified EOS ==!
! !
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 ) DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )
zt = pts (ji,jj,jk,jp_tem) - 10._wp ! pot. temperature anomaly (t-T0) zt = pts (ji,jj,jk,jp_tem) - rn_T0 ! pot. temperature anomaly (t-T0)
zs = pts (ji,jj,jk,jp_sal) - 35._wp ! abs. salinity anomaly (s-S0) zs = pts (ji,jj,jk,jp_sal) - rn_S0 ! abs. salinity anomaly (s-S0)
zh = gdept(ji,jj,jk,Kmm) ! depth in meters at t-point zh = gdept(ji,jj,jk,Kmm) ! depth in meters at t-point
ztm = tmask(ji,jj,jk) ! land/sea bottom mask = surf. mask ztm = tmask(ji,jj,jk) ! land/sea bottom mask = surf. mask
! !
zn = rn_a0 * ( 1._wp + rn_lambda1*zt + rn_mu1*zh ) + rn_nu*zs zn = rn_a0 * ( 1._wp + rn_lambda1*zt + rn_mu1*zh ) + rn_nu*zs
...@@ -973,9 +975,9 @@ CONTAINS ...@@ -973,9 +975,9 @@ CONTAINS
! !
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
! !
zt = pts (ji,jj,jp_tem) - 10._wp ! pot. temperature anomaly (t-T0) zt = pts (ji,jj,jp_tem) - rn_T0 ! pot. temperature anomaly (t-T0)
zs = pts (ji,jj,jp_sal) - 35._wp ! abs. salinity anomaly (s-S0) zs = pts (ji,jj,jp_sal) - rn_S0 ! abs. salinity anomaly (s-S0)
zh = pdep (ji,jj) ! depth at the partial step level zh = pdep (ji,jj) ! depth at the partial step level
! !
zn = rn_a0 * ( 1._wp + rn_lambda1*zt + rn_mu1*zh ) + rn_nu*zs zn = rn_a0 * ( 1._wp + rn_lambda1*zt + rn_mu1*zh ) + rn_nu*zs
pab(ji,jj,jp_tem) = zn * r1_rho0 ! alpha pab(ji,jj,jp_tem) = zn * r1_rho0 ! alpha
...@@ -1075,9 +1077,9 @@ CONTAINS ...@@ -1075,9 +1077,9 @@ CONTAINS
! !
CASE( np_seos ) !== simplified EOS ==! CASE( np_seos ) !== simplified EOS ==!
! !
zt = pts(jp_tem) - 10._wp ! pot. temperature anomaly (t-T0) zt = pts(jp_tem) - rn_T0 ! pot. temperature anomaly (t-T0)
zs = pts(jp_sal) - 35._wp ! abs. salinity anomaly (s-S0) zs = pts(jp_sal) - rn_S0 ! abs. salinity anomaly (s-S0)
zh = pdep ! depth at the partial step level zh = pdep ! depth at the partial step level
! !
zn = rn_a0 * ( 1._wp + rn_lambda1*zt + rn_mu1*zh ) + rn_nu*zs zn = rn_a0 * ( 1._wp + rn_lambda1*zt + rn_mu1*zh ) + rn_nu*zs
pab(jp_tem) = zn * r1_rho0 ! alpha pab(jp_tem) = zn * r1_rho0 ! alpha
...@@ -1164,28 +1166,37 @@ CONTAINS ...@@ -1164,28 +1166,37 @@ CONTAINS
!! Reference : TEOS-10, UNESCO !! Reference : TEOS-10, UNESCO
!! Rational approximation to TEOS10 algorithm (rms error on WOA13 values: 4.0e-5 degC) !! Rational approximation to TEOS10 algorithm (rms error on WOA13 values: 4.0e-5 degC)
!!---------------------------------------------------------------------- !!----------------------------------------------------------------------
REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: ctmp ! Cons. Temp [Celsius] REAL(wp), DIMENSION(:,:), INTENT(in ) :: ctmp ! Cons. Temp [Celsius]
REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: psal ! salinity [psu] REAL(wp), DIMENSION(:,:), INTENT(in ) :: psal ! salinity [psu]
! Leave result array automatic rather than making explicitly allocated ! Leave result array automatic rather than making explicitly allocated
REAL(wp), DIMENSION(jpi,jpj) :: ptmp ! potential temperature [Celsius] REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ptmp ! potential temperature [Celsius]
! !
INTEGER :: ji, jj ! dummy loop indices INTEGER :: ji, jj ! dummy loop indices
REAL(wp) :: zt , zs , ztm ! local scalars REAL(wp) :: zt , zs , ztm ! local scalars
REAL(wp) :: zn , zd ! local scalars REAL(wp) :: zn , zd ! local scalars
REAL(wp) :: zdeltaS , z1_S0 , z1_T0 REAL(wp) :: zdeltaS , z1_S0 , z1_T0
INTEGER :: ipi, ipj, iisht, ijsht ! dimensions and shift indices
!!---------------------------------------------------------------------- !!----------------------------------------------------------------------
! !
ipi = SIZE(psal,1) ! 1st dimension
ipj = SIZE(psal,2) ! 2nd dimension
!
iisht = ( jpi - ipi ) / 2
ijsht = ( jpj - ipj ) / 2 ! should be the same as iisht...
!
IF( .NOT.ALLOCATED(ptmp) ) ALLOCATE( ptmp(ipi,ipj) )
!
IF( ln_timing ) CALL timing_start('eos_pt_from_ct') IF( ln_timing ) CALL timing_start('eos_pt_from_ct')
! !
zdeltaS = 5._wp zdeltaS = 5._wp
z1_S0 = 0.875_wp/35.16504_wp z1_S0 = 0.875_wp/35.16504_wp
z1_T0 = 1._wp/40._wp z1_T0 = 1._wp/40._wp
! !
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) DO_2D( nn_hls-iisht, nn_hls-iisht, nn_hls-ijsht, nn_hls-ijsht )
! !
zt = ctmp (ji,jj) * z1_T0 zt = ctmp (ji-iisht,jj-ijsht) * z1_T0
zs = SQRT( ABS( psal(ji,jj) + zdeltaS ) * z1_S0 ) zs = SQRT( ABS( psal(ji-iisht,jj-ijsht) + zdeltaS ) * z1_S0 )
ztm = tmask(ji,jj,1) ztm = tmask(ji-iisht,jj-ijsht,1)
! !
zn = ((((-2.1385727895e-01_wp*zt & zn = ((((-2.1385727895e-01_wp*zt &
& - 2.7674419971e-01_wp*zs+1.0728094330_wp)*zt & & - 2.7674419971e-01_wp*zs+1.0728094330_wp)*zt &
...@@ -1200,7 +1211,7 @@ CONTAINS ...@@ -1200,7 +1211,7 @@ CONTAINS
& -3.4570358592e-01_wp*zs+5.6471810638_wp)*zt & & -3.4570358592e-01_wp*zs+5.6471810638_wp)*zt &
& + (1.5393993508_wp*zs-6.9394762624_wp)*zs+1.2750522650e+01_wp & + (1.5393993508_wp*zs-6.9394762624_wp)*zs+1.2750522650e+01_wp
! !
ptmp(ji,jj) = ( zt / z1_T0 + zn / zd ) * ztm ptmp(ji-iisht,jj-ijsht) = ( zt / z1_T0 + zn / zd ) * ztm
! !
END_2D END_2D
! !
...@@ -1211,9 +1222,9 @@ CONTAINS ...@@ -1211,9 +1222,9 @@ CONTAINS
SUBROUTINE eos_fzp_2d( psal, ptf, pdep ) SUBROUTINE eos_fzp_2d( psal, ptf, pdep )
!! !!
REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: psal ! salinity [psu] REAL(wp), DIMENSION(:,:), INTENT(in ) :: psal ! salinity [psu]
REAL(wp), DIMENSION(jpi,jpj), INTENT(in ), OPTIONAL :: pdep ! depth [m] REAL(wp), DIMENSION(:,:), INTENT(in ), OPTIONAL :: pdep ! depth [m]
REAL(wp), DIMENSION(:,:) , INTENT(out ) :: ptf ! freezing temperature [Celsius] REAL(wp), DIMENSION(:,:), INTENT(out ) :: ptf ! freezing temperature [Celsius]
!! !!
CALL eos_fzp_2d_t( psal, ptf, is_tile(ptf), pdep ) CALL eos_fzp_2d_t( psal, ptf, is_tile(ptf), pdep )
END SUBROUTINE eos_fzp_2d END SUBROUTINE eos_fzp_2d
...@@ -1231,35 +1242,52 @@ CONTAINS ...@@ -1231,35 +1242,52 @@ CONTAINS
!! !!
!! Reference : UNESCO tech. papers in the marine science no. 28. 1978 !! Reference : UNESCO tech. papers in the marine science no. 28. 1978
!!---------------------------------------------------------------------- !!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: kttf INTEGER , INTENT(in ) :: kttf
REAL(wp), DIMENSION(jpi,jpj) , INTENT(in ) :: psal ! salinity [psu] REAL(wp), DIMENSION(:,:), INTENT(in ) :: psal ! salinity [psu]
REAL(wp), DIMENSION(jpi,jpj) , INTENT(in ), OPTIONAL :: pdep ! depth [m] REAL(wp), DIMENSION(:,:), INTENT(in ), OPTIONAL :: pdep ! depth [m]
REAL(wp), DIMENSION(A2D_T(kttf)), INTENT(out ) :: ptf ! freezing temperature [Celsius] REAL(wp), DIMENSION(:,:), INTENT(out ) :: ptf ! freezing temperature [Celsius]
! !
INTEGER :: ji, jj ! dummy loop indices INTEGER :: ji, jj ! dummy loop indices
REAL(wp) :: zt, zs, z1_S0 ! local scalars REAL(wp) :: zt, zs, z1_S0 ! local scalars
INTEGER :: ipi, ipj, iisht, ijsht ! dimensions and shift indices
!!---------------------------------------------------------------------- !!----------------------------------------------------------------------
! !
ipi = SIZE(psal,1) ! 1st dimension
ipj = SIZE(psal,2) ! 2nd dimension
!
iisht = ( jpi - ipi ) / 2
ijsht = ( jpj - ipj ) / 2 ! should be the same as iisht...
!
SELECT CASE ( neos ) SELECT CASE ( neos )
! !
CASE ( np_teos10, np_seos ) !== CT,SA (TEOS-10 and S-EOS formulations) ==! CASE ( np_teos10, np_seos ) !== CT,SA (TEOS-10 and S-EOS formulations) ==!
! !
z1_S0 = 1._wp / 35.16504_wp z1_S0 = 1._wp / 35.16504_wp
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) DO_2D( nn_hls-iisht, nn_hls-iisht, nn_hls-ijsht, nn_hls-ijsht )
zs= SQRT( ABS( psal(ji,jj) ) * z1_S0 ) ! square root salinity zs= SQRT( ABS( psal(ji-iisht,jj-ijsht) ) * z1_S0 ) ! square root salinity
ptf(ji,jj) = ((((1.46873e-03_wp*zs-9.64972e-03_wp)*zs+2.28348e-02_wp)*zs & ptf(ji-iisht,jj-ijsht) = ((((1.46873e-03_wp*zs-9.64972e-03_wp)*zs+2.28348e-02_wp)*zs &
& - 3.12775e-02_wp)*zs+2.07679e-02_wp)*zs-5.87701e-02_wp & - 3.12775e-02_wp)*zs+2.07679e-02_wp)*zs-5.87701e-02_wp
ptf(ji-iisht,jj-ijsht) = ptf(ji-iisht,jj-ijsht) * psal(ji-iisht,jj-ijsht)
END_2D END_2D
ptf(:,:) = ptf(:,:) * psal(:,:)
! !
IF( PRESENT( pdep ) ) ptf(:,:) = ptf(:,:) - 7.53e-4 * pdep(:,:) IF( PRESENT( pdep ) ) THEN
DO_2D( nn_hls-iisht, nn_hls-iisht, nn_hls-ijsht, nn_hls-ijsht )
ptf(ji-iisht,jj-ijsht) = ptf(ji-iisht,jj-ijsht) - 7.53e-4 * pdep(ji-iisht,jj-ijsht)
END_2D
ENDIF
! !
CASE ( np_eos80 ) !== PT,SP (UNESCO formulation) ==! CASE ( np_eos80 ) !== PT,SP (UNESCO formulation) ==!
! !
ptf(:,:) = ( - 0.0575_wp + 1.710523e-3_wp * SQRT( psal(:,:) ) & DO_2D( nn_hls-iisht, nn_hls-iisht, nn_hls-ijsht, nn_hls-ijsht )
& - 2.154996e-4_wp * psal(:,:) ) * psal(:,:) ptf(ji-iisht,jj-ijsht) = ( - 0.0575_wp + 1.710523e-3_wp * SQRT( psal(ji-iisht,jj-ijsht) ) &
! & - 2.154996e-4_wp * psal(ji-iisht,jj-ijsht) ) * psal(ji-iisht,jj-ijsht)
IF( PRESENT( pdep ) ) ptf(:,:) = ptf(:,:) - 7.53e-4 * pdep(:,:) END_2D
!
IF( PRESENT( pdep ) ) THEN
DO_2D( nn_hls-iisht, nn_hls-iisht, nn_hls-ijsht, nn_hls-ijsht )
ptf(ji-iisht,jj-ijsht) = ptf(ji-iisht,jj-ijsht) - 7.53e-4 * pdep(ji-iisht,jj-ijsht)
END_2D
ENDIF
! !
CASE DEFAULT CASE DEFAULT
WRITE(ctmp1,*) ' bad flag value for neos = ', neos WRITE(ctmp1,*) ' bad flag value for neos = ', neos
...@@ -1336,10 +1364,10 @@ CONTAINS ...@@ -1336,10 +1364,10 @@ CONTAINS
!! pab_pe(:,:,:,jp_tem) is alpha_pe !! pab_pe(:,:,:,jp_tem) is alpha_pe
!! pab_pe(:,:,:,jp_sal) is beta_pe !! pab_pe(:,:,:,jp_sal) is beta_pe
!!---------------------------------------------------------------------- !!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: Kmm ! time level index INTEGER , INTENT(in ) :: Kmm ! time level index
REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: pts ! pot. temperature & salinity REAL(wp), DIMENSION(:,:,:,:), INTENT(in ) :: pts ! pot. temperature & salinity
REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT( out) :: pab_pe ! alpha_pe and beta_pe REAL(wp), DIMENSION(:,:,:,:), INTENT( out) :: pab_pe ! alpha_pe and beta_pe
REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT( out) :: ppen ! potential energy anomaly REAL(wp), DIMENSION(:,:,:) , INTENT( out) :: ppen ! potential energy anomaly
! !
INTEGER :: ji, jj, jk ! dummy loop indices INTEGER :: ji, jj, jk ! dummy loop indices
REAL(wp) :: zt , zh , zs , ztm ! local scalars REAL(wp) :: zt , zh , zs , ztm ! local scalars
...@@ -1352,7 +1380,7 @@ CONTAINS ...@@ -1352,7 +1380,7 @@ CONTAINS
! !
CASE( np_teos10, np_eos80 ) !== polynomial TEOS-10 / EOS-80 ==! CASE( np_teos10, np_eos80 ) !== polynomial TEOS-10 / EOS-80 ==!
! !
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 ) DO_3D( 0, 0, 0, 0, 1, jpkm1 )
! !
zh = gdept(ji,jj,jk,Kmm) * r1_Z0 ! depth zh = gdept(ji,jj,jk,Kmm) * r1_Z0 ! depth
zt = pts (ji,jj,jk,jp_tem) * r1_T0 ! temperature zt = pts (ji,jj,jk,jp_tem) * r1_T0 ! temperature
...@@ -1411,11 +1439,11 @@ CONTAINS ...@@ -1411,11 +1439,11 @@ CONTAINS
! !
CASE( np_seos ) !== Vallis (2006) simplified EOS ==! CASE( np_seos ) !== Vallis (2006) simplified EOS ==!
! !
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 ) DO_3D( 0, 0, 0, 0, 1, jpkm1 )
zt = pts(ji,jj,jk,jp_tem) - 10._wp ! temperature anomaly (t-T0) zt = pts (ji,jj,jk,jp_tem) - rn_T0 ! temperature anomaly (t-T0)
zs = pts (ji,jj,jk,jp_sal) - 35._wp ! abs. salinity anomaly (s-S0) zs = pts (ji,jj,jk,jp_sal) - rn_S0 ! abs. salinity anomaly (s-S0)
zh = gdept(ji,jj,jk,Kmm) ! depth in meters at t-point zh = gdept(ji,jj,jk,Kmm) ! depth in meters at t-point
ztm = tmask(ji,jj,jk) ! tmask ztm = tmask(ji,jj,jk) ! tmask
zn = 0.5_wp * zh * r1_rho0 * ztm zn = 0.5_wp * zh * r1_rho0 * ztm
! ! Potential Energy ! ! Potential Energy
ppen(ji,jj,jk) = ( rn_a0 * rn_mu1 * zt + rn_b0 * rn_mu2 * zs ) * zn ppen(ji,jj,jk) = ( rn_a0 * rn_mu1 * zt + rn_b0 * rn_mu2 * zs ) * zn
...@@ -1447,8 +1475,8 @@ CONTAINS ...@@ -1447,8 +1475,8 @@ CONTAINS
INTEGER :: ios ! local integer INTEGER :: ios ! local integer
INTEGER :: ioptio ! local integer INTEGER :: ioptio ! local integer
!! !!
NAMELIST/nameos/ ln_TEOS10, ln_EOS80, ln_SEOS, rn_a0, rn_b0, rn_lambda1, rn_mu1, & NAMELIST/nameos/ ln_TEOS10, ln_EOS80, ln_SEOS, rn_T0, rn_S0, rn_a0, rn_b0, rn_lambda1, rn_mu1, &
& rn_lambda2, rn_mu2, rn_nu & rn_lambda2, rn_mu2, rn_nu
!!---------------------------------------------------------------------- !!----------------------------------------------------------------------
! !
READ ( numnam_ref, nameos, IOSTAT = ios, ERR = 901 ) READ ( numnam_ref, nameos, IOSTAT = ios, ERR = 901 )
...@@ -1869,9 +1897,11 @@ CONTAINS ...@@ -1869,9 +1897,11 @@ CONTAINS
IF(lwp) THEN IF(lwp) THEN
WRITE(numout,*) WRITE(numout,*)
WRITE(numout,*) ' ==>>> use of simplified eos: ' WRITE(numout,*) ' ==>>> use of simplified eos: '
WRITE(numout,*) ' rhd(dT=T-10,dS=S-35,Z) = [-a0*(1+lambda1/2*dT+mu1*Z)*dT ' WRITE(numout,*) ' rhd(dT=T-rn_T0,dS=S-rn_S0,Z) = [-a0*(1+lambda1/2*dT+mu1*Z)*dT '
WRITE(numout,*) ' + b0*(1+lambda2/2*dT+mu2*Z)*dS - nu*dT*dS] / rho0' WRITE(numout,*) ' + b0*(1+lambda2/2*dT+mu2*Z)*dS - nu*dT*dS] / rho0'
WRITE(numout,*) ' with the following coefficients :' WRITE(numout,*) ' with the following coefficients :'
WRITE(numout,*) ' reference temperature rn_T0 = ', rn_T0
WRITE(numout,*) ' reference salinity rn_S0 = ', rn_S0
WRITE(numout,*) ' thermal exp. coef. rn_a0 = ', rn_a0 WRITE(numout,*) ' thermal exp. coef. rn_a0 = ', rn_a0
WRITE(numout,*) ' saline cont. coef. rn_b0 = ', rn_b0 WRITE(numout,*) ' saline cont. coef. rn_b0 = ', rn_b0
WRITE(numout,*) ' cabbeling coef. rn_lambda1 = ', rn_lambda1 WRITE(numout,*) ' cabbeling coef. rn_lambda1 = ', rn_lambda1
......
src/OCE/TRA/traadv.F90
View file @ 3efbcb17
...@@ -10,6 +10,7 @@ MODULE traadv ...@@ -10,6 +10,7 @@ MODULE traadv
!! - ! 2014-12 (G. Madec) suppression of cross land advection option !! - ! 2014-12 (G. Madec) suppression of cross land advection option
!! 3.6 ! 2015-06 (E. Clementi) Addition of Stokes drift in case of wave coupling !! 3.6 ! 2015-06 (E. Clementi) Addition of Stokes drift in case of wave coupling
!! 4.5 ! 2021-04 (G. Madec, S. Techene) add advective velocities as optional arguments !! 4.5 ! 2021-04 (G. Madec, S. Techene) add advective velocities as optional arguments
!! 4.5 ! 2022-06 (S. Techene, G, Madec) refactorization to reduce local memory usage
!!---------------------------------------------------------------------- !!----------------------------------------------------------------------
!!---------------------------------------------------------------------- !!----------------------------------------------------------------------
...@@ -178,7 +179,7 @@ CONTAINS ...@@ -178,7 +179,7 @@ CONTAINS
! !
!!gm ??? !!gm ???
! TEMP: [tiling] This copy-in not necessary after all lbc_lnks removed in the nn_hls = 2 case in tra_adv_fct ! TEMP: [tiling] This copy-in not necessary after all lbc_lnks removed in the nn_hls = 2 case in tra_adv_fct
IF( l_diaptr ) CALL dia_ptr( kt, Kmm, zvv(A2D(nn_hls),:) ) ! diagnose the effective MSF IF( l_diaptr ) CALL dia_ptr( kt, Kmm, zvv(A2D(0),:) ) ! diagnose the effective MSF
!!gm ??? !!gm ???
! !
...@@ -191,11 +192,19 @@ CONTAINS ...@@ -191,11 +192,19 @@ CONTAINS
SELECT CASE ( nadv ) !== compute advection trend and add it to general trend ==! SELECT CASE ( nadv ) !== compute advection trend and add it to general trend ==!
! !
CASE ( np_CEN ) ! Centered scheme : 2nd / 4th order CASE ( np_CEN ) ! Centered scheme : 2nd / 4th order
CALL tra_adv_cen ( kt, nit000, 'TRA', zuu, zvv, zww, Kmm, pts, jpts, Krhs, nn_cen_h, nn_cen_v ) IF( nn_hls == 1 ) THEN
CALL tra_adv_cen_hls1( kt, nit000, 'TRA', zuu, zvv, zww, Kmm, pts, jpts, Krhs, nn_cen_h, nn_cen_v )
ELSE
CALL tra_adv_cen ( kt, nit000, 'TRA', zuu, zvv, zww, Kmm, pts, jpts, Krhs, nn_cen_h, nn_cen_v )
ENDIF
CASE ( np_FCT ) ! FCT scheme : 2nd / 4th order CASE ( np_FCT ) ! FCT scheme : 2nd / 4th order
CALL tra_adv_fct ( kt, nit000, 'TRA', rDt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, nn_fct_h, nn_fct_v ) CALL tra_adv_fct ( kt, nit000, 'TRA', rDt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, nn_fct_h, nn_fct_v )
CASE ( np_MUS ) ! MUSCL CASE ( np_MUS ) ! MUSCL
IF( nn_hls == 1 ) THEN
CALL tra_adv_mus_hls1( kt, nit000, 'TRA', rDt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, ln_mus_ups )
ELSE
CALL tra_adv_mus( kt, nit000, 'TRA', rDt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, ln_mus_ups ) CALL tra_adv_mus( kt, nit000, 'TRA', rDt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, ln_mus_ups )
END IF
CASE ( np_UBS ) ! UBS CASE ( np_UBS ) ! UBS
CALL tra_adv_ubs ( kt, nit000, 'TRA', rDt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, nn_ubs_v ) CALL tra_adv_ubs ( kt, nit000, 'TRA', rDt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, nn_ubs_v )
CASE ( np_QCK ) ! QUICKEST CASE ( np_QCK ) ! QUICKEST
......
src/OCE/TRA/traadv_cen.F90
View file @ 3efbcb17
...@@ -4,6 +4,7 @@ MODULE traadv_cen ...@@ -4,6 +4,7 @@ MODULE traadv_cen
!! Ocean tracers: advective trend (2nd/4th order centered) !! Ocean tracers: advective trend (2nd/4th order centered)
!!====================================================================== !!======================================================================
!! History : 3.7 ! 2014-05 (G. Madec) original code !! History : 3.7 ! 2014-05 (G. Madec) original code
!! 4.5 ! 2022-06 (S. Techene, G, Madec) refactorization to reduce local memory usage
!!---------------------------------------------------------------------- !!----------------------------------------------------------------------
!!---------------------------------------------------------------------- !!----------------------------------------------------------------------
...@@ -29,7 +30,8 @@ MODULE traadv_cen ...@@ -29,7 +30,8 @@ MODULE traadv_cen
IMPLICIT NONE IMPLICIT NONE
PRIVATE PRIVATE
PUBLIC tra_adv_cen ! called by traadv.F90 PUBLIC tra_adv_cen ! called by traadv.F90
PUBLIC tra_adv_cen_hls1 ! called by traadv.F90
REAL(wp) :: r1_6 = 1._wp / 6._wp ! =1/6 REAL(wp) :: r1_6 = 1._wp / 6._wp ! =1/6
...@@ -47,7 +49,321 @@ MODULE traadv_cen ...@@ -47,7 +49,321 @@ MODULE traadv_cen
!!---------------------------------------------------------------------- !!----------------------------------------------------------------------
CONTAINS CONTAINS
SUBROUTINE tra_adv_cen_test( kt, kit000, cdtype, pU, pV, pW, &
& Kmm, pt, kjpt, Krhs, kn_cen_h, kn_cen_v )
!!----------------------------------------------------------------------
!! *** ROUTINE tra_adv_cen ***
!!
!! ** Purpose : Compute the now trend due to the advection of tracers
!! and add it to the general trend of passive tracer equations.
!!
!! ** Method : The advection is evaluated by a 2nd or 4th order scheme
!! using now fields (leap-frog scheme).
!! kn_cen_h = 2 ==>> 2nd order centered scheme on the horizontal
!! = 4 ==>> 4th order - - - -
!! kn_cen_v = 2 ==>> 2nd order centered scheme on the vertical
!! = 4 ==>> 4th order COMPACT scheme - -
!!
!! ** Action : - update pt(:,:,:,:,Krhs) with the now advective tracer trends
!! - send trends to trdtra module for further diagnostcs (l_trdtra=T)
!! - poleward advective heat and salt transport (l_diaptr=T)
!!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: kt ! ocean time-step index
INTEGER , INTENT(in ) :: Kmm, Krhs ! ocean time level indices
INTEGER , INTENT(in ) :: kit000 ! first time step index
CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)
INTEGER , INTENT(in ) :: kjpt ! number of tracers
INTEGER , INTENT(in ) :: kn_cen_h ! =2/4 (2nd or 4th order scheme)
INTEGER , INTENT(in ) :: kn_cen_v ! =2/4 (2nd or 4th order scheme)
! TEMP: [tiling] This can be A2D(nn_hls) after all lbc_lnks removed in the nn_hls = 2 case in tra_adv_fct
REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pU, pV, pW ! 3 ocean volume flux components
REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) :: pt ! tracers and RHS of tracer equation
!
INTEGER :: ji, jj, jk, jn ! dummy loop indices
INTEGER :: ierr ! local integer
REAL(wp) :: zC2t_u, zC4t_u ! local scalars
REAL(wp) :: zC2t_v, zC4t_v ! - -
REAL(wp) :: zftw_kp1
REAL(wp), DIMENSION(A2D(1)) :: zft_u, zft_v !, zft_w
REAL(wp), DIMENSION(:,:) , ALLOCATABLE :: zdt_u, zdt_v
REAL(wp), DIMENSION(:,:,:) , ALLOCATABLE :: ztw
!!----------------------------------------------------------------------
!
#if defined key_loop_fusion
CALL tra_adv_cen_lf ( kt, nit000, cdtype, pU, pV, pW, Kmm, pt, kjpt, Krhs, kn_cen_h, kn_cen_v )
#else
IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile
IF( kt == kit000 ) THEN
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) 'tra_adv_cen : centered advection scheme on ', cdtype, ' order h/v =', kn_cen_h,'/', kn_cen_v
IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~ '
ENDIF
! ! set local switches
l_trd = .FALSE.
l_hst = .FALSE.
l_ptr = .FALSE.
IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE.
IF( cdtype == 'TRA' .AND. ( iom_use( 'sophtadv' ) .OR. iom_use( 'sophtadv' ) ) ) l_ptr = .TRUE.
IF( cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. &
& iom_use("uadv_salttr") .OR. iom_use("vadv_salttr") ) ) l_hst = .TRUE.
ENDIF
!
! !* 2nd order centered
DO jn = 1, kjpt !== loop over the tracers ==!
!
DO jk = 1, jpkm1
!
DO_2D( 1, 0, 1, 0 ) ! Horizontal fluxes at layer jk
zft_u(ji,jj) = 0.5_wp * pU(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji+1,jj ,jk,jn,Kmm) )
zft_v(ji,jj) = 0.5_wp * pV(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji ,jj+1,jk,jn,Kmm) )
END_2D
!
DO_2D( 0, 0, 0, 0 ) ! Horizontal divergence of advective fluxes
pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( zft_u(ji,jj) - zft_u(ji-1,jj ) &
& + zft_v(ji,jj) - zft_v(ji ,jj-1) ) * r1_e1e2t(ji,jj) &
& / e3t(ji,jj,jk,Kmm)
END_2D
END DO
!
#define zft_w zft_u
IF( ln_linssh ) THEN !* top value (linear free surf. only as zwz is multiplied by wmask)
DO_2D( 0, 0, 0, 0 )
zft_w(ji,jj) = pW(ji,jj,1) * pt(ji,jj,1,jn,Kmm)
END_2D
ELSE
DO_2D( 0, 0, 0, 0 )
zft_w(ji,jj) = 0._wp
END_2D
ENDIF
DO jk = 1, jpk-2
DO_2D( 0, 0, 0, 0 ) ! Vertical fluxes
zftw_kp1 = 0.5 * pW(ji,jj,jk+1) * ( pt(ji,jj,jk+1,jn,Kmm) + pt(ji,jj,jk,jn,Kmm) ) * wmask(ji,jj,jk+1)
!
pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( zft_w(ji,jj) - zftw_kp1 ) * r1_e1e2t(ji,jj) &
& / e3t(ji,jj,jk,Kmm)
zft_w(ji,jj) = zftw_kp1
END_2D
END DO
jk = jpkm1 ! bottom vertical flux set to zero for all tracers
DO_2D( 0, 0, 0, 0 )
pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - zft_w(ji,jj) * r1_e1e2t(ji,jj) &
& / e3t(ji,jj,jk,Kmm)
END_2D
!
END DO
!
!
#undef zft_w
! ! trend diagnostics
!!gm + !!st to be done with the whole rewritting of trd
!! trd routine arguments MUST be changed adding jk and zwx, zwy in 2D
!! IF( l_trd ) THEN
!! CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_xad, zwx, pU, pt(:,:,:,jn,Kmm) )
!! CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_yad, zwy, pV, pt(:,:,:,jn,Kmm) )
!! CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_zad, zwz, pW, pt(:,:,:,jn,Kmm) )
!! ENDIF
!! ! ! "Poleward" heat and salt transports
!! IF( l_ptr ) CALL dia_ptr_hst( jn, 'adv', zwy(:,:,:) )
!! ! ! heat and salt transport
!! IF( l_hst ) CALL dia_ar5_hst( jn, 'adv', zwx(:,:,:), zwy(:,:,:) )
!
!
!
#endif
END SUBROUTINE tra_adv_cen_test
SUBROUTINE tra_adv_cen( kt, kit000, cdtype, pU, pV, pW, & SUBROUTINE tra_adv_cen( kt, kit000, cdtype, pU, pV, pW, &
& Kmm, pt, kjpt, Krhs, kn_cen_h, kn_cen_v )
!!----------------------------------------------------------------------
!! *** ROUTINE tra_adv_cen ***
!!
!! ** Purpose : Compute the now trend due to the advection of tracers
!! and add it to the general trend of passive tracer equations.
!!
!! ** Method : The advection is evaluated by a 2nd or 4th order scheme
!! using now fields (leap-frog scheme).
!! kn_cen_h = 2 ==>> 2nd order centered scheme on the horizontal
!! = 4 ==>> 4th order - - - -
!! kn_cen_v = 2 ==>> 2nd order centered scheme on the vertical
!! = 4 ==>> 4th order COMPACT scheme - -
!!
!! ** Action : - update pt(:,:,:,:,Krhs) with the now advective tracer trends
!! - send trends to trdtra module for further diagnostcs (l_trdtra=T)
!! - poleward advective heat and salt transport (l_diaptr=T)
!!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: kt ! ocean time-step index
INTEGER , INTENT(in ) :: Kmm, Krhs ! ocean time level indices
INTEGER , INTENT(in ) :: kit000 ! first time step index
CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)
INTEGER , INTENT(in ) :: kjpt ! number of tracers
INTEGER , INTENT(in ) :: kn_cen_h ! =2/4 (2nd or 4th order scheme)
INTEGER , INTENT(in ) :: kn_cen_v ! =2/4 (2nd or 4th order scheme)
! TEMP: [tiling] This can be A2D(nn_hls) after all lbc_lnks removed in the nn_hls = 2 case in tra_adv_fct
REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pU, pV, pW ! 3 ocean volume flux components
REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) :: pt ! tracers and RHS of tracer equation
!
INTEGER :: ji, jj, jk, jn ! dummy loop indices
INTEGER :: ierr ! local integer
REAL(wp) :: zC2t_u, zC4t_u ! local scalars
REAL(wp) :: zC2t_v, zC4t_v ! - -
REAL(wp) :: zftw_kp1
REAL(wp), DIMENSION(A2D(1)) :: zft_u, zft_v
REAL(wp), DIMENSION(:,:) , ALLOCATABLE :: zdt_u, zdt_v
REAL(wp), DIMENSION(:,:,:) , ALLOCATABLE :: ztw
!!----------------------------------------------------------------------
!
#if defined key_loop_fusion
CALL tra_adv_cen_lf ( kt, nit000, cdtype, pU, pV, pW, Kmm, pt, kjpt, Krhs, kn_cen_h, kn_cen_v )
#else
IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile
IF( kt == kit000 ) THEN
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) 'tra_adv_cen : centered advection scheme on ', cdtype, ' order h/v =', kn_cen_h,'/', kn_cen_v
IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~ '
ENDIF
! ! set local switches
l_trd = .FALSE.
l_hst = .FALSE.
l_ptr = .FALSE.
IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE.
IF( cdtype == 'TRA' .AND. ( iom_use( 'sophtadv' ) .OR. iom_use( 'sophtadv' ) ) ) l_ptr = .TRUE.
IF( cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. &
& iom_use("uadv_salttr") .OR. iom_use("vadv_salttr") ) ) l_hst = .TRUE.
ENDIF
!
IF( kn_cen_h == 4 ) ALLOCATE( zdt_u(A2D(2)) , zdt_v(A2D(2)) ) ! horizontal 4th order only
IF( kn_cen_v == 4 ) ALLOCATE( ztw(A2D(nn_hls),jpk) ) ! vertical 4th order only
!
DO jn = 1, kjpt !== loop over the tracers ==!
!
SELECT CASE( kn_cen_h ) !-- Horizontal divergence of advective fluxes --!
!
!!st limitation : does not take into acccount iceshelf specificity
!! in case of linssh
CASE( 2 ) !* 2nd order centered
DO jk = 1, jpkm1
!
DO_2D( 1, 0, 1, 0 ) ! Horizontal fluxes at layer jk
zft_u(ji,jj) = 0.5_wp * pU(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji+1,jj ,jk,jn,Kmm) )
zft_v(ji,jj) = 0.5_wp * pV(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji ,jj+1,jk,jn,Kmm) )
END_2D
!
DO_2D( 0, 0, 0, 0 ) ! Horizontal divergence of advective fluxes
pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( zft_u(ji,jj) - zft_u(ji-1,jj ) &
& + zft_v(ji,jj) - zft_v(ji ,jj-1) ) * r1_e1e2t(ji,jj) &
& / e3t(ji,jj,jk,Kmm)
END_2D
END DO
!
CASE( 4 ) !* 4th order centered
DO jk = 1, jpkm1
DO_2D( 2, 1, 2, 1 ) ! masked gradient
zdt_u(ji,jj) = ( pt(ji+1,jj ,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * umask(ji,jj,jk)
zdt_v(ji,jj) = ( pt(ji ,jj+1,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * vmask(ji,jj,jk)
END_2D
!
DO_2D( 1, 0, 1, 0 ) ! Horizontal advective fluxes
zC2t_u = pt(ji,jj,jk,jn,Kmm) + pt(ji+1,jj ,jk,jn,Kmm) ! C2 interpolation of T at u- & v-points (x2)
zC2t_v = pt(ji,jj,jk,jn,Kmm) + pt(ji ,jj+1,jk,jn,Kmm)
! ! C4 interpolation of T at u- & v-points (x2)
zC4t_u = zC2t_u + r1_6 * ( zdt_u(ji-1,jj ) - zdt_u(ji+1,jj ) )
zC4t_v = zC2t_v + r1_6 * ( zdt_v(ji ,jj-1) - zdt_v(ji ,jj+1) )
! ! C4 fluxes
zft_u(ji,jj) = 0.5_wp * pU(ji,jj,jk) * zC4t_u
zft_v(ji,jj) = 0.5_wp * pV(ji,jj,jk) * zC4t_v
END_2D
!
DO_2D( 0, 0, 0, 0 ) ! Horizontal divergence of advective fluxes
pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( zft_u(ji,jj) - zft_u(ji-1,jj ) &
& + zft_v(ji,jj) - zft_v(ji ,jj-1) ) * r1_e1e2t(ji,jj) &
& / e3t(ji,jj,jk,Kmm)
END_2D
END DO
!
CASE DEFAULT
CALL ctl_stop( 'traadv_cen: wrong value for nn_cen' )
END SELECT
!
#define zft_w zft_u
!
IF( ln_linssh ) THEN !* top value (linear free surf. only as zwz is multiplied by wmask)
DO_2D( 0, 0, 0, 0 )
zft_w(ji,jj) = pW(ji,jj,1) * pt(ji,jj,1,jn,Kmm)
END_2D
ELSE
DO_2D( 0, 0, 0, 0 )
zft_w(ji,jj) = 0._wp
END_2D
ENDIF
!
SELECT CASE( kn_cen_v ) !-- Vertical divergence of advective fluxes --! (interior)
!
CASE( 2 ) !* 2nd order centered
DO jk = 1, jpk-2
DO_2D( 0, 0, 0, 0 ) ! Vertical fluxes
zftw_kp1 = 0.5 * pW(ji,jj,jk+1) * ( pt(ji,jj,jk+1,jn,Kmm) + pt(ji,jj,jk,jn,Kmm) ) * wmask(ji,jj,jk+1)
!
pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( zft_w(ji,jj) - zftw_kp1 ) * r1_e1e2t(ji,jj) &
& / e3t(ji,jj,jk,Kmm)
zft_w(ji,jj) = zftw_kp1
END_2D
END DO
jk = jpkm1 ! bottom vertical flux set to zero for all tracers
DO_2D( 0, 0, 0, 0 )
pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - zft_w(ji,jj) * r1_e1e2t(ji,jj) &
& / e3t(ji,jj,jk,Kmm)
END_2D
!
CASE( 4 ) !* 4th order compact
CALL interp_4th_cpt( pt(:,:,:,jn,Kmm) , ztw ) ! ztw = interpolated value of T at w-point
!
DO jk = 1, jpk-2
!
DO_2D( 0, 0, 0, 0 )
zftw_kp1 = pW(ji,jj,jk+1) * ztw(ji,jj,jk+1) * wmask(ji,jj,jk+1)
! ! Divergence of advective fluxes
pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( zft_w(ji,jj) - zftw_kp1 ) * r1_e1e2t(ji,jj) &
& / e3t(ji,jj,jk,Kmm)
! ! update
zft_w(ji,jj) = zftw_kp1
END_2D
!
END DO
!
jk = jpkm1 ! bottom vertical flux set to zero for all tracers
DO_2D( 0, 0, 0, 0 )
pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - zft_w(ji,jj) * r1_e1e2t(ji,jj) &
& / e3t(ji,jj,jk,Kmm)
END_2D
!
END SELECT
!
#undef zft_w
! ! trend diagnostics
!!gm + !!st to be done with the whole rewritting of trd
!! trd routine arguments MUST be changed adding jk and zwx, zwy in 2D
!! IF( l_trd ) THEN
!! CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_xad, zwx, pU, pt(:,:,:,jn,Kmm) )
!! CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_yad, zwy, pV, pt(:,:,:,jn,Kmm) )
!! CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_zad, zwz, pW, pt(:,:,:,jn,Kmm) )
!! ENDIF
!! ! ! "Poleward" heat and salt transports
!! IF( l_ptr ) CALL dia_ptr_hst( jn, 'adv', zwy(:,:,:) )
!! ! ! heat and salt transport
!! IF( l_hst ) CALL dia_ar5_hst( jn, 'adv', zwx(:,:,:), zwy(:,:,:) )
!
END DO
!
IF( kn_cen_h == 4 ) DEALLOCATE( zdt_u , zdt_v ) ! horizontal 4th order only
IF( kn_cen_v == 4 ) DEALLOCATE( ztw ) ! vertical 4th order only
!
#endif
END SUBROUTINE tra_adv_cen
SUBROUTINE tra_adv_cen_hls1( kt, kit000, cdtype, pU, pV, pW, &
& Kmm, pt, kjpt, Krhs, kn_cen_h, kn_cen_v ) & Kmm, pt, kjpt, Krhs, kn_cen_h, kn_cen_v )
!!---------------------------------------------------------------------- !!----------------------------------------------------------------------
!! *** ROUTINE tra_adv_cen *** !! *** ROUTINE tra_adv_cen ***
...@@ -141,6 +457,12 @@ CONTAINS ...@@ -141,6 +457,12 @@ CONTAINS
CASE DEFAULT CASE DEFAULT
CALL ctl_stop( 'traadv_cen: wrong value for nn_cen' ) CALL ctl_stop( 'traadv_cen: wrong value for nn_cen' )
END SELECT END SELECT
DO_3D( 0, 0, 0, 0, 1, jpkm1 ) !-- Divergence of advective fluxes --!
pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) &
& - ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) &
& + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) ) &
& * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
END_3D
! !
SELECT CASE( kn_cen_v ) !-- Vertical fluxes --! (interior) SELECT CASE( kn_cen_v ) !-- Vertical fluxes --! (interior)
! !
...@@ -171,11 +493,17 @@ CONTAINS ...@@ -171,11 +493,17 @@ CONTAINS
! !
DO_3D( 0, 0, 0, 0, 1, jpkm1 ) !-- Divergence of advective fluxes --! DO_3D( 0, 0, 0, 0, 1, jpkm1 ) !-- Divergence of advective fluxes --!
pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) & pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) &
& - ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) & & - ( zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) &
& + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) &
& + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) &
& * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) & * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
END_3D END_3D
!
!!st DO_3D( 0, 0, 0, 0, 1, jpkm1 ) !-- Divergence of advective fluxes --!
!!st pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) &
!!st & - ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) &
!!st & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) &
!!st & + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) &
!!st & * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
!!st END_3D
! ! trend diagnostics ! ! trend diagnostics
IF( l_trd ) THEN IF( l_trd ) THEN
CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_xad, zwx, pU, pt(:,:,:,jn,Kmm) ) CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_xad, zwx, pU, pt(:,:,:,jn,Kmm) )
...@@ -183,14 +511,14 @@ CONTAINS ...@@ -183,14 +511,14 @@ CONTAINS
CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_zad, zwz, pW, pt(:,:,:,jn,Kmm) ) CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_zad, zwz, pW, pt(:,:,:,jn,Kmm) )
ENDIF ENDIF
! ! "Poleward" heat and salt transports ! ! "Poleward" heat and salt transports
IF( l_ptr ) CALL dia_ptr_hst( jn, 'adv', zwy(:,:,:) ) IF( l_ptr ) CALL dia_ptr_hst( jn, 'adv', zwy(A2D(0),:) )
! ! heat and salt transport ! ! heat and salt transport
IF( l_hst ) CALL dia_ar5_hst( jn, 'adv', zwx(:,:,:), zwy(:,:,:) ) IF( l_hst ) CALL dia_ar5_hst( jn, 'adv', zwx(:,:,:), zwy(:,:,:) )
! !
END DO END DO
! !
#endif #endif
END SUBROUTINE tra_adv_cen END SUBROUTINE tra_adv_cen_hls1
!!====================================================================== !!======================================================================
END MODULE traadv_cen END MODULE traadv_cen
src/OCE/TRA/traadv_cen_lf.F90
View file @ 3efbcb17
...@@ -175,7 +175,7 @@ CONTAINS ...@@ -175,7 +175,7 @@ CONTAINS
CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_zad, zwz, pW, pt(:,:,:,jn,Kmm) ) CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_zad, zwz, pW, pt(:,:,:,jn,Kmm) )
ENDIF ENDIF
! ! "Poleward" heat and salt transports ! ! "Poleward" heat and salt transports
IF( l_ptr ) CALL dia_ptr_hst( jn, 'adv', zwy(:,:,:) ) IF( l_ptr ) CALL dia_ptr_hst( jn, 'adv', zwy(A2D(0),:) )
! ! heat and salt transport ! ! heat and salt transport
IF( l_hst ) CALL dia_ar5_hst( jn, 'adv', zwx(:,:,:), zwy(:,:,:) ) IF( l_hst ) CALL dia_ar5_hst( jn, 'adv', zwx(:,:,:), zwy(:,:,:) )
! !
......
src/OCE/TRA/traadv_fct.F90
View file @ 3efbcb17
...@@ -130,7 +130,7 @@ CONTAINS ...@@ -130,7 +130,7 @@ CONTAINS
ENDIF ENDIF
! !
IF( l_ptr ) THEN IF( l_ptr ) THEN
ALLOCATE( zptry(A2D(nn_hls),jpk) ) ALLOCATE( zptry(A2D(0),jpk) )
zptry(:,:,:) = 0._wp zptry(:,:,:) = 0._wp
ENDIF ENDIF
! !
...@@ -213,7 +213,7 @@ CONTAINS ...@@ -213,7 +213,7 @@ CONTAINS
ztrdx(:,:,:) = zwx(:,:,:) ; ztrdy(:,:,:) = zwy(:,:,:) ; ztrdz(:,:,:) = zwz(:,:,:) ztrdx(:,:,:) = zwx(:,:,:) ; ztrdy(:,:,:) = zwy(:,:,:) ; ztrdz(:,:,:) = zwz(:,:,:)
END IF END IF
! ! "Poleward" heat and salt transports (contribution of upstream fluxes) ! ! "Poleward" heat and salt transports (contribution of upstream fluxes)
IF( l_ptr ) zptry(:,:,:) = zwy(:,:,:) IF( l_ptr ) zptry(:,:,:) = zwy(A2D(0),:)
! !
! !== anti-diffusive flux : high order minus low order ==! ! !== anti-diffusive flux : high order minus low order ==!
! !
...@@ -364,7 +364,7 @@ CONTAINS ...@@ -364,7 +364,7 @@ CONTAINS
! !
ENDIF ENDIF
IF( l_ptr ) THEN ! "Poleward" transports IF( l_ptr ) THEN ! "Poleward" transports
zptry(:,:,:) = zptry(:,:,:) + zwy(:,:,:) ! <<< add anti-diffusive fluxes zptry(:,:,:) = zptry(:,:,:) + zwy(A2D(0),:) ! <<< add anti-diffusive fluxes
CALL dia_ptr_hst( jn, 'adv', zptry(:,:,:) ) CALL dia_ptr_hst( jn, 'adv', zptry(:,:,:) )
ENDIF ENDIF
! !
...@@ -753,7 +753,7 @@ CONTAINS ...@@ -753,7 +753,7 @@ CONTAINS
ENDIF ENDIF
! !
IF( l_ptr ) THEN IF( l_ptr ) THEN
ALLOCATE( zptry(jpi,jpj,jpk) ) ALLOCATE( zptry(A2D(0),jpk) )
zptry(:,:,:) = 0._wp zptry(:,:,:) = 0._wp
ENDIF ENDIF
! !
...@@ -838,7 +838,7 @@ CONTAINS ...@@ -838,7 +838,7 @@ CONTAINS
ztrdx(:,:,:) = zwx_3d(:,:,:) ; ztrdy(:,:,:) = zwy_3d(:,:,:) ; ztrdz(:,:,:) = zwz(:,:,:) ztrdx(:,:,:) = zwx_3d(:,:,:) ; ztrdy(:,:,:) = zwy_3d(:,:,:) ; ztrdz(:,:,:) = zwz(:,:,:)
END IF END IF
! ! "Poleward" heat and salt transports (contribution of upstream fluxes) ! ! "Poleward" heat and salt transports (contribution of upstream fluxes)
IF( l_ptr ) zptry(:,:,:) = zwy_3d(:,:,:) IF( l_ptr ) zptry(:,:,:) = zwy_3d(A2D(0),:)
! !
! !== anti-diffusive flux : high order minus low order ==! ! !== anti-diffusive flux : high order minus low order ==!
! !
...@@ -986,7 +986,7 @@ CONTAINS ...@@ -986,7 +986,7 @@ CONTAINS
ENDIF ENDIF
! NOT TESTED - NEED l_ptr TRUE ! NOT TESTED - NEED l_ptr TRUE
IF( l_ptr ) THEN ! "Poleward" transports IF( l_ptr ) THEN ! "Poleward" transports
zptry(:,:,:) = zptry(:,:,:) + zwy_3d(:,:,:) ! <<< add anti-diffusive fluxes zptry(:,:,:) = zptry(:,:,:) + zwy_3d(A2D(0),:) ! <<< add anti-diffusive fluxes
CALL dia_ptr_hst( jn, 'adv', zptry(:,:,:) ) CALL dia_ptr_hst( jn, 'adv', zptry(:,:,:) )
ENDIF ENDIF
! !
......
src/OCE/TRA/traadv_mus.F90
View file @ 3efbcb17
...@@ -9,6 +9,7 @@ MODULE traadv_mus ...@@ -9,6 +9,7 @@ MODULE traadv_mus
!! 3.2 ! 2010-05 (C. Ethe, G. Madec) merge TRC-TRA + switch from velocity to transport !! 3.2 ! 2010-05 (C. Ethe, G. Madec) merge TRC-TRA + switch from velocity to transport
!! 3.4 ! 2012-06 (P. Oddo, M. Vichi) include the upstream where needed !! 3.4 ! 2012-06 (P. Oddo, M. Vichi) include the upstream where needed
!! 3.7 ! 2015-09 (G. Madec) add the ice-shelf cavities boundary condition !! 3.7 ! 2015-09 (G. Madec) add the ice-shelf cavities boundary condition
!! 4.5 ! 2022-06 (S. Techene, G, Madec) refactorization to reduce local memory usage
!!---------------------------------------------------------------------- !!----------------------------------------------------------------------
!!---------------------------------------------------------------------- !!----------------------------------------------------------------------
...@@ -34,7 +35,9 @@ MODULE traadv_mus ...@@ -34,7 +35,9 @@ MODULE traadv_mus
IMPLICIT NONE IMPLICIT NONE
PRIVATE PRIVATE
PUBLIC tra_adv_mus ! routine called by traadv.F90 PUBLIC tra_adv_mus ! routine called by traadv.F90
PUBLIC tra_adv_mus_hls1 ! routine called by traadv.F90
REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: upsmsk !: mixed upstream/centered scheme near some straits REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: upsmsk !: mixed upstream/centered scheme near some straits
! ! and in closed seas (orca 2 and 1 configurations) ! ! and in closed seas (orca 2 and 1 configurations)
...@@ -55,6 +58,217 @@ MODULE traadv_mus ...@@ -55,6 +58,217 @@ MODULE traadv_mus
CONTAINS CONTAINS
SUBROUTINE tra_adv_mus( kt, kit000, cdtype, p2dt, pU, pV, pW, & SUBROUTINE tra_adv_mus( kt, kit000, cdtype, p2dt, pU, pV, pW, &
& Kbb, Kmm, pt, kjpt, Krhs, ld_msc_ups )
!!----------------------------------------------------------------------
!! *** ROUTINE tra_adv_mus ***
!!
!! ** Purpose : Compute the now trend due to total advection of tracers
!! using a MUSCL scheme (Monotone Upstream-centered Scheme for
!! Conservation Laws) and add it to the general tracer trend.
!!
!! ** Method : MUSCL scheme plus centered scheme at ocean boundaries
!! ld_msc_ups=T :
!!
!! ** Action : - update pt(:,:,:,:,Krhs) with the now advective tracer trends
!! - send trends to trdtra module for further diagnostcs (l_trdtra=T)
!! - poleward advective heat and salt transport (ln_diaptr=T)
!!
!! References : Estubier, A., and M. Levy, Notes Techn. Pole de Modelisation
!! IPSL, Sept. 2000 (http://www.lodyc.jussieu.fr/opa)
!!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: kt ! ocean time-step index
INTEGER , INTENT(in ) :: Kbb, Kmm, Krhs ! ocean time level indices
INTEGER , INTENT(in ) :: kit000 ! first time step index
CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)
INTEGER , INTENT(in ) :: kjpt ! number of tracers
LOGICAL , INTENT(in ) :: ld_msc_ups ! use upstream scheme within muscl
REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step
! TEMP: [tiling] This can be A2D(nn_hls) after all lbc_lnks removed in the nn_hls = 2 case in tra_adv_fct
REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pU, pV, pW ! 3 ocean volume flux components
REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) :: pt ! tracers and RHS of tracer equation
!
INTEGER :: ji, jj, jk, jn ! dummy loop indices
INTEGER :: ierr ! local integer
INTEGER :: ik
REAL(wp) :: zu, z0u, zzslpx, zzwx, zw , zalpha ! local scalars
REAL(wp) :: zv, z0v, zzslpy, zzwy, z0w ! - -
REAL(wp) :: zdzt_kp2, zslpz_kp1, zfW_kp1
REAL(wp), DIMENSION(A2D(2)) :: zdxt, zslpx, zwx ! 2D workspace
REAL(wp), DIMENSION(A2D(2)) :: zdyt, zslpy, zwy ! - -
!!----------------------------------------------------------------------
!
IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile
IF( kt == kit000 ) THEN
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) 'tra_adv : MUSCL advection scheme on ', cdtype
IF(lwp) WRITE(numout,*) ' : mixed up-stream ', ld_msc_ups
IF(lwp) WRITE(numout,*) '~~~~~~~'
IF(lwp) WRITE(numout,*)
!
! Upstream / MUSCL scheme indicator
!
ALLOCATE( xind(jpi,jpj,jpk), STAT=ierr )
xind(:,:,:) = 1._wp ! set equal to 1 where up-stream is not needed
IF( ld_msc_ups ) THEN ! define the upstream indicator (if asked)
DO jk = 1, jpkm1
xind(:,:,jk) = 1._wp & ! =>1 where up-stream is not needed
& - rnfmsk(:,:) * rnfmsk_z(jk) * tmask(:,:,jk) ! =>0 near runoff mouths (& closed sea outflows)
END DO
ENDIF
!
ENDIF
!
l_trd = .FALSE.
l_hst = .FALSE.
l_ptr = .FALSE.
IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE.
IF( l_diaptr .AND. cdtype == 'TRA' .AND. ( iom_use( 'sophtadv' ) .OR. iom_use( 'sophtadv' ) ) ) l_ptr = .TRUE.
IF( l_iom .AND. cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. &
& iom_use("uadv_salttr") .OR. iom_use("vadv_salttr") ) ) l_hst = .TRUE.
ENDIF
!
!
DO jn = 1, kjpt !== loop over the tracers ==!
!
DO jk = 1, jpkm1
! !* Horizontal advective fluxes
!
! !-- first guess of the slopes
DO_2D( 2, 1, 2, 1 )
zdxt(ji,jj) = ( pt(ji+1,jj ,jk,jn,Kbb) - pt(ji,jj,jk,jn,Kbb) ) * umask(ji,jj,jk)
zdyt(ji,jj) = ( pt(ji ,jj+1,jk,jn,Kbb) - pt(ji,jj,jk,jn,Kbb) ) * vmask(ji,jj,jk)
END_2D
! !-- Slopes of tracer
DO_2D( 1, 1, 1, 1 )
! ! 1/2 Slopes at T-point (set to 0 if adjectent slopes are of opposite sign)
zzslpx = ( zdxt(ji,jj) + zdxt(ji-1,jj ) ) &
& * ( 0.25_wp + SIGN( 0.25_wp, zdxt(ji,jj) * zdxt(ji-1,jj ) ) )
zzslpy = ( zdyt(ji,jj) + zdyt(ji ,jj-1) ) &
& * ( 0.25_wp + SIGN( 0.25_wp, zdyt(ji,jj) * zdyt(ji ,jj-1) ) )
! ! Slopes limitation
zslpx(ji,jj) = SIGN( 1.0_wp, zzslpx ) * MIN( ABS( zzslpx ), &
& 2._wp*ABS( zdxt (ji-1,jj) ), &
& 2._wp*ABS( zdxt (ji ,jj) ) )
zslpy(ji,jj) = SIGN( 1.0_wp, zzslpy ) * MIN( ABS( zzslpy ), &
& 2._wp*ABS( zdyt (ji,jj-1) ), &
& 2._wp*ABS( zdyt (ji,jj ) ) )
END_2D
!!gm + !!st ticket ? comparaison pommes et carrottes ABS(zzslpx) et 2._wp*ABS( zdxt (ji-1,jj) )
!
DO_2D( 1, 0, 1, 0 ) !-- MUSCL horizontal advective fluxes
z0u = SIGN( 0.5_wp, pU(ji,jj,jk) )
zalpha = 0.5_wp - z0u
zu = z0u - 0.5_wp * pU(ji,jj,jk) * p2dt * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm)
zzwx = pt(ji+1,jj,jk,jn,Kbb) + xind(ji,jj,jk) * zu * zslpx(ji+1,jj)
zzwy = pt(ji ,jj,jk,jn,Kbb) + xind(ji,jj,jk) * zu * zslpx(ji ,jj)
zwx(ji,jj) = pU(ji,jj,jk) * ( zalpha * zzwx + (1._wp-zalpha) * zzwy )
!
z0v = SIGN( 0.5_wp, pV(ji,jj,jk) )
zalpha = 0.5_wp - z0v
zv = z0v - 0.5_wp * pV(ji,jj,jk) * p2dt * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm)
zzwx = pt(ji,jj+1,jk,jn,Kbb) + xind(ji,jj,jk) * zv * zslpy(ji,jj+1)
zzwy = pt(ji,jj ,jk,jn,Kbb) + xind(ji,jj,jk) * zv * zslpy(ji,jj )
zwy(ji,jj) = pV(ji,jj,jk) * ( zalpha * zzwx + (1._wp-zalpha) * zzwy )
END_2D
!
DO_2D( 0, 0, 0, 0 ) !-- Tracer advective trend
pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( zwx(ji,jj) - zwx(ji-1,jj ) &
& + zwy(ji,jj) - zwy(ji ,jj-1) ) &
& * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
END_2D
END DO
!!gm + !!st to be done with the whole rewritting of trd
!! trd routine arguments MUST be changed adding jk and zwx, zwy in 2D
!!
!! ! ! trend diagnostics
!! IF( l_trd ) THEN
!! CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jk, jptra_xad, zwx(:,:), pU, pt(:,:,:,jn,Kbb) )
!! CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jk, jptra_yad, zwy(:,:), pV, pt(:,:,:,jn,Kbb) )
!! END IF
!! ! ! "Poleward" heat and salt transports
!! IF( l_ptr ) CALL dia_ptr_hst( jn, jk, 'adv', zwy(:,:) )
!! ! ! heat transport
!! IF( l_hst ) CALL dia_ar5_hst( jn, jk, 'adv', zwx(:,:), zwy(:,:) )
!
! !* Vertical advective fluxes
!
#define zdzt_kp1 zdxt
#define zslpz zslpx
#define zfW zwx
zfW (A2D(0)) = 0._wp ! anciennement zwx at jk = 1
! ! anciennement zwx at jk = 2
DO_2D( 0, 0, 0, 0 )
zdzt_kp1(ji,jj) = tmask(ji,jj,2) * ( pt(ji,jj,1,jn,Kbb) - pt(ji,jj,2,jn,Kbb) )
END_2D
zslpz (A2D(0)) = 0._wp ! anciennement zslpx at jk = 1
!
IF( ln_linssh ) THEN !-- linear ssh : non zero top values
DO_2D( 0, 0, 0, 0 ) ! at the ocean surface
zfW(ji,jj) = pW(ji,jj,1) * pt(ji,jj,1,jn,Kbb) ! surface flux
END_2D
IF( ln_isfcav ) THEN ! ice-shelf cavities (top of the ocean)
DO_2D( 0, 0, 0, 0 ) ! update pt(Krhs) under the ice-shelf
ik = mikt(ji,jj) ! the flux at ik-1 is zero ( inside ice-shelf )
IF( ik > 1 ) THEN
pt(ji,jj,ik,jn,Krhs) = pt(ji,jj,ik,jn,Krhs) - pW(ji,jj,ik) * pt(ji,jj,ik,jn,Kbb) &
& * r1_e1e2t(ji,jj) / e3t(ji,jj,ik,Kmm)
ENDIF
END_2D
ENDIF
ENDIF
!
! wmask usage for computing zw and zwk is needed in isf case and linear ssh
!
!
DO jk = 1, jpkm1
IF( jk < jpkm1 ) THEN
DO_2D( 0, 0, 0, 0 )
! !-- Slopes of tracer
! ! masked vertical gradient at jk+2
zdzt_kp2 = ( pt(ji,jj,jk+1,jn,Kbb) - pt(ji,jj,jk+2,jn,Kbb) ) * tmask(ji,jj,jk+2) !!st wmask(ji,jj,jk+2)
! ! vertical slope at jk+1
zslpz_kp1 = ( zdzt_kp1(ji,jj) + zdzt_kp2 ) &
& * ( 0.25_wp + SIGN( 0.25_wp, zdzt_kp1(ji,jj) * zdzt_kp2 ) )
! ! slope limitation
zslpz_kp1 = SIGN( 1.0_wp, zslpz_kp1 ) * MIN( ABS( zslpz_kp1 ), &
& 2.*ABS( zdzt_kp2 ), &
& 2.*ABS( zdzt_kp1(ji,jj) ) )
! !-- vertical advective flux at jk+1
! ! caution: zfW_kp1 is masked for ice-shelf cavities
! ! since top fluxes already added to pt(Krhs) before the vertical loop
z0w = SIGN( 0.5_wp, pW(ji,jj,jk+1) )
zalpha = 0.5_wp + z0w
zw = z0w - 0.5_wp * pW(ji,jj,jk+1) * p2dt * r1_e1e2t(ji,jj) / e3w(ji,jj,jk+1,Kmm)
zzwx = pt(ji,jj,jk+1,jn,Kbb) + xind(ji,jj,jk) * zw * zslpz_kp1
zzwy = pt(ji,jj,jk ,jn,Kbb) + xind(ji,jj,jk) * zw * zslpz(ji,jj)
zfW_kp1 = pW(ji,jj,jk+1) * ( zalpha * zzwx + (1.-zalpha) * zzwy ) * wmask(ji,jj,jk)!!st * wmask(ji,jj,jk+1)
! !-- vertical advective trend at jk
pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( zfW(ji,jj) - zfW_kp1 ) &
& * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
! ! updates for next level
zdzt_kp1(ji,jj) = zdzt_kp2
zslpz (ji,jj) = zslpz_kp1
zfW (ji,jj) = zfW_kp1
END_2D
ELSE
DO_2D( 0, 0, 0, 0 ) !-- vertical advective trend at jpkm1
pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - zfW(ji,jj) &
& * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
END_2D
ENDIF
END DO ! end of jk loop
!
!!gm + !!st idem see above
!! ! ! send trends for diagnostic
!! IF( l_trd ) CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_zad, zwx, pW, pt(:,:,:,jn,Kbb) )
!
END DO ! end of tracer loop
!
END SUBROUTINE tra_adv_mus
SUBROUTINE tra_adv_mus_hls1( kt, kit000, cdtype, p2dt, pU, pV, pW, &
& Kbb, Kmm, pt, kjpt, Krhs, ld_msc_ups ) & Kbb, Kmm, pt, kjpt, Krhs, ld_msc_ups )
!!---------------------------------------------------------------------- !!----------------------------------------------------------------------
!! *** ROUTINE tra_adv_mus *** !! *** ROUTINE tra_adv_mus ***
...@@ -178,7 +392,7 @@ CONTAINS ...@@ -178,7 +392,7 @@ CONTAINS
! !
DO_3D( 0, 0, 0, 0, 1, jpkm1 ) !-- Tracer advective trend DO_3D( 0, 0, 0, 0, 1, jpkm1 ) !-- Tracer advective trend
pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) & pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) &
& + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) ) & & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) ) &
& * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) & * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
END_3D END_3D
! ! trend diagnostics ! ! trend diagnostics
...@@ -187,7 +401,7 @@ CONTAINS ...@@ -187,7 +401,7 @@ CONTAINS
CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_yad, zwy, pV, pt(:,:,:,jn,Kbb) ) CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_yad, zwy, pV, pt(:,:,:,jn,Kbb) )
END IF END IF
! ! "Poleward" heat and salt transports ! ! "Poleward" heat and salt transports
IF( l_ptr ) CALL dia_ptr_hst( jn, 'adv', zwy(:,:,:) ) IF( l_ptr ) CALL dia_ptr_hst( jn, 'adv', zwy(A2D(0),:) )
! ! heat transport ! ! heat transport
IF( l_hst ) CALL dia_ar5_hst( jn, 'adv', zwx(:,:,:), zwy(:,:,:) ) IF( l_hst ) CALL dia_ar5_hst( jn, 'adv', zwx(:,:,:), zwy(:,:,:) )
! !
...@@ -239,7 +453,7 @@ CONTAINS ...@@ -239,7 +453,7 @@ CONTAINS
! !
END DO ! end of tracer loop END DO ! end of tracer loop
! !
END SUBROUTINE tra_adv_mus END SUBROUTINE tra_adv_mus_hls1
!!====================================================================== !!======================================================================
END MODULE traadv_mus END MODULE traadv_mus
src/OCE/TRA/traadv_qck.F90
View file @ 3efbcb17
...@@ -301,7 +301,7 @@ CONTAINS ...@@ -301,7 +301,7 @@ CONTAINS
! ! trend diagnostics ! ! trend diagnostics
IF( l_trd ) CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_yad, zwy, pV, pt(:,:,:,jn,Kmm) ) IF( l_trd ) CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_yad, zwy, pV, pt(:,:,:,jn,Kmm) )
! ! "Poleward" heat and salt transports (contribution of upstream fluxes) ! ! "Poleward" heat and salt transports (contribution of upstream fluxes)
IF( l_ptr ) CALL dia_ptr_hst( jn, 'adv', zwy(:,:,:) ) IF( l_ptr ) CALL dia_ptr_hst( jn, 'adv', zwy(A2D(0),:) )
! !
END DO END DO
! !
......
src/OCE/TRA/traadv_qck_lf.F90
View file @ 3efbcb17
...@@ -268,7 +268,7 @@ CONTAINS ...@@ -268,7 +268,7 @@ CONTAINS
! ! trend diagnostics ! ! trend diagnostics
IF( l_trd ) CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_yad, zwy, pV, pt(:,:,:,jn,Kmm) ) IF( l_trd ) CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_yad, zwy, pV, pt(:,:,:,jn,Kmm) )
! ! "Poleward" heat and salt transports (contribution of upstream fluxes) ! ! "Poleward" heat and salt transports (contribution of upstream fluxes)
IF( l_ptr ) CALL dia_ptr_hst( jn, 'adv', zwy(:,:,:) ) IF( l_ptr ) CALL dia_ptr_hst( jn, 'adv', zwy(A2D(0),:) )
! !
END DO END DO
! !
......
src/OCE/TRA/traadv_ubs.F90
View file @ 3efbcb17
...@@ -184,7 +184,7 @@ CONTAINS ...@@ -184,7 +184,7 @@ CONTAINS
END IF END IF
! !
! ! "Poleward" heat and salt transports (contribution of upstream fluxes) ! ! "Poleward" heat and salt transports (contribution of upstream fluxes)
IF( l_ptr ) CALL dia_ptr_hst( jn, 'adv', ztv(:,:,:) ) IF( l_ptr ) CALL dia_ptr_hst( jn, 'adv', ztv(A2D(0),:) )
! ! heati/salt transport ! ! heati/salt transport
IF( l_hst ) CALL dia_ar5_hst( jn, 'adv', ztu(:,:,:), ztv(:,:,:) ) IF( l_hst ) CALL dia_ar5_hst( jn, 'adv', ztu(:,:,:), ztv(:,:,:) )
! !
......
src/OCE/TRA/traadv_ubs_lf.F90
View file @ 3efbcb17
...@@ -190,7 +190,7 @@ CONTAINS ...@@ -190,7 +190,7 @@ CONTAINS
END IF END IF
! !
! ! "Poleward" heat and salt transports (contribution of upstream fluxes) ! ! "Poleward" heat and salt transports (contribution of upstream fluxes)
IF( l_ptr ) CALL dia_ptr_hst( jn, 'adv', ztv(:,:,:) ) IF( l_ptr ) CALL dia_ptr_hst( jn, 'adv', ztv(A2D(0),:) )
! ! heati/salt transport ! ! heati/salt transport
IF( l_hst ) CALL dia_ar5_hst( jn, 'adv', ztu(:,:,:), ztv(:,:,:) ) IF( l_hst ) CALL dia_ar5_hst( jn, 'adv', ztu(:,:,:), ztv(:,:,:) )
! !
......
src/OCE/TRA/traatf.F90
View file @ 3efbcb17
...@@ -207,7 +207,7 @@ CONTAINS ...@@ -207,7 +207,7 @@ CONTAINS
! !
DO jn = 1, kjpt DO jn = 1, kjpt
! !
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 ) DO_3D( 0, 0, 0, 0, 1, jpkm1 )
ztn = pt(ji,jj,jk,jn,Kmm) ztn = pt(ji,jj,jk,jn,Kmm)
ztd = pt(ji,jj,jk,jn,Kaa) - 2._wp * ztn + pt(ji,jj,jk,jn,Kbb) ! time laplacian on tracers ztd = pt(ji,jj,jk,jn,Kaa) - 2._wp * ztn + pt(ji,jj,jk,jn,Kbb) ! time laplacian on tracers
! !
...@@ -238,8 +238,8 @@ CONTAINS ...@@ -238,8 +238,8 @@ CONTAINS
CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)
INTEGER , INTENT(in ) :: kjpt ! number of tracers INTEGER , INTENT(in ) :: kjpt ! number of tracers
REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) :: pt ! tracer fields REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) :: pt ! tracer fields
REAL(wp), DIMENSION(jpi,jpj ,kjpt) , INTENT(in ) :: psbc_tc ! surface tracer content REAL(wp), DIMENSION(A2D(0) ,kjpt) , INTENT(in ) :: psbc_tc ! surface tracer content
REAL(wp), DIMENSION(jpi,jpj ,kjpt) , INTENT(in ) :: psbc_tc_b ! before surface tracer content REAL(wp), DIMENSION(A2D(0) ,kjpt) , INTENT(in ) :: psbc_tc_b ! before surface tracer content
! !
LOGICAL :: ll_traqsr, ll_rnf, ll_isf ! local logical LOGICAL :: ll_traqsr, ll_rnf, ll_isf ! local logical
INTEGER :: ji, jj, jk, jn ! dummy loop indices INTEGER :: ji, jj, jk, jn ! dummy loop indices
...@@ -269,14 +269,11 @@ CONTAINS ...@@ -269,14 +269,11 @@ CONTAINS
ztrd_atf(:,:,:,:) = 0.0_wp ztrd_atf(:,:,:,:) = 0.0_wp
ENDIF ENDIF
! !
!!st variables only computed in the interior by traqsr
IF( ll_traqsr ) CALL lbc_lnk( 'traatf', qsr_hc_b(:,:,:) , 'T', 1.0_wp, qsr_hc(:,:,:) , 'T', 1.0_wp )
!
zfact = 1._wp / p2dt zfact = 1._wp / p2dt
zfact1 = rn_atfp * p2dt zfact1 = rn_atfp * p2dt
zfact2 = zfact1 * r1_rho0 zfact2 = zfact1 * r1_rho0
DO jn = 1, kjpt DO jn = 1, kjpt
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 ) DO_3D( 0, 0, 0, 0, 1, jpkm1 )
ze3t_b = e3t(ji,jj,jk,Kbb) ze3t_b = e3t(ji,jj,jk,Kbb)
ze3t_n = e3t(ji,jj,jk,Kmm) ze3t_n = e3t(ji,jj,jk,Kmm)
ze3t_a = e3t(ji,jj,jk,Kaa) ze3t_a = e3t(ji,jj,jk,Kaa)
......
src/OCE/TRA/traatf_qco.F90
View file @ 3efbcb17
...@@ -150,7 +150,7 @@ CONTAINS ...@@ -150,7 +150,7 @@ CONTAINS
ELSE ; CALL tra_atf_qco_lf( kt, Kbb, Kmm, Kaa, nit000, rn_Dt, 'TRA', pts, sbc_tsc, sbc_tsc_b, jpts ) ! non-linear free surface ELSE ; CALL tra_atf_qco_lf( kt, Kbb, Kmm, Kaa, nit000, rn_Dt, 'TRA', pts, sbc_tsc, sbc_tsc_b, jpts ) ! non-linear free surface
ENDIF ENDIF
! !
CALL lbc_lnk( 'traatfqco', pts(:,:,:,jp_tem,Kmm) , 'T', 1._wp, pts(:,:,:,jp_sal,Kmm) , 'T', 1._wp ) CALL lbc_lnk( 'traatf_qco', pts(:,:,:,jp_tem,Kmm) , 'T', 1._wp, pts(:,:,:,jp_sal,Kmm) , 'T', 1._wp )
! !
ENDIF ENDIF
! !
...@@ -203,7 +203,6 @@ CONTAINS ...@@ -203,7 +203,6 @@ CONTAINS
DO jn = 1, kjpt DO jn = 1, kjpt
! !
DO_3D( 0, 0, 0, 0, 1, jpkm1 ) DO_3D( 0, 0, 0, 0, 1, jpkm1 )
!!st DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )
ztn = pt(ji,jj,jk,jn,Kmm) ztn = pt(ji,jj,jk,jn,Kmm)
ztd = pt(ji,jj,jk,jn,Kaa) - 2._wp * ztn + pt(ji,jj,jk,jn,Kbb) ! time laplacian on tracers ztd = pt(ji,jj,jk,jn,Kaa) - 2._wp * ztn + pt(ji,jj,jk,jn,Kbb) ! time laplacian on tracers
! !
...@@ -234,8 +233,8 @@ CONTAINS ...@@ -234,8 +233,8 @@ CONTAINS
CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)
INTEGER , INTENT(in ) :: kjpt ! number of tracers INTEGER , INTENT(in ) :: kjpt ! number of tracers
REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) :: pt ! tracer fields REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) :: pt ! tracer fields
REAL(wp), DIMENSION(jpi,jpj ,kjpt) , INTENT(in ) :: psbc_tc ! surface tracer content REAL(wp), DIMENSION(A2D(0) ,kjpt) , INTENT(in ) :: psbc_tc ! surface tracer content
REAL(wp), DIMENSION(jpi,jpj ,kjpt) , INTENT(in ) :: psbc_tc_b ! before surface tracer content REAL(wp), DIMENSION(A2D(0) ,kjpt) , INTENT(in ) :: psbc_tc_b ! before surface tracer content
! !
LOGICAL :: ll_traqsr, ll_rnf, ll_isf ! local logical LOGICAL :: ll_traqsr, ll_rnf, ll_isf ! local logical
INTEGER :: ji, jj, jk, jn ! dummy loop indices INTEGER :: ji, jj, jk, jn ! dummy loop indices
...@@ -264,12 +263,12 @@ CONTAINS ...@@ -264,12 +263,12 @@ CONTAINS
ALLOCATE( ztrd_atf(jpi,jpj,jpk,kjpt) ) ALLOCATE( ztrd_atf(jpi,jpj,jpk,kjpt) )
ztrd_atf(:,:,:,:) = 0._wp ztrd_atf(:,:,:,:) = 0._wp
ENDIF ENDIF
!
zfact = 1._wp / p2dt zfact = 1._wp / p2dt
zfact1 = rn_atfp * p2dt zfact1 = rn_atfp * p2dt
zfact2 = zfact1 * r1_rho0 zfact2 = zfact1 * r1_rho0
DO jn = 1, kjpt DO jn = 1, kjpt
DO_3D( 0, 0, 0, 0, 1, jpkm1 ) DO_3D( 0, 0, 0, 0, 1, jpkm1 )
!!st DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )
ze3t_b = e3t(ji,jj,jk,Kbb) ze3t_b = e3t(ji,jj,jk,Kbb)
ze3t_n = e3t(ji,jj,jk,Kmm) ze3t_n = e3t(ji,jj,jk,Kmm)
ze3t_a = e3t(ji,jj,jk,Kaa) ze3t_a = e3t(ji,jj,jk,Kaa)
......
src/OCE/TRA/tradmp.F90
View file @ 3efbcb17
...@@ -94,7 +94,7 @@ CONTAINS ...@@ -94,7 +94,7 @@ CONTAINS
REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts ! active tracers and RHS of tracer equation REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts ! active tracers and RHS of tracer equation
! !
INTEGER :: ji, jj, jk, jn ! dummy loop indices INTEGER :: ji, jj, jk, jn ! dummy loop indices
REAL(wp), DIMENSION(A2D(nn_hls),jpk,jpts) :: zts_dta REAL(wp), DIMENSION(A2D(nn_hls),jpk,jpts) :: zts_dta
REAL(wp), DIMENSION(:,:,:) , ALLOCATABLE :: zwrk REAL(wp), DIMENSION(:,:,:) , ALLOCATABLE :: zwrk
REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: ztrdts REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: ztrdts
!!---------------------------------------------------------------------- !!----------------------------------------------------------------------
...@@ -146,7 +146,7 @@ CONTAINS ...@@ -146,7 +146,7 @@ CONTAINS
! !
! outputs (clem trunk) ! outputs (clem trunk)
IF( iom_use('hflx_dmp_cea') .OR. iom_use('sflx_dmp_cea') ) THEN IF( iom_use('hflx_dmp_cea') .OR. iom_use('sflx_dmp_cea') ) THEN
ALLOCATE( zwrk(A2D(nn_hls),jpk) ) ! Needed to handle expressions containing e3t when using key_qco or key_linssh ALLOCATE( zwrk(A2D(0),jpk) ) ! Needed to handle expressions containing e3t when using key_qco or key_linssh
zwrk(:,:,:) = 0._wp zwrk(:,:,:) = 0._wp
IF( iom_use('hflx_dmp_cea') ) THEN IF( iom_use('hflx_dmp_cea') ) THEN
......
src/OCE/TRA/traldf_iso.F90
View file @ 3efbcb17
...@@ -388,7 +388,7 @@ CONTAINS ...@@ -388,7 +388,7 @@ CONTAINS
! !
! ! "Poleward" diffusive heat or salt transports (T-S case only) ! ! "Poleward" diffusive heat or salt transports (T-S case only)
! note sign is reversed to give down-gradient diffusive transports ) ! note sign is reversed to give down-gradient diffusive transports )
IF( l_ptr ) CALL dia_ptr_hst( jn, 'ldf', -zftv(:,:,:) ) IF( l_ptr ) CALL dia_ptr_hst( jn, 'ldf', -zftv(A2D(0),:) )
! ! Diffusive heat transports ! ! Diffusive heat transports
IF( l_hst ) CALL dia_ar5_hst( jn, 'ldf', -zftu(:,:,:), -zftv(:,:,:) ) IF( l_hst ) CALL dia_ar5_hst( jn, 'ldf', -zftu(:,:,:), -zftv(:,:,:) )
! !
......
src/OCE/TRA/traldf_lap_blp.F90
View file @ 3efbcb17
...@@ -171,7 +171,7 @@ CONTAINS ...@@ -171,7 +171,7 @@ CONTAINS
IF( ( kpass == 1 .AND. .NOT.ln_traldf_blp ) .OR. & !== first pass only ( laplacian) ==! IF( ( kpass == 1 .AND. .NOT.ln_traldf_blp ) .OR. & !== first pass only ( laplacian) ==!
( kpass == 2 .AND. ln_traldf_blp ) ) THEN !== 2nd pass only (bilaplacian) ==! ( kpass == 2 .AND. ln_traldf_blp ) ) THEN !== 2nd pass only (bilaplacian) ==!
IF( l_ptr ) CALL dia_ptr_hst( jn, 'ldf', -ztv(:,:,:) ) IF( l_ptr ) CALL dia_ptr_hst( jn, 'ldf', -ztv(A2D(0),:) )
IF( l_hst ) CALL dia_ar5_hst( jn, 'ldf', -ztu(:,:,:), -ztv(:,:,:) ) IF( l_hst ) CALL dia_ar5_hst( jn, 'ldf', -ztu(:,:,:), -ztv(:,:,:) )
ENDIF ENDIF
! ! ================== ! ! ==================
......