source: codes/icosagcm/trunk/src/etat0_dcmip3.f90 @ 48

MODULE etat0_dcmip3_mod

! test cases DCMIP 2012, category 3 : Non-hydrostatic gravity waves

! Questions
! Replace ps0 by preff ??

  USE genmod
  USE dcmip_initial_conditions_test_1_2_3

  PRIVATE
! some of the following should be obtained via getin
!  REAL(rstd),PARAMETER :: ztop=10000.   ! Height position of the model top (m)
!  REAL(rstd),PARAMETER :: ptop=27391.9  ! Pressure at the model top (Pa)
!  REAL(rstd),PARAMETER :: Xscale=1.     ! Small-planet scaling factor (1=Earth)
!  REAL(rstd),PARAMETER :: Omega =0      ! Planetary rotation rate (s-1)
  REAL(rstd),PARAMETER :: u0=20.         ! Maximum amplitude of the zonal wind (m.s-1)
  REAL(rstd),PARAMETER :: zs=0.          ! Surface elevation (m)
  REAL(rstd),PARAMETER :: N=0.01         ! Brunt-Vaisala frequency (s-1)
  REAL(rstd),PARAMETER :: Teq=300.       ! Surface temperature at the equator (K)
  REAL(rstd),PARAMETER :: Peq=1e5        ! Reference surface pressure at the equator (hPa)
  REAL(rstd),PARAMETER :: d=5000.        ! Witdth parameter for theta
  REAL(rstd),PARAMETER :: lonc=2*pi/3    ! Longitudinal centerpoint of theta
  REAL(rstd),PARAMETER :: latc=0         ! Longitudinal centerpoint of theta
  REAL(rstd),PARAMETER :: dtheta=1.      ! Maximum amplitude of theta (K)
  REAL(rstd),PARAMETER :: Lz=20000.      ! Vertical wave lenght of the theta perturbation

  PUBLIC  etat0
CONTAINS


  SUBROUTINE etat0(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
  USE icosa
  IMPLICIT NONE
    TYPE(t_field),POINTER :: f_ps(:)
    TYPE(t_field),POINTER :: f_phis(:)
    TYPE(t_field),POINTER :: f_theta_rhodz(:)
    TYPE(t_field),POINTER :: f_u(:)
    TYPE(t_field),POINTER :: f_q(:)
    
    REAL(rstd),POINTER :: ps(:)
    REAL(rstd),POINTER :: phis(:)
    REAL(rstd),POINTER :: u(:,:)
    REAL(rstd),POINTER :: theta_rhodz(:,:)
    REAL(rstd),POINTER :: q(:,:,:)

    INTEGER :: ind

    DO ind=1,ndomain
      CALL swap_dimensions(ind)
      CALL swap_geometry(ind)

      ps=f_ps(ind)
      phis=f_phis(ind)
      u=f_u(ind)
      theta_rhodz=f_theta_rhodz(ind)
      q=f_q(ind)
      CALL compute_etat0_DCMIP3(ps,phis,u,theta_rhodz,q)
    ENDDO
            
  END SUBROUTINE etat0
  

  SUBROUTINE compute_etat0_DCMIP3(ps, phis, u, theta_rhodz,q)
  USE icosa
  USE pression_mod
  USE theta2theta_rhodz_mod
  USE wind_mod
  IMPLICIT NONE
  REAL(rstd), INTENT(OUT) :: ps(iim*jjm)
  REAL(rstd), INTENT(OUT) :: phis(iim*jjm)
  REAL(rstd), INTENT(OUT) :: u(3*iim*jjm,llm)
  REAL(rstd), INTENT(OUT) :: theta_rhodz(iim*jjm,llm)
  REAL(rstd), INTENT(OUT) :: q(iim*jjm,llm,nqtot)
  
  REAL(rstd) :: Ts(iim*jjm)
  REAL(rstd) :: s(iim*jjm)
  REAL(rstd) :: T(iim*jjm,llm)
  REAL(rstd) :: p(iim*jjm,llm+1)
  REAL(rstd) :: theta(iim*jjm,llm)
  REAL(rstd) :: ulon(3*iim*jjm,llm)
  REAL(rstd) :: ulat(3*iim*jjm,llm)
  
  
  INTEGER :: i,j,l,ij
  REAL(rstd) :: Rd        ! gas constant of dry air, P=rho.Rd.T
  REAL(rstd) :: alpha, rm, zs
  REAL(rstd) :: lon,lat, C0, C1, GG
  REAL(rstd) :: p0psk, pspsk,r,zz, thetab, thetap
  REAL(rstd) :: dummy, pp
  LOGICAL    :: use_dcmip_routine
  
  Rd=cpp*kappa
  
  GG=(g/N)**2/cpp
  C0=0.25*u0*(u0+2.*Omega*radius)
  
  q(:,:,:)=0
  
!  use_dcmip_routine=.TRUE.
  use_dcmip_routine=.FALSE.
  dummy=0.

  pp=peq
  DO j=jj_begin,jj_end
     DO i=ii_begin,ii_end
        ij=(j-1)*iim+i
        CALL xyz2lonlat(xyz_i(ij,:),lon,lat)

        IF(use_dcmip_routine) THEN
           CALL test3_gravity_wave(lon,lat,pp,dummy,0, dummy,dummy,dummy,dummy,phis(ij),ps(ij),dummy,dummy)
        ELSE
           C1=C0*(cos(2*lat)-1)
           
           !--- GROUND GEOPOTENTIAL
           phis(ij)=0.
           
           !--- GROUND TEMPERATURE
           Ts(ij) = GG+(Teq-GG)*EXP(-C1*(N/g)**2)
           
           !--- GROUND PRESSURE
           Ps(ij) = peq*EXP(C1/GG/Rd)*(Ts(ij)/Teq)**(1/kappa)
           
           
           r=radius*acos(sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc))
           s(ij)= d**2/(d**2+r**2)
        END IF
     END DO
  END DO
  
  CALL compute_pression(ps,p,0)
  
  DO l=1,llm
     DO j=jj_begin,jj_end
        DO i=ii_begin,ii_end
           ij=(j-1)*iim+i
           pp=0.5*(p(ij,l+1)+p(ij,l))  ! full-layer pressure
           IF(use_dcmip_routine) THEN
              CALL xyz2lonlat(xyz_i(ij,:),lon,lat)
              CALL test3_gravity_wave(lon,lat,pp,dummy,0,dummy,dummy,dummy,T(ij,l),dummy,dummy,dummy,dummy)
           ELSE
              pspsk=(pp/ps(ij))**kappa
              p0psk=(Peq/ps(ij))**kappa
              thetab = Ts(ij)*p0psk / ( Ts(ij) / GG * ( pspsk-1) +1)  ! background pot. temp.
              zz     = -g/N**2*log( Ts(ij)/GG * (pspsk -1)+1)         ! altitude
              thetap = dtheta *sin(2*Pi*zz/Lz) * s(ij)                ! perturbation pot. temp.
              theta(ij,l) = thetab + thetap
              T(ij,l) = theta(ij,l)* ((pp/peq)**kappa)
              ! T(ij,l) = Ts(ij)*pspsk / ( Ts(ij) / GG * ( pspsk-1) +1)  ! background temp.
           END IF
        ENDDO
     ENDDO
  ENDDO
  
  IF(use_dcmip_routine) THEN
     CALL compute_temperature2theta_rhodz(ps,T,theta_rhodz,0)
  ELSE
