vdif_ts.f90 File Reference

Go to the source code of this file.

Functions/Subroutines
subroutine	vdif_ts (NCON1, F)

Function/Subroutine Documentation

vdif_ts()

subroutine vdif_ts	(	integer	NCON1,
		real(sp), dimension(0:mt,kb)	F
	)

Definition at line 54 of file vdif_ts.f90.

  !------------------------------------------------------------------------------|

  USE all_vars
  USE mod_utils
  USE bcs
  USE mod_wd


  IMPLICIT NONE
  INTEGER :: I,K,NCON1,J,KI,hutemp
  !   REAL(SP) :: TMP,TMP1,TMP2,TMP3,QTMP,GW,ZDEP,FKH,UMOLPR,WETFAC
  !   REAL(SP), DIMENSION(0:MT,KB)  :: F
  !   REAL(SP), DIMENSION(M,KB)     :: FF,AF,CF,VHF,VHPF,RAD
  !   REAL(SP), DIMENSION(M)        :: KHBOTTOM,WFSURF,SWRADF
  REAL(DP) :: TMP,TMP1,TMP2,TMP3,QTMP,GW,ZDEP,FKH,UMOLPR,WETFAC
  REAL(SP), DIMENSION(0:MT,KB)  :: F
  REAL(DP), DIMENSION(M,KB)     :: FF,AF,CF,VHF,VHPF,RAD
  REAL(DP), DIMENSION(M)        :: KHBOTTOM,WFSURF,SWRADF,SASURF
  REAL(DP), DIMENSION(M)        :: COSGAMA1,COSGAMA2


  REAL(SP) :: WFTMP1, WFTMP2, WFTMP3

  IF(dbg_set(dbg_sbr)) WRITE(ipt,*)"Start: vdif_ts :", ncon1

  umolpr = umol*1.e0_sp
  sasurf = 0.0_sp
  gw = 0.0_dp
  !
  !------------------------------------------------------------------------------!
  !                                                                              !
  !        the following section solves the equation                             !
  !         dti*(kh*f')'-f=-fb                                                   !
  !                                                                              !
  !------------------------------------------------------------------------------!

  DO k = 2, kbm1
     DO i = 1, m
        IF(iswetn(i) == 1)THEN
           fkh = kh(i,k)

           IF(k == kbm1) THEN
              khbottom(i)=fkh
           END IF

           af(i,k-1)=-dti*(fkh+umolpr)/(dz(i,k-1)*dzz(i,k-1)*d(i)*d(i))
           cf(i,k)=-dti*(fkh+umolpr)/(dz(i,k)*dzz(i,k-1)*d(i)*d(i))
        END IF
     END DO
  END DO


  !------------------------------------------------------------------------------!
  !     the net heat flux input.                                                 !
  !     the method shown below can be used when we turn off the                  !
  !     body force in subroutine advt. be sure this method could                 !
  !     cause the surface overheated if the vertical resolution                  !
  !     is not high enough.                                                      !
  !------------------------------------------------------------------------------!

  SELECT CASE(ncon1)
  CASE(1)

     rad    = 0.0_sp
     wfsurf = 0.0_sp
     swradf = 0.0_sp
     IF(heating_type == 'flux')THEN
        !     DO I=1,M
        !       WFSURF(I)=WTSURF(I)
        !       SWRADF(I)=SWRAD(I)
        !       SASURF(I)=0.0_SP
        !     END DO
        wfsurf(1:m)=wtsurf(1:m)
        swradf(1:m)=swrad(1:m)
        sasurf=0.0_sp

!
!-----------------------------------------------------------------------
!  If net heat flux is cooling and SST is at freezing point or below
!  then suppress further cooling. Note: net heat flux sign convention is that 
!  positive (in FVCOM negitive - J. Qi) means heating the ocean (J Wilkin - ROMS).
!-----------------------------------------------------------------------
! 
!  Below the surface heat flux WFSURF is ZERO if cooling AND
!  the SST is cooler than the threshold.  The value is retained if
!  warming.
!
!    WFTMP3 = 0      if SST warmer than threshold (WFTMP1) - change nothing
!    WFTMP3 = 1      if SST colder than threshold (WFTMP1)
!
!    0.5*(WFTMP2+ABS(WFTMP2)) = 0                        if flux is warming
!                             = WFSURF                   if flux is cooling 
!
        wftmp1=-2.0_sp              ! nominal SST threshold to cease cooling
        DO i=1,m
          wftmp2=wfsurf(i)
          wftmp3=0.5_sp*(1.0_sp+sign(1.0_sp,wftmp1-f(i,1)))
          wfsurf(i) = wftmp2-wftmp3*0.5_sp*(wftmp2+abs(wftmp2))
        END DO

        DO k = 1, kb
           DO i = 1, m
              IF(iswetn(i) == 1)THEN
                 zdep = z(i,k)*d(i)
                 rad(i,k)=swradf(i)*(rheat*exp(zdep/zeta1)+(1-rheat)*exp(zdep/zeta2))
              END IF
           END DO
        END DO
     END IF

  CASE(2)

     swradf   = 0.0_sp
     wfsurf   =0.0_sp
     cosgama1 =1.0_sp
!---Considering the salinity conservation, the sea surface salinity flux-----!
!---can be set as zero, that is----------------------------------------------!
     sasurf   = 0.0_sp
     rad      =0.0_sp

  CASE DEFAULT
     CALL fatal_error('NCON NOT CORRECT IN VDIF_TS')
  END SELECT


  !------------------------------------------------------------------------------!
  !   surface bcs; wfsurf                                                        !
  !------------------------------------------------------------------------------!

  DO i = 1, m
     IF(iswetn(i) == 1)THEN
        vhf(i,1) = af(i,1) / (af(i,1)-1.)
        vhpf(i,1) = -dti *(sasurf(i)+wfsurf(i)-swradf(i) &
             +rad(i,1)-rad(i,2)) / (-dz(i,1)*d(i)) - f(i,1)
        vhpf(i,1) = vhpf(i,1) / (af(i,1)-1.)
     END IF
  END DO

  DO k = 2, kbm2
     DO i = 1, m
        IF(iswetn(i) == 1)THEN
           vhpf(i,k)=1./ (af(i,k)+cf(i,k)*(1.-vhf(i,k-1))-1.)
           vhf(i,k) = af(i,k) * vhpf(i,k)
           vhpf(i,k) = (cf(i,k)*vhpf(i,k-1)-dble(f(i,k)) &
                +dti*(rad(i,k)-rad(i,k+1))/(d(i)*dz(i,k)))*vhpf(i,k)
        END IF
     END DO
  END DO


!!$   DO  K = 1, KBM1
!!$     DO  I = 1, M
!!$#  if !defined (1)
!!$     IF (D(I) > 0.0_SP) THEN
!!$#  else
!!$     IF(ISWETN(I) == 1)THEN
!!$#  endif
!!$         FF(I,K) = F(I,K)
!!$       END IF
!!$     END DO
!!$   END DO


  DO  k = 1, kbm1
     DO  i = 1, m
        IF(iswetn(i) == 1)THEN
           ff(i,k) = f(i,k)
        END IF
     END DO
  END DO

  ! THIS PIECE OF CODE DESPERATELY NEEDS TO BE CLARIFIED
  ! AND STREAMLINED

  DO i = 1, m
     IF (isonb(i) /= 2) THEN
        IF(iswetn(i) == 1)THEN
           tmp1=pfpxb(i)*cos(sita_gd(i))+pfpyb(i)*sin(sita_gd(i))
           tmp2=ah_bottom(i)*phpn(i)
           tmp3=khbottom(i)+umolpr+ah_bottom(i)*phpn(i)*phpn(i)
           tmp=tmp1*tmp2/tmp3*(khbottom(i)+umolpr)

           ! -------------------------------------------------------------------
   IF(noflux_bot_condition)THEN
     SELECT CASE(ncon1)
     CASE(1)
       IF (tmp1 > 0.0_sp) tmp=0.0_sp
     CASE(2)
!!$       IF (TMP1 > 0.0_SP) TMP=0.0_SP
       IF (tmp1 < 0.0_sp) tmp=0.0_sp
       tmp = -tmp
     CASE DEFAULT
       CALL fatal_error('NCON NOT CORRECT IN VDIF_TS')
     END SELECT
   ELSE    
           tmp = 0.0_sp
   END IF      
           ! -------------------------------------------------------------------

!!$ THIS IS FOR AN OLDER VERSION OF GROUNDWATER
!!$           GW=0.0_SP
!!$           IF(NCON1 == 2) THEN ! MODIFIED FOR NEW BFWDIS
!!$              IF (BFWDIS(I) /= 0.0_SP) THEN
!!$
!!$                 QTMP=-(F(I,KBM1)*D(I)*DZ(I,KBM1)*BFWDIS(I))/ &
!!$                      (D(I)*DZ(I,KBM1)*ART1(I)+BFWDIS(I))
!!$                 GW=DTI/D(I)/DZ(I,KBM1)*QTMP
!!$                 
!!$                 TMP=0.0_SP
!!$
!!$              END IF
!!$           END IF


           ff(i,kbm1) = (cf(i,kbm1)*vhpf(i,kbm2)-ff(i,kbm1)-gw+dti*(rad(i,kbm1)-rad(i,kb)-tmp)/(d(i)*dz(i,kbm1))) &
                /(cf(i,kbm1)*(1._sp-vhf(i,kbm2))-1._sp)


        END IF
     END IF
  END DO

  DO  k = 2, kbm1
     ki = kb - k
     DO  i = 1, m
        IF(isonb(i) /= 2) THEN
           IF(iswetn(i) == 1)THEN
              ff(i,ki) = (vhf(i,ki)*ff(i,ki+1)+vhpf(i,ki))
           END IF
        END IF


     END DO
  END DO

  DO i = 1, m
     IF(iswetn(i)*iswetnt(i) == 1 )then
        DO k = 1, kbm1
           f(i,k) = ff(i,k)
        END DO
     END IF
  END DO


  IF(dbg_set(dbg_sbr)) WRITE(ipt,*)"End: vdif_ts"
  RETURN

Here is the call graph for this function:

Here is the caller graph for this function: