<p><b>laura@ucar.edu</b> 2013-02-28 15:35:04 -0700 (Thu, 28 Feb 2013)</p><p>merge the latest revision of the atmos_physics branch (revision 2516)<br>
</p><hr noshade><pre><font color="gray">Index: branches/atmphys_gfs_mmm
===================================================================
--- branches/atmphys_gfs_mmm 2013-02-28 22:21:22 UTC (rev 2525)
+++ branches/atmphys_gfs_mmm 2013-02-28 22:35:04 UTC (rev 2526)
Property changes on: branches/atmphys_gfs_mmm
___________________________________________________________________
Modified: svn:mergeinfo
## -1,3 +1,4 ##
+/branches/atmos_physics:2282-2525
/branches/cam_mpas_nh:1260-1270
/branches/ocean_projects/ale_split_exp:1437-1483
/branches/ocean_projects/ale_vert_coord:1225-1383
\ No newline at end of property
Modified: branches/atmphys_gfs_mmm/Makefile
===================================================================
--- branches/atmphys_gfs_mmm/Makefile 2013-02-28 22:21:22 UTC (rev 2525)
+++ branches/atmphys_gfs_mmm/Makefile 2013-02-28 22:35:04 UTC (rev 2526)
@@ -195,8 +195,15 @@
CPPINCLUDES = -I../inc -I$(NETCDF)/include -I$(PIO) -I$(PNETCDF)/include
FCINCLUDES = -I../inc -I$(NETCDF)/include -I$(PIO) -I$(PNETCDF)/include
-LIBS = -L$(PIO) -L$(PNETCDF)/lib -L$(NETCDF)/lib -lpio -lpnetcdf -lnetcdf
+LIBS = -L$(PIO) -L$(PNETCDF)/lib -L$(NETCDF)/lib -lpio -lpnetcdf
+NCLIB = -lnetcdf
+NCLIBF = -lnetcdff
+ifneq ($(wildcard $(NETCDF)/lib/libnetcdff.*), ) # CHECK FOR NETCDF4
+ LIBS += $(NCLIBF)
+endif # CHECK FOR NETCDF4
+LIBS += $(NCLIB)
+
RM = rm -f
CPP = cpp -C -P -traditional
RANLIB = ranlib
@@ -277,9 +284,6 @@
TAU_MESSAGE="TAU Hooks are off."
endif
-ifneq ($(wildcard $(NETCDF)/lib/libnetcdff.*), ) # CHECK FOR NETCDF4
- LIBS += -lnetcdff
-endif # CHECK FOR NETCDF4
####################################################
# Section for adding external libraries and includes
Modified: branches/atmphys_gfs_mmm/src/core_atmos_physics/Makefile
===================================================================
--- branches/atmphys_gfs_mmm/src/core_atmos_physics/Makefile 2013-02-28 22:21:22 UTC (rev 2525)
+++ branches/atmphys_gfs_mmm/src/core_atmos_physics/Makefile 2013-02-28 22:35:04 UTC (rev 2526)
@@ -19,16 +19,17 @@
OBJS = \
mpas_atmphys_driver_cloudiness.o \
mpas_atmphys_driver_convection_deep.o \
+ mpas_atmphys_driver_gwdo.o \
mpas_atmphys_driver_lsm.o \
mpas_atmphys_driver_microphysics.o \
mpas_atmphys_driver_radiation_lw.o \
mpas_atmphys_driver_radiation_sw.o \
mpas_atmphys_driver_sfclayer.o \
mpas_atmphys_driver_pbl.o \
+ mpas_atmphys_driver.o \
mpas_atmphys_camrad_init.o \
mpas_atmphys_control.o \
mpas_atmphys_date_time.o \
- mpas_atmphys_driver.o \
mpas_atmphys_init.o \
mpas_atmphys_landuse.o \
mpas_atmphys_lsm_noahinit.o \
@@ -75,6 +76,10 @@
./physics_wrf/module_cu_kfeta.o \
./physics_wrf/module_cu_tiedtke.o
+mpas_atmphys_driver_gwdo.o: \
+ mpas_atmphys_vars.o \
+ ./physics_wrf/module_bl_gwdo.o
+
mpas_atmphys_driver_lsm.o: \
mpas_atmphys_constants.o \
mpas_atmphys_landuse.o \
@@ -205,6 +210,7 @@
mpas_atmphys_driver.o: \
mpas_atmphys_driver_convection_deep.o \
+ mpas_atmphys_driver_gwdo.o \
mpas_atmphys_driver_pbl.o \
mpas_atmphys_driver_radiation_lw.o \
mpas_atmphys_driver_radiation_sw.o \
Modified: branches/atmphys_gfs_mmm/src/core_atmos_physics/mpas_atmphys_control.F
===================================================================
--- branches/atmphys_gfs_mmm/src/core_atmos_physics/mpas_atmphys_control.F 2013-02-28 22:21:22 UTC (rev 2525)
+++ branches/atmphys_gfs_mmm/src/core_atmos_physics/mpas_atmphys_control.F 2013-02-28 22:35:04 UTC (rev 2526)
@@ -41,6 +41,7 @@
write(0,*) ' config_eddy_scheme = ', trim(config_eddy_scheme)
write(0,*) ' config_lsm_scheme = ', trim(config_lsm_scheme)
write(0,*) ' config_pbl_scheme = ', trim(config_pbl_scheme)
+ write(0,*) ' config_gwdo_scheme = ', trim(config_gwdo_scheme)
write(0,*) ' config_radt_cld_scheme = ', trim(config_radt_cld_scheme)
write(0,*) ' config_radt_lw_scheme = ', trim(config_radt_lw_scheme)
write(0,*) ' config_radt_sw_scheme = ', trim(config_radt_sw_scheme)
@@ -93,6 +94,29 @@
endif
+!gravity wave drag over orography scheme:
+ if(.not. (config_gwdo_scheme .eq. 'off' .or. &
+ config_gwdo_scheme .eq. 'ysu_gwdo')) then
+
+ write(mpas_err_message,'(A,A10)') 'illegal value for gwdo_scheme: ', &
+ trim(config_gwdo_scheme)
+ call physics_error_fatal(mpas_err_message)
+
+ elseif(config_gwdo_scheme .eq. 'ysu_gwdo' .and. config_pbl_scheme .ne. 'ysu') then
+
+ write(mpas_err_message,'(A,A10)') 'turn YSU PBL scheme on with config_gwdo = ysu_gwdo:', &
+ trim(config_gwdo_scheme)
+ call physics_error_fatal(mpas_err_message)
+
+ endif
+ if(config_gwdo_scheme .eq. 'ysu_gwdo') then
+
+ write(mpas_err_message,'(A,A10)') 'gwdo_scheme being tested - do not use right now: ', &
+ trim(config_gwdo_scheme)
+ call physics_error_fatal(mpas_err_message)
+
+ endif
+
!diffusion scheme:
if(.not. (config_eddy_scheme .eq. 'off')) then
Modified: branches/atmphys_gfs_mmm/src/core_atmos_physics/mpas_atmphys_driver.F
===================================================================
--- branches/atmphys_gfs_mmm/src/core_atmos_physics/mpas_atmphys_driver.F 2013-02-28 22:21:22 UTC (rev 2525)
+++ branches/atmphys_gfs_mmm/src/core_atmos_physics/mpas_atmphys_driver.F 2013-02-28 22:35:04 UTC (rev 2526)
@@ -5,6 +5,7 @@
use mpas_atmphys_driver_cloudiness
use mpas_atmphys_driver_convection_deep
+ use mpas_atmphys_driver_gwdo
use mpas_atmphys_driver_pbl
use mpas_atmphys_driver_lsm
use mpas_atmphys_driver_radiation_sw
@@ -118,6 +119,14 @@
call deallocate_pbl
endif
+ !call to gravity wave drag over orography scheme:
+ if(config_gwdo_scheme .ne. 'off') then
+ call allocate_gwdo
+ call driver_gwdo(itimestep,block%mesh,block%sfc_input,block%diag_physics, &
+ block%tend_physics)
+ call deallocate_gwdo
+ endif
+
!call to convection scheme:
call update_convection_step1(block%mesh,block%diag_physics,block%tend_physics)
if(l_conv) then
Copied: branches/atmphys_gfs_mmm/src/core_atmos_physics/mpas_atmphys_driver_gwdo.F (from rev 2525, branches/atmos_physics/src/core_atmos_physics/mpas_atmphys_driver_gwdo.F)
===================================================================
--- branches/atmphys_gfs_mmm/src/core_atmos_physics/mpas_atmphys_driver_gwdo.F (rev 0)
+++ branches/atmphys_gfs_mmm/src/core_atmos_physics/mpas_atmphys_driver_gwdo.F 2013-02-28 22:35:04 UTC (rev 2526)
@@ -0,0 +1,213 @@
+!=============================================================================================
+ module mpas_atmphys_driver_gwdo
+ use mpas_configure, only: len_disp => config_len_disp
+ use mpas_grid_types
+
+ use mpas_atmphys_constants
+ use mpas_atmphys_vars
+
+!from wrf physics:
+ use module_bl_gwdo
+
+ implicit none
+ private
+ public:: allocate_gwdo, &
+ deallocate_gwdo, &
+ driver_gwdo
+
+ integer,private:: i,j,k
+
+ contains
+
+!=============================================================================================
+ subroutine allocate_gwdo
+!=============================================================================================
+
+ if(.not.allocated(area_p) ) allocate(area_p(ims:ime,jms:jme) )
+ if(.not.allocated(var2d_p) ) allocate(var2d_p(ims:ime,jms:jme) )
+ if(.not.allocated(con_p) ) allocate(con_p(ims:ime,jms:jme) )
+ if(.not.allocated(oa1_p) ) allocate(oa1_p(ims:ime,jms:jme) )
+ if(.not.allocated(oa2_p) ) allocate(oa2_p(ims:ime,jms:jme) )
+ if(.not.allocated(oa3_p) ) allocate(oa3_p(ims:ime,jms:jme) )
+ if(.not.allocated(oa4_p) ) allocate(oa4_p(ims:ime,jms:jme) )
+ if(.not.allocated(ol1_p) ) allocate(ol1_p(ims:ime,jms:jme) )
+ if(.not.allocated(ol2_p) ) allocate(ol2_p(ims:ime,jms:jme) )
+ if(.not.allocated(ol3_p) ) allocate(ol3_p(ims:ime,jms:jme) )
+ if(.not.allocated(ol4_p) ) allocate(ol4_p(ims:ime,jms:jme) )
+ if(.not.allocated(kpbl_p )) allocate(kpbl_p(ims:ime,jms:jme) )
+ if(.not.allocated(dusfcg_p)) allocate(dusfcg_p(ims:ime,jms:jme))
+ if(.not.allocated(dvsfcg_p)) allocate(dvsfcg_p(ims:ime,jms:jme))
+
+ if(.not.allocated(dtaux3d_p)) allocate(dtaux3d_p(ims:ime,kms:kme,jms:jme))
+ if(.not.allocated(dtauy3d_p)) allocate(dtauy3d_p(ims:ime,kms:kme,jms:jme))
+ if(.not.allocated(rublten_p)) allocate(rublten_p(ims:ime,kms:kme,jms:jme))
+ if(.not.allocated(rvblten_p)) allocate(rvblten_p(ims:ime,kms:kme,jms:jme))
+
+ end subroutine allocate_gwdo
+
+!=============================================================================================
+ subroutine deallocate_gwdo
+!=============================================================================================
+
+ if(allocated(area_p) ) deallocate(area_p )
+ if(allocated(var2d_p) ) deallocate(var2d_p )
+ if(allocated(con_p) ) deallocate(con_p )
+ if(allocated(oa1_p) ) deallocate(oa1_p )
+ if(allocated(oa2_p) ) deallocate(oa2_p )
+ if(allocated(oa3_p) ) deallocate(oa3_p )
+ if(allocated(oa4_p) ) deallocate(oa4_p )
+ if(allocated(ol1_p) ) deallocate(ol1_p )
+ if(allocated(ol2_p) ) deallocate(ol2_p )
+ if(allocated(ol3_p) ) deallocate(ol3_p )
+ if(allocated(ol4_p) ) deallocate(ol4_p )
+ if(allocated(kpbl_p )) deallocate(kpbl_p )
+ if(allocated(dusfcg_p)) deallocate(dusfcg_p)
+ if(allocated(dvsfcg_p)) deallocate(dvsfcg_p)
+
+ if(allocated(dtaux3d_p)) deallocate(dtaux3d_p)
+ if(allocated(dtauy3d_p)) deallocate(dtauy3d_p)
+ if(allocated(rublten_p)) deallocate(rublten_p)
+ if(allocated(rvblten_p)) deallocate(rvblten_p)
+
+ end subroutine deallocate_gwdo
+
+!=============================================================================================
+ subroutine gwdo_from_MPAS(mesh,sfc_input,diag_physics,tend_physics)
+!=============================================================================================
+
+!input arguments:
+ type(mesh_type),intent(in):: mesh
+ type(sfc_input_type),intent(in) :: sfc_input
+ type(diag_physics_type),intent(in):: diag_physics
+ type(tend_physics_type),intent(in):: tend_physics
+
+!---------------------------------------------------------------------------------------------
+
+ do j = jts,jte
+ do i = its,ite
+ area_p(i,j) = mesh % areaCell % array(i)
+ var2d_p(i,j) = sfc_input % var2d % array(i)
+ con_p(i,j) = sfc_input % con % array(i)
+ oa1_p(i,j) = sfc_input % oa1 % array(i)
+ oa2_p(i,j) = sfc_input % oa2 % array(i)
+ oa3_p(i,j) = sfc_input % oa3 % array(i)
+ oa4_p(i,j) = sfc_input % oa4 % array(i)
+ ol1_p(i,j) = sfc_input % ol1 % array(i)
+ ol2_p(i,j) = sfc_input % ol2 % array(i)
+ ol3_p(i,j) = sfc_input % ol3 % array(i)
+ ol4_p(i,j) = sfc_input % ol4 % array(i)
+ enddo
+ enddo
+
+ do j = jts,jte
+ do i = its,ite
+ kpbl_p(i,j) = diag_physics % kpbl % array(i)
+ dusfcg_p(i,j) = diag_physics % dusfcg % array(i)
+ dvsfcg_p(i,j) = diag_physics % dvsfcg % array(i)
+ enddo
+ enddo
+
+ do j = jts,jte
+ do k = kts,kte
+ do i = its,ite
+ dtaux3d_p(i,k,j) = diag_physics % dtaux3d % array(k,i)
+ dtauy3d_p(i,k,j) = diag_physics % dtauy3d % array(k,i)
+ rublten_p(i,k,j) = tend_physics % rublten % array(k,i)
+ rvblten_p(i,k,j) = tend_physics % rvblten % array(k,i)
+ enddo
+ enddo
+ enddo
+
+ end subroutine gwdo_from_MPAS
+
+!=============================================================================================
+ subroutine gwdo_to_MPAS(diag_physics,tend_physics)
+!=============================================================================================
+
+!inout arguments:
+ type(diag_physics_type),intent(inout):: diag_physics
+ type(tend_physics_type),intent(inout):: tend_physics
+
+!---------------------------------------------------------------------------------------------
+
+ do j = jts,jte
+ do i = its,ite
+ diag_physics % dusfcg % array(i) = dusfcg_p(i,j)
+ diag_physics % dvsfcg % array(i) = dvsfcg_p(i,j)
+ enddo
+ enddo
+
+ do j = jts,jte
+ do k = kts,kte
+ do i = its,ite
+ diag_physics % dtaux3d % array(k,i) = dtaux3d_p(i,k,j)
+ diag_physics % dtauy3d % array(k,i) = dtauy3d_p(i,k,j)
+ tend_physics % rublten % array(k,i) = rublten_p(i,k,j)
+ tend_physics % rvblten % array(k,i) = rvblten_p(i,k,j)
+ enddo
+ enddo
+ enddo
+
+ end subroutine gwdo_to_MPAS
+
+!=============================================================================================
+ subroutine driver_gwdo(itimestep,mesh,sfc_input,diag_physics,tend_physics)
+!=============================================================================================
+
+!input arguments:
+ type(mesh_type),intent(in):: mesh
+ type(sfc_input_type),intent(in):: sfc_input
+ integer,intent(in):: itimestep
+
+!inout arguments:
+ type(diag_physics_type),intent(inout):: diag_physics
+ type(tend_physics_type),intent(inout):: tend_physics
+
+!---------------------------------------------------------------------------------------------
+ write(0,*)
+ write(0,*) '--- enter subroutine driver_gwdo: dt_pbl=',dt_pbl
+
+!copy all MPAS arrays to rectanguler grid arrays:
+ call gwdo_from_MPAS(mesh,sfc_input,diag_physics,tend_physics)
+
+ gwdo_select: select case (trim(gwdo_scheme))
+
+ case("ysu_gwdo")
+#if defined(do_hydrostatic_pressure)
+!... REARRANGED CALL USING HYDROSTATIC PRESSURE:
+ call gwdo ( &
+ p3d = pres_hydd_p , p3di = pres2_hydd_p , pi3d = pi_p , &
+ u3d = u_p , v3d = v_p , t3d = t_p , &
+ qv3d = qv_p , z = z_p , rublten = rublten_p , &
+ rvblten = rvblten_p , dtaux3d = dtaux3d_p , dtauy3d = dtauy3d_p , &
+ dusfcg = dusfcg_p , dvsfcg = dvsfcg_p , kpbl2d = kpbl_p , &
+ areaCell = area_p , itimestep = itimestep , dt = dt_pbl , &
+ dx = len_disp , cp = cp , g = g , &
+ rd = R_d , rv = R_v , ep1 = ep_1 , &
+ pi = pii , var2d = var2d_p , oc12d = con_p , &
+ oa2d1 = oa1_p , oa2d2 = oa2_p , oa2d3 = oa3_p , &
+ oa2d4 = oa4_p , ol2d1 = ol1_p , ol2d2 = ol2_p , &
+ ol2d3 = ol3_p , ol2d4 = ol4_p , &
+ ids = ids , ide = ide , jds = jds , jde = jde , kds = kds , kde = kde , &
+ ims = ims , ime = ime , jms = jms , jme = jme , kms = kms , kme = kme , &
+ its = its , ite = ite , jts = jts , jte = jte , kts = kts , kte = kte &
+ )
+#else
+!... REARRANGED CALL:
+ call gwdo ( &
+ )
+#endif
+
+ case default
+
+ end select gwdo_select
+
+!copy all arrays back to the MPAS grid:
+ call gwdo_to_MPAS(diag_physics,tend_physics)
+ write(0,*) '--- end subroutine driver_gwdo'
+
+ end subroutine driver_gwdo
+
+!=============================================================================================
+ end module mpas_atmphys_driver_gwdo
+!=============================================================================================
Modified: branches/atmphys_gfs_mmm/src/core_atmos_physics/mpas_atmphys_driver_microphysics.F
===================================================================
--- branches/atmphys_gfs_mmm/src/core_atmos_physics/mpas_atmphys_driver_microphysics.F 2013-02-28 22:21:22 UTC (rev 2525)
+++ branches/atmphys_gfs_mmm/src/core_atmos_physics/mpas_atmphys_driver_microphysics.F 2013-02-28 22:35:04 UTC (rev 2526)
@@ -386,6 +386,7 @@
!local variables:
integer:: i,j,k
+ real(kind=RKIND):: rho_a
!---------------------------------------------------------------------------------------------
@@ -397,8 +398,9 @@
!precipitable water:
diag_physics % precipw % array(i) = 0._RKIND
do k = kts,kte
- diag_physics % precipw % array(i) = diag_physics % precipw % array(i) &
- + qv_p(i,k,j) * dz_p(i,k,j)
+ rho_a = rho_p(i,k,j) / (1._RKIND + qv_p(i,k,j))
+ diag_physics % precipw % array(i) = &
+ diag_physics % precipw % array(i) + qv_p(i,k,j) * rho_a * dz_p(i,k,j)
enddo
!time-step precipitation:
@@ -459,8 +461,8 @@
real(kind=RKIND),dimension(:),allocatable:: qv1d,qr1d,qs1d,qg1d,t1d,p1d,dBZ1d
!---------------------------------------------------------------------------------------------
- write(0,*)
- write(0,*) '--- enter subroutine COMPUTE_RADAR_REFLECTIVITY:'
+!write(0,*)
+!write(0,*) '--- enter subroutine COMPUTE_RADAR_REFLECTIVITY:'
microp_select: select case(microp_scheme)
@@ -501,8 +503,8 @@
! write(0,201) i,k,dBZ1d(k)
enddo
diag_physics % refl10cm_max % array(i) = maxval(dBZ1d(:))
- if(diag_physics % refl10cm_max % array(i) .gt. 0.) &
- write(0,201) j,i,diag_physics % refl10cm_max % array(i)
+! if(diag_physics % refl10cm_max % array(i) .gt. 0.) &
+! write(0,201) j,i,diag_physics % refl10cm_max % array(i)
enddo
enddo
@@ -517,7 +519,7 @@
case default
end select microp_select
- write(0,*) '--- end subroutine COMPUTE_RADAR_REFLECTIVITY'
+!write(0,*) '--- end subroutine COMPUTE_RADAR_REFLECTIVITY'
201 format(2i6,e15.8)
Modified: branches/atmphys_gfs_mmm/src/core_atmos_physics/mpas_atmphys_initialize_real.F
===================================================================
--- branches/atmphys_gfs_mmm/src/core_atmos_physics/mpas_atmphys_initialize_real.F 2013-02-28 22:21:22 UTC (rev 2525)
+++ branches/atmphys_gfs_mmm/src/core_atmos_physics/mpas_atmphys_initialize_real.F 2013-02-28 22:35:04 UTC (rev 2526)
@@ -80,9 +80,6 @@
call physics_init_sst(mesh,fg)
endif
-!initialize seaice points:
- call physics_init_seaice(mesh,fg)
-
!initialization of the surface background albedo: interpolation of the monthly values to the
!initial date:
initial_date = trim(config_start_time)
@@ -126,6 +123,9 @@
!initialization of soil layers properties:
call init_soil_layers(mesh,fg,dminfo)
+!initialize seaice points:
+ call physics_init_seaice(mesh,fg)
+
!define xland over land and ocean:
do iCell = 1, nCellsSolve
if(landmask(iCell) .eq. 1 .or. (landmask(iCell).eq.0 .and. seaice(iCell).eq.1._RKIND)) then
@@ -460,7 +460,7 @@
integer:: iCell,nCells
integer,dimension(:),pointer:: landmask
- real(kind=RKIND),dimension(:),pointer :: sst,tsk
+ real(kind=RKIND),dimension(:),pointer :: sst,tsk,xice
real(kind=RKIND),dimension(:,:),pointer:: tslb
!--------------------------------------------------------------------------------------------------
@@ -474,14 +474,23 @@
sst => input % sst % array
tsk => input % skintemp % array
tslb => input % tslb % array
+ xice => input % xice % array
!update the skin temperature and the soil temperature of the first soil layer with the updated
!sea-surface temperatures:
+!change made so that the SSTs read for the surface update file are the same as the skin temperature
+!over the oceans.
+!do iCell = 1, nCells
+! if(landmask(iCell) == 0 .and. xice(iCell) == 0) then
+! tsk(iCell) = sst(iCell)
+! endif
+!enddo
do iCell = 1, nCells
if(landmask(iCell) == 0) then
tsk(iCell) = sst(iCell)
endif
enddo
+
write(0,*) '--- end subroutine physics_update_sst:'
end subroutine physics_init_sst
Modified: branches/atmphys_gfs_mmm/src/core_atmos_physics/mpas_atmphys_interface_nhyd.F
===================================================================
--- branches/atmphys_gfs_mmm/src/core_atmos_physics/mpas_atmphys_interface_nhyd.F 2013-02-28 22:21:22 UTC (rev 2525)
+++ branches/atmphys_gfs_mmm/src/core_atmos_physics/mpas_atmphys_interface_nhyd.F 2013-02-28 22:35:04 UTC (rev 2526)
@@ -55,10 +55,13 @@
if(.not.allocated(qg_p) ) allocate(qg_p(ims:ime,kms:kme,jms:jme) )
!... arrays used for calculating the hydrostatic pressure and exner function:
- if(.not.allocated(psfc_hyd_p) ) allocate(psfc_hyd_p(ims:ime,jms:jme) )
- if(.not.allocated(pres_hyd_p) ) allocate(pres_hyd_p(ims:ime,kms:kme,jms:jme) )
- if(.not.allocated(pres2_hyd_p)) allocate(pres2_hyd_p(ims:ime,kms:kme,jms:jme))
- if(.not.allocated(znu_hyd_p) ) allocate(znu_hyd_p(ims:ime,kms:kme,jms:jme) )
+ if(.not.allocated(psfc_hyd_p) ) allocate(psfc_hyd_p(ims:ime,jms:jme) )
+ if(.not.allocated(psfc_hydd_p) ) allocate(psfc_hydd_p(ims:ime,jms:jme) )
+ if(.not.allocated(pres_hyd_p) ) allocate(pres_hyd_p(ims:ime,kms:kme,jms:jme) )
+ if(.not.allocated(pres_hydd_p) ) allocate(pres_hydd_p(ims:ime,kms:kme,jms:jme) )
+ if(.not.allocated(pres2_hyd_p) ) allocate(pres2_hyd_p(ims:ime,kms:kme,jms:jme) )
+ if(.not.allocated(pres2_hydd_p)) allocate(pres2_hydd_p(ims:ime,kms:kme,jms:jme))
+ if(.not.allocated(znu_hyd_p) ) allocate(znu_hyd_p(ims:ime,kms:kme,jms:jme) )
end subroutine allocate_forall_physics
@@ -97,10 +100,13 @@
if(allocated(qs_p) ) deallocate(qs_p )
if(allocated(qg_p) ) deallocate(qg_p )
- if(allocated(psfc_hyd_p) ) deallocate(psfc_hyd_p )
- if(allocated(pres_hyd_p) ) deallocate(pres_hyd_p )
- if(allocated(pres2_hyd_p)) deallocate(pres2_hyd_p )
- if(allocated(znu_hyd_p) ) deallocate(znu_hyd_p )
+ if(allocated(psfc_hyd_p) ) deallocate(psfc_hyd_p )
+ if(allocated(psfc_hydd_p) ) deallocate(psfc_hydd_p )
+ if(allocated(pres_hyd_p) ) deallocate(pres_hyd_p )
+ if(allocated(pres_hydd_p) ) deallocate(pres_hydd_p )
+ if(allocated(pres2_hyd_p) ) deallocate(pres2_hyd_p )
+ if(allocated(pres2_hydd_p)) deallocate(pres2_hydd_p)
+ if(allocated(znu_hyd_p) ) deallocate(znu_hyd_p )
end subroutine deallocate_forall_physics
@@ -125,7 +131,7 @@
real(kind=RKIND),dimension(:,:),pointer:: rho_zz,theta_m,qv,pressure_p,u,v,w
real(kind=RKIND),dimension(:,:),pointer:: qvs,rh
- real(kind=RKIND):: rho1,rho2,tem1,tem2
+ real(kind=RKIND):: rho_a,rho1,rho2,tem1,tem2
!---------------------------------------------------------------------------------------------
@@ -294,16 +300,21 @@
do i = its,ite
!pressure at w-points:
k = kte+1
- pres2_hyd_p(i,k,j) = pres2_p(i,k,j)
+ pres2_hyd_p(i,k,j) = pres2_p(i,k,j)
+ pres2_hydd_p(i,k,j) = pres2_p(i,k,j)
do k = kte,1,-1
- pres2_hyd_p(i,k,j) = pres2_hyd_p(i,k+1,j) + g*rho_p(i,k,j)*dz_p(i,k,j)
+ rho_a = rho_p(i,k,j) / (1.+qv_p(i,k,j))
+ pres2_hyd_p(i,k,j) = pres2_hyd_p(i,k+1,j) + g*rho_p(i,k,j)*dz_p(i,k,j)
+ pres2_hydd_p(i,k,j) = pres2_hydd_p(i,k+1,j) + g*rho_a*dz_p(i,k,j)
enddo
!pressure at theta-points:
do k = kte,1,-1
- pres_hyd_p(i,k,j) = 0.5*(pres2_hyd_p(i,k+1,j)+pres2_hyd_p(i,k,j))
+ pres_hyd_p(i,k,j) = 0.5*(pres2_hyd_p(i,k+1,j)+pres2_hyd_p(i,k,j))
+ pres_hydd_p(i,k,j) = 0.5*(pres2_hydd_p(i,k+1,j)+pres2_hydd_p(i,k,j))
enddo
!surface pressure:
psfc_hyd_p(i,j) = pres2_hyd_p(i,1,j)
+ psfc_hydd_p(i,j) = pres2_hydd_p(i,1,j)
!znu:
do k = kte,1,-1
znu_hyd_p(i,k,j) = pres_hyd_p(i,k,j) / psfc_hyd_p(i,j)
Modified: branches/atmphys_gfs_mmm/src/core_atmos_physics/mpas_atmphys_manager.F
===================================================================
--- branches/atmphys_gfs_mmm/src/core_atmos_physics/mpas_atmphys_manager.F 2013-02-28 22:21:22 UTC (rev 2525)
+++ branches/atmphys_gfs_mmm/src/core_atmos_physics/mpas_atmphys_manager.F 2013-02-28 22:35:04 UTC (rev 2526)
@@ -463,6 +463,7 @@
lsm_scheme = trim(config_lsm_scheme)
microp_scheme = trim(config_microp_scheme)
pbl_scheme = trim(config_pbl_scheme)
+ gwdo_scheme = trim(config_gwdo_scheme)
radt_cld_scheme = trim(config_radt_cld_scheme)
radt_lw_scheme = trim(config_radt_lw_scheme)
radt_sw_scheme = trim(config_radt_sw_scheme)
Modified: branches/atmphys_gfs_mmm/src/core_atmos_physics/mpas_atmphys_update_surface.F
===================================================================
--- branches/atmphys_gfs_mmm/src/core_atmos_physics/mpas_atmphys_update_surface.F 2013-02-28 22:21:22 UTC (rev 2525)
+++ branches/atmphys_gfs_mmm/src/core_atmos_physics/mpas_atmphys_update_surface.F 2013-02-28 22:35:04 UTC (rev 2526)
@@ -146,6 +146,9 @@
if(config_frac_seaice) then
if(xice(iCell).ne.xicem(iCell) .and. xicem(iCell).gt.xice_threshold) then
+ !Fractional values of sfc_albedo and sfc_emiss are valid according to the earlier
+ !fractional sea-ice fraction, xicem. We recompute them for the new sea-ice fraction,
+ !xice.
sfc_albedo(iCell) = 0.08 + (sfc_albedo(iCell) -0.08) * xice(iCell)/xicem(iCell)
sfc_emiss(iCell) = 0.98 + (sfc_emiss(iCell)-0.98) * xice(iCell)/xicem(iCell)
endif
@@ -178,12 +181,12 @@
sfc_emibck(iCell) = 0.98
elseif(xland(iCell).lt.1.5 .and. xice(iCell).lt.xice_threshold .and. &
- xicem(iCell).lt.xice_threshold) then
+ xicem(iCell).ge.xice_threshold) then
!sea-ice points turn to water points:
xicem(iCell) = xice(iCell)
xland(iCell) = 2.
- isltyp(iCell) = 16
+ isltyp(iCell) = 14
ivgtyp(iCell) = iswater
vegfra(iCell) = 0.
tmn(iCell) = sst(iCell)
Modified: branches/atmphys_gfs_mmm/src/core_atmos_physics/mpas_atmphys_vars.F
===================================================================
--- branches/atmphys_gfs_mmm/src/core_atmos_physics/mpas_atmphys_vars.F 2013-02-28 22:21:22 UTC (rev 2525)
+++ branches/atmphys_gfs_mmm/src/core_atmos_physics/mpas_atmphys_vars.F 2013-02-28 22:35:04 UTC (rev 2526)
@@ -14,6 +14,7 @@
character(len=StrKIND),public:: microp_scheme
character(len=StrKIND),public:: conv_deep_scheme
character(len=StrKIND),public:: conv_shallow_scheme
+ character(len=StrKIND),public:: gwdo_scheme
character(len=StrKIND),public:: lsm_scheme
character(len=StrKIND),public:: pbl_scheme
character(len=StrKIND),public:: radt_cld_scheme
@@ -100,10 +101,13 @@
!... arrays used for calculating the hydrostatic pressure and exner function:
real(kind=RKIND),dimension(:,:),allocatable:: &
- psfc_hyd_p !surface pressure [hPa]
+ psfc_hyd_p, &!surface pressure [hPa]
+ psfc_hydd_p !"dry" surface pressure [hPa]
real(kind=RKIND),dimension(:,:,:),allocatable:: &
pres_hyd_p, &!pressure located at theta levels [hPa]
- pres2_hyd_p, &!pressure located at w-velocity levels [hPa]
+ pres_hydd_p, &!"dry" pressure located at theta levels [hPa]
+ pres2_hyd_p, &!pressure located at w-velocity levels [hPa]
+ pres2_hydd_p, &!"dry" pressure located at w-velocity levels [hPa]
znu_hyd_p !(pres_hyd_p / P0) needed in the Tiedtke convection scheme [hPa]
!=============================================================================================
@@ -225,6 +229,30 @@
kzq_p !
!=============================================================================================
+!... variables and arrays related to parameterization of gravity wave drag over orography:
+!=============================================================================================
+
+ real(kind=RKIND),dimension(:,:),allocatable:: &
+ var2d_p, &!orographic variance (m2)
+ con_p, &!orographic convexity (m2)
+ oa1_p, &!orographic direction asymmetry function (-)
+ oa2_p, &!orographic direction asymmetry function (-)
+ oa3_p, &!orographic direction asymmetry function (-)
+ oa4_p, &!orographic direction asymmetry function (-)
+ ol1_p, &!orographic direction asymmetry function (-)
+ ol2_p, &!orographic direction asymmetry function (-)
+ ol3_p, &!orographic direction asymmetry function (-)
+ ol4_p !orographic direction asymmetry function (-)
+
+ real(kind=RKIND),dimension(:,:),allocatable:: &
+ dusfcg_p, &!vertically-integrated gwdo u-stress (Pa m s-1)
+ dvsfcg_p !vertically-integrated gwdo v -stress (Pa m s-1)
+
+ real(kind=RKIND),dimension(:,:,:),allocatable:: &
+ dtaux3d_p, &!gravity wave drag over orography u-stress (m s-1)
+ dtauy3d_p !gravity wave drag over orography u-stress (m s-1)
+
+!=============================================================================================
!... variables and arrays related to parameterization of surface layer:
!=============================================================================================
real(kind=RKIND),dimension(:,:),allocatable:: &
Modified: branches/atmphys_gfs_mmm/src/core_atmos_physics/physics_wrf/Makefile
===================================================================
--- branches/atmphys_gfs_mmm/src/core_atmos_physics/physics_wrf/Makefile 2013-02-28 22:21:22 UTC (rev 2525)
+++ branches/atmphys_gfs_mmm/src/core_atmos_physics/physics_wrf/Makefile 2013-02-28 22:35:04 UTC (rev 2526)
@@ -7,6 +7,7 @@
OBJS = \
libmassv.o \
+ module_bl_gwdo.o \
module_bl_ysu.o \
module_cam_shr_kind_mod.o \
module_cam_support.o \
Copied: branches/atmphys_gfs_mmm/src/core_atmos_physics/physics_wrf/module_bl_gwdo.F (from rev 2525, branches/atmos_physics/src/core_atmos_physics/physics_wrf/module_bl_gwdo.F)
===================================================================
--- branches/atmphys_gfs_mmm/src/core_atmos_physics/physics_wrf/module_bl_gwdo.F (rev 0)
+++ branches/atmphys_gfs_mmm/src/core_atmos_physics/physics_wrf/module_bl_gwdo.F 2013-02-28 22:35:04 UTC (rev 2526)
@@ -0,0 +1,728 @@
+! WRf:model_layer:physics
+!
+!
+!
+!
+!
+module module_bl_gwdo
+contains
+!
+!-------------------------------------------------------------------
+!
+ subroutine gwdo(u3d,v3d,t3d,qv3d,p3d,p3di,pi3d,z, &
+ rublten,rvblten, &
+ dtaux3d,dtauy3d,dusfcg,dvsfcg, &
+ var2d,oc12d,oa2d1,oa2d2,oa2d3,oa2d4,ol2d1,ol2d2,ol2d3,ol2d4, &
+ znu,znw,mut,p_top, &
+ cp,g,rd,rv,ep1,pi, &
+ dt,dx,kpbl2d,itimestep, &
+ areaCell, &
+ ids,ide, jds,jde, kds,kde, &
+ ims,ime, jms,jme, kms,kme, &
+ its,ite, jts,jte, kts,kte)
+!-------------------------------------------------------------------
+ implicit none
+!------------------------------------------------------------------------------
+!
+!-- u3d 3d u-velocity interpolated to theta points (m/s)
+!-- v3d 3d v-velocity interpolated to theta points (m/s)
+!-- t3d temperature (k)
+!-- qv3d 3d water vapor mixing ratio (kg/kg)
+!-- p3d 3d pressure (pa)
+!-- p3di 3d pressure (pa) at interface level
+!-- pi3d 3d exner function (dimensionless)
+!-- rublten u tendency due to
+! pbl parameterization (m/s/s)
+!-- rvblten v tendency due to
+!-- cp heat capacity at constant pressure for dry air (j/kg/k)
+!-- g acceleration due to gravity (m/s^2)
+!-- rd gas constant for dry air (j/kg/k)
+!-- z height above sea level (m)
+!-- rv gas constant for water vapor (j/kg/k)
+!-- dt time step (s)
+!-- dx model grid interval (m)
+!-- ep1 constant for virtual temperature (r_v/r_d - 1) (dimensionless)
+!-- ids start index for i in domain
+!-- ide end index for i in domain
+!-- jds start index for j in domain
+!-- jde end index for j in domain
+!-- kds start index for k in domain
+!-- kde end index for k in domain
+!-- ims start index for i in memory
+!-- ime end index for i in memory
+!-- jms start index for j in memory
+!-- jme end index for j in memory
+!-- kms start index for k in memory
+!-- kme end index for k in memory
+!-- its start index for i in tile
+!-- ite end index for i in tile
+!-- jts start index for j in tile
+!-- jte end index for j in tile
+!-- kts start index for k in tile
+!-- kte end index for k in tile
+!-------------------------------------------------------------------
+!
+ integer, intent(in ) :: ids,ide, jds,jde, kds,kde, &
+ ims,ime, jms,jme, kms,kme, &
+ its,ite, jts,jte, kts,kte
+ integer, intent(in ) :: itimestep
+!
+ real, intent(in ) :: dt,dx,cp,g,rd,rv,ep1,pi
+!
+ real, dimension( ims:ime, kms:kme, jms:jme ) , &
+ intent(in ) :: qv3d, &
+ p3d, &
+ pi3d, &
+ t3d, &
+ z
+ real, dimension( ims:ime, kms:kme, jms:jme ) , &
+ intent(in ) :: p3di
+!
+ real, dimension( ims:ime, kms:kme, jms:jme ) , &
+ intent(inout) :: rublten, &
+ rvblten
+ real, dimension( ims:ime, kms:kme, jms:jme ) , &
+ intent(inout) :: dtaux3d, &
+ dtauy3d
+!
+ real, dimension( ims:ime, kms:kme, jms:jme ) , &
+ intent(in ) :: u3d, &
+ v3d
+!
+ integer, dimension( ims:ime, jms:jme ) , &
+ intent(in ) :: kpbl2d
+ real, dimension( ims:ime, jms:jme ) , &
+ intent(inout ) :: dusfcg, &
+ dvsfcg
+!
+ real, dimension( ims:ime, jms:jme ) , &
+ intent(in ) :: var2d, &
+ oc12d, &
+ oa2d1,oa2d2,oa2d3,oa2d4, &
+ ol2d1,ol2d2,ol2d3,ol2d4
+!
+ real, dimension( ims:ime, jms:jme ) , &
+ optional , &
+ intent(in ) :: mut
+!
+ real, dimension( kms:kme ) , &
+ optional , &
+ intent(in ) :: znu, &
+ znw
+!
+ real, optional, intent(in ) :: p_top
+
+!MPAS specific (Laura D. Fowler):
+ real,intent(in),dimension(ims:ime,jms:jme),optional:: areaCell
+!MPAS specific end.
+
+!
+!local
+!
+ real, dimension( its:ite, kts:kte ) :: delprsi, &
+ pdh
+ real, dimension( its:ite, kts:kte+1 ) :: pdhi
+ real, dimension( its:ite, 4 ) :: oa4, &
+ ol4
+ integer :: i,j,k,kdt
+!
+ do j = jts,jte
+ if(present(mut))then
+! For ARW we will replace p and p8w with dry hydrostatic pressure
+ do k = kts,kte+1
+ do i = its,ite
+ if(k.le.kte)pdh(i,k) = mut(i,j)*znu(k) + p_top
+ pdhi(i,k) = mut(i,j)*znw(k) + p_top
+ enddo
+ enddo
+ else
+ do k = kts,kte+1
+ do i = its,ite
+ if(k.le.kte)pdh(i,k) = p3d(i,k,j)
+ pdhi(i,k) = p3di(i,k,j)
+ enddo
+ enddo
+ endif
+!
+ do k = kts,kte
+ do i = its,ite
+ delprsi(i,k) = pdhi(i,k)-pdhi(i,k+1)
+ enddo
+ enddo
+ do i = its,ite
+ oa4(i,1) = oa2d1(i,j)
+ oa4(i,2) = oa2d2(i,j)
+ oa4(i,3) = oa2d3(i,j)
+ oa4(i,4) = oa2d4(i,j)
+ ol4(i,1) = ol2d1(i,j)
+ ol4(i,2) = ol2d2(i,j)
+ ol4(i,3) = ol2d3(i,j)
+ ol4(i,4) = ol2d4(i,j)
+ enddo
+ call gwdo2d(dudt=rublten(ims,kms,j),dvdt=rvblten(ims,kms,j) &
+ ,dtaux2d=dtaux3d(ims,kms,j),dtauy2d=dtauy3d(ims,kms,j) &
+ ,u1=u3d(ims,kms,j),v1=v3d(ims,kms,j) &
+ ,t1=t3d(ims,kms,j),q1=qv3d(ims,kms,j) &
+ ,prsi=pdhi(its,kts),del=delprsi(its,kts) &
+ ,prsl=pdh(its,kts),prslk=pi3d(ims,kms,j) &
+ ,zl=z(ims,kms,j),rcl=1.0 &
+ ,dusfc=dusfcg(ims,j),dvsfc=dvsfcg(ims,j) &
+ ,var=var2d(ims,j),oc1=oc12d(ims,j) &
+ ,oa4=oa4,ol4=ol4 &
+ ,g=g,cp=cp,rd=rd,rv=rv,fv=ep1,pi=pi &
+ ,dxmeter=dx,deltim=dt &
+ ,kpbl=kpbl2d(ims,j),kdt=itimestep,lat=j &
+ ,areaCell=areaCell(ims,j) &
+ ,ids=ids,ide=ide, jds=jds,jde=jde, kds=kds,kde=kde &
+ ,ims=ims,ime=ime, jms=jms,jme=jme, kms=kms,kme=kme &
+ ,its=its,ite=ite, jts=jts,jte=jte, kts=kts,kte=kte )
+ enddo
+!
+!
+ end subroutine gwdo
+!
+!-------------------------------------------------------------------
+!
+!
+!
+!
+ subroutine gwdo2d(dudt,dvdt,dtaux2d,dtauy2d, &
+ u1,v1,t1,q1, &
+ prsi,del,prsl,prslk,zl,rcl, &
+ var,oc1,oa4,ol4,dusfc,dvsfc, &
+ g,cp,rd,rv,fv,pi,dxmeter,deltim,kpbl,kdt,lat, &
+ areaCell, &
+ ids,ide, jds,jde, kds,kde, &
+ ims,ime, jms,jme, kms,kme, &
+ its,ite, jts,jte, kts,kte)
+!-------------------------------------------------------------------
+!
+! this code handles the time tendencies of u v due to the effect of mountain
+! induced gravity wave drag from sub-grid scale orography. this routine
+! not only treats the traditional upper-level wave breaking due to mountain
+! variance (alpert 1988), but also the enhanced lower-tropospheric wave
+! breaking due to mountain convexity and asymmetry (kim and arakawa 1995).
+! thus, in addition to the terrain height data in a model grid gox,
+! additional 10-2d topographic statistics files are needed, including
+! orographic standard deviation (var), convexity (oc1), asymmetry (oa4)
+! and ol (ol4). these data sets are prepared based on the 30 sec usgs orography
+! hong (1999). the current scheme was implmented as in hong et al.(2008)
+!
+! coded by song-you hong and young-joon kim and implemented by song-you hong
+!
+! references:
+! hong et al. (2008), wea. and forecasting
+! kim and arakawa (1995), j. atmos. sci.
+! alpet et al. (1988), NWP conference.
+! hong (1999), NCEP office note 424.
+!
+! notice : comparible or lower resolution orography files than model resolution
+! are desirable in preprocess (wps) to prevent weakening of the drag
+!-------------------------------------------------------------------
+!
+! input
+! dudt (ims:ime,kms:kme) non-lin tendency for u wind component
+! dvdt (ims:ime,kms:kme) non-lin tendency for v wind component
+! u1(ims:ime,kms:kme) zonal wind / sqrt(rcl) m/sec at t0-dt
+! v1(ims:ime,kms:kme) meridional wind / sqrt(rcl) m/sec at t0-dt
+! t1(ims:ime,kms:kme) temperature deg k at t0-dt
+! q1(ims:ime,kms:kme) specific humidity at t0-dt
+!
+! rcl a scaling factor = reciprocal of square of cos(lat)
+! for mrf gsm. rcl=1 if u1 and v1 are wind components.
+! deltim time step secs
+! del(kts:kte) positive increment of pressure across layer (pa)
+!
+! output
+! dudt, dvdt wind tendency due to gwdo
+!
+!-------------------------------------------------------------------
+ implicit none
+!-------------------------------------------------------------------
+ integer :: kdt,lat,latd,lond, &
+ ids,ide, jds,jde, kds,kde, &
+ ims,ime, jms,jme, kms,kme, &
+ its,ite, jts,jte, kts,kte
+!
+ real :: g,rd,rv,fv,cp,pi,dxmeter,deltim,rcl
+ real :: dudt(ims:ime,kms:kme),dvdt(ims:ime,kms:kme), &
+ dtaux2d(ims:ime,kms:kme),dtauy2d(ims:ime,kms:kme), &
+ u1(ims:ime,kms:kme),v1(ims:ime,kms:kme), &
+ t1(ims:ime,kms:kme),q1(ims:ime,kms:kme), &
+ zl(ims:ime,kms:kme),prslk(ims:ime,kms:kme)
+ real :: prsl(its:ite,kts:kte),prsi(its:ite,kts:kte+1), &
+ del(its:ite,kts:kte)
+ real :: oa4(its:ite,4),ol4(its:ite,4)
+!
+ integer :: kpbl(ims:ime)
+ real :: var(ims:ime),oc1(ims:ime), &
+ dusfc(ims:ime),dvsfc(ims:ime)
+
+!MPAS specific (Laura D. Fowler): We take into accound the actual size of individual
+!grid-boxes:
+ real,intent(in),dimension(ims:ime),optional:: areaCell
+ real,dimension(its:ite):: cleff_area
+!MPAS specific end.
+
+! critical richardson number for wave breaking : ! larger drag with larger value
+!
+ real,parameter :: ric = 0.25
+!
+ real,parameter :: dw2min = 1.
+ real,parameter :: rimin = -100.
+ real,parameter :: bnv2min = 1.0e-5
+ real,parameter :: efmin = 0.0
+ real,parameter :: efmax = 10.0
+ real,parameter :: xl = 4.0e4
+ real,parameter :: critac = 1.0e-5
+ real,parameter :: gmax = 1.
+ real,parameter :: veleps = 1.0
+ real,parameter :: factop = 0.5
+ real,parameter :: frc = 1.0
+ real,parameter :: ce = 0.8
+ real,parameter :: cg = 0.5
+!
+! local variables
+!
+ integer :: i,k,lcap,lcapp1,nwd,idir,kpblmin,kpblmax, &
+ klcap,kp1,ikount,kk
+!
+ real :: rcs,rclcs,csg,fdir,cleff,cs,rcsks, &
+ wdir,ti,rdz,temp,tem2,dw2,shr2,bvf2,rdelks, &
+ wtkbj,coefm,tem,gfobnv,hd,fro,rim,temc,tem1,efact, &
+ temv,dtaux,dtauy
+!
+ logical :: ldrag(its:ite),icrilv(its:ite), &
+ flag(its:ite),kloop1(its:ite)
+!
+ real :: taub(its:ite),taup(its:ite,kts:kte+1), &
+ xn(its:ite),yn(its:ite), &
+ ubar(its:ite),vbar(its:ite), &
+ fr(its:ite),ulow(its:ite), &
+ rulow(its:ite),bnv(its:ite), &
+ oa(its:ite),ol(its:ite), &
+ roll(its:ite),dtfac(its:ite), &
+ brvf(its:ite),xlinv(its:ite), &
+ delks(its:ite),delks1(its:ite), &
+ bnv2(its:ite,kts:kte),usqj(its:ite,kts:kte), &
+ taud(its:ite,kts:kte),ro(its:ite,kts:kte), &
+ vtk(its:ite,kts:kte),vtj(its:ite,kts:kte), &
+ zlowtop(its:ite),velco(its:ite,kts:kte-1)
+!
+ integer :: kbl(its:ite),klowtop(its:ite), &
+ lowlv(its:ite)
+!
+ logical :: iope
+ integer,parameter :: mdir=8
+ integer :: nwdir(mdir)
+ data nwdir/6,7,5,8,2,3,1,4/
+!
+! initialize local variables
+!
+ kbl=0 ; klowtop=0 ; lowlv=0
+!
+!---- constants
+!
+ rcs = sqrt(rcl)
+ cs = 1. / sqrt(rcl)
+ csg = cs * g
+ lcap = kte
+ lcapp1 = lcap + 1
+ fdir = mdir / (2.0*pi)
+!
+!
+!!!!!!! cleff (subgrid mountain scale ) is highly tunable parameter
+!!!!!!! the bigger (smaller) value produce weaker (stronger) wave drag
+!
+ cleff = max(dxmeter,50.e3)
+!
+! initialize!!
+!
+ dtaux = 0.0
+ dtauy = 0.0
+ do k = kts,kte
+ do i = its,ite
+ usqj(i,k) = 0.0
+ bnv2(i,k) = 0.0
+ vtj(i,k) = 0.0
+ vtk(i,k) = 0.0
+ taup(i,k) = 0.0
+ taud(i,k) = 0.0
+ dtaux2d(i,k)= 0.0
+ dtauy2d(i,k)= 0.0
+ enddo
+ enddo
+ do i = its,ite
+ taup(i,kte+1) = 0.0
+ xlinv(i) = 1.0/xl
+ enddo
+!
+ do k = kts,kte
+ do i = its,ite
+ vtj(i,k) = t1(i,k) * (1.+fv*q1(i,k))
+ vtk(i,k) = vtj(i,k) / prslk(i,k)
+ ro(i,k) = 1./rd * prsl(i,k) / vtj(i,k) ! density kg/m**3
+ enddo
+ enddo
+!
+ do i = its,ite
+ zlowtop(i) = 2. * var(i)
+ enddo
+
+!MPAS specific (Laura D. Fowler)::
+ if(present(areaCell)) then
+ do i = its,ite
+ cleff_area(i) = max(sqrt(areaCell(i)),50.e3)
+ enddo
+ endif
+!MPAS specific end.
+
+!
+!--- determine new reference level > 2*var
+!
+ do i = its,ite
+ kloop1(i) = .true.
+ enddo
+ do k = kts+1,kte
+ do i = its,ite
+ if(kloop1(i).and.zl(i,k)-zl(i,1).ge.zlowtop(i)) then
+ klowtop(i) = k+1
+ kloop1(i) = .false.
+ endif
+ enddo
+ enddo
+!
+ kpblmax = 2
+ do i = its,ite
+ kbl(i) = max(2, kpbl(i))
+ kbl(i) = max(kbl(i), klowtop(i))
+ delks(i) = 1.0 / (prsi(i,1) - prsi(i,kbl(i)))
+ ubar (i) = 0.0
+ vbar (i) = 0.0
+ taup(i,1) = 0.0
+ oa(i) = 0.0
+ kpblmax = max(kpblmax,kbl(i))
+ flag(i) = .true.
+ lowlv(i) = 2
+ enddo
+ kpblmax = min(kpblmax+1,kte-1)
+!
+! compute low level averages within pbl
+!
+ do k = kts,kpblmax
+ do i = its,ite
+ if (k.lt.kbl(i)) then
+ rcsks = rcs * del(i,k) * delks(i)
+ ubar(i) = ubar(i) + rcsks * u1(i,k) ! pbl u mean
+ vbar(i) = vbar(i) + rcsks * v1(i,k) ! pbl v mean
+ endif
+ enddo
+ enddo
+!
+! figure out low-level horizontal wind direction
+!
+! nwd 1 2 3 4 5 6 7 8
+! wd w s sw nw e n ne se
+!
+ do i = its,ite
+ wdir = atan2(ubar(i),vbar(i)) + pi
+ idir = mod(nint(fdir*wdir),mdir) + 1
+ nwd = nwdir(idir)
+ oa(i) = (1-2*int( (nwd-1)/4 )) * oa4(i,mod(nwd-1,4)+1)
+ ol(i) = ol4(i,mod(nwd-1,4)+1)
+ enddo
+!
+ kpblmin = kte
+ do i = its,ite
+ kpblmin = min(kpblmin, kbl(i))
+ enddo
+!
+ do i = its,ite
+ if (oa(i).le.0.0) kbl(i) = kpbl(i) + 1
+ enddo
+!
+ do i = its,ite
+ delks(i) = 1.0 / (prsi(i,1) - prsi(i,kbl(i)))
+ delks1(i) = 1.0 / (prsl(i,1) - prsl(i,kbl(i)))
+ enddo
+!
+!--- saving richardson number in usqj for migwdi
+!
+ do k = kts,kte-1
+ do i = its,ite
+ ti = 2.0 / (t1(i,k)+t1(i,k+1))
+ rdz = 1./(zl(i,k+1) - zl(i,k))
+ tem1 = u1(i,k) - u1(i,k+1)
+ tem2 = v1(i,k) - v1(i,k+1)
+ dw2 = rcl*(tem1*tem1 + tem2*tem2)
+ shr2 = max(dw2,dw2min) * rdz * rdz
+ bvf2 = g*(g/cp+rdz*(vtj(i,k+1)-vtj(i,k))) * ti
+ usqj(i,k) = max(bvf2/shr2,rimin)
+ bnv2(i,k) = 2*g*rdz*(vtk(i,k+1)-vtk(i,k))/(vtk(i,k+1)+vtk(i,k))
+ bnv2(i,k) = max( bnv2(i,k), bnv2min )
+ enddo
+ enddo
+!
+!-----initialize arrays
+!
+ do i = its,ite
+ xn(i) = 0.0
+ yn(i) = 0.0
+ ubar (i) = 0.0
+ vbar (i) = 0.0
+ roll (i) = 0.0
+ taub (i) = 0.0
+ ulow (i) = 0.0
+ dtfac(i) = 1.0
+ ldrag(i) = .false.
+ icrilv(i) = .false. ! initialize critical level control vector
+ enddo
+!
+!---- compute low level averages
+!---- (u,v)*cos(lat) use uv=(u1,v1) which is wind at t0-1
+!---- use rcs=1/cos(lat) to get wind field
+!
+ do k = 1,kpblmax
+ do i = its,ite
+ if (k .lt. kbl(i)) then
+ rdelks = del(i,k) * delks(i)
+ rcsks = rcs * rdelks
+ ubar(i) = ubar(i) + rcsks * u1(i,k) ! u mean
+ vbar(i) = vbar(i) + rcsks * v1(i,k) ! v mean
+ roll(i) = roll(i) + rdelks * ro(i,k) ! ro mean
+ endif
+ enddo
+ enddo
+!
+!----compute the "low level" or 1/3 wind magnitude (m/s)
+!
+ do i = its,ite
+ ulow(i) = max(sqrt(ubar(i)*ubar(i) + vbar(i)*vbar(i)), 1.0)
+ rulow(i) = 1./ulow(i)
+ enddo
+!
+ do k = kts,kte-1
+ do i = its,ite
+ velco(i,k) = (0.5*rcs) * ((u1(i,k)+u1(i,k+1)) * ubar(i) &
+ + (v1(i,k)+v1(i,k+1)) * vbar(i))
+ velco(i,k) = velco(i,k) * rulow(i)
+ if ((velco(i,k).lt.veleps) .and. (velco(i,k).gt.0.)) then
+ velco(i,k) = veleps
+ endif
+ enddo
+ enddo
+!
+! no drag when critical level in the base layer
+!
+ do i = its,ite
+ ldrag(i) = velco(i,1).le.0.
+ enddo
+!
+ do k = kts+1,kpblmax-1
+ do i = its,ite
+ if (k .lt. kbl(i)) ldrag(i) = ldrag(i).or. velco(i,k).le.0.
+ enddo
+ enddo
+!
+! no drag when bnv2.lt.0
+!
+ do k = kts,kpblmax-1
+ do i = its,ite
+ if (k .lt. kbl(i)) ldrag(i) = ldrag(i).or. bnv2(i,k).lt.0.
+ enddo
+ enddo
+!
+!-----the low level weighted average ri is stored in usqj(1,1; im)
+!-----the low level weighted average n**2 is stored in bnv2(1,1; im)
+!---- this is called bnvl2 in phys_gwd_alpert_sub not bnv2
+!---- rdelks (del(k)/delks) vert ave factor so we can * instead of /
+!
+ do i = its,ite
+ wtkbj = (prsl(i,1)-prsl(i,2)) * delks1(i)
+ bnv2(i,1) = wtkbj * bnv2(i,1)
+ usqj(i,1) = wtkbj * usqj(i,1)
+ enddo
+!
+ do k = kts+1,kpblmax-1
+ do i = its,ite
+ if (k .lt. kbl(i)) then
+ rdelks = (prsl(i,k)-prsl(i,k+1)) * delks1(i)
+ bnv2(i,1) = bnv2(i,1) + bnv2(i,k) * rdelks
+ usqj(i,1) = usqj(i,1) + usqj(i,k) * rdelks
+ endif
+ enddo
+ enddo
+!
+ do i = its,ite
+ ldrag(i) = ldrag(i) .or. bnv2(i,1).le.0.0
+ ldrag(i) = ldrag(i) .or. ulow(i).eq.1.0
+ ldrag(i) = ldrag(i) .or. var(i) .le. 0.0
+ enddo
+!
+! ----- set all ri low level values to the low level value
+!
+ do k = kts+1,kpblmax-1
+ do i = its,ite
+ if (k .lt. kbl(i)) usqj(i,k) = usqj(i,1)
+ enddo
+ enddo
+!
+ do i = its,ite
+ if (.not.ldrag(i)) then
+ bnv(i) = sqrt( bnv2(i,1) )
+ fr(i) = bnv(i) * rulow(i) * var(i)
+ xn(i) = ubar(i) * rulow(i)
+ yn(i) = vbar(i) * rulow(i)
+ endif
+ enddo
+!
+! compute the base level stress and store it in taub
+! calculate enhancement factor, number of mountains & aspect
+! ratio const. use simplified relationship between standard
+! deviation & critical hgt
+!
+ do i = its,ite
+ if (.not. ldrag(i)) then
+ efact = (oa(i) + 2.) ** (ce*fr(i)/frc)
+ efact = min( max(efact,efmin), efmax )
+ coefm = (1. + ol(i)) ** (oa(i)+1.)
+ if(present(areaCell)) then
+ xlinv(i) = coefm / cleff_area(i)
+ else
+ xlinv(i) = coefm / cleff
+ endif
+ tem = fr(i) * fr(i) * oc1(i)
+ gfobnv = gmax * tem / ((tem + cg)*bnv(i))
+ taub(i) = xlinv(i) * roll(i) * ulow(i) * ulow(i) &
+ * ulow(i) * gfobnv * efact
+ else
+ taub(i) = 0.0
+ xn(i) = 0.0
+ yn(i) = 0.0
+ endif
+ enddo
+!
+! now compute vertical structure of the stress.
+!
+!----set up bottom values of stress
+!
+ do k = kts,kpblmax
+ do i = its,ite
+ if (k .le. kbl(i)) taup(i,k) = taub(i)
+ enddo
+ enddo
+!
+ do k = kpblmin, kte-1 ! vertical level k loop!
+ kp1 = k + 1
+ do i = its,ite
+!
+!-----unstablelayer if ri < ric
+!-----unstable layer if upper air vel comp along surf vel <=0 (crit lay)
+!---- at (u-c)=0. crit layer exists and bit vector should be set (.le.)
+!
+ if (k .ge. kbl(i)) then
+ icrilv(i) = icrilv(i) .or. ( usqj(i,k) .lt. ric) &
+ .or. (velco(i,k) .le. 0.0)
+ brvf(i) = max(bnv2(i,k),bnv2min) ! brunt-vaisala frequency squared
+ brvf(i) = sqrt(brvf(i)) ! brunt-vaisala frequency
+ endif
+ enddo
+!
+ do i = its,ite
+ if (k .ge. kbl(i) .and. (.not. ldrag(i))) then
+ if (.not.icrilv(i) .and. taup(i,k) .gt. 0.0 ) then
+ temv = 1.0 / velco(i,k)
+ tem1 = xlinv(i)*(ro(i,kp1)+ro(i,k))*brvf(i)*velco(i,k)*0.5
+ hd = sqrt(taup(i,k) / tem1)
+ fro = brvf(i) * hd * temv
+!
+! rim is the minimum-richardson number by shutts (1985)
+!
+ tem2 = sqrt(usqj(i,k))
+ tem = 1. + tem2 * fro
+ rim = usqj(i,k) * (1.-fro) / (tem * tem)
+!
+! check stability to employ the 'saturation hypothesis'
+! of lindzen (1981) except at tropospheric downstream regions
+!
+ if (rim .le. ric) then ! saturation hypothesis!
+ if ((oa(i) .le. 0. .or. kp1 .ge. lowlv(i) )) then
+ temc = 2.0 + 1.0 / tem2
+ hd = velco(i,k) * (2.*sqrt(temc)-temc) / brvf(i)
+ taup(i,kp1) = tem1 * hd * hd
+ endif
+ else ! no wavebreaking!
+ taup(i,kp1) = taup(i,k)
+ endif
+ endif
+ endif
+ enddo
+ enddo
+!
+ if(lcap.lt.kte) then
+ do klcap = lcapp1,kte
+ do i = its,ite
+ taup(i,klcap) = prsi(i,klcap) / prsi(i,lcap) * taup(i,lcap)
+ enddo
+ enddo
+ endif
+!
+! calculate - (g)*d(tau)/d(pressure) and deceleration terms dtaux, dtauy
+!
+ do k = kts,kte
+ do i = its,ite
+ taud(i,k) = 1. * (taup(i,k+1) - taup(i,k)) * csg / del(i,k)
+ enddo
+ enddo
+!
+!------limit de-acceleration (momentum deposition ) at top to 1/2 value
+!------the idea is some stuff must go out the 'top'
+!
+ do klcap = lcap,kte
+ do i = its,ite
+ taud(i,klcap) = taud(i,klcap) * factop
+ enddo
+ enddo
+!
+!------if the gravity wave drag would force a critical line
+!------in the lower ksmm1 layers during the next deltim timestep,
+!------then only apply drag until that critical line is reached.
+!
+ do k = kts,kpblmax-1
+ do i = its,ite
+ if (k .le. kbl(i)) then
+ if(taud(i,k).ne.0.) &
+ dtfac(i) = min(dtfac(i),abs(velco(i,k) &
+ /(deltim*rcs*taud(i,k))))
+ endif
+ enddo
+ enddo
+!
+ do i = its,ite
+ dusfc(i) = 0.
+ dvsfc(i) = 0.
+ enddo
+!
+ do k = kts,kte
+ do i = its,ite
+ taud(i,k) = taud(i,k) * dtfac(i)
+ dtaux = taud(i,k) * xn(i)
+ dtauy = taud(i,k) * yn(i)
+ dtaux2d(i,k) = dtaux
+ dtauy2d(i,k) = dtauy
+ dudt(i,k) = dtaux + dudt(i,k)
+ dvdt(i,k) = dtauy + dvdt(i,k)
+ dusfc(i) = dusfc(i) + dtaux * del(i,k)
+ dvsfc(i) = dvsfc(i) + dtauy * del(i,k)
+ enddo
+ enddo
+!
+ do i = its,ite
+ dusfc(i) = (-1./g*rcs) * dusfc(i)
+ dvsfc(i) = (-1./g*rcs) * dvsfc(i)
+ enddo
+!
+ return
+ end subroutine gwdo2d
+!-------------------------------------------------------------------
+end module module_bl_gwdo
Modified: branches/atmphys_gfs_mmm/src/core_atmos_physics/physics_wrf/module_bl_ysu.F
===================================================================
--- branches/atmphys_gfs_mmm/src/core_atmos_physics/physics_wrf/module_bl_ysu.F 2013-02-28 22:21:22 UTC (rev 2525)
+++ branches/atmphys_gfs_mmm/src/core_atmos_physics/physics_wrf/module_bl_ysu.F 2013-02-28 22:35:04 UTC (rev 2526)
@@ -639,6 +639,17 @@
delta(i) = 0.0
enddo
!
+!... ldf:
+ do k = kts,kte
+ do i = its,ite
+ xkzh(i,k) = 0.0
+ xkzm(i,k) = 0.0
+ xkzhl(i,k) = 0.0
+ xkzml(i,k) = 0.0
+ enddo
+ enddo
+!... end ldf.
+!
do k = kts,klpbl
do i = its,ite
wscalek(i,k) = 0.0
Modified: branches/atmphys_gfs_mmm/src/core_atmos_physics/physics_wrf/module_cu_tiedtke.F
===================================================================
--- branches/atmphys_gfs_mmm/src/core_atmos_physics/physics_wrf/module_cu_tiedtke.F 2013-02-28 22:21:22 UTC (rev 2525)
+++ branches/atmphys_gfs_mmm/src/core_atmos_physics/physics_wrf/module_cu_tiedtke.F 2013-02-28 22:35:04 UTC (rev 2526)
@@ -63,7 +63,7 @@
VTMPC1=RV/RD-1.0, &
VTMPC2=CPV/CPD-1.0, &
CVDIFTS=1.0, &
- CEVAPCU1=1.93E-6*261., &
+ CEVAPCU1=1.93E-6*261.0*0.5/G, & ! Correction from WRFV3.4.1 sourcecode.
CEVAPCU2=1.E3/(38.3*0.293) )
Modified: branches/atmphys_gfs_mmm/src/core_init_nhyd_atmos/Makefile
===================================================================
--- branches/atmphys_gfs_mmm/src/core_init_nhyd_atmos/Makefile 2013-02-28 22:21:22 UTC (rev 2525)
+++ branches/atmphys_gfs_mmm/src/core_init_nhyd_atmos/Makefile 2013-02-28 22:35:04 UTC (rev 2526)
@@ -8,6 +8,7 @@
mpas_init_atm_bitarray.o \
mpas_init_atm_queue.o \
mpas_init_atm_hinterp.o \
+ mpas_init_atm_static.o \
mpas_init_atm_surface.o \
read_geogrid.o \
mpas_atmphys_date_time.o \
@@ -19,8 +20,15 @@
core_hyd: $(OBJS)
ar -ru libdycore.a $(OBJS)
-mpas_init_atm_test_cases.o: mpas_atm_advection.o mpas_init_atm_read_met.o read_geogrid.o mpas_init_atm_llxy.o mpas_init_atm_hinterp.o \
- mpas_init_atm_surface.o mpas_atmphys_initialize_real.o
+mpas_init_atm_test_cases.o: \
+ read_geogrid.o \
+ mpas_atm_advection.o \
+ mpas_init_atm_read_met.o \
+ mpas_init_atm_llxy.o \
+ mpas_init_atm_hinterp.o \
+ mpas_init_atm_static.o \
+ mpas_init_atm_surface.o \
+ mpas_atmphys_initialize_real.o
mpas_init_atm_hinterp.o: mpas_init_atm_queue.o mpas_init_atm_bitarray.o
@@ -34,6 +42,12 @@
mpas_init_atm_mpas_core.o: mpas_advection.o mpas_init_atm_test_cases.o
+mpas_init_atm_static.o: \
+ mpas_atm_advection.o \
+ mpas_init_atm_hinterp.o \
+ mpas_init_atm_llxy.o \
+ mpas_atmphys_utilities.o
+
mpas_init_atm_surface.o: \
mpas_init_atm_hinterp.o \
mpas_init_atm_llxy.o \
Modified: branches/atmphys_gfs_mmm/src/core_init_nhyd_atmos/Registry
===================================================================
--- branches/atmphys_gfs_mmm/src/core_init_nhyd_atmos/Registry 2013-02-28 22:21:22 UTC (rev 2525)
+++ branches/atmphys_gfs_mmm/src/core_init_nhyd_atmos/Registry 2013-02-28 22:35:04 UTC (rev 2526)
@@ -108,7 +108,7 @@
var persistent real edgeNormalVectors ( R3 nEdges ) 0 io edgeNormalVectors mesh - -
var persistent real localVerticalUnitVectors ( R3 nCells ) 0 io localVerticalUnitVectors mesh - -
-var persistent real cellTangentPlane ( R3 TWO nEdges ) 0 io cellTangentPlane mesh - -
+var persistent real cellTangentPlane ( R3 TWO nCells ) 0 io cellTangentPlane mesh - -
var persistent integer cellsOnCell ( maxEdges nCells ) 0 io cellsOnCell mesh - -
var persistent integer verticesOnCell ( maxEdges nCells ) 0 io verticesOnCell mesh - -
@@ -140,6 +140,17 @@
var persistent real shdmin ( nCells ) 0 io shdmin mesh - -
var persistent real shdmax ( nCells ) 0 io shdmax mesh - -
var persistent real albedo12m ( nMonths nCells ) 0 io albedo12m mesh - -
+var persistent real varsso ( nCells ) 0 io varsso mesh - -
+var persistent real var2d ( nCells ) 0 io var2d mesh - -
+var persistent real con ( nCells ) 0 io con mesh - -
+var persistent real oa1 ( nCells ) 0 io oa1 mesh - -
+var persistent real oa2 ( nCells ) 0 io oa2 mesh - -
+var persistent real oa3 ( nCells ) 0 io oa3 mesh - -
+var persistent real oa4 ( nCells ) 0 io oa4 mesh - -
+var persistent real ol1 ( nCells ) 0 io ol1 mesh - -
+var persistent real ol2 ( nCells ) 0 io ol2 mesh - -
+var persistent real ol3 ( nCells ) 0 io ol3 mesh - -
+var persistent real ol4 ( nCells ) 0 io ol4 mesh - -
% description of the vertical grid structure
@@ -162,10 +173,10 @@
% Horizontally interpolated from first-guess data
var persistent real u_fg ( nFGLevels nEdges Time ) 1 - u fg - -
var persistent real v_fg ( nFGLevels nEdges Time ) 1 - v fg - -
-var persistent real t_fg ( nFGLevels nCells Time ) 1 - t fg - -
-var persistent real p_fg ( nFGLevels nCells Time ) 1 - p fg - -
-var persistent real z_fg ( nFGLevels nCells Time ) 1 - z fg - -
-var persistent real rh_fg ( nFGLevels nCells Time ) 1 - rh fg - -
+var persistent real t_fg ( nFGLevels nCells Time ) 1 o t fg - -
+var persistent real p_fg ( nFGLevels nCells Time ) 1 o p fg - -
+var persistent real z_fg ( nFGLevels nCells Time ) 1 o z fg - -
+var persistent real rh_fg ( nFGLevels nCells Time ) 1 o rh fg - -
var persistent real soilz_fg ( nCells Time ) 1 io soilz fg - -
var persistent real psfc_fg ( nCells Time ) 1 - psfc fg - -
var persistent real pmsl_fg ( nCells Time ) 1 - pmsl fg - -
@@ -216,6 +227,8 @@
var persistent real qv_init ( nVertLevels ) 0 io qv_init mesh - -
% Diagnostic fields: only written to output
+var persistent real precipw ( nCells Time ) 1 o precipw diag_physics - -
+var persistent real rh ( nVertLevels nCells Time ) 1 o rh diag - -
var persistent real rho ( nVertLevels nCells Time ) 1 o rho diag - -
var persistent real theta ( nVertLevels nCells Time ) 1 o theta diag - -
var persistent real v ( nVertLevels nEdges Time ) 1 o v diag - -
Copied: branches/atmphys_gfs_mmm/src/core_init_nhyd_atmos/mpas_init_atm_static.F (from rev 2525, branches/atmos_physics/src/core_init_nhyd_atmos/mpas_init_atm_static.F)
===================================================================
--- branches/atmphys_gfs_mmm/src/core_init_nhyd_atmos/mpas_init_atm_static.F (rev 0)
+++ branches/atmphys_gfs_mmm/src/core_init_nhyd_atmos/mpas_init_atm_static.F 2013-02-28 22:35:04 UTC (rev 2526)
@@ -0,0 +1,1892 @@
+!==================================================================================================
+ module mpas_init_atm_static
+!==================================================================================================
+ use atm_advection
+ use mpas_configure
+ use mpas_dmpar
+ use init_atm_hinterp
+ use init_atm_llxy
+
+ use mpas_atmphys_utilities
+
+ implicit none
+ private
+ public:: init_atm_static, &
+ init_atm_static_orogwd, &
+ init_atm_check_read_error, &
+ nearest_cell, &
+ sphere_distance
+
+ contains
+
+!==================================================================================================
+ subroutine init_atm_static(mesh)
+!==================================================================================================
+
+!inout arguments:
+ type(mesh_type),intent(inout):: mesh
+
+!local variables:
+ type(proj_info):: proj
+ type(dm_info),pointer :: dminfo
+
+ character(len=StrKIND):: fname
+
+ integer:: nx,ny,nz
+ integer:: endian,isigned,istatus,wordsize
+ integer:: i,j,k
+ integer:: iCell,iPoint,iTileStart,iTileEnd,jTileStart,jTileEnd
+ integer,dimension(5) :: interp_list
+ integer,dimension(:),allocatable :: nhs
+ integer,dimension(:,:),allocatable:: ncat
+
+ real(kind=4):: scalefactor
+ real(kind=4),dimension(:,:,:),allocatable:: rarray
+
+ real(kind=RKIND):: r_earth
+ real(kind=RKIND):: lat,lon,x,y
+ real(kind=RKIND):: lat_pt,lon_pt
+ real(kind=RKIND):: dx,dy,known_lat,known_lon,known_x,known_y
+ real(kind=RKIND),dimension(:,:),allocatable :: soiltemp_1deg
+ real(kind=RKIND),dimension(:,:),allocatable :: maxsnowalb
+ real(kind=RKIND),dimension(:,:,:),allocatable:: vegfra
+
+!--------------------------------------------------------------------------------------------------
+ write(0,*)
+ write(0,*) '--- enter subroutine init_atm_static:'
+
+!
+! Scale all distances and areas from a unit sphere to one with radius sphere_radius
+!
+
+ r_earth = mesh % sphere_radius
+
+ mesh % xCell % array = mesh % xCell % array * r_earth
+ mesh % yCell % array = mesh % yCell % array * r_earth
+ mesh % zCell % array = mesh % zCell % array * r_earth
+ mesh % xVertex % array = mesh % xVertex % array * r_earth
+ mesh % yVertex % array = mesh % yVertex % array * r_earth
+ mesh % zVertex % array = mesh % zVertex % array * r_earth
+ mesh % xEdge % array = mesh % xEdge % array * r_earth
+ mesh % yEdge % array = mesh % yEdge % array * r_earth
+ mesh % zEdge % array = mesh % zEdge % array * r_earth
+ mesh % dvEdge % array = mesh % dvEdge % array * r_earth
+ mesh % dcEdge % array = mesh % dcEdge % array * r_earth
+ mesh % areaCell % array = mesh % areaCell % array * r_earth**2.0
+ mesh % areaTriangle % array = mesh % areaTriangle % array * r_earth**2.0
+ mesh % kiteAreasOnVertex % array = mesh % kiteAreasOnVertex % array * r_earth**2.0
+
+ call atm_initialize_advection_rk(mesh)
+ call atm_initialize_deformation_weights(mesh)
+
+!
+! Interpolate HGT
+!
+!nx = 126
+!ny = 126
+ nx = 1206
+ ny = 1206
+ nz = 1
+ isigned = 1
+ endian = 0
+ wordsize = 2
+ scalefactor = 1.0
+ allocate(rarray(nx,ny,nz))
+ allocate(nhs(mesh%nCells))
+ nhs(:) = 0
+ mesh%ter%array(:) = 0.0
+
+ do jTileStart = 1,20401,ny-6
+ jTileEnd = jTileStart + ny - 1 - 6
+
+ do iTileStart=1,42001,nx-6
+ iTileEnd = iTileStart + nx - 1 - 6
+ write(fname,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') trim(config_geog_data_path)// &
+ 'topo_30s/',iTileStart,'-',iTileEnd,'.',jTileStart,'-',jTileEnd
+ write(0,*) trim(fname)
+
+ call read_geogrid(fname,len_trim(fname),rarray,nx,ny,nz,isigned,endian, &
+ scalefactor,wordsize,istatus)
+ call init_atm_check_read_error(istatus, fname, dminfo)
+
+ iPoint = 1
+ do j=4,ny-3
+ do i=4,nx-3
+ lat_pt = -89.99583 + (jTileStart + j - 5) * 0.0083333333
+ lon_pt = -179.99583 + (iTileStart + i - 5) * 0.0083333333
+ lat_pt = lat_pt * PI / 180.0
+ lon_pt = lon_pt * PI / 180.0
+
+ iPoint = nearest_cell(lat_pt,lon_pt,iPoint,mesh%nCells,mesh%maxEdges, &
+ mesh%nEdgesOnCell%array,mesh%cellsOnCell%array, &
+ mesh%latCell%array,mesh%lonCell%array)
+ mesh%ter%array(iPoint) = mesh%ter%array(iPoint) + rarray(i,j,1)
+ nhs(iPoint) = nhs(iPoint) + 1
+ end do
+ end do
+
+ end do
+ end do
+
+ do iCell = 1,mesh%nCells
+ mesh%ter%array(iCell) = mesh%ter%array(iCell) / real(nhs(iCell))
+ end do
+ deallocate(rarray)
+ deallocate(nhs)
+ write(0,*) '--- end interpolate TER'
+
+
+!
+! Interpolate LU_INDEX
+!
+ nx = 1200
+ ny = 1200
+ nz = 1
+ isigned = 1
+ endian = 0
+ wordsize = 1
+ scalefactor = 1.0
+ allocate(rarray(nx,ny,nz))
+ allocate(ncat(24,mesh%nCells))
+ ncat(:,:) = 0
+ mesh%lu_index%array(:) = 0.0
+
+ do jTileStart = 1,20401,ny
+ jTileEnd = jTileStart + ny - 1
+
+ do iTileStart = 1,42001,nx
+ iTileEnd = iTileStart + nx - 1
+ write(fname,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') trim(config_geog_data_path)// &
+ '/landuse_30s/',iTileStart,'-',iTileEnd,'.',jTileStart,'-',jTileEnd
+ write(0,*) trim(fname)
+
+ call read_geogrid(fname,len_trim(fname),rarray,nx,ny,nz,isigned,endian, &
+ scalefactor,wordsize,istatus)
+ call init_atm_check_read_error(istatus, fname, dminfo)
+
+ iPoint = 1
+ do j=1,ny
+ do i=1,nx
+ lat_pt = -89.99583 + (jTileStart + j - 2) * 0.0083333333
+ lon_pt = -179.99583 + (iTileStart + i - 2) * 0.0083333333
+ lat_pt = lat_pt * PI / 180.0
+ lon_pt = lon_pt * PI / 180.0
+
+ iPoint = nearest_cell(lat_pt,lon_pt,iPoint,mesh%nCells,mesh%maxEdges, &
+ mesh%nEdgesOnCell%array,mesh%cellsOnCell%array, &
+ mesh%latCell%array,mesh%lonCell%array)
+ ncat(int(rarray(i,j,1)),iPoint) = ncat(int(rarray(i,j,1)),iPoint) + 1
+ end do
+ end do
+
+ end do
+ end do
+
+ do iCell = 1,mesh%nCells
+ mesh%lu_index%array(iCell) = 1
+ do i = 2,24
+ if(ncat(i,iCell) > ncat(mesh%lu_index%array(iCell),iCell)) then
+ mesh%lu_index%array(iCell) = i
+ end if
+ end do
+ end do
+ deallocate(rarray)
+ deallocate(ncat)
+ write(0,*) '--- end interpolate LU_INDEX'
+
+
+!
+! Interpolate SOILCAT_TOP
+!
+ nx = 1200
+ ny = 1200
+ nz = 1
+ isigned = 1
+ endian = 0
+ wordsize = 1
+ scalefactor = 1.0
+ allocate(rarray(nx,ny,nz))
+ allocate(ncat(16,mesh%nCells))
+ ncat(:,:) = 0
+ mesh%soilcat_top%array(:) = 0.0
+
+ do jTileStart = 1,20401,ny
+ jTileEnd = jTileStart + ny - 1
+
+ do iTileStart = 1,42001,nx
+ iTileEnd = iTileStart + nx - 1
+ write(fname,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') trim(config_geog_data_path)// &
+ '/soiltype_top_30s/',iTileStart,'-',iTileEnd,'.',jTileStart,'-',jTileEnd
+ write(0,*) trim(fname)
+
+ call read_geogrid(fname,len_trim(fname),rarray,nx,ny,nz,isigned,endian, &
+ scalefactor,wordsize,istatus)
+ call init_atm_check_read_error(istatus, fname, dminfo)
+
+ iPoint = 1
+ do j=1,ny
+ do i=1,nx
+ lat_pt = -89.99583 + (jTileStart + j - 2) * 0.0083333333
+ lon_pt = -179.99583 + (iTileStart + i - 2) * 0.0083333333
+ lat_pt = lat_pt * PI / 180.0
+ lon_pt = lon_pt * PI / 180.0
+
+ iPoint = nearest_cell(lat_pt,lon_pt,iPoint,mesh%nCells,mesh%maxEdges, &
+ mesh%nEdgesOnCell%array,mesh%cellsOnCell%array, &
+ mesh%latCell%array,mesh%lonCell%array)
+ ncat(int(rarray(i,j,1)),iPoint) = ncat(int(rarray(i,j,1)),iPoint) + 1
+ end do
+ end do
+
+ end do
+ end do
+
+ do iCell = 1,mesh%nCells
+ mesh%soilcat_top%array(iCell) = 1
+ do i = 2,16
+ if(ncat(i,iCell) > ncat(mesh%soilcat_top%array(iCell),iCell)) then
+ mesh%soilcat_top%array(iCell) = i
+ end if
+ end do
+ end do
+ deallocate(rarray)
+ deallocate(ncat)
+ write(0,*) '--- end interpolate SOILCAT_TOP'
+
+
+!
+! Interpolate SOILCAT_BOT
+!
+ nx = 1200
+ ny = 1200
+ nz = 1
+ isigned = 1
+ endian = 0
+ wordsize = 1
+ scalefactor = 1.0
+ allocate(rarray(nx,ny,nz))
+ allocate(ncat(16,mesh%nCells))
+ ncat(:,:) = 0
+ mesh%soilcat_bot%array(:) = 0.0
+
+ do jTileStart = 1,20401,ny
+ jTileEnd = jTileStart + ny - 1
+
+ do iTileStart = 1,42001,nx
+ iTileEnd = iTileStart + nx - 1
+ write(fname,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') trim(config_geog_data_path)// &
+ '/soiltype_bot_30s/',iTileStart,'-',iTileEnd,'.',jTileStart,'-',jTileEnd
+ write(0,*) trim(fname)
+
+ call read_geogrid(fname,len_trim(fname),rarray,nx,ny,nz,isigned,endian, &
+ scalefactor,wordsize,istatus)
+ call init_atm_check_read_error(istatus, fname, dminfo)
+
+ iPoint = 1
+ do j=1,ny
+ do i=1,nx
+ lat_pt = -89.99583 + (jTileStart + j - 2) * 0.0083333333
+ lon_pt = -179.99583 + (iTileStart + i - 2) * 0.0083333333
+ lat_pt = lat_pt * PI / 180.0
+ lon_pt = lon_pt * PI / 180.0
+
+ iPoint = nearest_cell(lat_pt,lon_pt,iPoint,mesh%nCells,mesh%maxEdges, &
+ mesh%nEdgesOnCell%array,mesh%cellsOnCell%array, &
+ mesh%latCell%array,mesh%lonCell%array)
+ ncat(int(rarray(i,j,1)),iPoint) = ncat(int(rarray(i,j,1)),iPoint) + 1
+ end do
+ end do
+
+ end do
+ end do
+
+ do iCell =1, mesh%nCells
+ mesh%soilcat_bot%array(iCell) = 1
+ do i = 2,16
+ if(ncat(i,iCell) > ncat(mesh%soilcat_bot%array(iCell),iCell)) then
+ mesh%soilcat_bot%array(iCell) = i
+ end if
+ end do
+ end do
+ deallocate(rarray)
+ deallocate(ncat)
+ write(0,*) '--- end interpolate SOILCAT_BOT'
+
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! KLUDGE TO FIX SOIL TYPE OVER ANTARCTICA
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ where (mesh%lu_index%array == 24) mesh%soilcat_top%array = 16
+ where (mesh%lu_index%array == 24) mesh%soilcat_bot%array = 16
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! CORRECT INCONSISTENT SOIL AND LAND USE DATA
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ do iCell = 1,mesh%nCells
+ if (mesh%lu_index%array(iCell) == 16 .or. &
+ mesh%soilcat_top%array(iCell) == 14 .or. &
+ mesh%soilcat_bot%array(iCell) == 14) then
+ if (mesh%lu_index%array(iCell) /= 16) then
+ write(0,*) 'Turning lu_index into water at ', iCell
+ mesh%lu_index%array(iCell) = 16
+ end if
+ if (mesh%soilcat_top%array(iCell) /= 14) then
+ write(0,*) 'Turning soilcat_top into water at ', iCell
+ mesh%soilcat_top%array(iCell) = 14
+ end if
+ if (mesh%soilcat_bot%array(iCell) /= 14) then
+ write(0,*) 'Turning soilcat_bot into water at ', iCell
+ mesh%soilcat_bot%array(iCell) = 14
+ end if
+ end if
+ end do
+
+
+!
+! Derive LANDMASK
+!
+ mesh%landmask%array(:) = 0
+ do iCell=1, mesh%nCells
+ if (mesh%lu_index%array(iCell) /= 16) mesh%landmask%array(iCell) = 1
+ end do
+ write(0,*) '--- end interpolate LANDMASK'
+
+
+!
+! Interpolate SOILTEMP:
+!
+ nx = 186
+ ny = 186
+ nz = 1
+ isigned = 0
+ endian = 0
+ wordsize = 2
+ scalefactor = 0.01
+ allocate(rarray(nx,ny,nz))
+ allocate(soiltemp_1deg(-2:363,-2:183))
+ mesh%soiltemp%array(:) = 0.0
+
+ call map_set(PROJ_LATLON, proj, &
+ latinc = 1.0_RKIND, &
+ loninc = 1.0_RKIND, &
+ knowni = 1.0_RKIND, &
+ knownj = 1.0_RKIND, &
+ lat1 = -89.5_RKIND, &
+ lon1 = -179.5_RKIND)
+
+ write(fname,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') trim(config_geog_data_path)// &
+ 'soiltemp_1deg/',1,'-',180,'.',1,'-',180
+ write(0,*) trim(fname)
+
+ call read_geogrid(fname,len_trim(fname),rarray,nx,ny,nz,isigned, endian, &
+ scalefactor,wordsize,istatus)
+ call init_atm_check_read_error(istatus, fname, dminfo)
+ soiltemp_1deg(-2:180,-2:183) = rarray(1:183,1:186,1)
+
+ write(fname,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') trim(config_geog_data_path)// &
+ 'soiltemp_1deg/',181,'-',360,'.',1,'-',180
+ write(0,*) trim(fname)
+
+ call read_geogrid(fname, len_trim(fname),rarray,nx,ny,nz,isigned,endian, &
+ scalefactor,wordsize,istatus)
+ call init_atm_check_read_error(istatus,fname,dminfo)
+ soiltemp_1deg(181:363,-2:183) = rarray(4:186,1:186,1)
+
+ interp_list(1) = FOUR_POINT
+ interp_list(2) = W_AVERAGE4
+ interp_list(3) = W_AVERAGE16
+ interp_list(4) = SEARCH
+ interp_list(5) = 0
+
+ do iCell = 1,mesh%nCells
+
+ if(mesh%landmask%array(iCell) == 1) then
+ lat = mesh % latCell % array(iCell) * DEG_PER_RAD
+ lon = mesh % lonCell % array(iCell) * DEG_PER_RAD
+ call latlon_to_ij(proj, lat, lon, x, y)
+ if(x < 0.5) then
+ lon = lon + 360.0
+ call latlon_to_ij(proj, lat, lon, x, y)
+ else if (x >= 360.5) then
+ lon = lon - 360.0
+ call latlon_to_ij(proj, lat, lon, x, y)
+ end if
+ if (y < 1.0) y = 1.0
+ if (y > 179.0) y = 179.0
+ mesh%soiltemp%array(iCell) = interp_sequence(x,y,1,soiltemp_1deg,-2,363,-2,183, &
+ 1,1,0.0_RKIND,interp_list,1)
+ else
+ mesh%soiltemp%array(iCell) = 0.0
+ end if
+
+ end do
+ deallocate(rarray)
+ deallocate(soiltemp_1deg)
+ write(0,*) '--- end interpolate SOILTEMP'
+
+
+!
+! Interpolate SNOALB
+!
+ nx = 186
+ ny = 186
+ nz = 1
+ isigned = 0
+ endian = 0
+ wordsize = 1
+ scalefactor = 1.0
+ allocate(rarray(nx,ny,nz))
+ allocate(maxsnowalb(-2:363,-2:183))
+ mesh%snoalb%array(:) = 0.0
+
+ call map_set(PROJ_LATLON, proj, &
+ latinc = 1.0_RKIND, &
+ loninc = 1.0_RKIND, &
+ knowni = 1.0_RKIND, &
+ knownj = 1.0_RKIND, &
+ lat1 = -89.5_RKIND, &
+ lon1 = -179.5_RKIND)
+
+ write(fname,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') trim(config_geog_data_path)// &
+ 'maxsnowalb/',1,'-',180,'.',1,'-',180
+ write(0,*) trim(fname)
+
+ call read_geogrid(fname,len_trim(fname),rarray,nx,ny,nz,isigned,endian, &
+ scalefactor,wordsize,istatus)
+ call init_atm_check_read_error(istatus,fname,dminfo)
+ maxsnowalb(-2:180,-2:183) = rarray(1:183,1:186,1)
+
+ write(fname,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') trim(config_geog_data_path)// &
+ 'maxsnowalb/',181,'-',360,'.',1,'-',180
+ write(0,*) trim(fname)
+
+ call read_geogrid(fname,len_trim(fname),rarray,nx,ny,nz,isigned,endian, &
+ scalefactor,wordsize,istatus)
+ call init_atm_check_read_error(istatus, fname, dminfo)
+ maxsnowalb(181:363,-2:183) = rarray(4:186,1:186,1)
+
+ interp_list(1) = FOUR_POINT
+ interp_list(2) = W_AVERAGE4
+ interp_list(3) = W_AVERAGE16
+ interp_list(4) = SEARCH
+ interp_list(5) = 0
+
+ do iCell = 1,mesh%nCells
+
+ if(mesh%landmask%array(iCell) == 1) then
+ lat = mesh % latCell % array(iCell) * DEG_PER_RAD
+ lon = mesh % lonCell % array(iCell) * DEG_PER_RAD
+ call latlon_to_ij(proj, lat, lon, x, y)
+ if(x < 0.5) then
+ lon = lon + 360.0
+ call latlon_to_ij(proj, lat, lon, x, y)
+ else if (x >= 360.5) then
+ lon = lon - 360.0
+ call latlon_to_ij(proj, lat, lon, x, y)
+ end if
+ if (y < 1.0) y = 1.0
+ if (y > 179.0) y = 179.0
+ mesh%snoalb%array(iCell) = interp_sequence(x,y,1,maxsnowalb,-2,363,-2,183, &
+ 1,1,0.0_RKIND,interp_list,1)
+ else
+ mesh%snoalb%array(iCell) = 0.0
+ end if
+
+ end do
+ mesh%snoalb%array(:) = mesh%snoalb%array(:) / 100.0
+ deallocate(rarray)
+ deallocate(maxsnowalb)
+ write(0,*) '--- end interpolate SNOALB'
+
+
+!
+! Interpolate GREENFRAC
+!
+ nx = 1256
+ ny = 1256
+ nz = 12
+ isigned = 0
+ endian = 0
+ wordsize = 1
+ scalefactor = 1.0
+ allocate(rarray(nx,ny,nz))
+ allocate(vegfra(-2:2503,-2:1253,12))
+ mesh%greenfrac%array(:,:) = 0.0
+
+ call map_set(PROJ_LATLON, proj, &
+ latinc = 0.144_RKIND, &
+ loninc = 0.144_RKIND, &
+ knowni = 1.0_RKIND, &
+ knownj = 1.0_RKIND, &
+ lat1 = -89.928_RKIND, &
+ lon1 = -179.928_RKIND)
+
+ write(fname,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') trim(config_geog_data_path)// &
+ 'greenfrac/',1,'-',1250,'.',1,'-',1250
+ write(0,*) trim(fname)
+
+ call read_geogrid(fname,len_trim(fname),rarray,nx,ny,nz,isigned,endian, &
+ scalefactor,wordsize,istatus)
+ call init_atm_check_read_error(istatus,fname,dminfo)
+ vegfra(-2:1250,-2:1253,1:12) = rarray(1:1253,1:1256,1:12)
+
+ write(fname,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') trim(config_geog_data_path)// &
+ 'greenfrac/',1251,'-',2500,'.',1,'-',1250
+ write(0,*) trim(fname)
+
+ call read_geogrid(fname,len_trim(fname),rarray,nx,ny,nz,isigned,endian, &
+ scalefactor,wordsize,istatus)
+ call init_atm_check_read_error(istatus,fname,dminfo)
+ vegfra(1251:2503,-2:1253,1:12) = rarray(4:1256,1:1256,1:12)
+
+ do iCell = 1,mesh%nCells
+
+ if (mesh%landmask%array(iCell) == 1) then
+ lat = mesh % latCell % array(iCell) * DEG_PER_RAD
+ lon = mesh % lonCell % array(iCell) * DEG_PER_RAD
+ call latlon_to_ij(proj, lat, lon, x, y)
+ if(x < 0.5) then
+ lon = lon + 360.0
+ call latlon_to_ij(proj, lat, lon, x, y)
+ else if(x >= 2500.5) then
+ lon = lon - 360.0
+ call latlon_to_ij(proj, lat, lon, x, y)
+ end if
+ if (y < 1.0) y = 1.0
+ if (y > 1249.0) y = 1249.0
+ do k = 1,12
+ mesh%greenfrac%array(k,iCell) = interp_sequence(x,y,k,vegfra,-2,2503,-2,1253, &
+ 1,12,-1.e30_RKIND,interp_list,1)
+ end do
+ else
+ mesh%greenfrac%array(:,iCell) = 0.0
+ end if
+ mesh%shdmin%array(iCell) = minval(mesh%greenfrac%array(:,iCell))
+ mesh%shdmax%array(iCell) = maxval(mesh%greenfrac%array(:,iCell))
+
+ end do
+ deallocate(rarray)
+ deallocate(vegfra)
+ write(0,*) '--- end interpolate GREENFRAC'
+
+
+!
+! Interpolate ALBEDO12M
+!
+ nx = 1256
+ ny = 1256
+ nz = 12
+ isigned = 0
+ endian = 0
+ wordsize = 1
+ scalefactor = 1.0
+ allocate(rarray(nx,ny,nz))
+ allocate(vegfra(-2:2503,-2:1253,12))
+ mesh%albedo12m%array(:,:) = 0.0
+
+ call map_set(PROJ_LATLON, proj, &
+ latinc = 0.144_RKIND, &
+ loninc = 0.144_RKIND, &
+ knowni = 1.0_RKIND, &
+ knownj = 1.0_RKIND, &
+ lat1 = -89.928_RKIND, &
+ lon1 = -179.928_RKIND)
+
+ write(fname,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') trim(config_geog_data_path)// &
+ 'albedo_ncep/',1,'-',1250,'.',1,'-',1250
+ write(0,*) trim(fname)
+
+ call read_geogrid(fname,len_trim(fname),rarray,nx,ny,nz,isigned,endian, &
+ scalefactor, wordsize, istatus)
+ call init_atm_check_read_error(istatus,fname, dminfo)
+ vegfra(-2:1250,-2:1253,1:12) = rarray(1:1253,1:1256,1:12)
+
+ write(fname,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') trim(config_geog_data_path)// &
+ 'albedo_ncep/',1251,'-',2500,'.',1,'-',1250
+ write(0,*) trim(fname)
+
+ call read_geogrid(fname,len_trim(fname),rarray,nx,ny,nz,isigned,endian, &
+ scalefactor,wordsize,istatus)
+ call init_atm_check_read_error(istatus,fname,dminfo)
+ vegfra(1251:2503,-2:1253,1:12) = rarray(4:1256,1:1256,1:12)
+
+ do iCell = 1,mesh%nCells
+
+ if (mesh%landmask%array(iCell) == 1) then
+ lat = mesh % latCell % array(iCell) * DEG_PER_RAD
+ lon = mesh % lonCell % array(iCell) * DEG_PER_RAD
+ call latlon_to_ij(proj, lat, lon, x, y)
+ if(x < 0.5) then
+ lon = lon + 360.0
+ call latlon_to_ij(proj, lat, lon, x, y)
+ else if(x >= 2500.5) then
+ lon = lon - 360.0
+ call latlon_to_ij(proj, lat, lon, x, y)
+ end if
+ if (y < 1.0) y = 1.0
+ if (y > 1249.0) y = 1249.0
+ do k = 1,12
+ mesh%albedo12m%array(k,iCell) = interp_sequence(x,y,k,vegfra,-2,2503,-2,1253, &
+ 1,12,0.0_RKIND,interp_list,1)
+ end do
+ else
+ mesh%albedo12m%array(:,iCell) = 8.0
+ end if
+ end do
+ deallocate(rarray)
+ deallocate(vegfra)
+ write(0,*) '--- end interpolate ALBEDO12M'
+
+
+ end subroutine init_atm_static
+
+!==================================================================================================
+ subroutine init_atm_static_orogwd(mesh)
+!==================================================================================================
+
+!inout arguments:
+ type(mesh_type),intent(inout):: mesh
+
+!local variables:
+ type(proj_info):: proj
+ type(dm_info),pointer :: dminfo
+
+ character(len=StrKIND):: mess
+ character(len=StrKIND):: fname
+ character(len=StrKIND):: dir_gwdo,fname_g
+
+ integer:: nx,ny,nz
+ integer:: endian,isigned,istatus,wordsize
+ integer:: i,j
+ integer:: iCell,iPoint,iTileStart,iTileEnd,jTileStart,jTileEnd
+ integer,dimension(5) :: interp_list
+ integer,dimension(:),allocatable:: nhs
+
+ real(kind=4):: scalefactor
+ real(kind=4),dimension(:,:,:),allocatable:: rarray
+
+ real(kind=RKIND):: lat,lon,x,y
+ real(kind=RKIND):: lat_pt,lon_pt
+ real(kind=RKIND):: dx,dy,known_lat,known_lon,known_x,known_y
+ real(kind=RKIND):: minMeshD,maxMeshD
+ real(kind=RKIND):: mindcEdge,maxdcEdge
+ real(kind=RKIND),dimension(:,:),allocatable:: xarray
+
+!--------------------------------------------------------------------------------------------------
+ write(0,*)
+ write(0,*) '--- enter subroutine init_atm_static_orogwd:'
+
+!goto 100
+!
+! Interpolate VARSSO:
+ mesh%varsso%array(:) = 0.0_RKIND
+ nx = 600
+ ny = 600
+ nz = 1
+ isigned = 0
+ endian = 0
+ wordsize = 4
+ scalefactor = 1.0
+
+ dx = 0.00833333
+ dy = 0.00833333
+ known_x = 1.0
+ known_y = 1.0
+ known_lat = -59.99583
+ known_lon = -179.99583
+
+ allocate(rarray(nx,ny,nz))
+ allocate(nhs(mesh%nCells))
+ nhs(:) = 0
+ rarray(:,:,:) = 0._RKIND
+ do jTileStart = 1,13801,ny
+ jTileEnd = jTileStart + ny - 1
+
+ do iTileStart = 1,42601,nx
+ iTileEnd = iTileStart + nx -1
+ write(fname,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') trim(config_geog_data_path)//'varsso/', &
+ iTileStart,'-',iTileEnd,'.',jTileStart,'-',jTileEnd
+ write(0,*) trim(fname)
+
+ call read_geogrid(fname,len_trim(fname),rarray,nx,ny,nz,isigned,endian, &
+ scalefactor,wordsize,istatus)
+ call init_atm_check_read_error(istatus,fname,dminfo)
+
+ iPoint = 1
+ do j = 1,ny
+ do i = 1,nx
+ lat_pt = known_lat + (jTileStart + j - 2) * dy
+ lon_pt = known_lon + (iTileStart + i - 2) * dx
+ lat_pt = lat_pt * PI / 180.0
+ lon_pt = lon_pt * PI / 180.0
+
+ iPoint = nearest_cell(lat_pt,lon_pt,iPoint,mesh%nCells,mesh%maxEdges, &
+ mesh%nEdgesOnCell%array,mesh%cellsOnCell%array, &
+ mesh%latCell%array,mesh%lonCell%array)
+ mesh%varsso%array(iPoint) = mesh%varsso%array(iPoint) + rarray(i,j,1)
+ nhs(iPoint) = nhs(iPoint) + 1
+ enddo
+ enddo
+
+ enddo
+ enddo
+
+ do iCell = 1,mesh%nCells
+ if(nhs(iCell) .gt. 0) &
+ mesh%varsso%array(iCell) = mesh%varsso%array(iCell) / real(nhs(iCell))
+ enddo
+ deallocate(rarray)
+ deallocate(nhs)
+ write(0,*) '--- end interpolate VARSSO'
+
+ 100 continue
+!... statistic fields needed for the parameterization of gravity wavwe drag over orography. The
+!input directory depends on the mesh resolution, and the mesh must be a uniform mesh.
+ minMeshD = minval(mesh%meshDensity%array(1:mesh%nCells))
+ maxMeshD = maxval(mesh%meshDensity%array(1:mesh%nCells))
+ mindcEdge = minval(mesh%dcEdge%array(1:mesh%nEdges))
+ maxdcEdge = maxval(mesh%dcEdge%array(1:mesh%nEdges))
+
+ write(0,*)
+ write(0,*) 'BEGIN INTERPOLATION OF STATISTICAL FIELDS FOR GRAVITY WAVE DRAG OVER OROGRAPHY'
+ write(0,*) 'min MeshD =', minMeshD
+ write(0,*) 'max MeshD =', maxMeshD
+ write(0,*) 'min dcEdge =', mindcEdge
+ write(0,*) 'max dcEdge =', maxdcEdge
+
+ dir_gwdo = ' '
+ if(minMeshD == 1.0_RKIND .and. maxMeshD == 1.0_RKIND) then
+ !... uniform 10242 mesh:
+ if(mindcEdge .ge. 200000._RKIND .and. maxdcEdge .lt. 260000._RKIND) then
+ dir_gwdo = 'orogwd_2deg'
+ elseif(mindcEdge .ge. 90000._RKIND .and. maxdcEdge .lt. 150000_RKIND) then
+ dir_gwdo = 'orogwd_1deg'
+ elseif(mindcEdge .ge. 40000._RKIND .and. maxdcEdge .lt. 70000._RKIND) then
+ dir_gwdo = 'orogwd_30m'
+ else
+ write(0,*)
+! write(mess,*) 'GWDO: Interpolation not available. The initialization will abort'
+! call physics_error_fatal(mess)
+ write(mess,*) 'GWDO: Interpolation not available. Set config_gwdo_scheme = .false.'
+ return
+ endif
+ else
+ write(0,*)
+! write(mess,*) 'GWDO: The input mesh must be a uniform mesh. The initialization will abort'
+! call physics_error_fatal(mess)
+ write(mess,*) 'GWDO: The input mesh must be a uniform mesh. Set config_gwdo_scheme = .false.'
+ return
+ endif
+ write(0,*) 'dir_gwdo = ', trim(dir_gwdo)
+ write(0,*)
+
+!
+! Interpolate CON:
+!
+ mesh%con%array(:) = 0.0_RKIND
+
+ con_select: select case(dir_gwdo)
+ case("orogwd_2deg")
+ nx = 180
+ ny = 90
+ nz = 1
+ isigned = 0
+ endian = 0
+ wordsize = 2
+ scalefactor = 0.025
+ dx = 2.0
+ dy = 2.0
+ known_x = 1.0
+ known_y = 1.0
+ known_lat = -89.0
+ known_lon = 1.0
+ case("orogwd_1deg")
+ nx = 360
+ ny = 180
+ nz = 1
+ isigned = 0
+ endian = 0
+ wordsize = 2
+ scalefactor = 0.025
+ dx = 1.0
+ dy = 1.0
+ known_x = 1.0
+ known_y = 1.0
+ known_lat = -89.5
+ known_lon = 0.5
+ case("orogwd_30m")
+ nx = 720
+ ny = 360
+ nz = 1
+ isigned = 0
+ endian = 0
+ wordsize = 2
+ scalefactor = 0.025
+ dx = 0.5
+ dy = 0.5
+ known_x = 1.0
+ known_y = 1.0
+ known_lat = -89.75
+ known_lon = 0.25
+ case("orogwd_10m")
+ nx = 2160
+ ny = 1080
+ nz = 1
+ isigned = 0
+ endian = 0
+ wordsize = 2
+ scalefactor = 0.025
+ dx = 0.16666667
+ dy = 0.16666667
+ known_x = 1.0
+ known_y = 1.0
+ known_lat = -89.916667
+ known_lon = 0.0833333
+ case default
+ end select con_select
+
+ write(fname,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') &
+ trim(config_geog_data_path)//trim(dir_gwdo)//'/con/',1,'-',nx,'.',1,'-',ny
+ write(0,*) trim(fname)
+
+ allocate(xarray(nx,ny))
+ allocate(rarray(nx,ny,nz))
+ call read_geogrid(fname,len_trim(fname),rarray,nx,ny,nz,isigned,endian, &
+ scalefactor,wordsize,istatus)
+ call init_atm_check_read_error(istatus,fname,dminfo)
+ xarray(1:nx,1:ny) = rarray(1:nx,1:ny,1)
+
+ call map_set(PROJ_LATLON, proj, &
+ latinc = dy, &
+ loninc = dx, &
+ knowni = known_x, &
+ knownj = known_y, &
+ lat1 = known_lat, &
+ lon1 = known_lon)
+
+ interp_list(1) = AVERAGE4
+ interp_list(2) = AVERAGE4
+ interp_list(3) = AVERAGE4
+ interp_list(4) = AVERAGE4
+ interp_list(5) = 0
+
+ do iCell = 1,mesh%nCells
+ if(mesh % landmask % array(iCell) == 1) then
+ lat = mesh % latCell % array(iCell) * DEG_PER_RAD
+ lon = mesh % lonCell % array(iCell) * DEG_PER_RAD
+ call latlon_to_ij(proj, lat, lon, x, y)
+ mesh % con % array(iCell) = interp_sequence(x,y,1,xarray,1,nx,1,ny,1,1, &
+ 0.0_RKIND,interp_list,1)
+ endif
+ enddo
+ deallocate(rarray)
+ deallocate(xarray)
+ write(0,*) '--- end interpolate CON'
+
+!
+! Interpolate OA1:
+!
+ mesh%oa1%array(:) = 0.0_RKIND
+
+ oa1_select: select case(dir_gwdo)
+ case("orogwd_2deg")
+ nx = 180
+ ny = 90
+ nz = 1
+ isigned = 1
+ endian = 0
+ wordsize = 2
+ scalefactor = 0.0001
+ dx = 2.0
+ dy = 2.0
+ known_x = 1.0
+ known_y = 1.0
+ known_lat = -89.0
+ known_lon = 1.0
+ case("orogwd_1deg")
+ nx = 360
+ ny = 180
+ nz = 1
+ isigned = 1
+ endian = 0
+ wordsize = 2
+ scalefactor = 0.0001
+ dx = 1.0
+ dy = 1.0
+ known_x = 1.0
+ known_y = 1.0
+ known_lat = -89.5
+ known_lon = 0.5
+ case("orogwd_30m")
+ nx = 720
+ ny = 360
+ nz = 1
+ isigned = 1
+ endian = 0
+ wordsize = 2
+ scalefactor = 0.0001
+ dx = 0.5
+ dy = 0.5
+ known_x = 1.0
+ known_y = 1.0
+ known_lat = -89.75
+ known_lon = 0.25
+ case("orogwd_10m")
+ nx = 2160
+ ny = 1080
+ nz = 1
+ isigned = 1
+ endian = 0
+ wordsize = 2
+ scalefactor = 0.0001
+ dx = 0.16666667
+ dy = 0.16666667
+ known_x = 1.0
+ known_y = 1.0
+ known_lat = -89.916667
+ known_lon = 0.0833333
+ case default
+ end select oa1_select
+
+ write(fname,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') &
+ trim(config_geog_data_path)//trim(dir_gwdo)//'/oa1/',1,'-',nx,'.',1,'-',ny
+ write(0,*) trim(fname)
+
+ allocate(xarray(nx,ny))
+ allocate(rarray(nx,ny,nz))
+ call read_geogrid(fname,len_trim(fname),rarray,nx,ny,nz,isigned,endian, &
+ scalefactor,wordsize,istatus)
+ call init_atm_check_read_error(istatus,fname,dminfo)
+ xarray(1:nx,1:ny) = rarray(1:nx,1:ny,1)
+
+ call map_set(PROJ_LATLON, proj, &
+ latinc = dy, &
+ loninc = dx, &
+ knowni = known_x, &
+ knownj = known_y, &
+ lat1 = known_lat, &
+ lon1 = known_lon)
+
+ interp_list(1) = AVERAGE4
+ interp_list(2) = AVERAGE4
+ interp_list(3) = AVERAGE4
+ interp_list(4) = AVERAGE4
+ interp_list(5) = 0
+
+ do iCell = 1,mesh%nCells
+ if(mesh % landmask % array(iCell) == 1) then
+ lat = mesh % latCell % array(iCell) * DEG_PER_RAD
+ lon = mesh % lonCell % array(iCell) * DEG_PER_RAD
+ call latlon_to_ij(proj, lat, lon, x, y)
+ mesh % oa1 % array(iCell) = interp_sequence(x,y,1,xarray,1,nx,1,ny,1,1, &
+ 0.0_RKIND,interp_list,1)
+ endif
+ enddo
+ deallocate(rarray)
+ deallocate(xarray)
+ write(0,*) '--- end interpolate OA1'
+
+!
+! Interpolate OA2:
+ mesh%oa2%array(:) = 0.0_RKIND
+
+ oa2_select: select case(dir_gwdo)
+ case("orogwd_2deg")
+ nx = 180
+ ny = 90
+ nz = 1
+ isigned = 1
+ endian = 0
+ wordsize = 2
+ scalefactor = 0.0001
+ dx = 2.0
+ dy = 2.0
+ known_x = 1.0
+ known_y = 1.0
+ known_lat = -89.0
+ known_lon = 1.0
+ case("orogwd_1deg")
+ nx = 360
+ ny = 180
+ nz = 1
+ isigned = 1
+ endian = 0
+ wordsize = 2
+ scalefactor = 0.0001
+ dx = 1.0
+ dy = 1.0
+ known_x = 1.0
+ known_y = 1.0
+ known_lat = -89.5
+ known_lon = 0.5
+ case("orogwd_30m")
+ nx = 720
+ ny = 360
+ nz = 1
+ isigned = 1
+ endian = 0
+ wordsize = 2
+ scalefactor = 0.0001
+ dx = 0.5
+ dy = 0.5
+ known_x = 1.0
+ known_y = 1.0
+ known_lat = -89.75
+ known_lon = 0.25
+ case("orogwd_10m")
+ nx = 2160
+ ny = 1080
+ nz = 1
+ isigned = 1
+ endian = 0
+ wordsize = 2
+ scalefactor = 0.0001
+ dx = 0.16666667
+ dy = 0.16666667
+ known_x = 1.0
+ known_y = 1.0
+ known_lat = -89.916667
+ known_lon = 0.0833333
+ case default
+ end select oa2_select
+
+ write(fname,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') &
+ trim(config_geog_data_path)//trim(dir_gwdo)//'/oa2/',1,'-',nx,'.',1,'-',ny
+ write(0,*) trim(fname)
+
+ allocate(xarray(nx,ny))
+ allocate(rarray(nx,ny,nz))
+ call read_geogrid(fname,len_trim(fname),rarray,nx,ny,nz,isigned,endian, &
+ scalefactor,wordsize,istatus)
+ call init_atm_check_read_error(istatus,fname,dminfo)
+ xarray(1:nx,1:ny) = rarray(1:nx,1:ny,1)
+
+ call map_set(PROJ_LATLON, proj, &
+ latinc = dy, &
+ loninc = dx, &
+ knowni = known_x, &
+ knownj = known_y, &
+ lat1 = known_lat, &
+ lon1 = known_lon)
+
+ interp_list(1) = AVERAGE4
+ interp_list(2) = AVERAGE4
+ interp_list(3) = AVERAGE4
+ interp_list(4) = AVERAGE4
+ interp_list(5) = 0
+
+ do iCell = 1,mesh%nCells
+ if(mesh % landmask % array(iCell) == 1) then
+ lat = mesh % latCell % array(iCell) * DEG_PER_RAD
+ lon = mesh % lonCell % array(iCell) * DEG_PER_RAD
+ call latlon_to_ij(proj, lat, lon, x, y)
+ mesh % oa2 % array(iCell) = interp_sequence(x,y,1,xarray,1,nx,1,ny,1,1, &
+ 0.0_RKIND,interp_list,1)
+ endif
+ enddo
+ deallocate(rarray)
+ deallocate(xarray)
+ write(0,*) '--- end interpolate OA2'
+
+!
+! Interpolate OA3:
+!
+ mesh%oa3%array(:) = 0.0_RKIND
+
+ oa3_select: select case(dir_gwdo)
+ case("orogwd_2deg")
+ nx = 180
+ ny = 90
+ nz = 1
+ isigned = 1
+ endian = 0
+ wordsize = 2
+ scalefactor = 0.0001
+ dx = 2.0
+ dy = 2.0
+ known_x = 1.0
+ known_y = 1.0
+ known_lat = -89.0
+ known_lon = 1.0
+ case("orogwd_1deg")
+ nx = 360
+ ny = 180
+ nz = 1
+ isigned = 1
+ endian = 0
+ wordsize = 2
+ scalefactor = 0.0001
+ dx = 1.0
+ dy = 1.0
+ known_x = 1.0
+ known_y = 1.0
+ known_lat = -89.5
+ known_lon = 0.5
+ case("orogwd_30m")
+ nx = 720
+ ny = 360
+ nz = 1
+ isigned = 1
+ endian = 0
+ wordsize = 2
+ scalefactor = 0.0001
+ dx = 0.5
+ dy = 0.5
+ known_x = 1.0
+ known_y = 1.0
+ known_lat = -89.75
+ known_lon = 0.25
+ case("orogwd_10m")
+ nx = 2160
+ ny = 1080
+ nz = 1
+ isigned = 1
+ endian = 0
+ wordsize = 2
+ scalefactor = 0.0001
+ dx = 0.16666667
+ dy = 0.16666667
+ known_x = 1.0
+ known_y = 1.0
+ known_lat = -89.916667
+ known_lon = 0.0833333
+ case default
+ end select oa3_select
+
+ write(fname,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') &
+ trim(config_geog_data_path)//trim(dir_gwdo)//'/oa3/',1,'-',nx,'.',1,'-',ny
+ write(0,*) trim(fname)
+
+ allocate(xarray(nx,ny))
+ allocate(rarray(nx,ny,nz))
+ call read_geogrid(fname,len_trim(fname),rarray,nx,ny,nz,isigned,endian, &
+ scalefactor,wordsize,istatus)
+ call init_atm_check_read_error(istatus,fname,dminfo)
+ xarray(1:nx,1:ny) = rarray(1:nx,1:ny,1)
+
+ call map_set(PROJ_LATLON, proj, &
+ latinc = dy, &
+ loninc = dx, &
+ knowni = known_x, &
+ knownj = known_y, &
+ lat1 = known_lat, &
+ lon1 = known_lon)
+
+ interp_list(1) = AVERAGE4
+ interp_list(2) = AVERAGE4
+ interp_list(3) = AVERAGE4
+ interp_list(4) = AVERAGE4
+ interp_list(5) = 0
+
+ do iCell = 1,mesh%nCells
+ if(mesh % landmask % array(iCell) == 1) then
+ lat = mesh % latCell % array(iCell) * DEG_PER_RAD
+ lon = mesh % lonCell % array(iCell) * DEG_PER_RAD
+ call latlon_to_ij(proj, lat, lon, x, y)
+ mesh % oa3 % array(iCell) = interp_sequence(x,y,1,xarray,1,nx,1,ny,1,1, &
+ 0.0_RKIND,interp_list,1)
+ endif
+ enddo
+ deallocate(rarray)
+ deallocate(xarray)
+ write(0,*) '--- end interpolate OA3'
+
+!
+! Interpolate OA4:
+!
+ mesh%oa4%array(:) = 0.0_RKIND
+
+ oa4_select: select case(dir_gwdo)
+ case("orogwd_2deg")
+ nx = 180
+ ny = 90
+ nz = 1
+ isigned = 1
+ endian = 0
+ wordsize = 2
+ scalefactor = 0.0001
+ dx = 2.0
+ dy = 2.0
+ known_x = 1.0
+ known_y = 1.0
+ known_lat = -89.0
+ known_lon = 1.0
+ case("orogwd_1deg")
+ nx = 360
+ ny = 180
+ nz = 1
+ isigned = 1
+ endian = 0
+ wordsize = 2
+ scalefactor = 0.0001
+ dx = 1.0
+ dy = 1.0
+ known_x = 1.0
+ known_y = 1.0
+ known_lat = -89.5
+ known_lon = 0.5
+ case("orogwd_30m")
+ nx = 720
+ ny = 360
+ nz = 1
+ isigned = 1
+ endian = 0
+ wordsize = 2
+ scalefactor = 0.0001
+ dx = 0.5
+ dy = 0.5
+ known_x = 1.0
+ known_y = 1.0
+ known_lat = -89.75
+ known_lon = 0.25
+ case("orogwd_10m")
+ nx = 2160
+ ny = 1080
+ nz = 1
+ isigned = 1
+ endian = 0
+ wordsize = 2
+ scalefactor = 0.0001
+ dx = 0.16666667
+ dy = 0.16666667
+ known_x = 1.0
+ known_y = 1.0
+ known_lat = -89.916667
+ known_lon = 0.0833333
+ case default
+ end select oa4_select
+
+ write(fname,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') &
+ trim(config_geog_data_path)//trim(dir_gwdo)//'/oa4/',1,'-',nx,'.',1,'-',ny
+ write(0,*) trim(fname)
+
+ allocate(xarray(nx,ny))
+ allocate(rarray(nx,ny,nz))
+ call read_geogrid(fname,len_trim(fname),rarray,nx,ny,nz,isigned,endian, &
+ scalefactor,wordsize,istatus)
+ call init_atm_check_read_error(istatus,fname,dminfo)
+ xarray(1:nx,1:ny) = rarray(1:nx,1:ny,1)
+
+ call map_set(PROJ_LATLON, proj, &
+ latinc = dy, &
+ loninc = dx, &
+ knowni = known_x, &
+ knownj = known_y, &
+ lat1 = known_lat, &
+ lon1 = known_lon)
+
+ interp_list(1) = AVERAGE4
+ interp_list(2) = AVERAGE4
+ interp_list(3) = AVERAGE4
+ interp_list(4) = AVERAGE4
+ interp_list(5) = 0
+
+ do iCell = 1,mesh%nCells
+ if(mesh % landmask % array(iCell) == 1) then
+ lat = mesh % latCell % array(iCell) * DEG_PER_RAD
+ lon = mesh % lonCell % array(iCell) * DEG_PER_RAD
+ call latlon_to_ij(proj, lat, lon, x, y)
+ mesh % oa4 % array(iCell) = interp_sequence(x,y,1,xarray,1,nx,1,ny,1,1, &
+ 0.0_RKIND,interp_list,1)
+ endif
+ enddo
+ deallocate(rarray)
+ deallocate(xarray)
+ write(0,*) '--- end interpolate OA4'
+
+!
+! Interpolate OL1:
+!
+ mesh%ol1%array(:) = 0.0_RKIND
+
+ ol1_select: select case(dir_gwdo)
+ case("orogwd_2deg")
+ nx = 180
+ ny = 90
+ nz = 1
+ isigned = 0
+ endian = 0
+ wordsize = 2
+ scalefactor = 0.0001
+ dx = 2.0
+ dy = 2.0
+ known_x = 1.0
+ known_y = 1.0
+ known_lat = -89.0
+ known_lon = 1.0
+ case("orogwd_1deg")
+ nx = 360
+ ny = 180
+ nz = 1
+ isigned = 0
+ endian = 0
+ wordsize = 2
+ scalefactor = 0.0001
+ dx = 1.0
+ dy = 1.0
+ known_x = 1.0
+ known_y = 1.0
+ known_lat = -89.5
+ known_lon = 0.5
+ case("orogwd_30m")
+ nx = 720
+ ny = 360
+ nz = 1
+ isigned = 0
+ endian = 0
+ wordsize = 2
+ scalefactor = 0.0001
+ dx = 0.5
+ dy = 0.5
+ known_x = 1.0
+ known_y = 1.0
+ known_lat = -89.75
+ known_lon = 0.25
+ case("orogwd_10m")
+ nx = 2160
+ ny = 1080
+ nz = 1
+ isigned = 0
+ endian = 0
+ wordsize = 2
+ scalefactor = 0.0001
+ dx = 0.16666667
+ dy = 0.16666667
+ known_x = 1.0
+ known_y = 1.0
+ known_lat = -89.916667
+ known_lon = 0.0833333
+ case default
+ end select ol1_select
+
+ write(fname,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') &
+ trim(config_geog_data_path)//trim(dir_gwdo)//'/ol1/',1,'-',nx,'.',1,'-',ny
+ write(0,*) trim(fname)
+
+ allocate(xarray(nx,ny))
+ allocate(rarray(nx,ny,nz))
+ call read_geogrid(fname,len_trim(fname),rarray,nx,ny,nz,isigned,endian, &
+ scalefactor,wordsize,istatus)
+ call init_atm_check_read_error(istatus,fname,dminfo)
+ xarray(1:nx,1:ny) = rarray(1:nx,1:ny,1)
+
+ call map_set(PROJ_LATLON, proj, &
+ latinc = dy, &
+ loninc = dx, &
+ knowni = known_x, &
+ knownj = known_y, &
+ lat1 = known_lat, &
+ lon1 = known_lon)
+
+ interp_list(1) = AVERAGE4
+ interp_list(2) = AVERAGE4
+ interp_list(3) = AVERAGE4
+ interp_list(4) = AVERAGE4
+ interp_list(5) = 0
+
+ do iCell = 1,mesh%nCells
+ if(mesh % landmask % array(iCell) == 1) then
+ lat = mesh % latCell % array(iCell) * DEG_PER_RAD
+ lon = mesh % lonCell % array(iCell) * DEG_PER_RAD
+ call latlon_to_ij(proj, lat, lon, x, y)
+ mesh % ol1 % array(iCell) = interp_sequence(x,y,1,xarray,1,nx,1,ny,1,1, &
+ 0.0_RKIND,interp_list,1)
+ endif
+ enddo
+ deallocate(rarray)
+ deallocate(xarray)
+ write(0,*) '--- end interpolate OL1'
+
+!
+! Interpolate OL2:
+!
+ mesh%ol2%array(:) = 0.0_RKIND
+
+ ol2_select: select case(dir_gwdo)
+ case("orogwd_2deg")
+ nx = 180
+ ny = 90
+ nz = 1
+ isigned = 0
+ endian = 0
+ wordsize = 2
+ scalefactor = 0.0001
+ dx = 2.0
+ dy = 2.0
+ known_x = 1.0
+ known_y = 1.0
+ known_lat = -89.0
+ known_lon = 1.0
+ case("orogwd_1deg")
+ nx = 360
+ ny = 180
+ nz = 1
+ isigned = 0
+ endian = 0
+ wordsize = 2
+ scalefactor = 0.0001
+ dx = 1.0
+ dy = 1.0
+ known_x = 1.0
+ known_y = 1.0
+ known_lat = -89.5
+ known_lon = 0.5
+ case("orogwd_30m")
+ nx = 720
+ ny = 360
+ nz = 1
+ isigned = 0
+ endian = 0
+ wordsize = 2
+ scalefactor = 0.0001
+ dx = 0.5
+ dy = 0.5
+ known_x = 1.0
+ known_y = 1.0
+ known_lat = -89.75
+ known_lon = 0.25
+ case("orogwd_10m")
+ nx = 2160
+ ny = 1080
+ nz = 1
+ isigned = 0
+ endian = 0
+ wordsize = 2
+ scalefactor = 0.0001
+ dx = 0.16666667
+ dy = 0.16666667
+ known_x = 1.0
+ known_y = 1.0
+ known_lat = -89.916667
+ known_lon = 0.0833333
+ case default
+ end select ol2_select
+
+ write(fname,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') &
+ trim(config_geog_data_path)//trim(dir_gwdo)//'/ol2/',1,'-',nx,'.',1,'-',ny
+ write(0,*) trim(fname)
+
+ allocate(xarray(nx,ny))
+ allocate(rarray(nx,ny,nz))
+ call read_geogrid(fname,len_trim(fname),rarray,nx,ny,nz,isigned,endian, &
+ scalefactor,wordsize,istatus)
+ call init_atm_check_read_error(istatus,fname,dminfo)
+ xarray(1:nx,1:ny) = rarray(1:nx,1:ny,1)
+
+ call map_set(PROJ_LATLON, proj, &
+ latinc = dy, &
+ loninc = dx, &
+ knowni = known_x, &
+ knownj = known_y, &
+ lat1 = known_lat, &
+ lon1 = known_lon)
+
+ interp_list(1) = AVERAGE4
+ interp_list(2) = AVERAGE4
+ interp_list(3) = AVERAGE4
+ interp_list(4) = AVERAGE4
+ interp_list(5) = 0
+
+ do iCell = 1,mesh%nCells
+ if(mesh % landmask % array(iCell) == 1) then
+ lat = mesh % latCell % array(iCell) * DEG_PER_RAD
+ lon = mesh % lonCell % array(iCell) * DEG_PER_RAD
+ call latlon_to_ij(proj, lat, lon, x, y)
+ mesh % ol2 % array(iCell) = interp_sequence(x,y,1,xarray,1,nx,1,ny,1,1, &
+ 0.0_RKIND,interp_list,1)
+ endif
+ enddo
+ deallocate(rarray)
+ deallocate(xarray)
+ write(0,*) '--- end interpolate OL2'
+
+!
+! Interpolate OL3:
+!
+ mesh%ol3%array(:) = 0.0_RKIND
+
+ ol3_select: select case(dir_gwdo)
+ case("orogwd_2deg")
+ nx = 180
+ ny = 90
+ nz = 1
+ isigned = 0
+ endian = 0
+ wordsize = 2
+ scalefactor = 0.0001
+ dx = 2.0
+ dy = 2.0
+ known_x = 1.0
+ known_y = 1.0
+ known_lat = -89.0
+ known_lon = 1.0
+ case("orogwd_1deg")
+ nx = 360
+ ny = 180
+ nz = 1
+ isigned = 0
+ endian = 0
+ wordsize = 2
+ scalefactor = 0.0001
+ dx = 1.0
+ dy = 1.0
+ known_x = 1.0
+ known_y = 1.0
+ known_lat = -89.5
+ known_lon = 0.5
+ case("orogwd_30m")
+ nx = 720
+ ny = 360
+ nz = 1
+ isigned = 0
+ endian = 0
+ wordsize = 2
+ scalefactor = 0.0001
+ dx = 0.5
+ dy = 0.5
+ known_x = 1.0
+ known_y = 1.0
+ known_lat = -89.75
+ known_lon = 0.25
+ case("orogwd_10m")
+ nx = 2160
+ ny = 1080
+ nz = 1
+ isigned = 0
+ endian = 0
+ wordsize = 2
+ scalefactor = 0.0001
+ dx = 0.16666667
+ dy = 0.16666667
+ known_x = 1.0
+ known_y = 1.0
+ known_lat = -89.916667
+ known_lon = 0.0833333
+ case default
+ end select ol3_select
+
+ write(fname,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') &
+ trim(config_geog_data_path)//trim(dir_gwdo)//'/ol3/',1,'-',nx,'.',1,'-',ny
+ write(0,*) trim(fname)
+
+ allocate(xarray(nx,ny))
+ allocate(rarray(nx,ny,nz))
+ call read_geogrid(fname,len_trim(fname),rarray,nx,ny,nz,isigned,endian, &
+ scalefactor,wordsize,istatus)
+ call init_atm_check_read_error(istatus,fname,dminfo)
+ xarray(1:nx,1:ny) = rarray(1:nx,1:ny,1)
+
+ call map_set(PROJ_LATLON, proj, &
+ latinc = dy, &
+ loninc = dx, &
+ knowni = known_x, &
+ knownj = known_y, &
+ lat1 = known_lat, &
+ lon1 = known_lon)
+
+ interp_list(1) = AVERAGE4
+ interp_list(2) = AVERAGE4
+ interp_list(3) = AVERAGE4
+ interp_list(4) = AVERAGE4
+ interp_list(5) = 0
+
+ do iCell = 1,mesh%nCells
+ if(mesh % landmask % array(iCell) == 1) then
+ lat = mesh % latCell % array(iCell) * DEG_PER_RAD
+ lon = mesh % lonCell % array(iCell) * DEG_PER_RAD
+ call latlon_to_ij(proj, lat, lon, x, y)
+ mesh % ol3 % array(iCell) = interp_sequence(x,y,1,xarray,1,nx,1,ny,1,1, &
+ 0.0_RKIND,interp_list,1)
+ endif
+ enddo
+ deallocate(rarray)
+ deallocate(xarray)
+ write(0,*) '--- end interpolate OL3'
+
+!
+! Interpolate OL4:
+!
+ mesh%ol4%array(:) = 0.0_RKIND
+
+ ol4_select: select case(dir_gwdo)
+ case("orogwd_2deg")
+ nx = 180
+ ny = 90
+ nz = 1
+ isigned = 0
+ endian = 0
+ wordsize = 2
+ scalefactor = 0.0001
+ dx = 2.0
+ dy = 2.0
+ known_x = 1.0
+ known_y = 1.0
+ known_lat = -89.0
+ known_lon = 1.0
+ case("orogwd_1deg")
+ nx = 360
+ ny = 180
+ nz = 1
+ isigned = 0
+ endian = 0
+ wordsize = 2
+ scalefactor = 0.0001
+ dx = 1.0
+ dy = 1.0
+ known_x = 1.0
+ known_y = 1.0
+ known_lat = -89.5
+ known_lon = 0.5
+ case("orogwd_30m")
+ nx = 720
+ ny = 360
+ nz = 1
+ isigned = 0
+ endian = 0
+ wordsize = 2
+ scalefactor = 0.0001
+ dx = 0.5
+ dy = 0.5
+ known_x = 1.0
+ known_y = 1.0
+ known_lat = -89.75
+ known_lon = 0.25
+ case("orogwd_10m")
+ nx = 2160
+ ny = 1080
+ nz = 1
+ isigned = 0
+ endian = 0
+ wordsize = 2
+ scalefactor = 0.0001
+ dx = 0.16666667
+ dy = 0.16666667
+ known_x = 1.0
+ known_y = 1.0
+ known_lat = -89.916667
+ known_lon = 0.0833333
+ case default
+ end select ol4_select
+
+ write(fname,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') &
+ trim(config_geog_data_path)//trim(dir_gwdo)//'/ol4/',1,'-',nx,'.',1,'-',ny
+ write(0,*) trim(fname)
+
+ allocate(xarray(nx,ny))
+ allocate(rarray(nx,ny,nz))
+ call read_geogrid(fname,len_trim(fname),rarray,nx,ny,nz,isigned,endian, &
+ scalefactor,wordsize,istatus)
+ call init_atm_check_read_error(istatus,fname,dminfo)
+ xarray(1:nx,1:ny) = rarray(1:nx,1:ny,1)
+
+ call map_set(PROJ_LATLON, proj, &
+ latinc = dy, &
+ loninc = dx, &
+ knowni = known_x, &
+ knownj = known_y, &
+ lat1 = known_lat, &
+ lon1 = known_lon)
+
+ interp_list(1) = AVERAGE4
+ interp_list(2) = AVERAGE4
+ interp_list(3) = AVERAGE4
+ interp_list(4) = AVERAGE4
+ interp_list(5) = 0
+
+ do iCell = 1,mesh%nCells
+ if(mesh % landmask % array(iCell) == 1) then
+ lat = mesh % latCell % array(iCell) * DEG_PER_RAD
+ lon = mesh % lonCell % array(iCell) * DEG_PER_RAD
+ call latlon_to_ij(proj, lat, lon, x, y)
+ mesh % ol4 % array(iCell) = interp_sequence(x,y,1,xarray,1,nx,1,ny,1,1, &
+ 0.0_RKIND,interp_list,1)
+ endif
+ enddo
+ deallocate(rarray)
+ deallocate(xarray)
+ write(0,*) '--- end interpolate OL4'
+
+!
+! Interpolate VAR2D:
+!
+ mesh%var2d%array(:) = 0.0_RKIND
+
+ var2d_select: select case(dir_gwdo)
+ case("orogwd_2deg")
+ nx = 180
+ ny = 90
+ nz = 1
+ isigned = 0
+ endian = 0
+ wordsize = 4
+ scalefactor = 0.02
+ dx = 2.0
+ dy = 2.0
+ known_x = 1.0
+ known_y = 1.0
+ known_lat = -89.0
+ known_lon = 1.0
+ case("orogwd_1deg")
+ nx = 360
+ ny = 180
+ nz = 1
+ isigned = 0
+ endian = 0
+ wordsize = 4
+ scalefactor = 0.02
+ dx = 1.0
+ dy = 1.0
+ known_x = 1.0
+ known_y = 1.0
+ known_lat = -89.5
+ known_lon = 0.5
+ case("orogwd_30m")
+ nx = 720
+ ny = 360
+ nz = 1
+ isigned = 0
+ endian = 0
+ wordsize = 4
+ scalefactor = 0.02
+ dx = 0.5
+ dy = 0.5
+ known_x = 1.0
+ known_y = 1.0
+ known_lat = -89.75
+ known_lon = 0.25
+ case("orogwd_10m")
+ nx = 2160
+ ny = 1080
+ nz = 1
+ isigned = 0
+ endian = 0
+ wordsize = 2
+ scalefactor = 0.02
+ dx = 0.16666667
+ dy = 0.16666667
+ known_x = 1.0
+ known_y = 1.0
+ known_lat = -89.916667
+ known_lon = 0.0833333
+ case default
+ end select var2d_select
+
+ write(fname,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') &
+ trim(config_geog_data_path)//trim(dir_gwdo)//'/var/',1,'-',nx,'.',1,'-',ny
+ write(0,*) trim(fname)
+
+
+ allocate(xarray(nx,ny))
+ allocate(rarray(nx,ny,nz))
+ call read_geogrid(fname,len_trim(fname),rarray,nx,ny,nz,isigned,endian, &
+ scalefactor,wordsize,istatus)
+ call init_atm_check_read_error(istatus,fname,dminfo)
+ xarray(1:nx,1:ny) = rarray(1:nx,1:ny,1)
+
+ call map_set(PROJ_LATLON, proj, &
+ latinc = dy, &
+ loninc = dx, &
+ knowni = known_x, &
+ knownj = known_y, &
+ lat1 = known_lat, &
+ lon1 = known_lon)
+
+ interp_list(1) = AVERAGE4
+ interp_list(2) = AVERAGE4
+ interp_list(3) = AVERAGE4
+ interp_list(4) = AVERAGE4
+ interp_list(5) = 0
+
+ do iCell = 1,mesh%nCells
+ if(mesh % landmask % array(iCell) == 1) then
+ lat = mesh % latCell % array(iCell) * DEG_PER_RAD
+ lon = mesh % lonCell % array(iCell) * DEG_PER_RAD
+ call latlon_to_ij(proj, lat, lon, x, y)
+ mesh % var2d % array(iCell) = interp_sequence(x,y,1,xarray,1,nx,1,ny,1,1, &
+ 0.0_RKIND,interp_list,1)
+ endif
+ enddo
+ deallocate(rarray)
+ deallocate(xarray)
+ write(0,*) '--- end interpolate VAR2D'
+
+ end subroutine init_atm_static_orogwd
+
+!==================================================================================================
+ subroutine init_atm_check_read_error(istatus, fname, dminfo)
+!==================================================================================================
+ implicit none
+
+ integer, intent(in) :: istatus
+ character (len=*), intent(in) :: fname
+ type (dm_info), intent(in) :: dminfo
+
+ if (istatus /= 0) then
+ write(0,*) 'ERROR: Could not read file '//trim(fname)
+ call mpas_dmpar_abort(dminfo)
+ end if
+
+ end subroutine init_atm_check_read_error
+
+!==================================================================================================
+ integer function nearest_cell(target_lat, target_lon, start_cell, nCells, maxEdges, &
+ nEdgesOnCell, cellsOnCell, latCell, lonCell)
+!==================================================================================================
+ implicit none
+
+ real (kind=RKIND), intent(in) :: target_lat, target_lon
+ integer, intent(in) :: start_cell
+ integer, intent(in) :: nCells, maxEdges
+ integer, dimension(nCells), intent(in) :: nEdgesOnCell
+ integer, dimension(maxEdges,nCells), intent(in) :: cellsOnCell
+ real (kind=RKIND), dimension(nCells), intent(in) :: latCell, lonCell
+
+ integer :: i
+ integer :: iCell
+ integer :: current_cell
+ real (kind=RKIND) :: current_distance, d
+ real (kind=RKIND) :: nearest_distance
+
+ nearest_cell = start_cell
+ current_cell = -1
+
+ do while (nearest_cell /= current_cell)
+ current_cell = nearest_cell
+ current_distance = sphere_distance(latCell(current_cell), lonCell(current_cell), target_lat, &
+ target_lon, 1.0_RKIND)
+ nearest_cell = current_cell
+ nearest_distance = current_distance
+ do i = 1, nEdgesOnCell(current_cell)
+ iCell = cellsOnCell(i,current_cell)
+ if (iCell <= nCells) then
+ d = sphere_distance(latCell(iCell), lonCell(iCell), target_lat, target_lon, 1.0_RKIND)
+ if (d < nearest_distance) then
+ nearest_cell = iCell
+ nearest_distance = d
+ end if
+ end if
+ end do
+ end do
+
+ end function nearest_cell
+
+!==================================================================================================
+ real (kind=RKIND) function sphere_distance(lat1, lon1, lat2, lon2, radius)
+
+!Compute the great-circle distance between (lat1, lon1) and (lat2, lon2) on a
+!sphere with given radius.
+!==================================================================================================
+ implicit none
+
+ real (kind=RKIND), intent(in) :: lat1, lon1, lat2, lon2, radius
+ real (kind=RKIND) :: arg1
+
+ arg1 = sqrt( sin(0.5*(lat2-lat1))**2 + &
+ cos(lat1)*cos(lat2)*sin(0.5*(lon2-lon1))**2 )
+ sphere_distance = 2.*radius*asin(arg1)
+
+ end function sphere_distance
+
+!==================================================================================================
+ end module mpas_init_atm_static
+!==================================================================================================
Modified: branches/atmphys_gfs_mmm/src/core_init_nhyd_atmos/mpas_init_atm_surface.F
===================================================================
--- branches/atmphys_gfs_mmm/src/core_init_nhyd_atmos/mpas_init_atm_surface.F 2013-02-28 22:21:22 UTC (rev 2525)
+++ branches/atmphys_gfs_mmm/src/core_init_nhyd_atmos/mpas_init_atm_surface.F 2013-02-28 22:35:04 UTC (rev 2526)
@@ -149,7 +149,7 @@
!sea-surface data:
if(index(field % field, 'SKINTEMP') /= 0 .or. index(field % field, 'SST') /= 0) then
- write(0,*) '... Processing SST:'
+! write(0,*) '... Processing SST:'
fg % sst % array(1:mesh%nCells) = 0.0_RKIND
destField1d => fg % sst % array
@@ -171,7 +171,7 @@
!sea-ice data:
elseif(index(field % field, 'SEAICE') /= 0) then
- write(0,*) '... Processing SEAICE:'
+! write(0,*) '... Processing SEAICE:'
fg % xice % array(1:mesh%nCells) = 0.0_RKIND
destField1d => fg % xice % array
Modified: branches/atmphys_gfs_mmm/src/core_init_nhyd_atmos/mpas_init_atm_test_cases.F
===================================================================
--- branches/atmphys_gfs_mmm/src/core_init_nhyd_atmos/mpas_init_atm_test_cases.F 2013-02-28 22:21:22 UTC (rev 2525)
+++ branches/atmphys_gfs_mmm/src/core_init_nhyd_atmos/mpas_init_atm_test_cases.F 2013-02-28 22:35:04 UTC (rev 2526)
@@ -10,6 +10,7 @@
use mpas_RBF_interpolation
use mpas_vector_reconstruction
use mpas_timer
+ use mpas_init_atm_static
use mpas_init_atm_surface
! Added only clause to keep xlf90 from getting confused from the overloaded abs intrinsic in mpas_timekeeping
@@ -105,10 +106,15 @@
write(0,*) ' real-data GFS test case '
block_ptr => domain % blocklist
do while (associated(block_ptr))
+ if (config_static_interp) then
+ call init_atm_static(block_ptr % mesh)
+ call init_atm_static_orogwd(block_ptr % mesh)
+ endif
call init_atm_test_case_gfs(block_ptr % mesh, block_ptr % fg, &
block_ptr % state % time_levs(1) % state, block_ptr % diag, &
- config_test_case)
+ block_ptr % diag_physics, config_test_case)
if (config_met_interp) call physics_initialize_real(block_ptr % mesh, block_ptr % fg, domain % dminfo)
+
block_ptr => block_ptr % next
end do
@@ -2232,7 +2238,7 @@
end subroutine init_atm_test_case_mtn_wave
- subroutine init_atm_test_case_gfs(grid, fg, state, diag, test_case)
+ subroutine init_atm_test_case_gfs(grid, fg, state, diag, diag_physics, test_case)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Real-data test case using GFS data
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -2248,6 +2254,7 @@
type (fg_type), intent(inout) :: fg
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
type (block_type), pointer :: block
@@ -2292,10 +2299,9 @@
!This is temporary variable here. It just need when calculate tangential velocity v.
integer :: eoe, j
- integer, dimension(:), pointer :: nEdgesOnEdge, nEdgesOnCell
+ integer, dimension(:), pointer :: nEdgesOnCell
integer, dimension(:,:), pointer :: edgesOnEdge, cellsOnEdge, edgesOnCell, cellsOnCell
real (kind=RKIND), dimension(:), pointer :: dvEdge, dcEdge, AreaCell
- real (kind=RKIND), dimension(:,:), pointer :: weightsOnEdge
real (kind=RKIND), dimension(:,:), pointer :: v
real (kind=RKIND), dimension(:,:), pointer :: sorted_arr
@@ -2303,7 +2309,7 @@
type (field1DReal), target :: tempFieldTarget
real(kind=RKIND), dimension(:), pointer :: hs, hs1
- real(kind=RKIND) :: hm, zh, dzmin, dzmina, dzmina_global, dzminf, sm
+ real(kind=RKIND) :: hm, hm_global, zh, dzmin, dzmina, dzmina_global, dzminf, sm
integer :: nsmterrain, kz, sfc_k
logical :: hybrid, smooth
@@ -2311,19 +2317,10 @@
real (kind=RKIND) :: p_check
! For interpolating terrain and land use
- integer :: nx, ny, nzz, iPoint, subx, suby
- integer :: isigned, endian, wordsize, istatus
- integer :: iTileStart, iTileEnd
- integer :: jTileStart, jTileEnd
- integer, allocatable, dimension(:) :: nhs
- integer, allocatable, dimension(:,:) :: ncat
- real (kind=4) :: scalefactor ! NB: this should be a single-precision real
- real (kind=RKIND) :: lat_pt, lon_pt, lon_pt_o
- real (kind=4), allocatable, dimension(:,:,:) :: rarray ! NB: this should be a single-precision real array
+ integer :: nx, ny
+ integer :: istatus
+
real (kind=RKIND), allocatable, dimension(:,:) :: rslab, maskslab
- real (kind=RKIND), allocatable, dimension(:,:) :: maxsnowalb
- real (kind=RKIND), allocatable, dimension(:,:) :: soiltemp_1deg
- real (kind=RKIND), allocatable, dimension(:,:,:) :: vegfra
integer, dimension(:), pointer :: mask_array
integer, dimension(grid % nEdges), target :: edge_mask
character (len=StrKIND) :: fname
@@ -2365,8 +2362,6 @@
parinfo => block % parinfo
dminfo => block % domain % dminfo
- weightsOnEdge => grid % weightsOnEdge % array
- nEdgesOnEdge => grid % nEdgesOnEdge % array
nEdgesOnCell => grid % nEdgesOnCell % array
edgesOnEdge => grid % edgesOnEdge % array
edgesOnCell => grid % edgesOnCell % array
@@ -2422,622 +2417,9 @@
omega_e = omega
p0 = 1.e+05
- interp_list(1) = FOUR_POINT
- interp_list(2) = SEARCH
- interp_list(3) = 0
-
-
- !
- ! Scale all distances and areas from a unit sphere to one with radius sphere_radius
- !
-
- if (config_static_interp) then
-
- grid % xCell % array = grid % xCell % array * r_earth
- grid % yCell % array = grid % yCell % array * r_earth
- grid % zCell % array = grid % zCell % array * r_earth
- grid % xVertex % array = grid % xVertex % array * r_earth
- grid % yVertex % array = grid % yVertex % array * r_earth
- grid % zVertex % array = grid % zVertex % array * r_earth
- grid % xEdge % array = grid % xEdge % array * r_earth
- grid % yEdge % array = grid % yEdge % array * r_earth
- grid % zEdge % array = grid % zEdge % array * r_earth
- grid % dvEdge % array = grid % dvEdge % array * r_earth
- grid % dcEdge % array = grid % dcEdge % array * r_earth
- grid % areaCell % array = grid % areaCell % array * r_earth**2.0
- grid % areaTriangle % array = grid % areaTriangle % array * r_earth**2.0
- grid % kiteAreasOnVertex % array = grid % kiteAreasOnVertex % array * r_earth**2.0
-
scalars(:,:,:) = 0.
- call atm_initialize_advection_rk(grid)
- call atm_initialize_deformation_weights(grid)
-
-
- !
- ! Interpolate HGT
- !
-! nx = 126
-! ny = 126
- nx = 1206
- ny = 1206
- nzz = 1
- isigned = 1
- endian = 0
- wordsize = 2
- scalefactor = 1.0
- allocate(rarray(nx,ny,nzz))
- allocate(nhs(grid % nCells))
- nhs(:) = 0
- ter(:) = 0.0
-
- do jTileStart=1,20401,ny-6
-! do jTileStart=1,961,ny-6
- jTileEnd = jTileStart + ny - 1 - 6
- do iTileStart=1,42001,nx-6
-! do iTileStart=1,2041,nx-6
- iTileEnd = iTileStart + nx - 1 - 6
- write(fname,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') trim(config_geog_data_path)//'topo_30s/',iTileStart,'-',iTileEnd,'.',jTileStart,'-',jTileEnd
-! write(fname,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') trim(config_geog_data_path)//'topo_10m/',iTileStart,'-',iTileEnd,'.',jTileStart,'-',jTileEnd
-write(0,*) trim(fname)
-
- call read_geogrid(fname, len_trim(fname), &
- rarray, &
- nx, ny, nzz, &
- isigned, endian, scalefactor, wordsize, istatus)
- call init_atm_check_read_error(istatus, fname, dminfo)
-
- iPoint = 1
- do j=4,ny-3
- do i=4,nx-3
- lat_pt = -89.99583 + (jTileStart + j - 5) * 0.0083333333
- lon_pt = -179.99583 + (iTileStart + i - 5) * 0.0083333333
-! lat_pt = -89.91667 + (jTileStart + j - 5) * 0.166667
-! lon_pt = -179.91667 + (iTileStart + i - 5) * 0.166667
- lat_pt = lat_pt * pii / 180.0
- lon_pt = lon_pt * pii / 180.0
-
- iPoint = nearest_cell(lat_pt, lon_pt, &
- iPoint, &
- grid % nCells, grid % maxEdges, grid % nEdgesOnCell % array, grid % cellsOnCell % array, &
- grid % latCell % array, grid % lonCell % array)
-
- ter(iPoint) = ter(iPoint) + rarray(i,j,1)
- nhs(iPoint) = nhs(iPoint) + 1
-
- end do
- end do
-
- end do
- end do
-
- do iCell=1, grid % nCells
- ter(iCell) = ter(iCell) / real(nhs(iCell))
- end do
-
- deallocate(rarray)
- deallocate(nhs)
-
-
- !
- ! Interpolate LU_INDEX
- !
- nx = 1200
- ny = 1200
- nzz = 1
- isigned = 1
- endian = 0
- wordsize = 1
- scalefactor = 1.0
- allocate(rarray(nx,ny,nzz))
- allocate(ncat(24,grid % nCells))
- ncat(:,:) = 0
- grid % lu_index % array(:) = 0.0
-
- do jTileStart=1,20401,ny
- jTileEnd = jTileStart + ny - 1
- do iTileStart=1,42001,nx
- iTileEnd = iTileStart + nx - 1
- write(fname,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') trim(config_geog_data_path)//'/landuse_30s/',iTileStart,'-',iTileEnd,'.',jTileStart,'-',jTileEnd
-write(0,*) trim(fname)
-
- call read_geogrid(fname, len_trim(fname), &
- rarray, &
- nx, ny, nzz, &
- isigned, endian, scalefactor, wordsize, istatus)
- call init_atm_check_read_error(istatus, fname, dminfo)
-
- iPoint = 1
- do j=1,ny
- do i=1,nx
- lat_pt = -89.99583 + (jTileStart + j - 2) * 0.0083333333
- lon_pt = -179.99583 + (iTileStart + i - 2) * 0.0083333333
- lat_pt = lat_pt * pii / 180.0
- lon_pt = lon_pt * pii / 180.0
-
- iPoint = nearest_cell(lat_pt, lon_pt, &
- iPoint, &
- grid % nCells, grid % maxEdges, grid % nEdgesOnCell % array, grid % cellsOnCell % array, &
- grid % latCell % array, grid % lonCell % array)
-
- ncat(int(rarray(i,j,1)),iPoint) = ncat(int(rarray(i,j,1)),iPoint) + 1
-
- end do
- end do
-
- end do
- end do
-
- do iCell=1, grid % nCells
- grid % lu_index % array(iCell) = 1
- do i=2,24
- if (ncat(i,iCell) > ncat(grid % lu_index % array(iCell),iCell)) then
- grid % lu_index % array(iCell) = i
- end if
- end do
- end do
-
- deallocate(rarray)
- deallocate(ncat)
-
-
- !
- ! Interpolate SOILCAT_TOP
- !
- nx = 1200
- ny = 1200
- nzz = 1
- isigned = 1
- endian = 0
- wordsize = 1
- scalefactor = 1.0
- allocate(rarray(nx,ny,nzz))
- allocate(ncat(16,grid % nCells))
- ncat(:,:) = 0
- grid % soilcat_top % array(:) = 0.0
-
- do jTileStart=1,20401,ny
- jTileEnd = jTileStart + ny - 1
- do iTileStart=1,42001,nx
- iTileEnd = iTileStart + nx - 1
- write(fname,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') trim(config_geog_data_path)//'/soiltype_top_30s/',iTileStart,'-',iTileEnd,'.',jTileStart,'-',jTileEnd
-write(0,*) trim(fname)
-
- call read_geogrid(fname, len_trim(fname), &
- rarray, &
- nx, ny, nzz, &
- isigned, endian, scalefactor, wordsize, istatus)
- call init_atm_check_read_error(istatus, fname, dminfo)
-
- iPoint = 1
- do j=1,ny
- do i=1,nx
- lat_pt = -89.99583 + (jTileStart + j - 2) * 0.0083333333
- lon_pt = -179.99583 + (iTileStart + i - 2) * 0.0083333333
- lat_pt = lat_pt * pii / 180.0
- lon_pt = lon_pt * pii / 180.0
-
- iPoint = nearest_cell(lat_pt, lon_pt, &
- iPoint, &
- grid % nCells, grid % maxEdges, grid % nEdgesOnCell % array, grid % cellsOnCell % array, &
- grid % latCell % array, grid % lonCell % array)
-
- ncat(int(rarray(i,j,1)),iPoint) = ncat(int(rarray(i,j,1)),iPoint) + 1
-
- end do
- end do
-
- end do
- end do
-
- do iCell=1, grid % nCells
- grid % soilcat_top % array(iCell) = 1
- do i=2,16
- if (ncat(i,iCell) > ncat(grid % soilcat_top % array(iCell),iCell)) then
- grid % soilcat_top % array(iCell) = i
- end if
- end do
- end do
-
- deallocate(rarray)
- deallocate(ncat)
-
-
- !
- ! Interpolate SOILCAT_BOT
- !
- nx = 1200
- ny = 1200
- nzz = 1
- isigned = 1
- endian = 0
- wordsize = 1
- scalefactor = 1.0
- allocate(rarray(nx,ny,nzz))
- allocate(ncat(16,grid % nCells))
- ncat(:,:) = 0
- grid % soilcat_bot % array(:) = 0.0
-
- do jTileStart=1,20401,ny
- jTileEnd = jTileStart + ny - 1
- do iTileStart=1,42001,nx
- iTileEnd = iTileStart + nx - 1
- write(fname,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') trim(config_geog_data_path)//'/soiltype_bot_30s/',iTileStart,'-',iTileEnd,'.',jTileStart,'-',jTileEnd
-write(0,*) trim(fname)
-
- call read_geogrid(fname, len_trim(fname), &
- rarray, &
- nx, ny, nzz, &
- isigned, endian, scalefactor, wordsize, istatus)
- call init_atm_check_read_error(istatus, fname, dminfo)
-
- iPoint = 1
- do j=1,ny
- do i=1,nx
- lat_pt = -89.99583 + (jTileStart + j - 2) * 0.0083333333
- lon_pt = -179.99583 + (iTileStart + i - 2) * 0.0083333333
- lat_pt = lat_pt * pii / 180.0
- lon_pt = lon_pt * pii / 180.0
-
- iPoint = nearest_cell(lat_pt, lon_pt, &
- iPoint, &
- grid % nCells, grid % maxEdges, grid % nEdgesOnCell % array, grid % cellsOnCell % array, &
- grid % latCell % array, grid % lonCell % array)
-
- ncat(int(rarray(i,j,1)),iPoint) = ncat(int(rarray(i,j,1)),iPoint) + 1
-
- end do
- end do
-
- end do
- end do
-
- do iCell=1, grid % nCells
- grid % soilcat_bot % array(iCell) = 1
- do i=2,16
- if (ncat(i,iCell) > ncat(grid % soilcat_bot % array(iCell),iCell)) then
- grid % soilcat_bot % array(iCell) = i
- end if
- end do
- end do
-
- deallocate(rarray)
- deallocate(ncat)
-
-
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-! KLUDGE TO FIX SOIL TYPE OVER ANTARCTICA
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- where (grid % lu_index % array == 24) grid % soilcat_top % array = 16
- where (grid % lu_index % array == 24) grid % soilcat_bot % array = 16
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-! CORRECT INCONSISTENT SOIL AND LAND USE DATA
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- do iCell = 1,grid % nCells
- if (grid % lu_index % array(iCell) == 16 .or. &
- grid % soilcat_top % array(iCell) == 14 .or. &
- grid % soilcat_bot % array(iCell) == 14) then
- if (grid % lu_index % array(iCell) /= 16) then
- write(0,*) 'Turning lu_index into water at ', iCell
- grid % lu_index % array(iCell) = 16
- end if
- if (grid % soilcat_top % array(iCell) /= 14) then
- write(0,*) 'Turning soilcat_top into water at ', iCell
- grid % soilcat_top % array(iCell) = 14
- end if
- if (grid % soilcat_bot % array(iCell) /= 14) then
- write(0,*) 'Turning soilcat_bot into water at ', iCell
- grid % soilcat_bot % array(iCell) = 14
- end if
- end if
- end do
-
-
- !
- ! Derive LANDMASK
- !
- grid % landmask % array(:) = 0
- do iCell=1, grid % nCells
- if (grid % lu_index % array(iCell) /= 16) grid % landmask % array(iCell) = 1
- end do
-
-
- !
- ! Interpolate SOILTEMP:
- !
- nx = 186
- ny = 186
- nzz = 1
- isigned = 0
- endian = 0
- wordsize = 2
- scalefactor = 0.01
- allocate(rarray(nx,ny,nzz))
- allocate(soiltemp_1deg(-2:363,-2:183))
- grid % soiltemp % array(:) = 0.0
-
- call map_set(PROJ_LATLON, proj, &
- latinc = 1.0_RKIND, &
- loninc = 1.0_RKIND, &
- knowni = 1.0_RKIND, &
- knownj = 1.0_RKIND, &
- lat1 = -89.5_RKIND, &
- lon1 = -179.5_RKIND)
-
- write(fname,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') trim(config_geog_data_path)//'soiltemp_1deg/',1,'-',180,'.',1,'-',180
-write(0,*) trim(fname)
- call read_geogrid(fname, len_trim(fname), &
- rarray, &
- nx, ny, nzz, &
- isigned, endian, scalefactor, wordsize, istatus)
- call init_atm_check_read_error(istatus, fname, dminfo)
-
- soiltemp_1deg(-2:180,-2:183) = rarray(1:183,1:186,1)
-
- write(fname,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') trim(config_geog_data_path)//'soiltemp_1deg/',181,'-',360,'.',1,'-',180
-write(0,*) trim(fname)
- call read_geogrid(fname, len_trim(fname), &
- rarray, &
- nx, ny, nzz, &
- isigned, endian, scalefactor, wordsize, istatus)
- call init_atm_check_read_error(istatus, fname, dminfo)
-
- soiltemp_1deg(181:363,-2:183) = rarray(4:186,1:186,1)
-
- interp_list(1) = FOUR_POINT
- interp_list(2) = W_AVERAGE4
- interp_list(3) = W_AVERAGE16
- interp_list(4) = SEARCH
- interp_list(5) = 0
-
- do iCell=1,grid%nCells
-
- if (grid % landmask % array(iCell) == 1) then
- lat = grid % latCell % array(iCell)*DEG_PER_RAD
- lon = grid % lonCell % array(iCell)*DEG_PER_RAD
- call latlon_to_ij(proj, lat, lon, x, y)
- if (x < 0.5) then
- lon = lon + 360.0
- call latlon_to_ij(proj, lat, lon, x, y)
- else if (x >= 360.5) then
- lon = lon - 360.0
- call latlon_to_ij(proj, lat, lon, x, y)
- end if
-if (y < 1.0) y = 1.0
-if (y > 179.0) y = 179.0
-! grid % soiltemp % array(iCell) = interp_sequence(x, y, 1, soiltemp_1deg, 1, 360, 1, 180, 1, 1, -1.e30_RKIND, interp_list, 1)
- grid % soiltemp % array(iCell) = interp_sequence(x, y, 1, soiltemp_1deg, -2, 363, -2, 183, 1, 1, 0.0_RKIND, interp_list, 1)
- else
- grid % soiltemp % array(iCell) = 0.0
- end if
-
- end do
-
- deallocate(rarray)
- deallocate(soiltemp_1deg)
-
-
- !
- ! Interpolate SNOALB
- !
- nx = 186
- ny = 186
- nzz = 1
- isigned = 0
- endian = 0
- wordsize = 1
- scalefactor = 1.0
- allocate(rarray(nx,ny,nzz))
- allocate(maxsnowalb(-2:363,-2:183))
- grid % snoalb % array(:) = 0.0
-
- call map_set(PROJ_LATLON, proj, &
- latinc = 1.0_RKIND, &
- loninc = 1.0_RKIND, &
- knowni = 1.0_RKIND, &
- knownj = 1.0_RKIND, &
- lat1 = -89.5_RKIND, &
- lon1 = -179.5_RKIND)
-
- write(fname,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') trim(config_geog_data_path)//'maxsnowalb/',1,'-',180,'.',1,'-',180
-write(0,*) trim(fname)
- call read_geogrid(fname, len_trim(fname), &
- rarray, &
- nx, ny, nzz, &
- isigned, endian, scalefactor, wordsize, istatus)
- call init_atm_check_read_error(istatus, fname, dminfo)
-
- maxsnowalb(-2:180,-2:183) = rarray(1:183,1:186,1)
-
- write(fname,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') trim(config_geog_data_path)//'maxsnowalb/',181,'-',360,'.',1,'-',180
-write(0,*) trim(fname)
- call read_geogrid(fname, len_trim(fname), &
- rarray, &
- nx, ny, nzz, &
- isigned, endian, scalefactor, wordsize, istatus)
- call init_atm_check_read_error(istatus, fname, dminfo)
-
- maxsnowalb(181:363,-2:183) = rarray(4:186,1:186,1)
-
- interp_list(1) = FOUR_POINT
- interp_list(2) = W_AVERAGE4
- interp_list(3) = W_AVERAGE16
- interp_list(4) = SEARCH
- interp_list(5) = 0
-
- do iCell=1,grid%nCells
-
- if (grid % landmask % array(iCell) == 1) then
- lat = grid % latCell % array(iCell)*DEG_PER_RAD
- lon = grid % lonCell % array(iCell)*DEG_PER_RAD
- call latlon_to_ij(proj, lat, lon, x, y)
- if (x < 0.5) then
- lon = lon + 360.0
- call latlon_to_ij(proj, lat, lon, x, y)
- else if (x >= 360.5) then
- lon = lon - 360.0
- call latlon_to_ij(proj, lat, lon, x, y)
- end if
-if (y < 1.0) y = 1.0
-if (y > 179.0) y = 179.0
-! grid % snoalb % array(iCell) = interp_sequence(x, y, 1, maxsnowalb, 1, 360, 1, 180, 1, 1, -1.e30_RKIND, interp_list, 1)
- grid % snoalb % array(iCell) = interp_sequence(x, y, 1, maxsnowalb, -2, 363, -2, 183, 1, 1, 0.0_RKIND, interp_list, 1)
- else
- grid % snoalb % array(iCell) = 0.0
- end if
-
- end do
-
- grid % snoalb % array(:) = grid % snoalb % array(:) / 100.0
-
- deallocate(rarray)
- deallocate(maxsnowalb)
-
-
- !
- ! Interpolate GREENFRAC
- !
- nx = 1256
- ny = 1256
- nzz = 12
- isigned = 0
- endian = 0
- wordsize = 1
- scalefactor = 1.0
- allocate(rarray(nx,ny,nzz))
- allocate(vegfra(-2:2503,-2:1253,12))
-! grid % vegfra % array(:) = 0.0
- grid % greenfrac % array(:,:) = 0.0
-
- call map_set(PROJ_LATLON, proj, &
- latinc = 0.144_RKIND, &
- loninc = 0.144_RKIND, &
- knowni = 1.0_RKIND, &
- knownj = 1.0_RKIND, &
- lat1 = -89.928_RKIND, &
- lon1 = -179.928_RKIND)
-
- write(fname,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') trim(config_geog_data_path)//'greenfrac/',1,'-',1250,'.',1,'-',1250
-write(0,*) trim(fname)
- call read_geogrid(fname, len_trim(fname), &
- rarray, &
- nx, ny, nzz, &
- isigned, endian, scalefactor, wordsize, istatus)
- call init_atm_check_read_error(istatus, fname, dminfo)
-
- vegfra(-2:1250,-2:1253,1:12) = rarray(1:1253,1:1256,1:12)
-
- write(fname,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') trim(config_geog_data_path)//'greenfrac/',1251,'-',2500,'.',1,'-',1250
-write(0,*) trim(fname)
- call read_geogrid(fname, len_trim(fname), &
- rarray, &
- nx, ny, nzz, &
- isigned, endian, scalefactor, wordsize, istatus)
- call init_atm_check_read_error(istatus, fname, dminfo)
-
- vegfra(1251:2503,-2:1253,1:12) = rarray(4:1256,1:1256,1:12)
-
- do iCell=1,grid%nCells
- if (grid % landmask % array(iCell) == 1) then
- lat = grid % latCell % array(iCell)*DEG_PER_RAD
- lon = grid % lonCell % array(iCell)*DEG_PER_RAD
- call latlon_to_ij(proj, lat, lon, x, y)
- if (x < 0.5) then
- lon = lon + 360.0
- call latlon_to_ij(proj, lat, lon, x, y)
- else if (x >= 2500.5) then
- lon = lon - 360.0
- call latlon_to_ij(proj, lat, lon, x, y)
- end if
-if (y < 1.0) y = 1.0
-if (y > 1249.0) y = 1249.0
- do k=1,12
- grid % greenfrac % array(k,iCell) = interp_sequence(x, y, k, vegfra, -2, 2503, -2, 1253, 1, 12, -1.e30_RKIND, interp_list, 1)
- end do
- else
- grid % greenfrac % array(:,iCell) = 0.0
- end if
- grid % shdmin % array(iCell) = minval(grid % greenfrac % array(:,iCell))
- grid % shdmax % array(iCell) = maxval(grid % greenfrac % array(:,iCell))
-
- end do
-
- deallocate(rarray)
- deallocate(vegfra)
-
-
- !
- ! Interpolate ALBEDO12M
- !
- nx = 1256
- ny = 1256
- nzz = 12
- isigned = 0
- endian = 0
- wordsize = 1
- scalefactor = 1.0
- allocate(rarray(nx,ny,nzz))
- allocate(vegfra(-2:2503,-2:1253,12))
- grid % albedo12m % array(:,:) = 0.0
-
- call map_set(PROJ_LATLON, proj, &
- latinc = 0.144_RKIND, &
- loninc = 0.144_RKIND, &
- knowni = 1.0_RKIND, &
- knownj = 1.0_RKIND, &
- lat1 = -89.928_RKIND, &
- lon1 = -179.928_RKIND)
-
- write(fname,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') trim(config_geog_data_path)//'albedo_ncep/',1,'-',1250,'.',1,'-',1250
-write(0,*) trim(fname)
- call read_geogrid(fname, len_trim(fname), &
- rarray, &
- nx, ny, nzz, &
- isigned, endian, scalefactor, wordsize, istatus)
- call init_atm_check_read_error(istatus, fname, dminfo)
-
- vegfra(-2:1250,-2:1253,1:12) = rarray(1:1253,1:1256,1:12)
-
- write(fname,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') trim(config_geog_data_path)//'albedo_ncep/',1251,'-',2500,'.',1,'-',1250
-write(0,*) trim(fname)
- call read_geogrid(fname, len_trim(fname), &
- rarray, &
- nx, ny, nzz, &
- isigned, endian, scalefactor, wordsize, istatus)
- call init_atm_check_read_error(istatus, fname, dminfo)
-
- vegfra(1251:2503,-2:1253,1:12) = rarray(4:1256,1:1256,1:12)
-
- do iCell=1,grid%nCells
- if (grid % landmask % array(iCell) == 1) then
- lat = grid % latCell % array(iCell)*DEG_PER_RAD
- lon = grid % lonCell % array(iCell)*DEG_PER_RAD
- call latlon_to_ij(proj, lat, lon, x, y)
- if (x < 0.5) then
- lon = lon + 360.0
- call latlon_to_ij(proj, lat, lon, x, y)
- else if (x >= 2500.5) then
- lon = lon - 360.0
- call latlon_to_ij(proj, lat, lon, x, y)
- end if
-if (y < 1.0) y = 1.0
-if (y > 1249.0) y = 1249.0
- do k=1,12
- grid % albedo12m % array(k,iCell) = interp_sequence(x, y, k, vegfra, -2, 2503, -2, 1253, 1, 12, 0.0_RKIND, interp_list, 1)
- end do
- else
- grid % albedo12m % array(:,iCell) = 8.0
- end if
- end do
-
- deallocate(rarray)
- deallocate(vegfra)
-
-
- end if ! config_static_interp
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! BEGIN ADOPT GFS TERRAIN HEIGHT
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -3153,7 +2535,9 @@
hx(:,iCell) = ter(iCell)
end do
- hm = maxval(ter(:))
+ hm = maxval(ter(1:nCellsSolve))
+ call mpas_dmpar_max_real(dminfo, hm, hm_global)
+ hm = hm_global
write(0,*) "max ter = ", hm
! Metrics for hybrid coordinate and vertical stretching
@@ -3293,7 +2677,7 @@
call mpas_dmpar_exch_halo_field(tempField)
! dzmina = minval(hs(:)-hx(k-1,:))
- dzmina = minval(zw(k)+ah(k)*hs(1:grid%nCellsSolve)-zw(k-1)-ah(k-1)*hx(k-1,1:grid%nCellsSolve))
+ dzmina = minval(zw(k)+ah(k)*hs(1:nCellsSolve)-zw(k-1)-ah(k-1)*hx(k-1,1:nCellsSolve))
call mpas_dmpar_min_real(dminfo, dzmina, dzmina_global)
! write(0,*) ' k,i, dzmina, dzmin, zw(k)-zw(k-1) ', k,i, dzmina, dzmin, zw(k)-zw(k-1)
if (dzmina_global >= dzmin*(zw(k)-zw(k-1))) then
@@ -4130,6 +3514,7 @@
do k=1,grid%nVertLevels
target_z = 0.5 * (grid % zgrid % array(k,iCell) + grid % zgrid % array(k+1,iCell))
state % scalars % array(state % index_qv,k,iCell) = vertical_interp(target_z, config_nfglevels, sorted_arr, order=1, extrap=0)
+ diag % rh % array(k,iCell) = state % scalars % array(state % index_qv,k,iCell)
end do
@@ -4256,8 +3641,8 @@
! QV
es = 6.112 * exp((17.27*(state % theta_m % array(k,iCell) - 273.16))/(state % theta_m % array(k,iCell) - 35.86))
- rs = 0.622 * es / (diag % pressure % array(k,iCell) - es)
- scalars(state % index_qv,k,iCell) = rs * scalars(state % index_qv,k,iCell)
+ rs = 0.622 * es * 100. / (diag % pressure % array(k,iCell) - es * 100.)
+ scalars(state % index_qv,k,iCell) = 0.01 * rs * scalars(state % index_qv,k,iCell)
! PI
p(k,iCell) = (diag % pressure % array(k,iCell) / p0) ** (rgas / cp)
@@ -4274,6 +3659,17 @@
!
+ ! Calculation of the initial precipitable water:
+ !
+ do iCell = 1,grid%nCells
+ diag_physics%precipw%array(iCell) = 0.0
+ do k = 1,grid%nVertLevels
+ diag_physics%precipw%array(iCell) = diag_physics%precipw%array(iCell) &
+ + rho_zz(k,iCell)*scalars(state%index_qv,k,iCell)*(zgrid(k+1,iCell)-zgrid(k,iCell))
+ enddo
+ enddo
+
+ !
! Reference state based on a dry isothermal atmosphere
!
do iCell=1,grid % nCells
@@ -5977,48 +5373,6 @@
end subroutine init_atm_test_case_resting_atmosphere
- integer function nearest_cell(target_lat, target_lon, &
- start_cell, &
- nCells, maxEdges, nEdgesOnCell, cellsOnCell, latCell, lonCell)
-
- implicit none
-
- real (kind=RKIND), intent(in) :: target_lat, target_lon
- integer, intent(in) :: start_cell
- integer, intent(in) :: nCells, maxEdges
- integer, dimension(nCells), intent(in) :: nEdgesOnCell
- integer, dimension(maxEdges,nCells), intent(in) :: cellsOnCell
- real (kind=RKIND), dimension(nCells), intent(in) :: latCell, lonCell
-
- integer :: i
- integer :: iCell
- integer :: current_cell
- real (kind=RKIND) :: current_distance, d
- real (kind=RKIND) :: nearest_distance
-
- nearest_cell = start_cell
- current_cell = -1
-
- do while (nearest_cell /= current_cell)
- current_cell = nearest_cell
- current_distance = sphere_distance(latCell(current_cell), lonCell(current_cell), target_lat, target_lon, 1.0_RKIND)
- nearest_cell = current_cell
- nearest_distance = current_distance
- do i = 1, nEdgesOnCell(current_cell)
- iCell = cellsOnCell(i,current_cell)
- if (iCell <= nCells) then
- d = sphere_distance(latCell(iCell), lonCell(iCell), target_lat, target_lon, 1.0_RKIND)
- if (d < nearest_distance) then
- nearest_cell = iCell
- nearest_distance = d
- end if
- end if
- end do
- end do
-
- end function nearest_cell
-
-
integer function nearest_edge(target_lat, target_lon, &
start_edge, &
nCells, nEdges, maxEdges, nEdgesOnCell, edgesOnCell, cellsOnEdge, latCell, lonCell, latEdge, lonEdge)
@@ -6160,44 +5514,8 @@
end function vertical_interp
- subroutine init_atm_check_read_error(istatus, fname, dminfo)
-
- implicit none
-
- integer, intent(in) :: istatus
- character (len=*), intent(in) :: fname
- type (dm_info), intent(in) :: dminfo
-
- if (istatus /= 0) then
- write(0,*) 'ERROR: Could not read file '//trim(fname)
- call mpas_dmpar_abort(dminfo)
- end if
-
- end subroutine init_atm_check_read_error
-
-
!----------------------------------------------------------------------------------------------------------
- real (kind=RKIND) function sphere_distance(lat1, lon1, lat2, lon2, radius)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Compute the great-circle distance between (lat1, lon1) and (lat2, lon2) on a
- ! sphere with given radius.
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- implicit none
-
- real (kind=RKIND), intent(in) :: lat1, lon1, lat2, lon2, radius
-
- real (kind=RKIND) :: arg1
-
- arg1 = sqrt( sin(0.5*(lat2-lat1))**2 + &
- cos(lat1)*cos(lat2)*sin(0.5*(lon2-lon1))**2 )
- sphere_distance = 2.*radius*asin(arg1)
-
- end function sphere_distance
-
-!--------------------------------------------------------------------
-
real (kind=RKIND) function env_qv( z, temperature, pressure, rh_max )
implicit none
Modified: branches/atmphys_gfs_mmm/src/core_nhyd_atmos/Registry
===================================================================
--- branches/atmphys_gfs_mmm/src/core_nhyd_atmos/Registry 2013-02-28 22:21:22 UTC (rev 2525)
+++ branches/atmphys_gfs_mmm/src/core_nhyd_atmos/Registry 2013-02-28 22:35:04 UTC (rev 2526)
@@ -115,7 +115,7 @@
var persistent real edgeNormalVectors ( R3 nEdges ) 0 iro edgeNormalVectors mesh - -
var persistent real localVerticalUnitVectors ( R3 nCells ) 0 iro localVerticalUnitVectors mesh - -
-var persistent real cellTangentPlane ( R3 TWO nEdges ) 0 iro cellTangentPlane mesh - -
+var persistent real cellTangentPlane ( R3 TWO nCells ) 0 iro cellTangentPlane mesh - -
var persistent integer cellsOnCell ( maxEdges nCells ) 0 iro cellsOnCell mesh - -
var persistent integer verticesOnCell ( maxEdges nCells ) 0 iro verticesOnCell mesh - -
@@ -381,6 +381,7 @@
namelist character physics config_eddy_scheme off
namelist character physics config_lsm_scheme off
namelist character physics config_pbl_scheme off
+namelist character physics config_gwdo_scheme off
namelist character physics config_radt_cld_scheme off
namelist character physics config_radt_lw_scheme off
namelist character physics config_radt_sw_scheme off
@@ -520,6 +521,41 @@
var persistent real kzq ( nVertLevels nCells Time ) 1 o kzq diag_physics - -
%--------------------------------------------------------------------------------------------------
+%... PARAMETERIZATION OF GRAVITY WAVE DRAG OVER OROGRAPHY:
+%--------------------------------------------------------------------------------------------------
+
+% var2d : orographic variance (m2)
+% con : orographic convexity (m2)
+% oa1 : orographic direction asymmetry function (-)
+% oa2 : orographic direction asymmetry function (-)
+% oa3 : orographic direction asymmetry function (-)
+% oa4 : orographic direction asymmetry function (-)
+% ol1 : orographic direction asymmetry function (-)
+% ol2 : orographic direction asymmetry function (-)
+% ol3 : orographic direction asymmetry function (-)
+% ol4 : orographic direction asymmetry function (-)
+% dusfcg : vertically-integrated gravity wave drag over orography u-stress (Pa m s-1)
+% dvsfcg : vertically-integrated gravity wave drag over orography v-stress (Pa m s-1)
+% dtaux3d : gravity wave drag over orography u-stress (m s-1)
+% dtauy3d : gravity wave drag over orography v-stress (m s-1)
+
+var persistent real var2d ( nCells ) 0 ro var2d sfc_input - -
+var persistent real con ( nCells ) 0 ro con sfc_input - -
+var persistent real oa1 ( nCells ) 0 ro oa1 sfc_input - -
+var persistent real oa2 ( nCells ) 0 ro oa2 sfc_input - -
+var persistent real oa3 ( nCells ) 0 ro oa3 sfc_input - -
+var persistent real oa4 ( nCells ) 0 ro oa4 sfc_input - -
+var persistent real ol1 ( nCells ) 0 ro ol1 sfc_input - -
+var persistent real ol2 ( nCells ) 0 ro ol2 sfc_input - -
+var persistent real ol3 ( nCells ) 0 ro ol3 sfc_input - -
+var persistent real ol4 ( nCells ) 0 ro ol4 sfc_input - -
+
+var persistent real dusfcg ( nCells Time ) 1 ro dusfcg diag_physics - -
+var persistent real dvsfcg ( nCells Time ) 1 ro dvsfcg diag_physics - -
+var persistent real dtaux3d ( nVertLevels nCells Time ) 1 ro dtaux3d diag_physics - -
+var persistent real dtauy3d ( nVertLevels nCells Time ) 1 ro dtauy3d diag_physics - -
+
+%--------------------------------------------------------------------------------------------------
%... PARAMETERIZATION OF SURFACE LAYER PROCESSES:
%--------------------------------------------------------------------------------------------------
Modified: branches/atmphys_gfs_mmm/src/core_nhyd_atmos/mpas_atm_interp_diagnostics.F
===================================================================
--- branches/atmphys_gfs_mmm/src/core_nhyd_atmos/mpas_atm_interp_diagnostics.F 2013-02-28 22:21:22 UTC (rev 2525)
+++ branches/atmphys_gfs_mmm/src/core_nhyd_atmos/mpas_atm_interp_diagnostics.F 2013-02-28 22:35:04 UTC (rev 2526)
@@ -187,7 +187,7 @@
do iCell = 1, nCells
do k = 1, nVertLevels
kk = nVertLevels+1-k
- field_in(iCell,kk) = uzonal(k,iCell)
+ field_in(iCell,kk) = umeridional(k,iCell)
enddo
enddo
call interp_tofixed_pressure(nCells,nVertLevels,nIntP,press_in,field_in,press_interp,field_interp)
Index: branches/atmphys_gfs_mmm/src/core_ocean
===================================================================
--- branches/atmphys_gfs_mmm/src/core_ocean 2013-02-28 22:21:22 UTC (rev 2525)
+++ branches/atmphys_gfs_mmm/src/core_ocean 2013-02-28 22:35:04 UTC (rev 2526)
Property changes on: branches/atmphys_gfs_mmm/src/core_ocean
___________________________________________________________________
Modified: svn:mergeinfo
## -1,3 +1,4 ##
+/branches/atmos_physics/src/core_ocean:2282-2525
/branches/cam_mpas_nh/src/core_ocean:1260-1270
/branches/ocean_projects/ale_split_exp/src/core_ocean:1437-1483
/branches/ocean_projects/ale_vert_coord/src/core_ocean:1225-1383
\ No newline at end of property
Modified: branches/atmphys_gfs_mmm/src/framework/mpas_io.F
===================================================================
--- branches/atmphys_gfs_mmm/src/framework/mpas_io.F 2013-02-28 22:21:22 UTC (rev 2525)
+++ branches/atmphys_gfs_mmm/src/framework/mpas_io.F 2013-02-28 22:35:04 UTC (rev 2526)
@@ -913,9 +913,12 @@
type (fieldlist_type), pointer :: field_cursor
integer :: pio_type
- integer :: ndims, pd
- integer :: i, i1, i2, i3, i4, i5, indx
- integer, dimension(:), pointer :: dimlist, compdof
+ integer :: ndims
+ integer (kind=PIO_OFFSET) :: pd, indx
+ integer :: i
+ integer (kind=PIO_OFFSET) :: i1, i2, i3, i4, i5
+ integer, dimension(:), pointer :: dimlist
+ integer (kind=PIO_OFFSET), dimension(:), pointer :: compdof
type (decomplist_type), pointer :: decomp_cursor, new_decomp
! write(0,*) 'Called MPAS_io_set_var_indices()'
@@ -1053,11 +1056,11 @@
do i=1,ndims-1
dimlist(i) = field_cursor % fieldhandle % dims(i) % dimsize
new_decomp % decomphandle % dims(i) = dimlist(i)
- pd = pd * dimlist(i)
+ pd = pd * int(dimlist(i),PIO_OFFSET)
end do
new_decomp % decomphandle % dims(ndims) = field_cursor % fieldhandle % dims(ndims) % dimsize
dimlist(ndims) = size(indices)
- pd = pd * dimlist(ndims)
+ pd = pd * int(dimlist(ndims),PIO_OFFSET)
allocate(compdof(pd))
@@ -1069,10 +1072,10 @@
do i2=1,dimlist(2)
do i1=1,dimlist(1)
compdof(indx) = i1 &
- + (i2-1)*dimlist(1) &
- + (i3-1)*dimlist(2)*dimlist(1) &
- + (i4-1)*dimlist(3)*dimlist(2)*dimlist(1) &
- + (indices(i5)-1)*dimlist(4)*dimlist(3)*dimlist(2)*dimlist(1)
+ + (i2-1)*int(dimlist(1),PIO_OFFSET) &
+ + (i3-1)*int(dimlist(2),PIO_OFFSET)*int(dimlist(1),PIO_OFFSET) &
+ + (i4-1)*int(dimlist(3),PIO_OFFSET)*int(dimlist(2),PIO_OFFSET)*int(dimlist(1),PIO_OFFSET) &
+ + int(indices(i5)-1,PIO_OFFSET)*int(dimlist(4),PIO_OFFSET)*int(dimlist(3),PIO_OFFSET)*int(dimlist(2),PIO_OFFSET)*int(dimlist(1),PIO_OFFSET)
indx = indx + 1
end do
end do
@@ -1085,9 +1088,9 @@
do i2=1,dimlist(2)
do i1=1,dimlist(1)
compdof(indx) = i1 &
- + (i2-1)*dimlist(1) &
- + (i3-1)*dimlist(2)*dimlist(1) &
- + (indices(i4)-1)*dimlist(3)*dimlist(2)*dimlist(1)
+ + (i2-1)*int(dimlist(1),PIO_OFFSET) &
+ + (i3-1)*int(dimlist(2),PIO_OFFSET)*int(dimlist(1),PIO_OFFSET) &
+ + int(indices(i4)-1,PIO_OFFSET)*int(dimlist(3),PIO_OFFSET)*int(dimlist(2),PIO_OFFSET)*int(dimlist(1),PIO_OFFSET)
indx = indx + 1
end do
end do
@@ -1097,7 +1100,7 @@
do i3=1,dimlist(3)
do i2=1,dimlist(2)
do i1=1,dimlist(1)
- compdof(indx) = i1 + (i2-1)*dimlist(1) + (indices(i3)-1)*dimlist(2)*dimlist(1)
+ compdof(indx) = i1 + (i2-1)*int(dimlist(1),PIO_OFFSET) + int(indices(i3)-1,PIO_OFFSET)*int(dimlist(2),PIO_OFFSET)*int(dimlist(1),PIO_OFFSET)
indx = indx + 1
end do
end do
@@ -1105,13 +1108,13 @@
else if (ndims == 2) then
do i2=1,dimlist(2)
do i1=1,dimlist(1)
- compdof(indx) = i1 + (indices(i2)-1)*dimlist(1)
+ compdof(indx) = i1 + int(indices(i2)-1,PIO_OFFSET)*int(dimlist(1),PIO_OFFSET)
indx = indx + 1
end do
end do
else if (ndims == 1) then
do i1=1,dimlist(1)
- compdof(indx) = indices(i1)
+ compdof(indx) = int(indices(i1),PIO_OFFSET)
indx = indx + 1
end do
end if
Modified: branches/atmphys_gfs_mmm/src/framework/mpas_io_input.F
===================================================================
--- branches/atmphys_gfs_mmm/src/framework/mpas_io_input.F 2013-02-28 22:21:22 UTC (rev 2525)
+++ branches/atmphys_gfs_mmm/src/framework/mpas_io_input.F 2013-02-28 22:35:04 UTC (rev 2526)
@@ -516,7 +516,7 @@
integer :: ierr
- call MPAS_readStream(input_obj % io_stream, 1, ierr)
+ call MPAS_readStream(input_obj % io_stream, input_obj % time, ierr)
end subroutine mpas_read_and_distribute_fields!}}}
</font>
</pre>