!     CALL compute_theta2theta_rhodz(ps,theta,theta_rhodz,0)
     CALL compute_temperature2theta_rhodz(ps,T,theta_rhodz,0)
  END IF
  
  pp=peq
  DO l=1,llm
     DO j=jj_begin-1,jj_end+1
        DO i=ii_begin-1,ii_end+1
           ij=(j-1)*iim+i
           IF(use_dcmip_routine) THEN
              CALL xyz2lonlat(xyz_e(ij+u_right,:),lon,lat)
              CALL test3_gravity_wave(lon,lat,pp,0.,0, &
                   ulon(ij+u_right,l),ulat(ij+u_right,l),&
                   dummy,dummy,dummy,dummy,dummy,dummy)
              CALL xyz2lonlat(xyz_e(ij+u_lup,:),lon,lat)
              CALL test3_gravity_wave(lon,lat,pp,0.,0,&
                   ulon(ij+u_lup,l),ulat(ij+u_lup,l),&
                   dummy,dummy,dummy,dummy,dummy,dummy)
              CALL xyz2lonlat(xyz_e(ij+u_ldown,:),lon,lat)
              CALL test3_gravity_wave(lon,lat,pp,0.,0,&
                   ulon(ij+u_ldown,l),ulat(ij+u_ldown,l),&
                   dummy,dummy,dummy,dummy,dummy,dummy)
           ELSE
              CALL xyz2lonlat(xyz_e(ij+u_right,:),lon,lat)
              ulon(ij+u_right,l)=u0*cos(lat)
              ulat(ij+u_right,l)=0
              
              CALL xyz2lonlat(xyz_e(ij+u_lup,:),lon,lat)
              ulon(ij+u_lup,l)=u0*cos(lat)
              ulat(ij+u_lup,l)=0
              
              CALL xyz2lonlat(xyz_e(ij+u_ldown,:),lon,lat)
              ulon(ij+u_ldown,l)=u0*cos(lat)
              ulat(ij+u_ldown,l)=0
           END IF
        ENDDO
     ENDDO
  ENDDO
  
  CALL compute_wind_perp_from_lonlat_compound(ulon,ulat,u)    
  
END SUBROUTINE compute_etat0_DCMIP3


END MODULE etat0_DCMIP3_mod