vdif_ts_gom.f90

Go to the documentation of this file.

!/===========================================================================/
! Copyright (c) 2007, The University of Massachusetts Dartmouth 
! Produced at the School of Marine Science & Technology 
! Marine Ecosystem Dynamics Modeling group
! All rights reserved.
!
! FVCOM has been developed by the joint UMASSD-WHOI research team. For 
! details of authorship and attribution of credit please see the FVCOM
! technical manual or contact the MEDM group.
!
! 
! This file is part of FVCOM. For details, see http://fvcom.smast.umassd.edu 
! The full copyright notice is contained in the file COPYRIGHT located in the 
! root directory of the FVCOM code. This original header must be maintained
! in all distributed versions.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 
! AND ANY EXPRESS OR  IMPLIED WARRANTIES, INCLUDING,  BUT NOT  LIMITED TO,
! THE IMPLIED WARRANTIES OF MERCHANTABILITY AND  FITNESS FOR A PARTICULAR
! PURPOSE ARE DISCLAIMED.  
!
!/---------------------------------------------------------------------------/
! CVS VERSION INFORMATION
! $Id$
! $Name$
! $Revision$
!/===========================================================================/

!==============================================================================|
!     this subroutine is used to calculate the true temperature                !
!     and salinity by including vertical diffusion implicitly.                 !
!==============================================================================|
!     See NOTES in VDIF_TS.F about possible optimization!

! REMOVED DEPTH CHECK FOR NONE WET DRY CASE. THIS IS A LEGACY FROM
! ECOM-SI WHICH HAS LAND IN THE DOMAIN. IT IS NOT NEEDED IN FVCOM

SUBROUTINE vdif_ts_gom(NCON1,F)                

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

  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


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

  umolpr = umol*1.e0_sp
  sasurf = 0.0_sp
  !
  !------------------------------------------------------------------------------!
  !                                                                              !
  !        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)
           !     END DO
           wfsurf(1:m)=wtsurf(1:m)
           swradf(1:m)=swrad(1:m)


           DO k = 1, kb
              DO i = 1, m
                 IF(iswetn(i) == 1)THEN
                    zdep = z(i,k)*d(i)

                    !       SPECIFIED PARAMETER FOR GOM CASE (JUNE 18, 2003)
                    IF(d(i) < 100.0_sp) THEN
                       rheat=0.78_sp
                       zeta1=1.4_sp
                       zeta2=6.3_sp
                    ELSE IF(d(i) > 150.0_sp) THEN
                       rheat=0.58_sp
                       zeta1=0.35_sp
                       zeta2=23.0_sp
                    ELSE
                       tmp1=(d(i)-100.0_sp)/(150.0_sp-100.0_sp)
                       rheat=0.78_sp+(0.58_sp-0.78_sp)*tmp1
                       zeta1=1.4_sp+ (0.35_sp-1.4_sp )*tmp1
                       zeta2=6.3_sp+ (23.0_sp-6.3_sp )*tmp1
                    END IF

                    rad(i,k)=swradf(i)*(rheat*exp(zdep/zeta1)+(1-rheat)*exp(zdep/zeta2))
                 END IF

              END DO
           END DO

        ELSE IF(heating_type == 'body') THEN !! H_TYPE = 'body_h'   
           rad    = 0.0_sp
           wfsurf = 0.0_sp
           swradf = 0.0_sp
        ELSE IF(heating_on) THEN
           CALL fatal_error("HEATING TYPE IS SET INCORRECTLY:",&
                & trim(heating_type) )
        END IF

     CASE(2)

        swradf= 0.0_sp
        wfsurf=0.0_sp
        cosgama1=1.0_sp
           !       SASURF(I)=-(QEVAP3(I)-QPREC3(I))*F(I,1)*ROFVROS*COSGAMA1(I)
           !---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_GOM')
     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)
!!$#  if defined (1)
!!$                    END IF
!!$# endif
!!$                 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 .AND.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
       tmp = -tmp
     CASE DEFAULT
       CALL fatal_error('NCON NOT CORRECT IN VDIF_TS_GOM')
     END SELECT
   ELSE    
           tmp = 0.0_sp
   END IF      
!!$              IF (TMP1 > 0.0_SP) TMP=0.0_SP
!!$              !       !also try TMP = 0.0_SP
              !----------------------------------------------------------------------
!!$ THIS IS FOR AN OLDER VERSION OF GROUNDWATER
!!$              GW=0.0_SP
!!$              IF(NCON1 == 2) 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
!!$                 IF (BFWDIS(I) /= 0.0_SP) TMP=0.0_SP
!!$
!!$              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 .AND.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_gom"
     
   END SUBROUTINE vdif_ts_gom
!==============================================================================|