source: trunk/00FlowSolve_PL/PROJECTS/NEW_SRC/set_ambient_stratification.f90

!============================================================================
! SET_AMBIENT_STRATIFICATION
! FUNCTION TO DEFINE s1_bar [K] & s2_bar [ppt]
subroutine set_ambient_stratification

  ! routine to set the ambient profiles for the scalars s1 and s2
  ! where total_s1 = s1_bar(x,z) + s1(x,y,z,t)
  ! and   total_s2 = s2_bar(x,z) + s2(x,y,z,t)

  use grid_info
  use user_parameters       ! for z_N,z_N_scale,Smax,Smin,Tmax,Tmin,dTdz
  use dependent_variables   ! s1_bar, s2_bar,rho_bar
  use dimensional_scales    ! s1_scale, s2_scale
  use counters_flags_etc
  use mpi_parameters
  implicit none
  include '../input/problem_size.h'
  include 'mpif.h'

!!! local variables
  real        :: xpos,zpos,Temp,Salinity
  integer     :: kglobal
  real,dimension(:,:),allocatable :: tabforcing


  print*,'Ambient_stratification, Lz',Lz 

  if ( pomme_data_flag == 'no') then


     !!! these are just values pulled out of a hat for testing
     Smax =  35.0          ! [psu]
     Smin =  35.0          ! [psu]
     Tmax =  15.0          ! [deg C]
     Tmin =  10.0          ! [deg C]   !!! must be the same as T0 in EOS.f90

     !!! defined in params.h
     !!! dTdz=1.4133829e-3*4

     !!! we'll use s1 for temperature (deg C) and s2 for salinity (g/kg or ppt)
     ! both of these test profiles are stable
     do k=1,locnz
        do i=1,nx

           zpos = Grid(0)%z(i,k)*L_scale     ! [m]
           Salinity= Smin + (Smax-Smin)/Lz * zpos
           Temp    = Tmin + dTdz * zpos             !!!!!(Tmax-Tmin)/Lz * zpos

           s1_bar(i,k) = Temp / s1_scale       
           s2_bar(i,k) = Salinity / s2_scale

        enddo
     enddo

  else

     !!! xlv lecture profil vertical [T,S] POMME
     !============================================
     allocate (tabforcing(nz,2))

     if(myid==0) write(0,*) '... Reading POMME data into file : '
     if(myid==0) print*,'   --> ',trim(pomme_filename)
     open(911,file=pomme_filename,status='old')
     do k=1,nz
        read(911,*) tabforcing(k,1),tabforcing(k,2)
     enddo
     close(911)

     do k=1,locnz
        do i=1,nx
           kglobal = k + myid*locnz
           s1_bar(i,k)=tabforcing(kglobal,1)
           s2_bar(i,k)=tabforcing(kglobal,2)
        enddo
     enddo

     s1_bar = s1_bar / s1_scale  
     s2_bar = s2_bar / s2_scale 

     if(myid==0) write(0,*) '... POMME data saved....'

     deallocate( tabforcing )
     !===================================

  endif


  print*,'MINVAL(T)=',minval(s1_bar),'MAXVAL(T)=',maxval(s1_bar)
  print*,'MINVAL(S)=',minval(s2_bar),'MAXVAL(S)=',maxval(s2_bar)



  return
end subroutine set_ambient_stratification
!============================================================================



!=============================================================================
!  TEMPD
!  FUNCTION TO GET INTERPOLATED VALUE OF TEMP FROM GRIDDED DATA
!
      REAL FUNCTION TEMPD(ZK)
!
      PARAMETER(NT=23)
!
      REAL ZK
!
      REAL ZD(NT),TD(NT)
!
      DATA IZTOPT/NT/
