<p><b>duda</b> 2012-01-13 17:47:18 -0700 (Fri, 13 Jan 2012)</p><p>Merging atmos_physics branch to trunk.<br>
<br>
<br>
M src/core_init_nhyd_atmos/mpas_init_atm_test_cases.F<br>
M src/core_init_nhyd_atmos/Registry<br>
M src/core_atmos_physics/mpas_atmphys_date_time.F<br>
M src/core_atmos_physics/mpas_atmphys_driver_radiation_sw.F<br>
M src/core_atmos_physics/mpas_atmphys_todynamics.F<br>
M src/core_atmos_physics/mpas_atmphys_camrad_init.F<br>
M src/core_atmos_physics/mpas_atmphys_vars.F<br>
M src/core_atmos_physics/mpas_atmphys_interface_nhyd.F<br>
M src/core_atmos_physics/mpas_atmphys_lsm_noahinit.F<br>
M src/core_atmos_physics/mpas_atmphys_driver_radiation_lw.F<br>
M src/core_atmos_physics/mpas_atmphys_driver.F<br>
M src/core_atmos_physics/mpas_atmphys_initialize_real.F<br>
M src/core_nhyd_atmos/mpas_atm_time_integration.F<br>
M src/core_nhyd_atmos/Registry<br>
</p><hr noshade><pre><font color="gray">Modified: trunk/mpas/src/core_atmos_physics/mpas_atmphys_camrad_init.F
===================================================================
--- trunk/mpas/src/core_atmos_physics/mpas_atmphys_camrad_init.F 2012-01-13 22:13:27 UTC (rev 1368)
+++ trunk/mpas/src/core_atmos_physics/mpas_atmphys_camrad_init.F 2012-01-14 00:47:18 UTC (rev 1369)
@@ -327,36 +327,6 @@
DM_BCAST_MACRO(ln_ah2ow)
DM_BCAST_MACRO(ln_eh2ow)
-!101 format(i4,7(1x,e19.12))
-!write(0,*) 'ah2onw'
-!do i_te = 1,21
-! write(0,101) i_te,(ah2onw(1,1,1,i_te,i_rh),i_rh=1,n_rh)
-!enddo
-!write(0,*) 'eh2onw'
-!do i_te = 1,21
-! write(0,101) i_te,(eh2onw(1,1,1,i_te,i_rh),i_rh=1,n_rh)
-!enddo
-!write(0,*) 'ah2ow'
-!do i_te = 1,21
-! write(0,101) i_te,(ah2ow(1,1,1,i_te,i_rh),i_rh=1,n_rh)
-!enddo
-!write(0,*) 'cn_ah2ow'
-!do i_te = 1,21
-! write(0,101) i_te,(cn_ah2ow(1,1,1,i_te,i_rh),i_rh=1,n_rh)
-!enddo
-!write(0,*) 'cn_eh2ow'
-!do i_te = 1,21
-! write(0,101) i_te,(cn_eh2ow(1,1,1,i_te,i_rh),i_rh=1,n_rh)
-!enddo
-!write(0,*) 'ln_ah2ow'
-!do i_te = 1,21
-! write(0,101) i_te,(ln_ah2ow(1,1,1,i_te,i_rh),i_rh=1,n_rh)
-!enddo
-!write(0,*) 'ln_eh2ow'
-!do i_te = 1,21
-! write(0,101) i_te,(ln_eh2ow(1,1,1,i_te,i_rh),i_rh=1,n_rh)
-!enddo
-
if(dminfo % my_proc_id == IO_NODE) close(cam_abs_unit)
! Set up table of H2O saturation vapor pressures for use in calculation effective path RH.
@@ -612,14 +582,15 @@
DM_BCAST_MACRO(kcb)
DM_BCAST_MACRO(wcb)
DM_BCAST_MACRO(gcb)
- DM_BCAST_MACRO(kvolc)
- DM_BCAST_MACRO(wvolc)
DM_BCAST_MACRO(kdst)
DM_BCAST_MACRO(wdst)
DM_BCAST_MACRO(gdst)
DM_BCAST_MACRO(kbg)
DM_BCAST_MACRO(wbg)
DM_BCAST_MACRO(gbg)
+ DM_BCAST_MACRO(kvolc)
+ DM_BCAST_MACRO(wvolc)
+ DM_BCAST_MACRO(gvolc)
if(dminfo % my_proc_id == IO_NODE) close(cam_aer_unit)
@@ -676,7 +647,7 @@
enddo
enddo
-!write(0,*) '--- end subroutine aer_optics_initialize:'
+ write(0,*) '--- end subroutine aer_optics_initialize:'
end subroutine aer_optics_initialize
Modified: trunk/mpas/src/core_atmos_physics/mpas_atmphys_date_time.F
===================================================================
--- trunk/mpas/src/core_atmos_physics/mpas_atmphys_date_time.F 2012-01-13 22:13:27 UTC (rev 1368)
+++ trunk/mpas/src/core_atmos_physics/mpas_atmphys_date_time.F 2012-01-14 00:47:18 UTC (rev 1369)
@@ -98,7 +98,7 @@
write(0,*)
write(0,*) '--- enter subroutine monthly_interp_to_date:'
- write(0,*) '--- current_date =',date_str
+ write(0,*) '--- current_date = ',date_str
write(day15,fmt='(I2.2)') 15
do l = 1 , 12
@@ -131,8 +131,8 @@
month2 = month1 + 1
endif
if(n == 1) then
- write(0,*) '--- month1=',month1
- write(0,*) '--- month2=',month2
+ write(0,*) '--- month1 =',month1
+ write(0,*) '--- month2 =',month2
endif
field_out(n) = ( field_in(month2,n) * (target_date - middle(l)) &
+ field_in(month1,n) * (middle(l+1) - target_date)) &
Modified: trunk/mpas/src/core_atmos_physics/mpas_atmphys_driver.F
===================================================================
--- trunk/mpas/src/core_atmos_physics/mpas_atmphys_driver.F 2012-01-13 22:13:27 UTC (rev 1368)
+++ trunk/mpas/src/core_atmos_physics/mpas_atmphys_driver.F 2012-01-14 00:47:18 UTC (rev 1369)
@@ -71,7 +71,7 @@
!call to short wave radiation scheme:
if(l_radtsw) then
- call allocate_radiation_sw
+ call allocate_radiation_sw(xtime_s)
call driver_radiation_sw(itimestep,block%mesh,block%state%time_levs(1)%state, &
block%diag_physics,block%sfc_input,block%tend_physics, &
xtime_s)
@@ -79,7 +79,7 @@
!call to long wave radiation scheme:
if(l_radtlw) then
- call allocate_radiation_lw
+ call allocate_radiation_lw(xtime_s)
call driver_radiation_lw(itimestep,block%mesh,block%state%time_levs(1)%state, &
block%diag_physics,block%sfc_input,block%tend_physics, &
xtime_s)
@@ -89,7 +89,7 @@
if(config_radt_sw_scheme.ne.'off' .or. config_radt_lw_scheme.ne.'off') &
call deallocate_cloudiness
if(config_radt_sw_scheme.ne.'off') call deallocate_radiation_sw
- if(config_radt_lw_scheme.ne.'off') call deallocate_radiation_lw
+ if(config_radt_lw_scheme.ne.'off') call deallocate_radiation_lw(xtime_s)
!call to surface-layer scheme:
if(config_sfclayer_scheme .ne. 'off') then
Modified: trunk/mpas/src/core_atmos_physics/mpas_atmphys_driver_radiation_lw.F
===================================================================
--- trunk/mpas/src/core_atmos_physics/mpas_atmphys_driver_radiation_lw.F 2012-01-13 22:13:27 UTC (rev 1368)
+++ trunk/mpas/src/core_atmos_physics/mpas_atmphys_driver_radiation_lw.F 2012-01-14 00:47:18 UTC (rev 1369)
@@ -1,5 +1,6 @@
!=============================================================================================
module mpas_atmphys_driver_radiation_lw
+ use mpas_configure, only: config_do_restart
use mpas_grid_types
use mpas_timer
@@ -26,12 +27,17 @@
contains
!=============================================================================================
- subroutine allocate_radiation_lw
+ subroutine allocate_radiation_lw(xtime_s)
!=============================================================================================
- if(.not.allocated(f_ice) ) allocate(f_ice(ims:ime,kms:kme,jms:jme) )
- if(.not.allocated(f_rain) ) allocate(f_rain(ims:ime,kms:kme,jms:jme) )
+!input arguments:
+ real(kind=RKIND),intent(in):: xtime_s
+!---------------------------------------------------------------------------------------------
+
+ if(.not.allocated(f_ice) ) allocate(f_ice(ims:ime,kms:kme,jms:jme) )
+ if(.not.allocated(f_rain) ) allocate(f_rain(ims:ime,kms:kme,jms:jme) )
+
if(.not.allocated(sfc_emiss_p) ) allocate(sfc_emiss_p(ims:ime,jms:jme) )
if(.not.allocated(snow_p) ) allocate(snow_p(ims:ime,jms:jme) )
if(.not.allocated(tsk_p) ) allocate(tsk_p(ims:ime,jms:jme) )
@@ -93,12 +99,21 @@
if(.not.allocated(aerosolcp_p) ) &
allocate(aerosolcp_p(ims:ime,1:num_aerlevels,jms:jme,num_aerosols) )
- if(.not.allocated(emstot_p) ) allocate(emstot_p(ims:ime,kms:kme,jms:jme) )
- if(.not.allocated(abstot_p) ) &
- allocate(abstot_p(ims:ime,kms:kme,cam_abs_dim2,jms:jme) )
- if(.not.allocated(absnxt_p) ) &
- allocate(absnxt_p(ims:ime,kms:kme,cam_abs_dim1,jms:jme) )
+ !allocate these arrays on the first time step, only:
+ if(xtime_s .lt. 1.e-12) then
+ write(0,*)
+ write(0,*) '--- end subroutine allocate_radiation_lw:'
+ write(0,*) '--- allocate emstot_p,abstot_p,absnxt_p'
+
+ if(.not.allocated(emstot_p) ) allocate(emstot_p(ims:ime,kms:kme,jms:jme) )
+ if(.not.allocated(abstot_p) ) &
+ allocate(abstot_p(ims:ime,kms:kme,cam_abs_dim2,jms:jme) )
+ if(.not.allocated(absnxt_p) ) &
+ allocate(absnxt_p(ims:ime,kms:kme,cam_abs_dim1,jms:jme) )
+
+ endif
+
case default
end select radiation_lw_select
@@ -106,9 +121,14 @@
end subroutine allocate_radiation_lw
!=============================================================================================
- subroutine deallocate_radiation_lw
+ subroutine deallocate_radiation_lw(xtime_s)
!=============================================================================================
+!input arguments:
+ real(kind=RKIND),intent(in):: xtime_s
+
+!---------------------------------------------------------------------------------------------
+
if(allocated(f_ice) ) deallocate(f_ice )
if(allocated(f_rain) ) deallocate(f_rain )
if(allocated(sfc_emiss_p) ) deallocate(sfc_emiss_p )
@@ -169,10 +189,6 @@
if(allocated(aerosolcn_p) ) deallocate(aerosolcn_p )
if(allocated(aerosolcp_p) ) deallocate(aerosolcp_p )
- if(allocated(abstot_p) ) deallocate(abstot_p )
- if(allocated(absnxt_p) ) deallocate(absnxt_p )
- if(allocated(emstot_p) ) deallocate(emstot_p )
-
case default
end select radiation_lw_select
@@ -180,7 +196,7 @@
end subroutine deallocate_radiation_lw
!=============================================================================================
- subroutine radiation_lw_from_MPAS(mesh,state,diag_physics,sfc_input)
+ subroutine radiation_lw_from_MPAS(mesh,state,diag_physics,sfc_input,xtime_s)
!=============================================================================================
!input arguments:
@@ -189,6 +205,8 @@
type(sfc_input_type) ,intent(in):: sfc_input
type(diag_physics_type),intent(in):: diag_physics
+ real(kind=RKIND),intent(in):: xtime_s
+
!---------------------------------------------------------------------------------------------
do j = jts,jte
@@ -283,42 +301,34 @@
enddo
enddo
- call mpas_timer_start("CAM lw: read arrays for infrared absorption")
- !infrared absorption:
- do j = jts,jte
- do n = 1,cam_abs_dim1
- do k = kts,kte
- do i = its,ite
- absnxt_p(i,k,n,j) = diag_physics % absnxt % array(k,n,i)
- enddo
- enddo
- enddo
- do n = 1,cam_abs_dim2
- do k = kts,kte+1
- do i = its,ite
- abstot_p(i,k,n,j) = diag_physics % abstot % array(k,n,i)
- enddo
- enddo
- enddo
- do k = kts,kte+1
- do i = its,ite
- emstot_p(i,k,j) = diag_physics % emstot % array(k,i)
- enddo
- enddo
- enddo
- call mpas_timer_stop("CAM lw: read arrays for infrared absorption")
-! write(0,*) '--- end radiation_lw_from_MPAS: doabsems=',doabsems
-! do k = kts,kte+1
-! write(0,102) k,(maxval(abstot_p(its:ite,k,n,jts:jte)),n=1,10)
-! enddo
-! write(0,*)
-! do k = kts,kte
-! write(0,102) k,(maxval(absnxt_p(its:ite,k,n,jts:jte)),n=1,cam_abs_dim1)
-! enddo
-! write(0,*)
-! do k = kts,kte+1
-! write(0,102) k,maxval(emstot_p(its:ite,k,jts:jte))
-! enddo
+ call mpas_timer_start("CAM lw: fill arrays for infrared absorption")
+ if(xtime_s .lt. 1.e-12) then
+ write(0,*)
+ write(0,*) '--- radiation_lw_from_MPAS:'
+ write(0,*) '--- initialize emstot_p,abstot_p,absnxt_p'
+ !infrared absorption:
+ do j = jts,jte
+ do n = 1,cam_abs_dim1
+ do k = kts,kte
+ do i = its,ite
+ absnxt_p(i,k,n,j) = 0.
+ enddo
+ enddo
+ enddo
+ do n = 1,cam_abs_dim2
+ do k = kts,kte+1
+ do i = its,ite
+ abstot_p(i,k,n,j) = 0.
+ enddo
+ enddo
+ enddo
+ do k = kts,kte+1
+ do i = its,ite
+ emstot_p(i,k,j) = 0.
+ enddo
+ enddo
+ enddo
+ endif
call mpas_timer_start("CAM lw: ozone and aerosols")
!ozone mixing ratio:
@@ -367,13 +377,15 @@
end subroutine radiation_lw_from_MPAS
!=============================================================================================
- subroutine radiation_lw_to_MPAS(diag_physics,tend_physics)
+ subroutine radiation_lw_to_MPAS(diag_physics,tend_physics,xtime_s)
!=============================================================================================
!input arguments:
type(diag_physics_type),intent(inout):: diag_physics
type(tend_physics_type),intent(inout):: tend_physics
+ real(kind=RKIND),intent(in):: xtime_s
+
!---------------------------------------------------------------------------------------------
do j = jts,jte
@@ -407,48 +419,6 @@
enddo
enddo
- radiation_lw_select: select case (trim(radt_lw_scheme))
-
- case("cam_lw")
- !infrared absorption:
- do j = jts,jte
- do n = 1,cam_abs_dim1
- do k = kts,kte
- do i = its,ite
- diag_physics % absnxt % array(k,n,i) = absnxt_p(i,k,n,j)
- enddo
- enddo
- enddo
- do n = 1,cam_abs_dim2
- do k = kts,kte+1
- do i = its,ite
- diag_physics % abstot % array(k,n,i) = abstot_p(i,k,n,j)
- enddo
- enddo
- enddo
- do k = kts,kte+1
- do i = its,ite
- diag_physics % emstot % array(k,i) = emstot_p(i,k,j)
- enddo
- enddo
- enddo
-! write(0,*) '--- end radiation_lw_to_MPAS: doabsems=',doabsems
-! do k = kts,kte+1
-! write(0,102) k,(maxval(abstot_p(its:ite,k,n,jts:jte)),n=1,10)
-! enddo
-! write(0,*)
-! do k = kts,kte
-! write(0,102) k,(maxval(absnxt_p(its:ite,k,n,jts:jte)),n=1,cam_abs_dim1)
-! enddo
-! write(0,*)
-! do k = kts,kte+1
-! write(0,102) k,maxval(emstot_p(its:ite,k,jts:jte))
-! enddo
-
- case default
-
- end select radiation_lw_select
-
!format:
101 format(i3,2i6,12(1x,e15.8))
102 format(i6,12(1x,e15.8))
@@ -517,7 +487,7 @@
101 format(i8,12(1x,e15.8))
!copy all MPAS arrays to rectangular grid:
- call radiation_lw_from_MPAS(mesh,state,diag_physics,sfc_input)
+ call radiation_lw_from_MPAS(mesh,state,diag_physics,sfc_input,xtime_s)
!call to longwave radiation scheme:
radiation_lw_select: select case (trim(radt_lw_scheme))
@@ -560,6 +530,8 @@
radt = dt_radtlw/60.
call mpas_timer_start("camrad")
+ write(0,*) '--- enter subroutine camrad_lw: doabsems=',doabsems
+ call mpas_timer_start("camrad")
call camrad( dolw = .true. , dosw = .false. , &
rthratenlw = rthratenlw_p , rthratensw = rthratensw_p , &
swupt = swupt_p , swuptc = swuptc_p , &
@@ -579,7 +551,7 @@
xlat = xlat_p , xlong = xlon_p , &
t_phy = t_p , pi_phy = pi_p , &
p_phy = pres_p , p8w = pres2_p , &
- z = z_p , dz8w = dz_p , &
+ z = zmid_p , dz8w = dz_p , &
rho_phy = rho_p , qv3d = qv_p , &
qc3d = qc_p , qr3d = qr_p , &
qi3d = qi_p , qs3d = qs_p , &
@@ -609,28 +581,27 @@
its = its , ite = ite , jts = jts , jte = jte , kts = kts , kte = kte &
)
call mpas_timer_stop("camrad")
- write(0,*) '--- end subroutince camrad lw: doabsems=',doabsems
write(0,*) 'max lwupb =',maxval(lwupb_p(its:ite,jms:jme))
- write(0,*) 'min lwupb =',minval(lwupb_p(its:ite,jms:jme))
write(0,*) 'max lwupbc =',maxval(lwupbc_p(its:ite,jms:jme))
- write(0,*) 'min lwupbc =',minval(lwupbc_p(its:ite,jms:jme))
write(0,*) 'max lwupt =',maxval(lwupt_p(its:ite,jms:jme))
- write(0,*) 'min lwupt =',minval(lwupt_p(its:ite,jms:jme))
write(0,*) 'max lwuptc =',maxval(lwuptc_p(its:ite,jms:jme))
- write(0,*) 'min lwuptc =',minval(lwuptc_p(its:ite,jms:jme))
write(0,*) 'max rthratenlw =',maxval(rthratenlw_p(its:ite,kts:kte,jms:jme))
write(0,*) 'min rthratenlw =',minval(rthratenlw_p(its:ite,kts:kte,jms:jme))
+ write(0,*) '--- end subroutine camrad lw: doabsems=',doabsems
case default
end select radiation_lw_select
!copy all arrays back to MPAS geodesic grid:
- call radiation_lw_to_MPAS(diag_physics,tend_physics)
+ call radiation_lw_to_MPAS(diag_physics,tend_physics,xtime_s)
- write(0,*) '--- end subroutine driver_radiation_lw:'
+ write(0,*) '--- end subroutine driver_radiation_lw'
call mpas_timer_stop("radiation_lw")
+!formats:
+ 200 format(i3,i3,8(1x,e15.8))
+
end subroutine driver_radiation_lw
!=============================================================================================
Modified: trunk/mpas/src/core_atmos_physics/mpas_atmphys_driver_radiation_sw.F
===================================================================
--- trunk/mpas/src/core_atmos_physics/mpas_atmphys_driver_radiation_sw.F 2012-01-13 22:13:27 UTC (rev 1368)
+++ trunk/mpas/src/core_atmos_physics/mpas_atmphys_driver_radiation_sw.F 2012-01-14 00:47:18 UTC (rev 1369)
@@ -27,12 +27,17 @@
contains
!=============================================================================================
- subroutine allocate_radiation_sw
+ subroutine allocate_radiation_sw(xtime_s)
!=============================================================================================
- if(.not.allocated(f_ice) ) allocate(f_ice(ims:ime,kms:kme,jms:jme) )
- if(.not.allocated(f_rain) ) allocate(f_rain(ims:ime,kms:kme,jms:jme) )
+!input arguments:
+ real(kind=RKIND),intent(in):: xtime_s
+!---------------------------------------------------------------------------------------------
+
+ if(.not.allocated(f_ice) ) allocate(f_ice(ims:ime,kms:kme,jms:jme) )
+ if(.not.allocated(f_rain) ) allocate(f_rain(ims:ime,kms:kme,jms:jme) )
+
if(.not.allocated(xlat_p) ) allocate(xlat_p(ims:ime,jms:jme) )
if(.not.allocated(xlon_p) ) allocate(xlon_p(ims:ime,jms:jme) )
@@ -65,9 +70,6 @@
if(.not.allocated(swupflxc_p) ) allocate(swupflxc_p(ims:ime,kms:kme+1,jms:jme) )
case("cam_sw")
- if(.not.allocated(xlat_p) ) allocate(xlat_p(ims:ime,jms:jme) )
- if(.not.allocated(xlon_p) ) allocate(xlon_p(ims:ime,jms:jme) )
-
if(.not.allocated(glw_p) ) allocate(glw_p(ims:ime,jms:jme) )
if(.not.allocated(lwcf_p) ) allocate(lwcf_p(ims:ime,jms:jme) )
if(.not.allocated(lwdnb_p) ) allocate(lwdnb_p(ims:ime,jms:jme) )
@@ -98,12 +100,17 @@
if(.not.allocated(aerosolcp_p) ) &
allocate(aerosolcp_p(ims:ime,1:num_aerlevels,jms:jme,num_aerosols) )
- if(.not.allocated(emstot_p) ) allocate(emstot_p(ims:ime,kms:kme,jms:jme) )
- if(.not.allocated(abstot_p) ) &
- allocate(abstot_p(ims:ime,kms:kme,cam_abs_dim2,jms:jme) )
- if(.not.allocated(absnxt_p) ) &
- allocate(absnxt_p(ims:ime,kms:kme,cam_abs_dim1,jms:jme) )
+ !allocate these arrays on the first time step, only:
+ if(xtime_s .lt. 1.e-12) then
+ if(.not.allocated(emstot_p) ) allocate(emstot_p(ims:ime,kms:kme,jms:jme) )
+ if(.not.allocated(abstot_p) ) &
+ allocate(abstot_p(ims:ime,kms:kme,cam_abs_dim2,jms:jme) )
+ if(.not.allocated(absnxt_p) ) &
+ allocate(absnxt_p(ims:ime,kms:kme,cam_abs_dim1,jms:jme) )
+
+ endif
+
case default
end select radiation_sw_select
@@ -176,10 +183,6 @@
if(allocated(aerosolcn_p) ) deallocate(aerosolcn_p )
if(allocated(aerosolcp_p) ) deallocate(aerosolcp_p )
- if(allocated(abstot_p) ) deallocate(abstot_p )
- if(allocated(absnxt_p) ) deallocate(absnxt_p )
- if(allocated(emstot_p) ) deallocate(emstot_p )
-
case default
end select radiation_sw_select
@@ -187,7 +190,7 @@
end subroutine deallocate_radiation_sw
!=============================================================================================
- subroutine radiation_sw_from_MPAS(mesh,state,diag_physics,sfc_input)
+ subroutine radiation_sw_from_MPAS(mesh,state,diag_physics,sfc_input,xtime_s)
!=============================================================================================
!input arguments:
@@ -196,6 +199,8 @@
type(diag_physics_type),intent(inout):: diag_physics
type(sfc_input_type) ,intent(inout):: sfc_input
+ real(kind=RKIND),intent(in):: xtime_s
+
!---------------------------------------------------------------------------------------------
do j = jts,jte
@@ -265,7 +270,6 @@
enddo
case("cam_sw")
-
do j = jts,jte
do i = its,ite
sfc_emiss_p(i,j) = diag_physics % sfc_emiss % array(i)
@@ -292,39 +296,29 @@
enddo
enddo
!infrared absorption:
- do j = jts,jte
- do n = 1,cam_abs_dim1
- do k = kts,kte
- do i = its,ite
- absnxt_p(i,k,n,j) = diag_physics % absnxt % array(k,n,i)
- enddo
- enddo
- enddo
- do n = 1,cam_abs_dim2
- do k = kts,kte+1
- do i = its,ite
- abstot_p(i,k,n,j) = diag_physics % abstot % array(k,n,i)
- enddo
- enddo
- enddo
- do k = kts,kte+1
- do i = its,ite
- emstot_p(i,k,j) = diag_physics % emstot % array(k,i)
- enddo
- enddo
- enddo
-! write(0,*) '--- end radiation_sw_from_MPAS: doabsems=',doabsems
-! do k = kts,kte+1
-! write(0,102) k,(maxval(abstot_p(its:ite,k,n,jts:jte)),n=1,10)
-! enddo
-! write(0,*)
-! do k = kts,kte
-! write(0,102) k,(maxval(absnxt_p(its:ite,k,n,jts:jte)),n=1,cam_abs_dim1)
-! enddo
-! write(0,*)
-! do k = kts,kte+1
-! write(0,102) k,maxval(emstot_p(its:ite,k,jts:jte))
-! enddo
+ if(xtime_s .lt. 1.e-12) then
+ do j = jts,jte
+ do n = 1,cam_abs_dim1
+ do k = kts,kte
+ do i = its,ite
+ absnxt_p(i,k,n,j) = 0.
+ enddo
+ enddo
+ enddo
+ do n = 1,cam_abs_dim2
+ do k = kts,kte+1
+ do i = its,ite
+ abstot_p(i,k,n,j) = 0.
+ enddo
+ enddo
+ enddo
+ do k = kts,kte+1
+ do i = its,ite
+ emstot_p(i,k,j) = 0.
+ enddo
+ enddo
+ enddo
+ endif
!ozone mixing ratio:
do k = 1, num_oznlevels
pin_p(k) = mesh % pin % array(k,1)
@@ -412,35 +406,6 @@
enddo
- radiation_sw_select: select case (trim(radt_sw_scheme))
- case("cam_sw")
- !infrared absorption:
- do j = jts,jte
- do n = 1,cam_abs_dim1
- do k = kts,kte
- do i = its,ite
- diag_physics % absnxt % array(k,n,i) = absnxt_p(i,k,n,j)
- enddo
- enddo
- enddo
- do n = 1,cam_abs_dim2
- do k = kts,kte+1
- do i = its,ite
- diag_physics % abstot % array(k,n,i) = abstot_p(i,k,n,j)
- enddo
- enddo
- enddo
- do k = kts,kte+1
- do i = its,ite
- diag_physics % emstot % array(k,i) = emstot_p(i,k,j)
- enddo
- enddo
- enddo
-
- case default
-
- end select radiation_sw_select
-
end subroutine radiation_sw_to_MPAS
!=============================================================================================
@@ -456,7 +421,7 @@
!---------------------------------------------------------------------------------------------
write(0,*)
- write(0,*) '--- begin radiation_sw initialization:'
+ write(0,*) '--- enter radiation_sw initialization:'
!call to shortwave radiation scheme:
radiation_sw_select: select case (trim(radt_sw_scheme))
@@ -475,7 +440,7 @@
end select radiation_sw_select
- write(0,*) '--- end radiation_sw initialization:'
+ write(0,*) '--- end radiation_sw initialization'
end subroutine init_radiation_sw
@@ -517,7 +482,7 @@
xtime_m = xtime_s/60.
!copy all MPAS arrays to rectangular grid:
- call radiation_sw_from_MPAS(mesh,state,diag_physics,sfc_input)
+ call radiation_sw_from_MPAS(mesh,state,diag_physics,sfc_input,xtime_s)
!... calculates solar declination:
!call radconst(declin,solcon,julday,degrad,dpd)
@@ -531,6 +496,7 @@
radiation_sw_select: select case (trim(radt_sw_scheme))
case ("rrtmg_sw")
+
write(0,*) '--- enter subroutine rrtmg_swrad:'
call rrtmg_swrad( &
rthratensw = rthratensw_p , swupt = swupt_p , swuptc = swuptc_p , &
@@ -562,6 +528,8 @@
write(0,*) '--- exit subroutine rrtmg_swrad'
case ("cam_sw")
+
+ write(0,*) '--- enter subroutine camrad_sw:'
call camrad( dolw = .false. , dosw = .true. , &
rthratenlw = rthratenlw_p , rthratensw = rthratensw_p , &
swupt = swupt_p , swuptc = swuptc_p , &
@@ -581,7 +549,7 @@
xlat = xlat_p , xlong = xlon_p , &
t_phy = t_p , pi_phy = pi_p , &
p_phy = pres_p , p8w = pres2_p , &
- z = z_p , dz8w = dz_p , &
+ z = zmid_p , dz8w = dz_p , &
rho_phy = rho_p , qv3d = qv_p , &
qc3d = qc_p , qr3d = qr_p , &
qi3d = qi_p , qs3d = qs_p , &
@@ -610,19 +578,14 @@
ims = ims , ime = ime , jms = jms , jme = jme , kms = kms , kme = kme , &
its = its , ite = ite , jts = jts , jte = jte , kts = kts , kte = kte &
)
-
- write(0,*) '--- end subroutine camrad sw'
write(0,*) 'doabsems =',doabsems
write(0,*) 'max swupb =',maxval(swupb_p(its:ite,jms:jme))
- write(0,*) 'min swupb =',minval(swupb_p(its:ite,jms:jme))
write(0,*) 'max swupbc =',maxval(swupbc_p(its:ite,jms:jme))
- write(0,*) 'min swupbc =',minval(swupbc_p(its:ite,jms:jme))
write(0,*) 'max swupt =',maxval(swupt_p(its:ite,jms:jme))
- write(0,*) 'min swupt =',minval(swupt_p(its:ite,jms:jme))
write(0,*) 'max swuptc =',maxval(swuptc_p(its:ite,jms:jme))
- write(0,*) 'min swuptc =',minval(swuptc_p(its:ite,jms:jme))
write(0,*) 'max rthratensw =',maxval(rthratensw_p(its:ite,kts:kte,jms:jme))
write(0,*) 'min rthratensw =',minval(rthratensw_p(its:ite,kts:kte,jms:jme))
+ write(0,*) '--- end subroutine camrad sw'
case default
@@ -631,8 +594,11 @@
!copy all arrays back to MPAS geodesic grid:
call radiation_sw_to_MPAS(diag_physics,tend_physics)
- write(0,*) '--- end subroutine driver_radiation_sw:'
+ write(0,*) '--- end subroutine driver_radiation_sw'
+!formats:
+ 200 format(i3,i6,8(1x,e15.8))
+
end subroutine driver_radiation_sw
!=============================================================================================
Modified: trunk/mpas/src/core_atmos_physics/mpas_atmphys_initialize_real.F
===================================================================
--- trunk/mpas/src/core_atmos_physics/mpas_atmphys_initialize_real.F 2012-01-13 22:13:27 UTC (rev 1368)
+++ trunk/mpas/src/core_atmos_physics/mpas_atmphys_initialize_real.F 2012-01-14 00:47:18 UTC (rev 1369)
@@ -426,8 +426,8 @@
nCells = mesh % nCells
nSoilLevels = mesh % nSoilLevels
nFGSoilLevels = mesh % nFGSoilLevels
- write(0,*) 'nSoilLevels =',nSoilLevels
- write(0,*) 'nFGSoilLevels=',nFGSoilLevels
+ write(0,*) 'nSoilLevels =',nSoilLevels
+ write(0,*) 'nFGSoilLevels =',nFGSoilLevels
landmask => mesh % landmask % array
@@ -607,7 +607,7 @@
endif
!make sure that all the cells flagged as sea-ice cells are defined as ocean cells:
-!num_seaice_changes = 0
+ num_seaice_changes = 0
do iCell = 1, nCellsSolve
if((landmask(iCell).eq.1 .and. xice(iCell).gt.0.) .or. xice(iCell).gt.200.) then
num_seaice_changes = num_seaice_changes + 1
Modified: trunk/mpas/src/core_atmos_physics/mpas_atmphys_interface_nhyd.F
===================================================================
--- trunk/mpas/src/core_atmos_physics/mpas_atmphys_interface_nhyd.F 2012-01-13 22:13:27 UTC (rev 1368)
+++ trunk/mpas/src/core_atmos_physics/mpas_atmphys_interface_nhyd.F 2012-01-14 00:47:18 UTC (rev 1369)
@@ -23,15 +23,18 @@
subroutine allocate_forall_physics
!=============================================================================================
- if(.not.allocated(psfc_p)) allocate(psfc_p(ims:ime,jms:jme) )
- if(.not.allocated(ptop_p)) allocate(ptop_p(ims:ime,jms:jme) )
+ if(.not.allocated(psfc_p) ) allocate(psfc_p(ims:ime,jms:jme) )
+ if(.not.allocated(ptop_p) ) allocate(ptop_p(ims:ime,jms:jme) )
if(.not.allocated(u_p) ) allocate(u_p(ims:ime,kms:kme,jms:jme) )
if(.not.allocated(v_p) ) allocate(v_p(ims:ime,kms:kme,jms:jme) )
+ if(.not.allocated(fzm_p) ) allocate(fzm_p(ims:ime,kms:kme,jms:jme) )
+ if(.not.allocated(fzp_p) ) allocate(fzp_p(ims:ime,kms:kme,jms:jme) )
if(.not.allocated(zz_p) ) allocate(zz_p(ims:ime,kms:kme,jms:jme) )
if(.not.allocated(pres_p) ) allocate(pres_p(ims:ime,kms:kme,jms:jme) )
if(.not.allocated(pi_p) ) allocate(pi_p(ims:ime,kms:kme,jms:jme) )
if(.not.allocated(z_p) ) allocate(z_p(ims:ime,kms:kme,jms:jme) )
+ if(.not.allocated(zmid_p) ) allocate(zmid_p(ims:ime,kms:kme,jms:jme) )
if(.not.allocated(dz_p) ) allocate(dz_p(ims:ime,kms:kme,jms:jme) )
if(.not.allocated(t_p) ) allocate(t_p(ims:ime,kms:kme,jms:jme) )
if(.not.allocated(th_p) ) allocate(th_p(ims:ime,kms:kme,jms:jme) )
@@ -42,44 +45,7 @@
if(.not.allocated(w_p) ) allocate(w_p(ims:ime,kms:kme,jms:jme) )
if(.not.allocated(pres2_p)) allocate(pres2_p(ims:ime,kms:kme,jms:jme) )
if(.not.allocated(t2_p) ) allocate(t2_p(ims:ime,kms:kme,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) )
- do j = jms,jme
- do i = ims,ime
- psfc_p(i,j) = 0.
- ptop_p(i,j) = 0.
- enddo
- enddo
-
- do j = jms,jme
- do k = kms,kme
- do i = ims,ime
- u_p(i,k,j) = 0.
- v_p(i,k,j) = 0.
- w_p(i,k,j) = 0.
- pres_p(i,k,j) = 0.
- pi_p(i,k,j) = 0.
- z_p(i,k,j) = 0.
- dz_p(i,k,j) = 0.
- t_p(i,k,j) = 0.
- th_p(i,k,j) = 0.
- al_p(i,k,j) = 0.
- rho_p(i,k,j) = 0.
- rh_p(i,k,j) = 0.
-
- w_p(i,k,j) = 0.
- pres2_p(i,k,j) = 0.
- t2_p(i,k,j) = 0.
-
- pres_hyd_p(i,k,j) = 0.
- pres2_hyd_p(i,k,j) = 0.
- enddo
- enddo
- enddo
-
-!allocate moist species (to be revisited!):
if(.not.allocated(qv_p) ) allocate(qv_p(ims:ime,kms:kme,jms:jme) )
if(.not.allocated(qc_p) ) allocate(qc_p(ims:ime,kms:kme,jms:jme) )
if(.not.allocated(qr_p) ) allocate(qr_p(ims:ime,kms:kme,jms:jme) )
@@ -93,16 +59,18 @@
subroutine deallocate_forall_physics
!=============================================================================================
-!de-allocation of all physics arrays:
if(allocated(psfc_p) ) deallocate(psfc_p )
if(allocated(ptop_p) ) deallocate(ptop_p )
if(allocated(u_p) ) deallocate(u_p )
if(allocated(v_p) ) deallocate(v_p )
+ if(allocated(fzm_p) ) deallocate(fzm_p )
+ if(allocated(fzp_p) ) deallocate(fzp_p )
if(allocated(zz_p) ) deallocate(zz_p )
if(allocated(pres_p) ) deallocate(pres_p )
if(allocated(pi_p) ) deallocate(pi_p )
if(allocated(z_p) ) deallocate(z_p )
+ if(allocated(zmid_p) ) deallocate(zmid_p )
if(allocated(dz_p) ) deallocate(dz_p )
if(allocated(t_p) ) deallocate(t_p )
if(allocated(th_p) ) deallocate(th_p )
@@ -114,9 +82,6 @@
if(allocated(pres2_p) ) deallocate(pres2_p )
if(allocated(t2_p) ) deallocate(t2_p )
- if(allocated(pres_hyd_p) ) deallocate(pres_hyd_p )
- if(allocated(pres2_hyd_p)) deallocate(pres2_hyd_p )
-
if(allocated(qv_p) ) deallocate(qv_p )
if(allocated(qc_p) ) deallocate(qc_p )
if(allocated(qr_p) ) deallocate(qr_p )
@@ -141,12 +106,13 @@
real(kind=RKIND),dimension(:),pointer:: latCell,lonCell
real(kind=RKIND),dimension(:),pointer :: fzm,fzp,rdzw
+ real(kind=RKIND),dimension(:),pointer :: sfc_pressure
real(kind=RKIND),dimension(:,:),pointer:: zgrid
real(kind=RKIND),dimension(:,:),pointer:: zz,exner,pressure_b,rtheta_p,rtheta_b
real(kind=RKIND),dimension(:,:),pointer:: rho_zz,theta_m,qv,pressure_p,u,v,w
real(kind=RKIND),dimension(:,:),pointer:: qvs,rh
- integer:: ip,iEdg
+ real(kind=RKIND):: rho1,rho2,tem1,tem2
!---------------------------------------------------------------------------------------------
@@ -164,25 +130,50 @@
latCell => mesh % latCell % array
lonCell => mesh % lonCell % array
- fzm => mesh % fzm % array
- fzp => mesh % fzp % array
- rdzw => mesh % rdzw % array
- zgrid => mesh % zgrid % array
- zz => mesh % zz % array
- exner => diag % exner % array
- pressure_b => diag % pressure_base % array
- pressure_p => diag % pressure_p % array
- rtheta_p => diag % rtheta_p % array
- rtheta_b => diag % rtheta_base % array
+ fzm => mesh % fzm % array
+ fzp => mesh % fzp % array
+ rdzw => mesh % rdzw % array
+ zgrid => mesh % zgrid % array
+ zz => mesh % zz % array
+ sfc_pressure => diag % surface_pressure % array
+ exner => diag % exner % array
+ pressure_b => diag % pressure_base % array
+ pressure_p => diag % pressure_p % array
+ rtheta_p => diag % rtheta_p % array
+ rtheta_b => diag % rtheta_base % array
- rho_zz => state % rho_zz % array
- theta_m => state % theta_m % array
- qv => state % scalars % array(state%index_qv,:,:)
+ rho_zz => state % rho_zz % array
+ theta_m => state % theta_m % array
+ qv => state % scalars % array(state%index_qv,:,:)
- w => state % w % array
- u => diag % uReconstructZonal % array
- v => diag % uReconstructMeridional % array
+ w => state % w % array
+ u => diag % uReconstructZonal % array
+ v => diag % uReconstructMeridional % array
+!ldf (2012-01-06): updates the surface pressure as is done in subroutine microphysics_to_MPAS.
+!do j = jts,jte
+!do i = its,ite
+! sfc_pressure(i) = 0.5*g*(zgrid(2,i)-zgrid(1,i)) &
+! * (1.25 * rho_zz(1,i) * zz(1,i) * (1. + qv(1,i)) &
+! - 0.25 * rho_zz(2,i) * zz(2,i) * (1. + qv(1,i)))
+! sfc_pressure(i) = sfc_pressure(i) + pressure_p(1,i) + pressure_b(1,i)
+!enddo
+!enddo
+!ldf end.
+!ldf (2012-01-09): updates the surface pressure using zgrid.
+ do j = jts,jte
+ do i = its,ite
+ tem1 = zgrid(2,i)-zgrid(1,i)
+ tem2 = zgrid(3,i)-zgrid(2,i)
+ rho1 = rho_zz(1,i) * zz(1,i) * (1. + qv(1,i))
+ rho2 = rho_zz(2,i) * zz(2,i) * (1. + qv(2,i))
+ sfc_pressure(i) = 0.5*g*(zgrid(2,i)-zgrid(1,i)) &
+ * (rho1 + 0.5*(rho2-rho1)*tem1/(tem1+tem2))
+ sfc_pressure(i) = sfc_pressure(i) + pressure_p(1,i) + pressure_b(1,i)
+ enddo
+ enddo
+!ldf end.
+
!copy sounding variables from the geodesic grid to the rectangular grid:
do j = jts, jte
do i = its, ite
@@ -215,11 +206,18 @@
pi_p(i,k,j) = exner(k,i)
pres_p(i,k,j) = pressure_p(k,i) + pressure_b(k,i)
- dz_p(i,k,j) = zgrid(k+1,i)-zgrid(k,i)
+ zmid_p(i,k,j) = 0.5*(zgrid(k+1,i)+zgrid(k,i))
+ dz_p(i,k,j) = zgrid(k+1,i)-zgrid(k,i)
- !arrays located at w points:
+ enddo
+ enddo
+ enddo
+
+ do j = jts, jte
+ do k = kts,kte+1
+ do i = its,ite
w_p(i,k,j) = w(k,i)
-
+ z_p(i,k,j) = zgrid(k,i)
enddo
enddo
enddo
@@ -247,11 +245,23 @@
enddo
!interpolation of pressure and temperature from theta points to w points:
+!do j = jts,jte
+!do k = kts+1,kte
+!do i = its,ite
+! t2_p(i,k,j) = fzm(k)*t_p(i,k,j) + fzp(k)*t_p(i,k-1,j)
+! pres2_p(i,k,j) = fzm(k)*pres_p(i,k,j) + fzp(k)*pres_p(i,k-1,j)
+!enddo
+!enddo
+!enddo
+!ldf(2011-01-10):
do j = jts,jte
do k = kts+1,kte
do i = its,ite
- t2_p(i,k,j) = fzm(k)*t_p(i,k,j) + fzp(k)*t_p(i,k-1,j)
- pres2_p(i,k,j) = fzm(k)*pres_p(i,k,j) + fzp(k)*pres_p(i,k-1,j)
+ tem1 = 1./(zgrid(k+1,i)-zgrid(k-1,i))
+ fzm_p(i,k,j) = (zgrid(k,i)-zgrid(k-1,i)) * tem1
+ fzp_p(i,k,j) = (zgrid(k+1,i)-zgrid(k,i)) * tem1
+ t2_p(i,k,j) = fzm_p(i,k,j)*t_p(i,k,j) + fzp_p(i,k,j)*t_p(i,k-1,j)
+ pres2_p(i,k,j) = fzm_p(i,k,j)*pres_p(i,k,j) + fzp_p(i,k,j)*pres_p(i,k-1,j)
enddo
enddo
enddo
@@ -311,17 +321,6 @@
enddo
enddo
-!calculation of the hydrostatic pressure at w points:
- do j = jts,jte
- do i = its,ite
- pres2_hyd_p(i,1,j) = psfc_p(i,j)
- do k = kts+1,kte+1
- pres2_hyd_p(i,k,j) = pres2_hyd_p(i,k-1,j) &
- - rho_p(i,k-1,j)*(1+qv_p(i,k-1,j))*g*dz_p(i,k-1,j)
- enddo
- enddo
- enddo
-
!formats:
201 format(3i8,10(1x,e15.8))
202 format(2i6,10(1x,e15.8))
@@ -480,6 +479,13 @@
real(kind=RKIND),dimension(:,:),pointer:: rho_zz,theta_m,pressure_p
real(kind=RKIND),dimension(:,:),pointer:: rt_diabatic_tend
+!ldf(2011-11-12): surface pressure.
+ real(kind=RKIND):: rho1,rho2,tem1,tem2
+ real(kind=RKIND),dimension(:),pointer:: rdzw
+ real(kind=RKIND),dimension(:),pointer:: sfc_pressure
+ real(kind=RKIND),dimension(:,:),pointer:: zgrid
+!ldf end.
+
!---------------------------------------------------------------------------------------------
write(0,*)
@@ -487,6 +493,7 @@
!initialization:
zz => mesh % zz % array
+ zgrid => mesh % zgrid % array
exner => diag % exner % array
exner_b => diag % exner_base % array
pressure_b => diag % pressure_base % array
@@ -499,6 +506,11 @@
rt_diabatic_tend => tend % rt_diabatic_tend % array
+!ldf (2011-11-12): update surface pressure.
+ rdzw => mesh % rdzw % array
+ sfc_pressure => diag % surface_pressure % array
+!ldf end.
+
!variables common to all cloud microphysics schemes:
do j = jts, jte
@@ -529,7 +541,30 @@
enddo
enddo
enddo
-
+
+!updates the surface pressure.
+!do j = jts,jte
+!do i = its,ite
+! sfc_pressure(i) = 0.5*g*(zgrid(2,i)-zgrid(1,i)) &
+! * (1.25 * rho_zz(1,i) * zz(1,i) * (1. + qv_p(i,1,j)) &
+! - 0.25 * rho_zz(2,i) * zz(2,i) * (1. + qv_p(i,2,j)))
+! sfc_pressure(i) = sfc_pressure(i) + pressure_p(1,i) + pressure_b(1,i)
+!enddo
+!enddo
+!ldf (2012-01-09):
+ do j = jts,jte
+ do i = its,ite
+ tem1 = zgrid(2,i)-zgrid(1,i)
+ tem2 = zgrid(3,i)-zgrid(2,i)
+ rho1 = rho_zz(1,i) * zz(1,i) * (1. + qv_p(i,1,j))
+ rho2 = rho_zz(2,i) * zz(2,i) * (1. + qv_P(i,2,j))
+ sfc_pressure(i) = 0.5*g*(zgrid(2,i)-zgrid(1,i)) &
+ * (rho1 + 0.5*(rho2-rho1)*tem1/(tem1+tem2))
+ sfc_pressure(i) = sfc_pressure(i) + pressure_p(1,i) + pressure_b(1,i)
+ enddo
+ enddo
+!ldf end.
+
!variables specific to different cloud microphysics schemes:
microp_select_init: select case(microp_scheme)
Modified: trunk/mpas/src/core_atmos_physics/mpas_atmphys_lsm_noahinit.F
===================================================================
--- trunk/mpas/src/core_atmos_physics/mpas_atmphys_lsm_noahinit.F 2012-01-13 22:13:27 UTC (rev 1368)
+++ trunk/mpas/src/core_atmos_physics/mpas_atmphys_lsm_noahinit.F 2012-01-14 00:47:18 UTC (rev 1369)
@@ -9,7 +9,7 @@
use mpas_configure, only: restart => config_do_restart, &
mminlu => input_landuse_data, &
mminsl => input_soil_data , &
- input_sfc_albedo => config_sfc_albedo
+ input_sfc_albedo => config_sfc_snowalbedo
use mpas_dmpar
use mpas_grid_types
Modified: trunk/mpas/src/core_atmos_physics/mpas_atmphys_todynamics.F
===================================================================
--- trunk/mpas/src/core_atmos_physics/mpas_atmphys_todynamics.F 2012-01-13 22:13:27 UTC (rev 1368)
+++ trunk/mpas/src/core_atmos_physics/mpas_atmphys_todynamics.F 2012-01-14 00:47:18 UTC (rev 1369)
@@ -3,7 +3,7 @@
use mpas_configure
use mpas_grid_types
- use mpas_atmphys_constants, only: R_d,R_v
+ use mpas_atmphys_constants, only: R_d,R_v,degrad
implicit none
private
@@ -34,9 +34,9 @@
!local variables:
!----------------
- integer:: i,k,nCells,nCellsSolve,nEdges,nEdgesSolve,nVertLevels
+ integer:: i,iCell,k,n,nCells,nCellsSolve,nEdges,nEdgesSolve,nVertLevels
- real(kind=RKIND),dimension(:,:),pointer:: theta,theta_m,qv
+ real(kind=RKIND),dimension(:,:),pointer:: theta_m,qv
real(kind=RKIND),dimension(:,:),pointer:: rthblten,rqvblten,rqcblten, &
rqiblten,rublten,rvblten
real(kind=RKIND),dimension(:,:),pointer:: rthcuten,rqvcuten,rqccuten, &
@@ -49,6 +49,10 @@
real(kind=RKIND):: tem
real(kind=RKIND),dimension(:,:),allocatable:: rublten_Edge
+!ldf (2011-12-16):
+ real(kind=RKIND),dimension(:,:),allocatable:: theta,tend_th
+!ldf end.
+
!=============================================================================================
!write(0,*)
!write(0,*) '--- enter subroutine physics_add_tend:'
@@ -59,7 +63,7 @@
nEdgesSolve = mesh % nEdgesSolve
nVertLevels = mesh % nVertLevels
- theta => diag % theta % array
+!theta => diag % theta % array
theta_m => state % theta_m % array
qv => state % scalars % array(state%index_qv,:,:)
@@ -84,6 +88,11 @@
tend_theta => tend % theta_m % array
tend_scalars => tend % scalars % array
+!initialize the tendency for the potential temperature and all scalars due to PBL, convection,
+!and longwave and shortwave radiation:
+ allocate(theta(nVertLevels,nCellsSolve) )
+ allocate(tend_th(nVertLevels,nCellsSolve))
+ tend_th = 0.
tend_scalars = 0.
!add coupled tendencies due to PBL processes:
@@ -96,24 +105,25 @@
tend_u(k,i)=tend_u(k,i)+rublten_Edge(k,i)*mass_edge(k,i)
enddo
enddo
-
deallocate(rublten_Edge)
do i = 1, nCellsSolve
do k = 1, nVertLevels
- tend_theta(k,i)=tend_theta(k,i)+rthblten(k,i)*mass(k,i)
+ tend_th(k,i) = tend_th(k,i)+rthblten(k,i)*mass(k,i)
tend_scalars(tend%index_qv,k,i)=tend_scalars(tend%index_qv,k,i)+rqvblten(k,i)*mass(k,i)
tend_scalars(tend%index_qc,k,i)=tend_scalars(tend%index_qc,k,i)+rqcblten(k,i)*mass(k,i)
tend_scalars(tend%index_qi,k,i)=tend_scalars(tend%index_qi,k,i)+rqiblten(k,i)*mass(k,i)
enddo
enddo
endif
+ write(0,*) 'max rthblten = ',maxval(rthblten(:,1:nCellsSolve))
+ write(0,*) 'min rthblten = ',minval(rthblten(:,1:nCellsSolve))
!add coupled tendencies due to convection:
if(config_conv_deep_scheme .ne. 'off') then
do i = 1, nCellsSolve
do k = 1, nVertLevels
- tend_theta(k,i)=tend_theta(k,i)+rthcuten(k,i)*mass(k,i)
+ tend_th(k,i)=tend_th(k,i)+rthcuten(k,i)*mass(k,i)
tend_scalars(tend%index_qv,k,i)=tend_scalars(tend%index_qv,k,i)+rqvcuten(k,i)*mass(k,i)
tend_scalars(tend%index_qc,k,i)=tend_scalars(tend%index_qc,k,i)+rqccuten(k,i)*mass(k,i)
tend_scalars(tend%index_qr,k,i)=tend_scalars(tend%index_qr,k,i)+rqrcuten(k,i)*mass(k,i)
@@ -122,64 +132,55 @@
enddo
enddo
endif
+ write(0,*) 'max rthcuten = ',maxval(rthcuten(:,1:nCellsSolve))
+ write(0,*) 'min rthcuten = ',minval(rthcuten(:,1:nCellsSolve))
!add coupled tendencies due to longwave radiation:
if(config_radt_lw_scheme .ne. 'off') then
do i = 1, nCellsSolve
do k = 1, nVertLevels
- tend_theta(k,i)=tend_theta(k,i)+rthratenlw(k,i)*mass(k,i)
+ tend_th(k,i)=tend_th(k,i)+rthratenlw(k,i)*mass(k,i)
enddo
enddo
endif
+ write(0,*) 'max rthratenlw = ',maxval(rthratenlw(:,1:nCellsSolve))
+ write(0,*) 'min rthratenlw = ',minval(rthratenlw(:,1:nCellsSolve))
!add coupled tendencies due to shortwave radiation:
if(config_radt_sw_scheme .ne. 'off') then
do i = 1, nCellsSolve
do k = 1, nVertLevels
- tend_theta(k,i)=tend_theta(k,i)+rthratensw(k,i)*mass(k,i)
+ tend_th(k,i)=tend_th(k,i)+rthratensw(k,i)*mass(k,i)
enddo
enddo
endif
+ write(0,*) 'max rthratensw = ',maxval(rthratensw(:,1:nCellsSolve))
+ write(0,*) 'min rthratensw = ',minval(rthratensw(:,1:nCellsSolve))
!if non-hydrostatic core, convert the tendency for the potential temperature to a
!tendency for the modified potential temperature:
#ifdef non_hydrostatic_core
do i = 1, nCellsSolve
do k = 1, nVertLevels
-
theta(k,i) = theta_m(k,i) / (1. + R_v/R_d * qv(k,i))
- tend_theta(k,i) = (1. + R_v/R_d * qv(k,i)) * tend_theta(k,i) &
+ tend_th(k,i) = (1. + R_v/R_d * qv(k,i)) * tend_th(k,i) &
+ R_v/R_d * theta(k,i) * tend_scalars(tend%index_qv,k,i)
+ tend_theta(k,i) = tend_theta(k,i) + tend_th(k,i)
enddo
- enddo
+ enddo
+#elif hydrostatic_core
+ do i = 1, nCellsSolve
+ do k = 1, nVertLevels
+ tend_theta(k,i) = tend_theta(k,i) + tend_th(k,i)
+ enddo
+ enddo
#endif
+ deallocate(theta)
+ deallocate(tend_th)
-!write(0,*) 'max PBL tendencies:'
-!write(0,*) 'max rthblten=',maxval(rthblten(:,:))
-!write(0,*) 'max rqvblten=',maxval(rqvblten(:,:))
-!write(0,*) 'max rqcblten=',maxval(rqcblten(:,:))
-!write(0,*) 'max rqiblten=',maxval(rqiblten(:,:))
-!write(0,*) 'max rublten =',maxval(rublten(:,:))
-!write(0,*) 'max rvblten =',maxval(rvblten(:,:))
-!write(0,*)
-!write(0,*) 'max CU tendencies:'
-!write(0,*) 'max rthcuten=',maxval(rthcuten(:,:))
-!write(0,*) 'max rqvcuten=',maxval(rqvcuten(:,:))
-!write(0,*) 'max rqccuten=',maxval(rqccuten(:,:))
-!write(0,*) 'max rqrcuten=',maxval(rqrcuten(:,:))
-!write(0,*) 'max rqicuten=',maxval(rqicuten(:,:))
-!write(0,*) 'max rqscuten=',maxval(rqscuten(:,:))
-!write(0,*)
-!write(0,*) 'max tend_scalars:'
-!write(0,*) 'max tend qv=',maxval(tend_scalars(tend%index_qv,:,:))
-!write(0,*) 'max tend qc=',maxval(tend_scalars(tend%index_qc,:,:))
-!write(0,*) 'max tend qr=',maxval(tend_scalars(tend%index_qr,:,:))
-!write(0,*) 'max tend qi=',maxval(tend_scalars(tend%index_qi,:,:))
-!write(0,*) 'max tend qs=',maxval(tend_scalars(tend%index_qs,:,:))
-!write(0,*)
-
!formats:
- 201 format(2i6,8(1x,e15.8))
+ 201 format(2i6,10(1x,e15.8))
+ 202 format(3i6,10(1x,e15.8))
end subroutine physics_addtend
Modified: trunk/mpas/src/core_atmos_physics/mpas_atmphys_vars.F
===================================================================
--- trunk/mpas/src/core_atmos_physics/mpas_atmphys_vars.F 2012-01-13 22:13:27 UTC (rev 1368)
+++ trunk/mpas/src/core_atmos_physics/mpas_atmphys_vars.F 2012-01-14 00:47:18 UTC (rev 1369)
@@ -34,7 +34,7 @@
integer,public:: its,ite,jts,jte,kts,kte
integer,public:: n_microp
- integer,public:: num_months !number of months [-]
+ integer,public:: num_months !number of months [-]
real(kind=RKIND),public:: dt_dyn !time-step for dynamics
real(kind=RKIND),public:: dt_microp !time-step for cloud microphysics parameterization.
@@ -44,15 +44,19 @@
real(kind=RKIND),public:: xice_threshold
real(kind=RKIND),dimension(:,:),allocatable:: &
- area_p !grid cell area [m2]
+ area_p !grid cell area [m2]
!... arrays related to surface:
real(kind=RKIND),dimension(:,:),allocatable:: &
psfc_p, &!surface pressure [Pa]
ptop_p !model-top pressure [Pa]
+ real(kind=RKIND),dimension(:,:,:),allocatable:: &
+ fzm_p, &!weight for interpolation to w points [-]
+ fzp_p !weight for interpolation to w points [-]
+
+ real(kind=RKIND),dimension(:,:,:),allocatable:: &
!... arrays related to u- and v-velocities interpolated to theta points:
- real(kind=RKIND),dimension(:,:,:),allocatable:: &
u_p, &!u-velocity interpolated to theta points [m/s]
v_p !v-velocity interpolated to theta points [m/s]
@@ -62,6 +66,7 @@
pres_p, &!pressure [Pa]
pi_p, &!(p_phy/p0)**(r_d/cp) [-]
z_p, &!height of layer [m]
+ zmid_p, &!height of middle of layer [m]
dz_p, &!layer thickness [m]
t_p, &!temperature [K]
th_p, &!potential temperature [K]
@@ -70,10 +75,6 @@
rh_p !relative humidity [-]
real(kind=RKIND),dimension(:,:,:),allocatable:: &
- pres_hyd_p, &!hydrostatic pressure at theta points [Pa]
- pres2_hyd_p !hydrostatic pressure at w points [Pa]
-
- real(kind=RKIND),dimension(:,:,:),allocatable:: &
qv_p, &!water vapor mixing ratio [kg/kg]
qc_p, &!cloud water mixing ratio [kg/kg]
qr_p, &!rain mixing ratio [kg/kg]
@@ -325,6 +326,7 @@
integer,public:: &
num_soils !number of soil layers [-]
+
integer,dimension(:,:),allocatable:: &
isltyp_p, &!dominant soil type category [-]
ivgtyp_p !dominant vegetation category [-]
Modified: trunk/mpas/src/core_init_nhyd_atmos/Registry
===================================================================
--- trunk/mpas/src/core_init_nhyd_atmos/Registry 2012-01-13 22:13:27 UTC (rev 1368)
+++ trunk/mpas/src/core_init_nhyd_atmos/Registry 2012-01-14 00:47:18 UTC (rev 1369)
@@ -229,7 +229,7 @@
var persistent real cqw ( nVertLevels nCells Time ) 1 - cqw diag - -
-var persistent real surface_pressure ( nCells Time ) 1 o surface_pressure diag - -
+var persistent real surface_pressure ( nCells Time ) 1 io surface_pressure diag - -
% coupled variables needed by the solver, but not output...
var persistent real ru ( nVertLevels nEdges Time ) 1 - ru diag - -
Modified: trunk/mpas/src/core_init_nhyd_atmos/mpas_init_atm_test_cases.F
===================================================================
--- trunk/mpas/src/core_init_nhyd_atmos/mpas_init_atm_test_cases.F 2012-01-13 22:13:27 UTC (rev 1368)
+++ trunk/mpas/src/core_init_nhyd_atmos/mpas_init_atm_test_cases.F 2012-01-14 00:47:18 UTC (rev 1369)
@@ -6,6 +6,8 @@
use mpas_dmpar
use atm_advection
use mpas_atmphys_initialize_real
+ use mpas_RBF_interpolation
+ use mpas_vector_reconstruction
! Added only clause to keep xlf90 from getting confused from the overloaded abs intrinsic in mpas_timekeeping
use mpas_timekeeping !, only: MPAS_Time_type, MPAS_TimeInterval_type, MPAS_Clock_type, &
@@ -3306,6 +3308,13 @@
if (config_met_interp) then
+ !ldf (2011-11-19): added initialization of the sea-surface temperature, seaice fraction, and
+ !seaice flag:
+ fg % sst % array = 0.
+ fg % xice % array = 0.
+ fg % seaice % array = 0.
+ !ldf end.
+
!
! First, try to locate the LANDSEA field for use as an interpolation mask
!
@@ -3823,7 +3832,8 @@
! Set SST based on SKINTEMP field if it wasn't found in input data
if (minval(fg % sst % array) == 0.0 .and. maxval(fg % sst % array) == 0.0) then
write(0,*) 'Setting SST from SKINTEMP'
- where (grid % landmask % array == 0) fg % sst % array = fg % skintemp % array
+ !where (grid % landmask % array == 0) fg % sst % array = fg % skintemp % array
+ fg % sst % array = fg % skintemp % array
end if
! Set SNOWC (snow-cover flag) based on SNOW
@@ -4085,6 +4095,20 @@
!
+ ! Reconstruct zonal and meridional winds for diagnostic puposes:
+ !
+ call mpas_rbf_interp_initialize(grid)
+ call mpas_init_reconstruct(grid)
+ call mpas_reconstruct(grid, state % u % array, &
+ diag % uReconstructX % array, &
+ diag % uReconstructY % array, &
+ diag % uReconstructZ % array, &
+ diag % uReconstructZonal % array, &
+ diag % uReconstructMeridional % array &
+ )
+
+
+ !
! Adjust surface pressure for difference in topography
!
do sfc_k=1,config_nfglevels
@@ -4243,7 +4267,8 @@
end if ! config_met_interp
- ! Calculate surface pressure
+ ! Calculate surface pressure (This is an ad-hoc calculation. The actual surface pressure is actually re-calculated at
+ !the top of the subroutine MPAS_to_physics in ../core_atmos_physics/mpas_atmphys_interface_nhyd.F
do iCell=1,grid%nCells
diag % surface_pressure % array(iCell) = 0.5*gravity/rdzw(1) &
* (1.25* rho_zz(1,iCell) * (1. + scalars(state % index_qv, 1, iCell)) &
Modified: trunk/mpas/src/core_nhyd_atmos/Registry
===================================================================
--- trunk/mpas/src/core_nhyd_atmos/Registry 2012-01-13 22:13:27 UTC (rev 1368)
+++ trunk/mpas/src/core_nhyd_atmos/Registry 2012-01-14 00:47:18 UTC (rev 1369)
@@ -329,6 +329,7 @@
%... NAMELIST VARIABLES ADDED FOR PHYSICS CONFIGURATION:
namelist logical physics config_frac_seaice false
namelist logical physics config_sfc_albedo false
+namelist logical physics config_sfc_snowalbedo false
namelist logical physics config_sst_update false
namelist logical physics config_sstdiurn_update false
namelist logical physics config_deepsoiltemp_update false
Modified: trunk/mpas/src/core_nhyd_atmos/mpas_atm_time_integration.F
===================================================================
--- trunk/mpas/src/core_nhyd_atmos/mpas_atm_time_integration.F 2012-01-13 22:13:27 UTC (rev 1368)
+++ trunk/mpas/src/core_nhyd_atmos/mpas_atm_time_integration.F 2012-01-14 00:47:18 UTC (rev 1369)
@@ -145,10 +145,6 @@
call mpas_dmpar_exch_halo_field2d_real(domain % dminfo, block % diag % rtheta_p % array(:,:), &
block % mesh % nVertLevels, block % mesh % nCells, &
block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
-!surface_pressure
- call mpas_dmpar_exch_halo_field1d_real(domain % dminfo, block % diag % surface_pressure % array(:), &
- block % mesh % nCells, &
- block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
block => block % next
end do
@@ -192,7 +188,7 @@
block => domain % blocklist
do while (associated(block))
call physics_addtend( domain % dminfo , block % parinfo % cellsToSend, block % parinfo % cellsToRecv, &
- block % mesh , block % state % time_levs(2) % state, block % diag, block % tend, &
+ block % mesh , block % state % time_levs(1) % state, block % diag, block % tend, &
block % tend_physics , block % state % time_levs(2) % state % rho_zz % array(:,:), &
block % diag % rho_edge % array(:,:) )
block => block % next
@@ -1010,7 +1006,7 @@
rho_pp, rho_zz, rho_base, ruAvg, ru_save, ru_p, u, ru, &
exner, exner_base, rtheta_base, pressure_p, &
zz, theta_m, pressure_b, qvapor
- real (kind=RKIND), dimension(:), pointer :: fzm, fzp, dvEdge, AreaCell, rdzw, surface_pressure
+ real (kind=RKIND), dimension(:), pointer :: fzm, fzp, dvEdge, AreaCell
real (kind=RKIND), dimension(:,:,:), pointer :: zb, zb3
integer, dimension(:,:), pointer :: cellsOnEdge
@@ -1053,9 +1049,7 @@
pressure_p => diag % pressure_p % array
pressure_b => diag % pressure_base % array
- surface_pressure => diag % surface_pressure % array
- rdzw => grid % rdzw % array
zz => grid % zz % array
zb => grid % zb % array
zb3 => grid % zb3 % array
@@ -1130,13 +1124,6 @@
* (exner(k,iCell)-exner_base(k,iCell)))
end do
- !calculation of the surface pressure:
- surface_pressure(iCell) = 0.5*gravity/rdzw(1) &
- * (1.25* rho_zz(1,iCell) * (1. + qvapor(1,iCell)) &
- - 0.25* rho_zz(2,iCell) * (1. + qvapor(2,iCell)))
- surface_pressure(iCell) = surface_pressure(iCell) + pressure_p(1,iCell) + pressure_b(1,iCell)
-
-
end do
! recover time-averaged ruAvg on all edges of owned cells (for upcoming scalar transport).
@@ -1872,7 +1859,8 @@
real (kind=RKIND) :: r_earth
real (kind=RKIND), dimension(:,:), pointer :: ur_cell, vr_cell
- real (kind=RKIND), parameter :: c_s = 0.25
+ real (kind=RKIND), parameter :: c_s = 0.125
+! real (kind=RKIND), parameter :: c_s = 0.25
real (kind=RKIND), dimension( grid % nVertLevels ) :: d_diag, d_off_diag, flux_arr
real (kind=RKIND), dimension(:,:), pointer :: defc_a, defc_b
logical :: delsq_horiz_mixing, newpx
</font>
</pre>