~felix-schueller/+junk/trunk

« back to all changes in this revision

Viewing changes to upslope.f90

Committer: Felix Schueller
Date: 2008-04-22 12:24:27 UTC
Revision ID: felix.schueller@gmail.com-20080422122427-60f0eikhqbxipces

FSS: moved calculation of the mterm to a subroutine (as proposed by Idar)

files modified:
upslope.f90

Show diffs side-by-side

added added

removed removed

upslope.f90

!! -OUTPUT:

!! -DESCRIPTION: This module contains the routine containing the model itself.

!! -TODO: write idealized topo also to netcdf

! -Last modified: Mon Apr 21, 2008 16:56

! -Last modified: Tue Apr 22, 2008 14:22

!----------------------------------------------------------

use nc_io

use utils_mod

public :: upslope !FSS: the main routine

private :: m_term

private :: sum_write

private :: update

private :: bilinear

!! -OUTPUT:

!! -DESCRIPTION:

!! -TODO:

!! -Last modified: Mon Apr 21, 2008 16:56

!! -Last modified: Tue Apr 22, 2008 14:22

!----------------------------------------------------------

implicit none

!** Variables and constants used in computations

356

354

!* DYNAMICS:

357

355

!** this section take into account the dynamical effect of upstream lifting

358

356

359

select case (dynamics)

360

361

362

!** hydrostatic version:

363

364

case('hydro')

365

366

denom=sigma**2 -coriolis**2

367

if(abs(real(denom)).lt.1E-18)denom=sign(1.,denom)*1E-18

368

mtermsq = (nterm**2/denom) * (kx**2 + ky**2)

369

370

!** propagating case:

371

if(mtermsq.ge.0.)then

372

mterm = cmplx(sign(sqrt(mtermsq),sigma),0.)

373

!** envanecent case:

374

else

375

mterm = cmplx(0.,sqrt(-mtermsq ))

376

endif

377

378

case('nonhydro')

379

380

!** here we need two branches since NH-version doesn't ensure real sqrt(m);

381

!** non-hydrostatic vertical wave number:

382

383

denom=sigma**2 -coriolis**2

384

385

if(abs(real(denom)).lt.1E-18)then

386

denom=sign(1.,denom)*1E-18

387

if(nterm.eq.0. .and. coriolis.eq.0.)then !pot.flow case

388

mtermsq = - (kx**2 + ky**2)

389

goto 233

390

endif

391

endif

392

mtermsq = ((nterm**2-sigma**2)/denom) * (kx**2 + ky**2)

393

394

233 continue

395

396

!** propagating case:

397

if(mtermsq.ge.0.)then

398

mterm = cmplx(sign(sqrt(mtermsq),sigma),0.)

399

!** envanecent case:

400

else

401

mterm = cmplx(0.,sqrt(-mtermsq ))

402

endif

403

404

case('nodyn')

405

406

mterm = 0.

407

408

end select

357

call m_term(mterm,sigma,kx,ky,nterm,coriolis,dynamics)

409

358

410

359

!* MICRO-PHYSICAL TIME SCALE

411

360

706

655

707

656

708

657

658

subroutine m_term(mterm,sigma,kx,ky,nterm,coriolis,dynamics)

659

!---------written by Idar Barstad-----------------

660

!! -INPUT: sigma :

661

!! kx,ky :

662

!! nterm :

663

!! coriolis :

664

!! dynamics :

665

!! -OUTPUT: mterm :

666

!! -DESCRIPTION: This routine will calculate the m-term (which is complex)

667

!! -TODO:

668

!! -Last modified:

669

!! modified by FSS

670

!----------------------------------------------------------

671

implicit none

672

673

!** Variables and constants used in computations

674

complex, intent(out) :: mterm

675

real, intent(in) :: sigma

676

real, intent(in) :: kx,ky

677

real, intent(in) :: nterm

678

real, intent(in) :: coriolis

679

character(len=8), intent(in) :: dynamics

680

real :: denom

681

real :: mtermsq

682

complex :: ci

683

684

685

select case (dynamics)

686

687

!** HOW TO CHOOSE THE ROOT OF m-term:

688

! when interpreting what root of m and what sign it should have,

689

! we look at the expression for displacement in Fourier space:

690

! eta = terrain * exp(imz) and interpret this physically.

691

! If this expression is evanescent, we need exp(-|m| z),

692

! if this is trignometric, we need exp(i |m| z), using the sign(sigma)

693

! to make sure that the wave propagate upward.

694

695

696

!** hydrostatic version:

697

698

case('hydro')

699

700

denom = sigma**2 -coriolis**2

701

702

if(abs(denom).lt.1E-18) denom = sign(1.,denom)*1E-18

703

704

mtermsq = (nterm**2/denom) * (kx**2 + ky**2)

705

706

707

! ** propagating case:

708

! ** ( you need the real part of the root and you will take the

709

! ** root that has the same sign as sigma, sign(sigma), to get

710

! ** the correct direction of the phase lines (upstream tilt) )

711

712

if(mtermsq.ge.0.)then

713

714

mterm = cmplx(sign(sqrt(mtermsq),sigma),0.)

715

716

else

717

718

! ** envanecent case:

719

! ** ( you need the imaginary part, and the root having the

720

! ** positive sign, which will become negative in eta-equation )

721

722

mterm = cmplx(0.,sqrt(-mtermsq ))

723

724

endif

725

726

727

728

729

!** Non-hydrostatic case:

730

731

case('nonhydro')

732

733

!** here we need two branches since NH-version doesn't ensure real sqrt(m);

734

!** non-hydrostatic vertical wave number:

735

736

denom=sigma**2 -coriolis**2

737

738

if(abs(denom).lt.1E-18)then

739

denom=sign(1.,denom)*1E-18

740

741

! ** potential flow case:

742

if(nterm.eq.0. .and. coriolis.eq.0.)then

743

mtermsq = - (kx**2 + ky**2)

744

goto 233

745

endif

746

747

endif

748

749

mtermsq = (( nterm**2 - sigma**2 ) / denom) * (kx**2 + ky**2)

750

751

233 continue

752

753

! ** propagating case:

754

! (non-dispersive waves)

755

! (you need the real part of the positive root)

756

757

if(mtermsq.ge.0.)then

758

mterm = cmplx(sign(sqrt(mtermsq),sigma),0.)

759

760

761

! ** envanecent case:

762

! ( you need imaginary part of the positive root )

763

else

764

765

mterm = cmplx(0.,sqrt(-mtermsq ))

766

767

endif

768

769

case('nodyn')

770

771

mterm = 0.

772

773

end select

774

end subroutine m_term

775

709

776

!***************************************

710

777

!*********** update you case-variables:

711

778

!***************************************

Older »