![(please configure the [header_logo] section in trac.ini)](/dynamico/chrome/site/your_project_logo.png){kind=logo}
Changeset 345
- Timestamp:
- 07/31/15 19:17:42 (9 years ago)
- Location:
- codes/icosagcm/trunk/src
- Files:
-
- 2 edited
-
etat0.f90 (modified) (8 diffs)
-
etat0_dcmip3.f90 (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
codes/icosagcm/trunk/src/etat0.f90
r344 r345 16 16 USE mpipara, ONLY : is_mpi_root 17 17 USE disvert_mod 18 ! Newinterface18 ! Generic interface 19 19 USE etat0_dcmip1_mod, ONLY : getin_etat0_dcmip1=>getin_etat0 20 20 USE etat0_dcmip2_mod, ONLY : getin_etat0_dcmip2=>getin_etat0 … … 22 22 USE etat0_williamson_mod, ONLY : getin_etat0_williamson=>getin_etat0 23 23 USE etat0_temperature_mod, ONLY: getin_etat0_temperature=>getin_etat0 24 ! Old interface24 ! Ad hoc interfaces 25 25 USE etat0_academic_mod, ONLY : etat0_academic=>etat0 26 USE etat0_dcmip3_mod, ONLY : etat0_dcmip3=>etat027 26 USE etat0_dcmip4_mod, ONLY : etat0_dcmip4=>etat0 28 27 USE etat0_heldsz_mod, ONLY : etat0_heldsz=>etat0 … … 53 52 CALL getin("etat0",etat0_type) 54 53 55 !------------------- Newinterface ---------------------54 !------------------- Generic interface --------------------- 56 55 collocated=.TRUE. 57 56 SELECT CASE (TRIM(etat0_type)) … … 65 64 CASE ('dcmip2_mountain','dcmip2_schaer_noshear','dcmip2_schaer_shear') 66 65 CALL getin_etat0_dcmip2 66 CASE ('dcmip3') 67 67 CASE ('dcmip5') 68 68 CALL getin_etat0_dcmip5 … … 74 74 END SELECT 75 75 76 !------------------- Old interface--------------------76 !------------------- Ad hoc interfaces -------------------- 77 77 SELECT CASE (TRIM(etat0_type)) 78 78 CASE ('start_file') … … 86 86 CALL etat0_venus(f_ps, f_phis, f_theta_rhodz, f_u, f_q) 87 87 PRINT *, "Venus (Lebonnois et al., 2012) test case" 88 CASE ('dcmip3')89 CALL etat0_dcmip3(f_ps,f_phis,f_theta_rhodz,f_u, f_q)90 88 CASE ('dcmip4') 91 89 IF(nqtot<2) THEN … … 171 169 USE etat0_dcmip1_mod, ONLY : compute_dcmip1 => compute_etat0 172 170 USE etat0_dcmip2_mod, ONLY : compute_dcmip2 => compute_etat0 171 USE etat0_dcmip3_mod, ONLY : compute_dcmip3 => compute_etat0 173 172 USE etat0_dcmip5_mod, ONLY : compute_dcmip5 => compute_etat0 174 173 USE etat0_williamson_mod, ONLY : compute_w91_6 => compute_etat0 … … 211 210 CALL compute_dcmip2(iim*jjm,lon_i,lat_i, phis, ps, temp_i, ulon_i, ulat_i) 212 211 CALL compute_dcmip2(3*iim*jjm,lon_e,lat_e, phis_e, ps_e, temp_e, ulon_e, ulat_e) 212 CASE('dcmip3') 213 CALL compute_dcmip3(iim*jjm,lon_i,lat_i, phis, ps, temp_i, ulon_i, ulat_i, q) 214 CALL compute_dcmip3(3*iim*jjm,lon_e,lat_e, phis_e, ps_e, temp_e, ulon_e, ulat_e, q_e) 213 215 CASE('dcmip5') 214 216 CALL compute_dcmip5(iim*jjm,lon_i,lat_i, phis, ps, temp_i, ulon_i, ulat_i, q) -
codes/icosagcm/trunk/src/etat0_dcmip3.f90
r295 r345 1 1 MODULE etat0_dcmip3_mod 2 3 ! test cases DCMIP 2012, category 3 : Non-hydrostatic gravity waves 4 5 ! Questions 6 ! Replace ps0 by preff ?? 7 8 USE genmod 9 USE dcmip_initial_conditions_test_1_2_3 10 2 ! test cases DCMIP 2012, category 3 : Non-hydrostatic gravity waves 3 IMPLICIT NONE 11 4 PRIVATE 12 13 PUBLIC etat0 14 5 PUBLIC :: compute_etat0 6 15 7 CONTAINS 16 17 18 SUBROUTINE etat0(f_ps,f_phis,f_theta_rhodz,f_u, f_q)19 USE icosa20 USE theta2theta_rhodz_mod21 IMPLICIT NONE22 TYPE(t_field),POINTER :: f_ps(:)23 TYPE(t_field),POINTER :: f_phis(:)24 TYPE(t_field),POINTER :: f_theta_rhodz(:)25 TYPE(t_field),POINTER :: f_u(:)26 TYPE(t_field),POINTER :: f_q(:)27 TYPE(t_field),POINTER,SAVE :: f_temp(:)28 29 REAL(rstd),POINTER :: ps(:)30 REAL(rstd),POINTER :: phis(:)31 REAL(rstd),POINTER :: u(:,:)32 REAL(rstd),POINTER :: Temp(:,:)33 REAL(rstd),POINTER :: q(:,:,:)34 35 INTEGER :: ind36 37 CALL allocate_field(f_temp,field_t,type_real,llm,name='temp')38 39 DO ind=1,ndomain40 IF (.NOT. assigned_domain(ind)) CYCLE41 CALL swap_dimensions(ind)42 CALL swap_geometry(ind)43 44 ps=f_ps(ind)45 phis=f_phis(ind)46 u=f_u(ind)47 q=f_q(ind)48 temp=f_temp(ind)49 CALL compute_etat0_DCMIP3(ps,phis,u,Temp,q)50 ENDDO51 52 CALL temperature2theta_rhodz(f_ps,f_temp,f_theta_rhodz)53 CALL deallocate_field(f_temp)54 55 END SUBROUTINE etat056 8 57 58 SUBROUTINE compute_etat0_DCMIP3(ps, phis, u, temp,q) 59 USE icosa 60 USE pression_mod 61 USE theta2theta_rhodz_mod 62 USE wind_mod 63 IMPLICIT NONE 64 REAL(rstd),PARAMETER :: u0=20. ! Maximum amplitude of the zonal wind (m.s-1) 65 REAL(rstd),PARAMETER :: N=0.01 ! Brunt-Vaisala frequency (s-1) 66 REAL(rstd),PARAMETER :: Teq=300. ! Surface temperature at the equator (K) 67 REAL(rstd),PARAMETER :: Peq=1e5 ! Reference surface pressure at the equator (hPa) 68 REAL(rstd),PARAMETER :: d=5000. ! Witdth parameter for theta 69 REAL(rstd),PARAMETER :: lonc=2*pi/3 ! Longitudinal centerpoint of theta 70 REAL(rstd),PARAMETER :: latc=0 ! Longitudinal centerpoint of theta 71 REAL(rstd),PARAMETER :: dtheta=1. ! Maximum amplitude of theta (K) 72 REAL(rstd),PARAMETER :: Lz=20000. ! Vertical wave lenght of the theta perturbation 73 74 REAL(rstd), INTENT(OUT) :: ps(iim*jjm) 75 REAL(rstd), INTENT(OUT) :: phis(iim*jjm) 76 REAL(rstd), INTENT(OUT) :: u(3*iim*jjm,llm) 77 REAL(rstd), INTENT(OUT) :: Temp(iim*jjm,llm) 78 REAL(rstd), INTENT(OUT) :: q(iim*jjm,llm,nqtot) 79 80 REAL(rstd) :: Ts(iim*jjm) 81 REAL(rstd) :: s(iim*jjm) 82 REAL(rstd) :: p(iim*jjm,llm+1) 83 REAL(rstd) :: theta(iim*jjm,llm) 84 REAL(rstd) :: ulon(3*iim*jjm,llm) 85 REAL(rstd) :: ulat(3*iim*jjm,llm) 86 87 88 INTEGER :: i,j,l,ij 89 REAL(rstd) :: Rd ! gas constant of dry air, P=rho.Rd.T 90 REAL(rstd) :: alpha, rm 91 REAL(rstd) :: C0, C1, GG 92 REAL(rstd) :: p0psk, pspsk,r,zz, thetab, thetap 93 REAL(rstd) :: dummy, pp 94 LOGICAL :: use_dcmip_routine 95 96 Rd=cpp*kappa 97 98 GG=(g/N)**2/cpp 99 C0=0.25*u0*(u0+2.*Omega*radius) 100 101 q(:,:,:)=0 102 103 ! use_dcmip_routine=.TRUE. 104 use_dcmip_routine=.FALSE. 105 dummy=0. 106 107 pp=peq 108 DO j=jj_begin,jj_end 109 DO i=ii_begin,ii_end 110 ij=(j-1)*iim+i 111 112 IF(use_dcmip_routine) THEN 113 CALL test3_gravity_wave(lon_i(ij),lat_i(ij),pp,dummy,0, dummy,dummy,dummy,dummy,phis(ij),ps(ij),dummy,dummy) 114 ELSE 115 C1=C0*(cos(2*lat_i(ij))-1) 116 117 !--- GROUND GEOPOTENTIAL 118 phis(ij)=0. 119 120 !--- GROUND TEMPERATURE 121 Ts(ij) = GG+(Teq-GG)*EXP(-C1*(N/g)**2) 122 123 !--- GROUND PRESSURE 124 Ps(ij) = peq*EXP(C1/GG/Rd)*(Ts(ij)/Teq)**(1/kappa) 125 126 127 r=radius*acos(sin(latc)*sin(lat_i(ij))+cos(latc)*cos(lat_i(ij))*cos(lon_i(ij)-lonc)) 128 s(ij)= d**2/(d**2+r**2) 129 END IF 130 END DO 131 END DO 132 133 !$OMP BARRIER 134 CALL compute_pression(ps,p,0) 135 !$OMP BARRIER 136 137 DO l=1,llm 138 DO j=jj_begin,jj_end 139 DO i=ii_begin,ii_end 140 ij=(j-1)*iim+i 141 pp=0.5*(p(ij,l+1)+p(ij,l)) ! full-layer pressure 142 IF(use_dcmip_routine) THEN 143 CALL test3_gravity_wave(lon_i(ij),lat_i(ij),pp,dummy,0, & 144 dummy,dummy,dummy,Temp(ij,l),dummy,dummy,dummy,dummy) 145 ELSE 146 pspsk=(pp/ps(ij))**kappa 147 p0psk=(Peq/ps(ij))**kappa 148 thetab = Ts(ij)*p0psk / ( Ts(ij) / GG * ( pspsk-1) +1) ! background pot. temp. 149 zz = -g/N**2*log( Ts(ij)/GG * (pspsk -1)+1) ! altitude 150 thetap = dtheta *sin(2*Pi*zz/Lz) * s(ij) ! perturbation pot. temp. 151 theta(ij,l) = thetab + thetap 152 Temp(ij,l) = theta(ij,l)* ((pp/peq)**kappa) 153 ! T(ij,l) = Ts(ij)*pspsk / ( Ts(ij) / GG * ( pspsk-1) +1) ! background temp. 154 END IF 155 ENDDO 156 ENDDO 157 ENDDO 158 159 ! IF(use_dcmip_routine) THEN 160 ! CALL compute_temperature2theta_rhodz(ps,T,theta_rhodz,0) 161 ! ELSE 162 ! CALL compute_temperature2theta_rhodz(ps,T,theta_rhodz,0) 163 ! END IF 164 165 pp=peq 166 DO l=1,llm 167 DO j=jj_begin-1,jj_end+1 168 DO i=ii_begin-1,ii_end+1 169 ij=(j-1)*iim+i 170 IF(use_dcmip_routine) THEN 171 CALL test3_gravity_wave(lon_e(ij+u_right),lat_e(ij+u_right), & 172 pp,0.,0, ulon(ij+u_right,l),ulat(ij+u_right,l),& 173 dummy,dummy,dummy,dummy,dummy,dummy) 174 CALL test3_gravity_wave(lon_e(ij+u_lup),lat_e(ij+u_lup), & 175 pp,0.,0, ulon(ij+u_lup,l),ulat(ij+u_lup,l),& 176 dummy,dummy,dummy,dummy,dummy,dummy) 177 CALL test3_gravity_wave(lon_e(ij+u_ldown),lat_e(ij+u_ldown), & 178 pp,0.,0, ulon(ij+u_ldown,l),ulat(ij+u_ldown,l),& 179 dummy,dummy,dummy,dummy,dummy,dummy) 180 ELSE 181 ulon(ij+u_right,l) = u0*cos(lat_e(ij+u_right)) 182 ulat(ij+u_right,l) = 0 183 ulon(ij+u_lup,l) = u0*cos(lat_e(ij+u_lup)) 184 ulat(ij+u_lup,l) = 0 185 ulon(ij+u_ldown,l) = u0*cos(lat_e(ij+u_ldown)) 186 ulat(ij+u_ldown,l) = 0 187 END IF 188 ENDDO 189 ENDDO 190 ENDDO 191 192 CALL compute_wind_perp_from_lonlat_compound(ulon,ulat,u) 193 194 END SUBROUTINE compute_etat0_DCMIP3 195 9 SUBROUTINE compute_etat0(ngrid,lon,lat, phis,ps,temp,ulon,ulat,q) 10 USE genmod 11 USE dcmip_initial_conditions_test_1_2_3 12 USE disvert_mod 13 INTEGER, INTENT(IN) :: ngrid 14 REAL(rstd), INTENT(IN) :: lon(ngrid) 15 REAL(rstd), INTENT(IN) :: lat(ngrid) 16 REAL(rstd), INTENT(OUT) :: phis(ngrid) 17 REAL(rstd), INTENT(OUT) :: ps(ngrid) 18 REAL(rstd), INTENT(OUT) :: ulon(ngrid,llm) 19 REAL(rstd), INTENT(OUT) :: ulat(ngrid,llm) 20 REAL(rstd), INTENT(OUT) :: temp(ngrid,llm) 21 REAL(rstd), INTENT(OUT) :: q(ngrid,llm,nqtot) 22 REAL(rstd),PARAMETER :: Peq=1e5 ! Reference surface pressure at the equator (hPa) 23 REAL(rstd) :: dummy, pp 24 INTEGER :: l,ij 25 pp=peq 26 DO ij=1,ngrid 27 CALL test3_gravity_wave(lon(ij),lat(ij),pp,dummy,0, & 28 dummy,dummy,dummy,dummy,phis(ij),ps(ij),dummy,dummy) 29 END DO 30 DO l=1,llm 31 DO ij=1,ngrid 32 pp = .5*(ap(l)+ap(l+1)) + .5*(bp(l)+bp(l+1))*ps(ij) ! full-layer pressure 33 CALL test3_gravity_wave(lon(ij),lat(ij),pp,dummy,0, & 34 ulon(ij,l),ulat(ij,l),dummy,Temp(ij,l),dummy,dummy,dummy,dummy) 35 END DO 36 END DO 37 q=0. 38 END SUBROUTINE compute_etat0 196 39 197 40 END MODULE etat0_DCMIP3_mod
Note: See TracChangeset
for help on using the changeset viewer.
