Don't like ads? PRO users don't see any ads ;-)
Public Pastes
Guest

Untitled

By: a guest on Jun 22nd, 2012  |  syntax: Fortran  |  size: 8.20 KB  |  hits: 16  |  expires: Never
download  |  raw  |  embed  |  report abuse  |  print
Text below is selected. Please press Ctrl+C to copy to your clipboard. (⌘+C on Mac)
  1. PROGRAM uebung
  2. !*********************************************************************
  3. !
  4. ! Description:
  5. !  Numerics exercise extracted from V. Wirth "Das Eliassen-Problem:
  6. !   Reaktion einer balancierten Stroemung auf eine schwache externe
  7. !   Stoerung. Ein numerischer Praktikumsversuch", Universitaet Mainz  
  8. !
  9. !   Full exercise:
  10. !   www.staff.uni-mainz.de/vwirth/ps/projekt_eliassen.ps
  11. !
  12. !   There are several publications that are based on the original code (and I think that
  13. !   Prof. Wirth would be glad to provide it in case you are interested in using
  14. !   it in your future research):
  15. !
  16. !   Wirth, V. Thermal versus dynamical tropopause in upper tropospheric
  17. !    balanced flow anomalies Quart. J.~Roy. Met. Soc., 2000, 126, 299-317
  18. !
  19. !   Wirth, V. Cyclone-Anticyclone Asymmetry Concerning the Height of the Thermal
  20. !    and the Dynamical Tropopause J.~Atmos. Sci., 2001, 58, 26-37
  21. !
  22. !   Wirth, V. A Dynamical Mechanism for Tropopause Sharpening,
  23. !     Meteorol. Zeitschrift, N.F. , 2004, 13, 477-484
  24. !
  25. ! Method: define idealized jet as base state                                      
  26. !         calculate ageostophic mass stream function due to weak Rayleigh damping    
  27. !
  28. ! Output Files:
  29. !     ../run/out.nc
  30. !
  31. ! History:
  32. !     Student version
  33. !    
  34. !**********************************************************************
  35.  
  36. ! mo_output provides a simplified interface for writing netCDF data
  37. !     see supplement.pdf for details          
  38.   USE mo_output,   ONLY:  open_out_file,   &
  39.                           init_out,        &
  40.                           write_out,       &
  41.                           close_out_file
  42.  
  43. ! subroutine to initialize the zonal wind field  
  44.   USE mo_jet, ONLY : init_jet
  45.  
  46.   IMPLICIT NONE
  47.  
  48. ! 1. variable declaration (parameters) -----------------------------------
  49.     ! number of grid-points in y-direction
  50.    INTEGER, PARAMETER :: jmax = 80
  51.     ! number of grid-points in z-direction
  52.    INTEGER, PARAMETER ::  kmax = 80
  53.     ! horizontal domain length in km
  54.    DOUBLE PRECISION, PARAMETER :: lykm = 1000.D0
  55.     ! vertical domain length in km
  56.    DOUBLE PRECISION,  PARAMETER :: lzkm = 10.D0
  57.     ! relaxation time in days
  58.    DOUBLE PRECISION, PARAMETER :: tau  = 5.D0
  59.     ! Coriolis parameter
  60.    DOUBLE PRECISION, PARAMETER :: fcori = 1.D-4
  61.     ! Bouyancy frequency
  62.    DOUBLE PRECISION, PARAMETER :: Nbouy = 1.D-2
  63.     ! density
  64.     DOUBLE PRECISION, PARAMETER :: rho0=1.0D0
  65.  
  66. ! 2. other variables -----------------------------------
  67.    DOUBLE PRECISION, POINTER :: deltay, deltaz ! grid spacings in y- and z-direction
  68.                                                ! (pointer attribute needed for netCDF interface)
  69.     ! stream function
  70.    DOUBLE PRECISION, POINTER, DIMENSION(:,:) :: phi
  71.  
  72.     ! zonal wind
  73.    DOUBLE PRECISION, POINTER, DIMENSION(:,:) :: uwnd
  74.  
  75.     ! damping term  
  76.    DOUBLE PRECISION, POINTER, DIMENSION(:,:) :: xdamp
  77.    
  78.     ! partial derivative w.r.t z of x
  79.    DOUBLE PRECISION, POINTER, DIMENSION(:,:) :: dxodz
  80.  
  81.     ! right hand side of Poisson Eq.
  82.    DOUBLE PRECISION, POINTER, DIMENSION(:,:) :: rhs
  83.  
  84.     ! coefficients in Poisson Eq.
  85.     DOUBLE PRECISION :: acoeff, bcoeff
  86.  
  87.     ! grid spacings in m
  88.    DOUBLE PRECISION :: ly, lz
  89.  
  90.     ! netcdf output unit                          
  91.    INTEGER :: ncid  
  92.  
  93.     ! for choice of iteration style
  94.    CHARACTER (len=1) :: choice  
  95. !-end variable declaration --------------------------------
  96.  
  97.     ! write to screen:
  98.    WRITE(6,*) 'This is uebung.exe'
  99.  
  100.    ALLOCATE(deltay)
  101.    ALLOCATE(deltaz)
  102.    ALLOCATE(phi(jmax,kmax))
  103.    ALLOCATE(uwnd(jmax,kmax))
  104.    ALLOCATE(xdamp(jmax,kmax))
  105.    ALLOCATE(dxodz(jmax,kmax))
  106.    ALLOCATE(rhs(jmax,kmax))
  107.  
  108.     ! convert to km
  109.    ly=1.d3*lykm
  110.    lz=1.d3*lzkm
  111.            
  112.     ! grid spacings (m)
  113.    deltay = ly/DBLE((jmax-1))
  114.    deltaz = lz/DBLE((kmax-1))
  115.  
  116.     ! coefficients in Poisson eq.
  117.    acoeff = Nbouy**2 / rho0  
  118.    bcoeff = fcori**2 / rho0
  119.  
  120.  
  121.  
  122.    ! initialize phi
  123.     phi=0.D0
  124.  
  125.    ! initalize the zonal wind field
  126.    CALL init_jet(uwnd, jmax, kmax, deltay, deltaz)
  127.  
  128.    ! the Rayleigh damping term
  129.    xdamp=-uwnd/(tau*86400.D0)
  130.  
  131.    ! calculate dxodz (partial derivative of X w.r.t. z)
  132.    CALL zderiv(xdamp, dxodz, deltaz, jmax,kmax)
  133.  
  134.    ! right hand side of the Poisson equation
  135.    rhs=fcori*dxodz
  136.  
  137. !===================================================================================!
  138.    ! choosing iterative solver
  139. print*, '========================================================================='
  140. print*, 'To continue, choose and type either:'
  141. print*, '"c" for Jacobi, or "d" for Gauss-Seidel"'
  142. READ(5,*) choice
  143.  
  144. IF (choice .EQ. 'c') THEN
  145.    CALL solve_with_Jacobi(phi, rhs, deltay, deltaz, acoeff, bcoeff, jmax, kmax)
  146. END IF
  147.  
  148. IF (choice .EQ. 'd') THEN
  149.    CALL solve_with_Gauss_Seidel(phi, rhs, deltay, deltaz, acoeff, bcoeff, jmax, kmax)
  150. END IF
  151. !===================================================================================!
  152.  
  153.  
  154.     ! write netCDF output
  155.    CALL open_out_file(ncid, "out.nc")
  156.    CALL init_out(ncid, uwnd, "uwnd", "m/s","y","z")
  157.    CALL init_out(ncid, phi, "phi", "kg -1 m-1)","y","z")
  158.    CALL init_out(ncid, deltay, "dy","m")
  159.    CALL init_out(ncid, deltaz, "dz","m")
  160.    CALL init_out(ncid, dxodz, "dxodz", "s-3","y","z")
  161.    CALL write_out(ncid)
  162.    CALL close_out_file(ncid)
  163.  
  164.    ! deallocate memory
  165.   DEALLOCATE(phi, uwnd, xdamp, dxodz, rhs)
  166.  
  167.   WRITE(6,*) 'uebung.exe: DONE.'
  168.  
  169. END PROGRAM uebung
  170.  
  171. SUBROUTINE zderiv( f, der, deltaz, jmax, kmax)
  172.  
  173.   IMPLICIT NONE
  174.  
  175.   DOUBLE PRECISION, INTENT(IN ) :: deltaz
  176.   INTEGER, INTENT(IN ) :: jmax, kmax
  177.   DOUBLE PRECISION, DIMENSION( jmax , kmax ) :: f, der
  178.  
  179.   INTEGER :: j,k
  180.  
  181.   DO j = 1,jmax
  182.     DO k = 2, kmax-1
  183.       der(j,k) = (f(j,k+1) - f(j,k-1))/(2.D0*deltaz)
  184.     END DO
  185.       der(j,1) = (f(j,2) - 0.D0)/(2.D0*deltaz)
  186.       der(j,kmax) = (0.D0 - f(j,kmax-1))/(2.D0*deltaz)
  187.   END DO
  188.  
  189. END SUBROUTINE zderiv
  190.  
  191. !================================MY CODE===============================================!
  192.  
  193. SUBROUTINE solve_with_Jacobi(ufct, rhs, deltax, deltay, acoeff, bcoeff, jmax, kmax )
  194.  
  195.   IMPLICIT NONE
  196.  
  197.   INTEGER, INTENT(IN ) :: jmax, kmax
  198.   DOUBLE PRECISION, INTENT(IN ) :: deltax, deltay, acoeff, bcoeff
  199.   DOUBLE PRECISION, DIMENSION( jmax, kmax ) :: ufct, rhs, pfct
  200.   INTEGER :: j, k, m, n
  201.   DOUBLE PRECISION, PARAMETER :: c=0.5d0
  202.  
  203.  
  204.  
  205. DO WHILE (1336 .lt. 1337)
  206.   pfct=ufct
  207.  
  208.   DO k=2, kmax-1
  209.     ufct(1,k)=(acoeff/(deltax**2.d0)*0.D0+pfct(2,k))+&
  210.                bcoeff/(deltay**2.d0)*(pfct(1,k-1)+pfct(1,k+1))-&
  211.                rhs(1,k))/&
  212.                (2.d0*(acoeff/deltax**2+bcoeff/deltay**2))
  213.   END DO
  214.  
  215.   ufct(
  216.  
  217.   DO j=2,jmax-1
  218.  
  219.     DO k=2, kmax-1
  220.       ufct(j,k)=(acoeff/(deltax**2.d0)*(pfct(j-1,k)+pfct(j+1,k))+&
  221.                  bcoeff/(deltay**2.d0)*(pfct(j,k-1)+pfct(j,k+1))-&
  222.                  rhs(j,k))/&
  223.                  (2.d0*(acoeff/deltax**2+bcoeff/deltay**2))
  224.     END DO
  225.   END DO
  226. n=n+1
  227. IF (c .GT. MAXVAL(ufct-pfct)) EXIT
  228. END DO
  229.  
  230. print*,'Anzahl der Durchlaeufe mit Jacobi:', n
  231. print*, '========================================================================='
  232.  
  233. END SUBROUTINE solve_with_Jacobi
  234.  
  235. SUBROUTINE solve_with_Gauss_Seidel(ufct, rhs, deltax, deltay, acoeff, bcoeff, jmax, kmax )
  236.   IMPLICIT NONE
  237.  
  238.   INTEGER, INTENT(IN ) :: jmax, kmax
  239.   DOUBLE PRECISION, INTENT(IN ) :: deltax, deltay, acoeff, bcoeff
  240.   DOUBLE PRECISION, DIMENSION( jmax, kmax ) :: ufct, rhs, pfct
  241.   INTEGER :: j, k, m, n
  242.   DOUBLE PRECISION, PARAMETER :: c=1.d-10
  243.  
  244.  
  245.  
  246. DO WHILE (1336 .lt. 1337)
  247.     pfct=ufct
  248.   ufct(1,1)=(acoeff/(deltax**2.d0)*(pfct(1-1,1)+pfct(1+1,1))+&
  249.              bcoeff/(deltay**2.d0)*(pfct(1,1-1)+pfct(1,1+1))-&
  250.              rhs(1,1))/&
  251.             (2.d0*(acoeff/deltax**2+bcoeff/deltay**2))
  252.  
  253.   DO j=2,jmax
  254.     DO k=2, kmax
  255.       ufct(j,k)=(acoeff/(deltax**2.d0)*(ufct(j-1,k)+pfct(j+1,k))+&
  256.                  bcoeff/(deltay**2.d0)*(ufct(j,k-1)+pfct(j,k+1))-&
  257.                  rhs(j,k))/&
  258.                  (2.d0*(acoeff/deltax**2+bcoeff/deltay**2))
  259.     END DO
  260.   END DO
  261.  
  262. n=n+1
  263. IF (c .GT. MAXVAL(ufct-pfct)) EXIT
  264. END DO
  265.  
  266. print*,'Anzahl der Durchlaeufe mit Gauss-Seidel:', n
  267. print*, '========================================================================='
  268.  
  269.  
  270. END SUBROUTINE solve_with_Gauss_Seidel
create a new version of this paste RAW Paste Data