!
!!!      DATA ZD/-60.,-55.,-50.,-47.2,-45.,-44.3,-42.2,-40.,-38.5,-35., &
!!!      -34.5,-32.5,-30.,-28.9,-27.4,-25.,-24.1,-20.,-17.6,-15.,-10., &
!!!      -5.,0./
      DATA ZD/0.,5.,10.,12.8,15.,15.7,17.8,20.,21.5,25.,25.5,27.5,30.,31.1, &
              32.6,35.,35.9,40.,42.4,45.,50.,55.,60./
      DATA TD/-1.15,-1.15,-1.15,-1.05,-1.16,-1.18,-0.92,-0.73,-0.55, &
      -0.05,-0.02,-0.10,-0.28,-0.45,-0.70,-1.26,-1.43,-1.58,-1.61, &
      -1.61,-1.62,-1.62,-1.62/
!
!!!      ZD(:) = 60. + ZD(:)
      TEMPD = FK(ZK,ZD,TD,NT,IZTOPT)
!
      RETURN
      END
!==============================================================================



!============================================================================
!  SALTD
!  FUNCTION TO GET INTERPOLATED VALUE OF SALINITY FROM GRIDDED DATA
!
      REAL FUNCTION SALTD(ZK)
!
      PARAMETER(NS=20)
!
      REAL ZK
!
      REAL ZD(NS),SD(NS)
!
      DATA IZTOPS/NS/
!
!!!      DATA ZD/-60.,-55.,-50.,-47.1,-45.,-43.3,-40.,-37.9,-35.,-32.7, &
!!!      -30.,-28.9,-27.8,-25.,-24.3,-20.,-15.,-10.,-5.,0./
      DATA ZD/0.,5.,10.,12.9,15.,16.7,20.,22.1,25.,27.3,30.,31.1,32.2, &
              35.,35.7,40.,45.,50.,55.,60./
      DATA SD/31.83,31.83,31.83,31.82,31.78,31.75,31.72,31.68,31.58, &
      31.53,31.44,31.38,31.27,30.96,30.87,30.80,30.79,3*30.78/
!
!!!      ZD(:) = 60. + ZD(:)
      SALTD = FK(ZK,ZD,SD,NS,IZTOPS)
!
      RETURN
      END
!=============================================================================



!==============================================================================
!  FK
!  FK IS A FUNCTION THAT USES LINEAR INTERPOLATION TO FIND  FK,
!  GIVEN  ZK  AND THE ARRAYS  Z  AND  F. FOR EXAMPLE, IN ROUTINE
!  DERIVS,  BVF = FK(ALT,BVFZ,BVFF)  WILL BE USED TO GET  BVF  AT  ALT
!  FROM THE ARRAYS  BVFZ  AND  BVFF.
      REAL FUNCTION FK(ZK,Z,F,NZ,IZTOP)
      REAL Z(*), F(*)
!
!  MAKE SURE Z(IZTOP) .GE. ZK
!!      PRINT 99, zk,nz,iztop
!!   99 FORMAT('zk,nz,iztop = ',1pe15.5,2i15)
    5 CONTINUE
      IF(ZK.GT.Z(IZTOP)) THEN
         IZTOP = IZTOP + 1
         IF(IZTOP.GT.NZ) THEN
            WRITE(2,9) ZK,Z(NZ)
    9       FORMAT(' ZK TOO HIGH. ZK,Z(NZ) = ',2(1PE15.5))
            CLOSE(2)
            CLOSE(3)
            STOP
         ENDIF
         GO TO 5
      ENDIF
!
      IF(ZK.EQ.Z(IZTOP)) THEN
         FK = F(IZTOP)
         RETURN
      ENDIF
!
!  FIND IZK SUCH THAT  Z(IZK) .GT. ZK .GE. Z(IZK-1)
      DO 100 IZ=IZTOP,2,-1
      IZM1 = IZ - 1
      IF(Z(IZ).GT.ZK .AND. ZK.GE.Z(IZM1)) THEN
         ZRATIO = (ZK - Z(IZM1)) / (Z(IZ) - Z(IZM1))
         FK = F(IZM1) + ZRATIO * (F(IZ) - F(IZM1))
         IZTOP = IZM1 + 1
         RETURN
      ENDIF
  100 CONTINUE
!
!  ZK TOO LOW
      WRITE(2,109) ZK,Z(1)
  109 FORMAT(' ZK TOO LOW. ZK,Z(1) = ',2(1PE15.5))
      CLOSE(2)
      CLOSE(3)
      STOP
!
      END
!