#include"cppdefs.h" !----------------------------------------------------------------------- !----------------------------------------------------------------------- !BOP ! ! !ROUTINE: Calculate steady-state Richardson number from c3\label{sec:Rist} ! ! !INTERFACE: REALTYPE function compute_rist(c1,c2,c3) ! ! !DESCRIPTION: ! Numerically computes the steady-state Richardson-number $Ri_{st}$ ! for two-equations models from the given ! $c_{\psi 3}$ and the parameters ! $c_{\psi 1}$ and $c_{\psi 2}$ according to \eq{Ri_st}. ! A (very tricky) double Newton-iteration is used to solve the resulting ! implicit non-linear equation. ! ! !USES: use turbulence, only: as,an,cmue1,cmue2 use turbulence, only: cm0 use turbulence, only: turb_method,stab_method use turbulence, only: cm0_fix,Prandtl0_fix use turbulence, only: Constant use turbulence, only: MunkAnderson use turbulence, only: SchumGerz use turbulence, only: EiflerSchrimpf IMPLICIT NONE ! ! !INPUT PARAMETERS: REALTYPE, intent(in) :: c1,c2,c3 ! ! !REVISION HISTORY: ! Original author(s): Hans Burchard, Lars Umlauf ! !EOP !----------------------------------------------------------------------- ! !LOCAL VARIABLES: integer :: i,j integer,parameter :: imax=100 REALTYPE :: cc3,fc,fp,step,ann REALTYPE :: ffc,ffp,Ri,Rii REALTYPE :: NN(0:2),SS(0:2) logical :: converged REALTYPE,parameter :: e=1.e-9,ee=1.e-4 ! !----------------------------------------------------------------------- !BOC NN=0. SS=0. Ri=0.18 ! Initial guess for Rist converged = .true. do j=0,imax Rii=Ri ann=0.1 do i=0,imax an(1)=ann as(1)=an(1)/Rii if (turb_method.eq.2) then select case(stab_method) case(Constant) cmue1=cm0_fix cmue2=cm0_fix/Prandtl0_fix case(MunkAnderson) call cmue_ma(2) case(SchumGerz) call cmue_sg(2) case(EiflerSchrimpf) call cmue_rf(2) end select else call cmue_d(2) end if fc=cmue1(1)*an(1)/Rii-cmue2(1)*an(1)-cm0**(-3) an(1)=ann+e as(1)=an(1)/Rii if (turb_method.eq.2) then select case(stab_method) case(Constant) cmue1=cm0_fix cmue2=cm0_fix/Prandtl0_fix case(MunkAnderson) call cmue_ma(2) case(SchumGerz) call cmue_sg(2) case(EiflerSchrimpf) call cmue_rf(2) end select else call cmue_d(2) end if fp=cmue1(1)*an(1)/Rii-cmue2(1)*an(1)-cm0**(-3) step=-fc/((fp-fc)/e) ann=ann+0.5*step if (abs(step).gt.100.) then STDERR ' ' STDERR 'Method for calculating the steady-state Richardson' STDERR 'number Rist does not converge.' STDERR 'Probably, the prescribed value for c3' STDERR 'is outside the range of the chosen stability function.' STDERR 'Please change gotmturb.nml accordingly.' STDERR 'You have chosen the stability function no.',stab_method STDERR 'If the problem persists, then use another' STDERR 'stability function.' STDERR 'I will continue anyway. Good luck!' STDERR ' ' converged=.false. goto 333 endif if (abs(step).lt.1.e-10) goto 111 end do 111 an(1)=ann as(1)=an(1)/Rii if (turb_method.eq.2) then select case(stab_method) case(Constant) cmue1=cm0_fix cmue2=cm0_fix/Prandtl0_fix case(MunkAnderson) call cmue_ma(2) case(SchumGerz) call cmue_sg(2) case(EiflerSchrimpf) call cmue_rf(2) end select else call cmue_d(2) end if cc3=c2+(c1-c2)/Rii*cmue1(1)/cmue2(1) ffc=cc3-c3 Rii=Ri+ee ann=0.1 do i=0,imax an(1)=ann as(1)=an(1)/Rii if (turb_method.eq.2) then select case(stab_method) case(Constant) cmue1=cm0_fix cmue2=cm0_fix/Prandtl0_fix case(MunkAnderson) call cmue_ma(2) case(SchumGerz) call cmue_sg(2) case(EiflerSchrimpf) call cmue_rf(2) end select else call cmue_d(2) end if fc=cmue1(1)*an(1)/Rii-cmue2(1)*an(1)-cm0**(-3) an(1)=ann+e as(1)=an(1)/Rii if (turb_method.eq.2) then select case(stab_method) case(Constant) cmue1=cm0_fix cmue2=cm0_fix/Prandtl0_fix case(MunkAnderson) call cmue_ma(2) case(SchumGerz) call cmue_sg(2) case(EiflerSchrimpf) call cmue_rf(2) end select else call cmue_d(2) end if fp=cmue1(1)*an(1)/Rii-cmue2(1)*an(1)-cm0**(-3) step=-fc/((fp-fc)/e) ann=ann+0.5*step if (abs(step).gt.100.) then STDERR 'Method for calculating the steady-state Richardson' STDERR 'number Rist does not converge.' STDERR 'Probably, the prescribed value for c3' STDERR 'is outside the range of the chosen stability function.' STDERR 'Please change gotmturb.nml accordingly.' STDERR 'You have chosen the stability function no.',stab_method STDERR 'If the problem persists, then use another' STDERR 'stability function.' STDERR 'I will continue anyway. Good luck!' STDERR ' ' converged=.false. goto 333 endif if (abs(step).lt.1.e-10) goto 222 end do 222 an(1)=ann as(1)=an(1)/Rii if (turb_method.eq.2) then select case(stab_method) case(Constant) cmue1=cm0_fix cmue2=cm0_fix/Prandtl0_fix case(MunkAnderson) call cmue_ma(2) case(SchumGerz) call cmue_sg(2) case(EiflerSchrimpf) call cmue_rf(2) end select else call cmue_d(2) end if cc3=c2+(c1-c2)/Rii*cmue1(1)/cmue2(1) ffp=cc3-c3 step=-ffc/((ffp-ffc)/ee) Ri=Ri+0.25*step if (abs(step).gt.100.) then STDERR ' ' STDERR 'Method for calculating the steady-state Richardson' STDERR 'number Rist does not converge.' STDERR 'Probably, the prescribed value for c3' STDERR 'is outside the range of the chosen stability function.' STDERR 'Please change gotmturb.nml accordingly.' STDERR 'You have chosen the stability function no.',stab_method STDERR 'If the problem persists, then use another' STDERR 'stability function.' STDERR 'I will continue anyway. Good luck!' STDERR ' ' converged=.false. goto 333 endif if (abs(step).lt.1.e-10) goto 333 end do 333 compute_rist=Ri if (.not.converged) compute_rist=-999. return end function compute_rist !EOC !----------------------------------------------------------------------- ! Copyright by the GOTM-team under the GNU Public License - www.gnu.org !-----------------------------------------------------------------------