<p><b>laura@ucar.edu</b> 2010-11-18 15:59:27 -0700 (Thu, 18 Nov 2010)</p><p>corrected the calculation of the relative humidity,saturation mixing ratio,and water vapor mixing ratio so that the initial moist soundings are the same in the hydrostatic and non-hydrostatic cores<br>
</p><hr noshade><pre><font color="gray">Modified: branches/atmos_physics/src/core_hyd_atmos/module_test_cases.F
===================================================================
--- branches/atmos_physics/src/core_hyd_atmos/module_test_cases.F 2010-11-18 22:55:46 UTC (rev 625)
+++ branches/atmos_physics/src/core_hyd_atmos/module_test_cases.F 2010-11-18 22:59:27 UTC (rev 626)
@@ -35,7 +35,8 @@
if (config_test_case == 3) write(0,*) ' normal-mode perturbation included '
block_ptr => domain % blocklist
do while (associated(block_ptr))
- call hyd_test_case_1(block_ptr % mesh, block_ptr % state % time_levs(1) % state, config_test_case)
+ call hyd_test_case_1(block_ptr % mesh, block_ptr % state % time_levs(1) % state, &
+ block_ptr % diag_physics, config_test_case)
do i=2,nTimeLevs
call copy_state(block_ptr % state % time_levs(i) % state, block_ptr % state % time_levs(1) % state)
end do
@@ -52,7 +53,7 @@
!----------------------------------------------------------------------------------------------------------
- subroutine hyd_test_case_1(grid, state, test_case)
+ subroutine hyd_test_case_1(grid, state, diag_physics, test_case)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Setup baroclinic wave test case from Jablonowski and Williamson 2008 (QJRMS)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -61,6 +62,7 @@
type (mesh_type), intent(inout) :: grid
type (state_type), intent(inout) :: state
+ type (diag_physics_type), intent(inout) :: diag_physics
integer, intent(in) :: test_case
real (kind=RKIND), parameter :: u0 = 35.0
@@ -70,24 +72,31 @@
real (kind=RKIND), parameter :: u_perturbation = 1., pert_radius = 0.1, latitude_pert = 40., longitude_pert = 20.
real (kind=RKIND), parameter :: theta_c = pii/4.0
real (kind=RKIND), parameter :: lambda_c = 3.0*pii/2.0
- real (kind=RKIND), parameter :: rh_max = 0.4 ! Maximum relative humidity
real (kind=RKIND), parameter :: k_x = 9. ! Normal mode wave number
real (kind=RKIND), dimension(:), pointer :: rdnu, rdnw, fnm, fnp, dbn, dnu, dnw
real (kind=RKIND), dimension(:), pointer :: surface_pressure
real (kind=RKIND), dimension(:,:), pointer :: pressure, theta, alpha, geopotential, h
- real (kind=RKIND), dimension(:,:,:), pointer :: scalars
- integer :: iCell, iEdge, vtx1, vtx2, ivtx, k, nz, nz1, index_qv
+ integer :: iCell, iEdge, vtx1, vtx2, ivtx, k, nz, nz1
real (kind=RKIND) :: u, v, flux, fluxk, lat1, lat2, eta_v, r_pert, u_pert, lat_pert, lon_pert, r
real (kind=RKIND) :: ptop, p0, phi
real (kind=RKIND) :: lon_Edge
- real (kind=RKIND), dimension(grid % nVertLevels, grid % nCells) :: rel_hum, temperature
+ real (kind=RKIND), dimension(grid % nVertLevels, grid % nCells) :: temperature
real (kind=RKIND) :: ptmp, es, qvs
integer :: iter
+!.. initialization of moisture:
+ integer:: index_qv
+! real (kind=RKIND),parameter :: rh_max = 0.40 ! Maximum relative humidity
+ real (kind=RKIND),parameter :: rh_max = 0.70 ! Maximum relative humidity
+ real (kind=RKIND):: ttemp
+ real (kind=RKIND),dimension(:,:), pointer:: qsat,relhum
+ real (kind=RKIND),dimension(:,:,:),pointer:: scalars
+!.. end initialization of moisture.
+
! real (kind=RKIND), dimension(grid % nVertLevelsP1 ) :: hyai, hybi, znu, znw, znwc, znwv, hyam, hybm
! real (kind=RKIND), dimension(grid % nVertLevelsP1 ) :: znuc, znuv, bn, divh, dpn, teta, phi
real (kind=RKIND), dimension(grid % nVertLevels + 1 ) :: hyai, hybi, znu, znw, znwc, znwv, hyam, hybm
@@ -149,10 +158,16 @@
alpha => state % alpha % array
geopotential => state % geopotential % array
h => state % h % array
+
+!.. initialization of moisture:
scalars => state % scalars % array
+ qsat => diag_physics % qsat % array
+ relhum => diag_physics % relhum % array
+ scalars(:,:,:) = 0.0
+ qsat(:,:) = 0.0
+ relhum(:,:) = 0.0
+!.. end initialization of moisture.
- scalars(:,:,:) = 0.
-
p0 = 100000.
bn (1) = 1.
znw(1) = 1.
@@ -366,27 +381,24 @@
theta(k,1)*(pressure(k,1)/p0)**(rgas/cp)
end do
-!
-! initialization for moisture
-!
- if (config_mp_physics /= 0) then
+!.. calculation of relative humidity:
+! if(config_mp_physics /= 0) then
+! do iCell = 1, grid % nCells
+! do k = 1, nz1
+! ptmp = 0.5*(pressure(k,iCell) + pressure(k+1,iCell))
+! if(ptmp < 50000.) then
+! relhum(k,iCell) = 0.0
+! else
+! relhum(k,iCell) = (1.-((p0-ptmp)/50000.)**1.25)
+! end if
+! relhum(k,iCell) = min(rh_max,relhum(k,iCell))
+! enddo
+! enddo
+! else
+! relhum(:,:) = 0.
+! endif
+!.. end calculation of relative humidity.
- do iCell=1,grid % nCells
- do k=1,nz1
- ptmp = 0.5*(pressure(k,iCell) + pressure(k+1,iCell))
- if (ptmp < 50000.) then
- rel_hum(k,iCell) = 0.0
- else
- rel_hum(k,iCell) = (1.-((p0-ptmp)/50000.)**1.25)
- end if
- rel_hum(k,iCell) = min(rh_max,rel_hum(k,iCell))
- end do
- end do
- else
- rel_hum(:,:) = 0.
- end if
-
-
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! iteration
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -411,25 +423,38 @@
*2.*u0*cos(znuv(k))**1.5 &
+(1.6*cos(phi)**3 &
*(sin(phi)**2+2./3.)-pii/4.)*a*omega)
+
+ if(config_mp_physics /= 0) then
+ !.. ldf (11-15-2010): move calculation of the relative humidity inside iteration loop to
+ !.. calculate the saturation mixing ratio and initialize the water vapor mixing ratio.
+ !.. We assume the temperature computed above is the virtual temperature when moisture is
+ !.. present:
+ if(ptmp < 50000.) then
+ relhum(k,iCell) = 0.0
+ else
+ relhum(k,iCell) = (1.-((p0-ptmp)/50000.)**1.25)
+ endif
+ relhum(k,iCell) = min(rh_max,relhum(k,iCell))
+
+ !.. conversion from virtual temperature to temperature:
+ ttemp = temperature(k,iCell)/(1+0.608*scalars(index_qv,k,iCell))
+
+ !.. calculation of water vapor mixing ratio:
+ if (ttemp > 273.15) then
+ es = 1000.*0.6112*exp(17.67*(ttemp-273.15)/(ttemp-29.65))
+ else
+ es = 1000.*0.6112*exp(21.8745584*(ttemp-273.15)/(ttemp-7.66))
+ end if
+ qsat(k,iCell) = (287.04/461.6)*es/(ptmp-es)
+ if(relhum(k,iCell) .eq. 0.0) qsat(k,iCell) = 0.0
+ scalars(index_qv,k,iCell) = relhum(k,iCell)*qsat(k,iCell)
+ endif
+
+ !conversion of virtual temperature to potential temperature:
+ theta (k,iCell) = temperature(k,iCell)/(1+0.608*scalars(index_qv,k,iCell))*(ptmp/p0)**(-rgas/cp)
+ alpha (k,iCell) = (rgas/p0)*theta(k,iCell)*(1.+1.61*scalars(index_qv,k,iCell))*(ptmp/p0)**cvpm
+ enddo
- temperature(k,iCell) = temperature(k,iCell)/(1.+0.61*scalars(index_qv,k,iCell))
-
- theta (k,iCell) = temperature(k,iCell)* &
- (0.5*(pressure(k,iCell)+pressure(k+1,iCell))/p0)**(-rgas/cp)
- alpha (k,iCell) = (rgas/p0)*theta(k,iCell)*(1.+1.61*scalars(index_qv,k,iCell))* &
- (0.5*(pressure(k,iCell)+pressure(k+1,iCell))/p0)**cvpm
-
- if (temperature(k,iCell) > 273.15) then
- es = 1000.*0.6112*exp(17.67*(temperature(k,iCell)-273.15)/(temperature(k,iCell)-29.65))
- else
- es = 1000.*0.6112*exp(21.8745584*(temperature(k,iCell)-273.16)/(temperature(k,iCell)-7.66))
- end if
- qvs = (287.04/461.6)*es/(ptmp-es)
-! qvs = 380.*exp(17.27*(temperature(k,iCell)-273.)/(temperature(k,iCell)-36.))/ptmp
-
- scalars(index_qv,k,iCell) = rel_hum(k,iCell)*qvs
- end do
-
do k=nz1,1,-1
pressure(k,iCell) = pressure(k+1,iCell)-dnw(k)*h(k,iCell)*(1.+scalars(index_qv,k,iCell))
geopotential(k,iCell) = geopotential(k+1,iCell)+dnw(k)*h(k,iCell)*alpha(k,iCell)
@@ -438,13 +463,23 @@
end do
end do
- write(6,*) 'ptop = ',ptop,' zt = ',geopotential(nz,1)/gravity
+ 202 format(2i6,10(1x,e15.8))
+ write(0,*) '--- initialization of water vapor:'
+ do iCell = 1, grid % nCells
+ if(iCell == 1 .or. iCell == grid % nCells) then
+ write(0,*)
+ do k = nz1, 1, -1
+ write(0,202) iCell,k,temperature(k,iCell),relhum(k,iCell),qsat(k,iCell),scalars(index_qv,k,iCell)
+ enddo
+ endif
+ enddo
- write(6,*) ' full sounding with moisture'
- do k=1,nz1
- write(6,*) k, geopotential(k,1)/gravity, 0.01*pressure(k,1), theta(k,1), &
- theta(k,1)*(pressure(k,1)/p0)**(rgas/cp)
- end do
+! write(6,*) 'ptop = ',ptop,' zt = ',geopotential(nz,1)/gravity
+! write(6,*) ' full sounding with moisture'
+! do k=1,nz1
+! write(6,*) k, geopotential(k,1)/gravity, 0.01*pressure(k,1), theta(k,1), &
+! theta(k,1)*(pressure(k,1)/p0)**(rgas/cp)
+! end do
! When initializing a scalar, be sure not to put unreasonably large values
! into indices in the moist class
Modified: branches/atmos_physics/src/core_nhyd_atmos/module_test_cases.F
===================================================================
--- branches/atmos_physics/src/core_nhyd_atmos/module_test_cases.F 2010-11-18 22:55:46 UTC (rev 625)
+++ branches/atmos_physics/src/core_nhyd_atmos/module_test_cases.F 2010-11-18 22:59:27 UTC (rev 626)
@@ -44,7 +44,8 @@
block_ptr => domain % blocklist
do while (associated(block_ptr))
write(0,*) ' calling test case setup '
- call nhyd_test_case_jw(block_ptr % mesh, block_ptr % state % time_levs(1) % state, block_ptr % diag, config_test_case)
+ call nhyd_test_case_jw(block_ptr % mesh, block_ptr % state % time_levs(1) % state, block_ptr % diag, &
+ block_ptr % diag_physics ,config_test_case)
write(0,*) ' returned from test case setup '
do i=2,nTimeLevs
call copy_state(block_ptr % state % time_levs(i) % state, block_ptr % state % time_levs(1) % state)
@@ -96,7 +97,7 @@
!----------------------------------------------------------------------------------------------------------
- subroutine nhyd_test_case_jw(grid, state, diag, test_case)
+ subroutine nhyd_test_case_jw(grid, state, diag, diag_physics, test_case)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Setup baroclinic wave test case from Jablonowski and Williamson 2008 (QJRMS)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -106,6 +107,7 @@
type (mesh_type), intent(inout) :: grid
type (state_type), intent(inout) :: state
type (diag_type), intent(inout) :: diag
+ type (diag_physics_type), intent(inout) :: diag_physics
integer, intent(in) :: test_case
real (kind=RKIND), parameter :: u0 = 35.0
@@ -115,22 +117,26 @@
real (kind=RKIND), parameter :: u_perturbation = 1., pert_radius = 0.1, latitude_pert = 40., longitude_pert = 20.
real (kind=RKIND), parameter :: theta_c = pii/4.0
real (kind=RKIND), parameter :: lambda_c = 3.0*pii/2.0
-! real (kind=RKIND), parameter :: rh_max = 0.4 ! Maximum relative humidity
- real (kind=RKIND), parameter :: rh_max = 0.65 ! Maximum relative humidity
real (kind=RKIND), parameter :: k_x = 9. ! Normal mode wave number
real (kind=RKIND), dimension(:), pointer :: rdzw, dzu, rdzu, fzm, fzp
+ real (kind=RKIND), dimension(:), pointer :: surface_pressure
real (kind=RKIND), dimension(:,:), pointer :: zgrid, zx, zz, hx
real (kind=RKIND), dimension(:,:), pointer :: pressure, ppb, pb, rho, rb, rr, tb, rtb, p, pp, dss, t, rt
real (kind=RKIND), dimension(:,:,:), pointer :: zf, zf3, zb, zb3
- real (kind=RKIND), dimension(:,:), pointer :: rh
- real (kind=RKIND), dimension(:,:,:), pointer :: scalars
real (kind=RKIND), dimension(:,:,:), pointer :: deriv_two
+
+!.. initialization of moisture:
+ integer:: index_qv
+! real (kind=RKIND),parameter :: rh_max = 0.40 ! Maximum relative humidity
+ real (kind=RKIND),parameter :: rh_max = 0.70 ! Maximum relative humidity
+ real (kind=RKIND),dimension(:,:), pointer:: qsat,relhum
+ real (kind=RKIND),dimension(:,:,:),pointer:: scalars
+ real (kind=RKIND),dimension(grid%nVertLevels,grid%nCells):: qv
+!.. end initialization of moisture.
integer :: iCell, iCell1, iCell2 , iEdge, vtx1, vtx2, ivtx, i, k, nz, nz1, itr, itrp, cell1, cell2, nCellsSolve
- integer :: index_qv
-
!This is temporary variable here. It just need when calculate tangential velocity v.
integer :: eoe, j
integer, dimension(:), pointer :: nEdgesOnEdge
@@ -145,8 +151,8 @@
real (kind=RKIND) :: r_earth, etavs, ztemp, zd, zt, dz, gam, delt, str
- real (kind=RKIND), dimension(grid % nVertLevels, grid % nCells) :: rel_hum, temperature, qv
- real (kind=RKIND) :: ptmp, es, qvs, xnutr, znut, ptemp
+ real (kind=RKIND), dimension(grid % nVertLevels) :: temperature
+ real (kind=RKIND) :: ptmp, es, qvs, xnutr, znut, ptemp, ttemp
integer :: iter
real (kind=RKIND), dimension(grid % nVertLevels + 1 ) :: hyai, hybi, znu, znw, znwc, znwv, hyam, hybm
@@ -227,18 +233,26 @@
rr => diag % rho_p % array
t => state % theta % array
rt => diag % rtheta_p % array
- rh => diag % rh % array
+ surface_pressure => diag % surface_pressure % array
+
+!.. initialization of moisture:
scalars => state % scalars % array
+ qsat => diag_physics % qsat % array
+ relhum => diag_physics % relhum % array
+ scalars(:,:,:) = 0.0
+ qsat(:,:) = 0.0
+ relhum(:,:) = 0.0
+ qv(:,:) = 0.0
+!.. end initialization of moisture.
+ surface_pressure(:) = 0.0
+
call initialize_advection_rk(grid)
call initialize_deformation_weights(grid)
index_qv = state % index_qv
- rh(:,:) = 0.
- scalars(:,:,:) = 0.
-
xnutr = 0.
zd = 12000.
znut = eta_t
@@ -502,39 +516,20 @@
rr (k,i) = 0.
end do
- if(i == 1) then
- do k=1,nz1
- write(0,*) ' k, ppb, pb, rb, tb (k,1) ',k,ppb(k,1),pb(k,1),rb(k,1)*zz(k,1),tb(k,1)
- enddo
- end if
+! if(i == 1) then
+! do k=1,nz1
+! write(0,*) ' k, ppb, pb, rb, tb (k,1) ',k,ppb(k,1),pb(k,1),rb(k,1)*zz(k,1),tb(k,1)
+! enddo
+! end if
+ 200 format(4i6,8(1x,e15.8))
201 format(3i6,8(1x,e15.8))
202 format(2i6,10(1x,e15.8))
-#ifdef DO_PHYSICS
-!... initialization for moisture (following initialization in hydrostatic core):
-! note that the initial profile of relative humidity is different from that
-! in the hydrostatic core since the grid is different.
-! write(0,*)
-! write(0,*) 'RELATIVE HUMIDITY:'
- do k = 1, nz1
- ztemp = .5*(zgrid(k+1,i)+zgrid(k,i))
- if(ztemp .gt. 5000.) then
- rel_hum(k,i) = 0.0
- elseif(ztemp .lt. 0.0) then
- rel_hum(k,i) = 1.0
- else
- rel_hum(k,i) = (1.-0.75*(ztemp/zt)**1.25)
- endif
- rel_hum(k,i) = min(rh_max,rel_hum(k,i))
-! if(ztemp .lt. 0.0) write(0,202) i,k,ztemp,rel_hum(k,i)
- enddo
-#else
- rel_hum(:,i) = 0.0
-#endif
-!
+ 203 format(i6,10(1x,e15.8))
+
! iterations to converge temperature as a function of pressure
!
- do itr = 1,10
+ do itr = 1,30
do k=1,nz1
eta (k) = (ppb(k,i)+pp(k,i))/p0
@@ -547,8 +542,8 @@
end do
phi = grid % latCell % array (i)
do k=1,nz1
- tt(k) = 0.
- tt(k) = teta(k)+.75*eta(k)*pii*u0/rgas*sin(etav(k)) &
+! tt(k) = teta(k)+.75*eta(k)*pii*u0/rgas*sin(etav(k)) &
+ temperature(k) = teta(k)+.75*eta(k)*pii*u0/rgas*sin(etav(k)) &
*sqrt(cos(etav(k)))* &
((-2.*sin(phi)**6 &
*(cos(phi)**2+1./3.)+10./63.) &
@@ -556,29 +551,40 @@
+(1.6*cos(phi)**3 &
*(sin(phi)**2+2./3.)-pii/4.)*r_earth*omega_e)
-
- !write(0,*) ' k, tt(k) ',k,tt(k)
ztemp = .5*(zgrid(k,i)+zgrid(k+1,i))
ptemp = ppb(k,i) + pp(k,i)
-! qv(k,i) = env_qv( ztemp, tt(k), ptemp, 0 )
#ifdef DO_PHYSICS
-!... computes the water vapor mixing ratio (following initialization in hydrostatic core):
-! note that the initial profile of water vapor is different from that in the hydrostatic
-! core since temperature and relative humidity are different.
- if (tt(k) > 273.15) then
- es = 1000.*0.6112*exp(17.67*(tt(k)-273.15)/(tt(k)-29.65))
+ !.. ldf (11-15-2010): move calculation of the relative humidity inside iteration loop to
+ !.. calculate the saturation mixing ratio and initialize the water vapor mixing ratio.
+ !.. We assume the temperature computed above is the virtual temperature when moisture is
+ !.. present:
+ if(ptemp < 50000.) then
+ relhum(k,i) = 0.0
+ elseif(ptemp > p0) then
+ relhum(k,i) = 1.0
else
- es = 1000.*0.6112*exp(21.8745584*(tt(k)-273.16)/(tt(k)-7.66))
- end if
- qvs = (287.04/461.6)*es/(ptemp-es)
- rh(k,i) = rel_hum(k,i)
- qv(k,i) = amin1(0.014,rel_hum(k,i)*qvs)
- scalars(index_qv,k,i) = amin1(0.014,rel_hum(k,i)*qvs)
-#else
- qv(k,i) = 0.
- scalars(index_qv,k,i) = 0.
+ relhum(k,i) = (1.-((p0-ptemp)/50000.)**1.25)
+ endif
+ relhum(k,i) = min(rh_max,relhum(k,i))
+
+ !.. conversion from virtual temperature to temperature:
+ ttemp = temperature(k)/(1.+0.608*scalars(index_qv,k,i))
+
+ !.. calculation of water vapor mixing ratio:
+ if (ttemp > 273.15) then
+ es = 1000.*0.6112*exp(17.67*(ttemp-273.15)/(ttemp-29.65))
+ else
+ es = 1000.*0.6112*exp(21.8745584*(ttemp-273.15)/(ttemp-7.66))
+ end if
+ qsat(k,i) = (287.04/461.6)*es/(ptemp-es)
+ if(relhum(k,i) .eq. 0.0) qsat(k,i) = 0.0
+ scalars(index_qv,k,i) = relhum(k,i)*qsat(k,i)
+ qv(k,i) = scalars(index_qv,k,i)
#endif
+ !.. conversion from virtual temperature to "modified" temperature:
+ tt(k) = temperature(k)*(1+scalars(index_qv,k,i))
+
end do
! do k=2,nz1
! cqw(k,i) = 1./(1.+.5*(qv(k,i)+qv(k-1,i)))
@@ -612,19 +618,35 @@
do k=1,nz1
p (k,i) = ((ppb(k,i)+pp(k,i))/p0)**(rgas/cp)
t (k,i) = tt(k)/p(k,i)
-! t (k,i) = tt(k)*(1.+1.61*scalars(index_qv,k,i))/p(k,i)
rt (k,i) = t(k,i)*rr(k,i)+rb(k,i)*(t(k,i)-tb(k,i))
rho (k,i) = rb(k,i) + rr(k,i)
end do
- if(i == 1) then
- do k=1,nz1
- write(0,*) ' k, p, t, rt ',k,p(k,1),t(k,1),rt(k,1)
- enddo
- end if
+ !calculation of surface pressure:
+ surface_pressure(i) = 0.5*dzw(1)*gravity &
+ * (1.25*(rr(1,i) + rb(1,i)) * (1. + scalars(index_qv,1,i)) &
+ - 0.25*(rr(2,i) + rb(2,i)) * (1. + scalars(index_qv,2,i)))
+ surface_pressure(i) = surface_pressure(i) + pp(1,i) + ppb(1,i)
+! if(i == 1) then
+! do k=1,nz1
+! write(0,*) ' k, p, t, rt ',k,p(k,1),t(k,1),rt(k,1)
+! enddo
+! end if
+
end do ! end loop over cells
+ write(0,*)
+ write(0,*) '--- initialization of water vapor:'
+ do iCell = 1, grid % nCells
+ if(iCell == 1 .or. iCell == grid % nCells) then
+ do k = nz1, 1, -1
+ write(0,202) iCell,k,t(k,iCell),relhum(k,iCell),qsat(k,iCell),scalars(index_qv,k,iCell)
+ enddo
+ write(0,*)
+ endif
+ enddo
+
lat_pert = latitude_pert*pii/180.
lon_pert = longitude_pert*pii/180.
</font>
</pre>