Changeset 48 for codes/icosagcm/trunk


Ignore:
Timestamp:
07/29/12 03:14:26 (12 years ago)
Author:
dubos
Message:

Introduced DCMIP-supplied routines for optional use in test cases 1,2,3.

Location:
codes/icosagcm/trunk/src
Files:
1 added
1 edited

Legend:

Unmodified
Added
Removed

  • codes/icosagcm/trunk/src/etat0_dcmip3.f90

    r37 r48  
    77 
    88  USE genmod 
     9  USE dcmip_initial_conditions_test_1_2_3 
     10 
    911  PRIVATE 
    1012! some of the following should be obtained via getin 
    … …  
    1517  REAL(rstd),PARAMETER :: u0=20.         ! Maximum amplitude of the zonal wind (m.s-1) 
    1618  REAL(rstd),PARAMETER :: zs=0.          ! Surface elevation (m) 
    17   REAL(rstd),PARAMETER :: N=0.01         ! Brunt-Vaissala frequency (s-1) 
     19  REAL(rstd),PARAMETER :: N=0.01         ! Brunt-Vaisala frequency (s-1) 
    1820  REAL(rstd),PARAMETER :: Teq=300.       ! Surface temperature at the equator (K) 
    1921  REAL(rstd),PARAMETER :: Peq=1e5        ! Reference surface pressure at the equator (hPa) 
    … …  
    7476  REAL(rstd) :: Ts(iim*jjm) 
    7577  REAL(rstd) :: s(iim*jjm) 
    76   REAL(rstd) :: z(iim*jjm,llm) 
    7778  REAL(rstd) :: T(iim*jjm,llm) 
    7879  REAL(rstd) :: p(iim*jjm,llm+1) 
    79   REAL(rstd) :: thetap(iim*jjm,llm) 
    80   REAL(rstd) :: thetap_rhodz(iim*jjm,llm) 
     80  REAL(rstd) :: theta(iim*jjm,llm) 
    8181  REAL(rstd) :: ulon(3*iim*jjm,llm) 
    8282  REAL(rstd) :: ulat(3*iim*jjm,llm) 
    … …  
    8686  REAL(rstd) :: Rd        ! gas constant of dry air, P=rho.Rd.T 
    8787  REAL(rstd) :: alpha, rm, zs 
    88   REAL(rstd) :: lon,lat 
    89   REAL(rstd) :: C0,C1 
    90   REAL(rstd) :: GG 
    91   REAL(rstd) :: pspsk,r 
    92    
    93   Rd=cpp/kappa 
    94      
    95 !    alpha=g/Rd/lapse_rate              ! g/(Rd*Gamma) 
    96  
    97     GG=(g/N)**2/cpp 
    98     C0=0.25*u0*(u0+2.*Omega*radius) 
    99  
    100     q(:,:,:)=0 
    101      
    102      
    103     DO j=jj_begin,jj_end 
    104       DO i=ii_begin,ii_end 
     88  REAL(rstd) :: lon,lat, C0, C1, GG 
     89  REAL(rstd) :: p0psk, pspsk,r,zz, thetab, thetap 
     90  REAL(rstd) :: dummy, pp 
     91  LOGICAL    :: use_dcmip_routine 
     92   
     93  Rd=cpp*kappa 
     94   
     95  GG=(g/N)**2/cpp 
     96  C0=0.25*u0*(u0+2.*Omega*radius) 
     97   
     98  q(:,:,:)=0 
     99   
     100!  use_dcmip_routine=.TRUE. 
     101  use_dcmip_routine=.FALSE. 
     102  dummy=0. 
     103 
     104  pp=peq 
     105  DO j=jj_begin,jj_end 
     106     DO i=ii_begin,ii_end 
    105107        ij=(j-1)*iim+i 
    106108        CALL xyz2lonlat(xyz_i(ij,:),lon,lat) 
    107         C1=C0*(cos(2*lat)-1) 
    108          
    109         !--- GROUND GEOPOTENTIAL 
    110         phis(ij)=0. 
    111  
    112         !--- GROUND TEMPERATURE 
    113         Ts(ij) = GG+(Teq-GG)*EXP(-C1*(N/g)**2) 
    114  
    115         !--- GROUND PRESSURE 
    116         Ps(ij) = peq*EXP(C1/GG/Rd)*(Ts(ij)/Teq)**(1/kappa) 
    117  
    118  
    119         r=radius*acos(sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc)) 
    120         s(ij)= d**2/(d**2+r**2) 
    121       END DO 
    122     END DO 
    123  
    124     CALL compute_pression(ps,p,0) 
    125      
    126     DO l=1,llm 
    127       DO j=jj_begin,jj_end 
     109 
     110        IF(use_dcmip_routine) THEN 
     111           CALL test3_gravity_wave(lon,lat,pp,dummy,0, dummy,dummy,dummy,dummy,phis(ij),ps(ij),dummy,dummy) 
     112        ELSE 
     113           C1=C0*(cos(2*lat)-1) 
     114            
     115           !--- GROUND GEOPOTENTIAL 
     116           phis(ij)=0. 
     117            
     118           !--- GROUND TEMPERATURE 
     119           Ts(ij) = GG+(Teq-GG)*EXP(-C1*(N/g)**2) 
     120            
     121           !--- GROUND PRESSURE 
     122           Ps(ij) = peq*EXP(C1/GG/Rd)*(Ts(ij)/Teq)**(1/kappa) 
     123            
     124            
     125           r=radius*acos(sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc)) 
     126           s(ij)= d**2/(d**2+r**2) 
     127        END IF 
     128     END DO 
     129  END DO 
     130   
     131  CALL compute_pression(ps,p,0) 
     132   
     133  DO l=1,llm 
     134     DO j=jj_begin,jj_end 
    128135        DO i=ii_begin,ii_end 
    129           ij=(j-1)*iim+i 
    130  
    131           pspsk=(0.5*(p(ij,l+1)+p(ij,l))/ps(ij) )**kappa 
    132           T(ij,l)=Ts(ij)*pspsk / ( Ts(ij) / GG * ( pspsk-1) +1) 
    133           z(ij,l)=-g/N**2*log( Ts(ij)/GG * (pspsk -1)+1)  
     136           ij=(j-1)*iim+i 
     137           pp=0.5*(p(ij,l+1)+p(ij,l))  ! full-layer pressure 
     138           IF(use_dcmip_routine) THEN 
     139              CALL xyz2lonlat(xyz_i(ij,:),lon,lat) 
     140              CALL test3_gravity_wave(lon,lat,pp,dummy,0,dummy,dummy,dummy,T(ij,l),dummy,dummy,dummy,dummy) 
     141           ELSE 
     142              pspsk=(pp/ps(ij))**kappa 
     143              p0psk=(Peq/ps(ij))**kappa 
     144              thetab = Ts(ij)*p0psk / ( Ts(ij) / GG * ( pspsk-1) +1)  ! background pot. temp. 
     145              zz     = -g/N**2*log( Ts(ij)/GG * (pspsk -1)+1)         ! altitude 
     146              thetap = dtheta *sin(2*Pi*zz/Lz) * s(ij)                ! perturbation pot. temp. 
     147              theta(ij,l) = thetab + thetap 
     148              T(ij,l) = theta(ij,l)* ((pp/peq)**kappa) 
     149              ! T(ij,l) = Ts(ij)*pspsk / ( Ts(ij) / GG * ( pspsk-1) +1)  ! background temp. 
     150           END IF 
    134151        ENDDO 
    135       ENDDO 
    136     ENDDO 
    137      
    138     CALL compute_temperature2theta_rhodz(ps,T,theta_rhodz,0) 
    139  
    140     DO l=1,llm 
    141       DO j=jj_begin,jj_end 
    142         DO i=ii_begin,ii_end 
    143           ij=(j-1)*iim+i 
    144           thetap(ij,l)=dtheta *sin(2*Pi*z(ij,l)/Lz) * s(ij) 
     152     ENDDO 
     153  ENDDO 
     154   
     155  IF(use_dcmip_routine) THEN 
     156     CALL compute_temperature2theta_rhodz(ps,T,theta_rhodz,0) 
     157  ELSE 
     158!     CALL compute_theta2theta_rhodz(ps,theta,theta_rhodz,0) 
     159     CALL compute_temperature2theta_rhodz(ps,T,theta_rhodz,0) 
     160  END IF 
     161   
     162  pp=peq 
     163  DO l=1,llm 
     164     DO j=jj_begin-1,jj_end+1 
     165        DO i=ii_begin-1,ii_end+1 
     166           ij=(j-1)*iim+i 
     167           IF(use_dcmip_routine) THEN 
     168              CALL xyz2lonlat(xyz_e(ij+u_right,:),lon,lat) 
     169              CALL test3_gravity_wave(lon,lat,pp,0.,0, & 
     170                   ulon(ij+u_right,l),ulat(ij+u_right,l),& 
     171                   dummy,dummy,dummy,dummy,dummy,dummy) 
     172              CALL xyz2lonlat(xyz_e(ij+u_lup,:),lon,lat) 
     173              CALL test3_gravity_wave(lon,lat,pp,0.,0,& 
     174                   ulon(ij+u_lup,l),ulat(ij+u_lup,l),& 
     175                   dummy,dummy,dummy,dummy,dummy,dummy) 
     176              CALL xyz2lonlat(xyz_e(ij+u_ldown,:),lon,lat) 
     177              CALL test3_gravity_wave(lon,lat,pp,0.,0,& 
     178                   ulon(ij+u_ldown,l),ulat(ij+u_ldown,l),& 
     179                   dummy,dummy,dummy,dummy,dummy,dummy) 
     180           ELSE 
     181              CALL xyz2lonlat(xyz_e(ij+u_right,:),lon,lat) 
     182              ulon(ij+u_right,l)=u0*cos(lat) 
     183              ulat(ij+u_right,l)=0 
     184               
     185              CALL xyz2lonlat(xyz_e(ij+u_lup,:),lon,lat) 
     186              ulon(ij+u_lup,l)=u0*cos(lat) 
     187              ulat(ij+u_lup,l)=0 
     188               
     189              CALL xyz2lonlat(xyz_e(ij+u_ldown,:),lon,lat) 
     190              ulon(ij+u_ldown,l)=u0*cos(lat) 
     191              ulat(ij+u_ldown,l)=0 
     192           END IF 
    145193        ENDDO 
    146       ENDDO 
    147     ENDDO   
    148      
    149     CALL compute_theta2theta_rhodz(ps,thetap,thetap_rhodz,0) 
    150      
    151      
    152     DO l=1,llm 
    153       DO j=jj_begin,jj_end 
    154         DO i=ii_begin,ii_end 
    155           ij=(j-1)*iim+i 
    156           theta_rhodz(ij,l)=theta_rhodz(ij,l)+thetap_rhodz(ij,l) 
    157         ENDDO 
    158       ENDDO 
    159     ENDDO   
    160  
    161     DO l=1,llm 
    162       DO j=jj_begin-1,jj_end+1 
    163         DO i=ii_begin-1,ii_end+1 
    164           ij=(j-1)*iim+i 
    165           CALL xyz2lonlat(xyz_e(ij+u_right,:),lon,lat) 
    166           ulon(ij+u_right,l)=u0*cos(lat) 
    167           ulat(ij+u_right,l)=0 
    168  
    169           CALL xyz2lonlat(xyz_e(ij+u_lup,:),lon,lat) 
    170           ulon(ij+u_lup,l)=u0*cos(lat) 
    171           ulat(ij+u_lup,l)=0 
    172  
    173           CALL xyz2lonlat(xyz_e(ij+u_ldown,:),lon,lat) 
    174           ulon(ij+u_ldown,l)=u0*cos(lat) 
    175           ulat(ij+u_ldown,l)=0 
    176         ENDDO 
    177       ENDDO 
    178     ENDDO 
    179      
    180   
    181     CALL compute_wind_perp_from_lonlat_compound(ulon,ulat,u) 
    182  
    183      
    184           
    185   END SUBROUTINE compute_etat0_DCMIP3 
    186  
    187  
    188 END MODULE etat0_DCMIP3_mod  
     194     ENDDO 
     195  ENDDO 
     196   
     197  CALL compute_wind_perp_from_lonlat_compound(ulon,ulat,u)     
     198   
     199END SUBROUTINE compute_etat0_DCMIP3 
     200 
     201 
     202END MODULE etat0_DCMIP3_mod 
Note: See TracChangeset for help on using the changeset viewer.