<p><b>laura@ucar.edu</b> 2011-01-13 16:36:43 -0700 (Thu, 13 Jan 2011)</p><p>Initialization of Noah land-surface model<br>
</p><hr noshade><pre><font color="gray">Added: branches/atmos_physics/src/core_physics/module_physics_lsm_noahinit.F
===================================================================
--- branches/atmos_physics/src/core_physics/module_physics_lsm_noahinit.F (rev 0)
+++ branches/atmos_physics/src/core_physics/module_physics_lsm_noahinit.F 2011-01-13 23:36:43 UTC (rev 689)
@@ -0,0 +1,476 @@
+#define DM_BCAST_CHAR(A) call dmpar_bcast_char(dminfo,A)
+#define DM_BCAST_INTEGER(A) call dmpar_bcast_int(dminfo,A)
+#define DM_BCAST_INTEGERS(A) call dmpar_bcast_ints(dminfo,size(A),A)
+#define DM_BCAST_REAL(A) call dmpar_bcast_real(dminfo,A)
+#define DM_BCAST_REALS(A) call dmpar_bcast_reals(dminfo,size(A),A)
+
+!=============================================================================================
+ module module_physics_lsm_noahinit
+ use configure, only: restart => config_do_restart, &
+ mminlu => input_landuse_data, &
+ mminsl => input_soil_data , &
+ input_sfc_albedo
+ use dmpar
+ use grid_types
+
+ use module_physics_constants, grav => g
+ use module_physics_utilities
+!wrf physics
+ use module_sf_noahlsm
+
+ implicit none
+ private
+ public:: noah_init_forMPAS
+
+ contains
+
+!=============================================================================================
+ subroutine noah_init_forMPAS(dminfo,mesh,diag_physics,sfc_input)
+!=============================================================================================
+
+!input arguments:
+ type(dm_info):: dminfo
+ type(mesh_type):: mesh
+
+!inout arguments:
+ type(diag_physics_type),intent(inout):: diag_physics
+ type(sfc_input_type),intent(inout):: sfc_input
+
+!---------------------------------------------------------------------------------------------
+
+!read formatted files needed for land-surface model:
+ call lsminit(dminfo,mesh,diag_physics,sfc_input)
+
+ end subroutine noah_init_forMPAS
+
+!=============================================================================================
+ subroutine lsminit(dminfo,mesh,diag_physics,sfc_input)
+!=============================================================================================
+
+!input arguments:
+ type(dm_info),intent(in):: dminfo
+ type(mesh_type),intent(in):: mesh
+
+!inout arguments:
+ type(diag_physics_type),intent(inout):: diag_physics
+ type(sfc_input_type),intent(inout):: sfc_input
+
+!local variables:
+ logical,parameter:: fndsnowh = .true.
+
+ integer:: iCell,nCells,nSoilLevels
+ integer:: errflag,ns
+
+ real(kind=RKIND):: bx,fk,smcmax,psisat,free
+ real(kind=RKIND),parameter:: blim = 5.5
+ real(kind=RKIND),parameter:: hlice = 3.335e5
+ real(kind=RKIND),parameter:: t0 = 273.15
+
+ integer,dimension(:),pointer:: ivgtyp,isltyp
+ real(kind=RKIND),dimension(:),pointer:: snoalb,snow,snowh
+ real(kind=RKIND),dimension(:,:),pointer:: tslb,smois,sh2o
+
+!---------------------------------------------------------------------------------------------
+
+ nCells = mesh % nCells
+ nSoilLevels = mesh % nSoilLevels
+
+ sh2o => diag_physics % sh2o % array
+ snow => diag_physics % snow % array
+ snowh => diag_physics % snowh % array
+
+ isltyp => sfc_input % isltyp % array
+ ivgtyp => sfc_input % ivgtyp % array
+ smois => sfc_input % smois % array
+ tslb => sfc_input % tslb % array
+ snoalb => sfc_input % snoalb % array
+
+!reads the NOAH LSM tables:
+ call physics_message( 'initialize NOAH LSM tables' )
+ call soil_veg_gen_parm(dminfo,mminlu,mminsl)
+
+ if(.not.restart) then
+
+ errflag = 0
+ do iCell = 1, nCells
+ if(isltyp(iCell) .lt. 1) then
+ errflag = 1
+ write(err_message,*) "module_sf_noahlsm.F: lsminit: out of range ISLTYP ", &
+ iCell,isltyp(iCell)
+ call physics_message(err_message)
+ endif
+ if(.not. input_sfc_albedo) snoalb(iCell) = maxalb(ivgtyp(iCell))*0.01
+ enddo
+! if(errflag .eq. 1) &
+! call physics_error_fatal("module_sf_noahlsm.F: lsminit: out of range value "// &
+! "of ISLTYP. Is this field in the input?" )
+
+!initializes soil liquid water content SH2O:
+ do iCell = 1, nCells
+ bx = bb(isltyp(iCell))
+ smcmax = maxsmc(isltyp(iCell))
+ psisat = satpsi(isltyp(iCell))
+ if((bx > 0.0).and.(smcmax > 0.0).and.(psisat > 0.0)) then
+ do ns = 1, nSoilLevels
+! ----------------------------------------------------------------------
+!SH2O <= SMOIS for T < 273.149K (-0.001C)
+ if(tslb(iCell,ns) < 273.149) then
+! ----------------------------------------------------------------------
+! first guess following explicit solution for Flerchinger Eqn from Koren
+! et al, JGR, 1999, Eqn 17 (KCOUNT=0 in FUNCTION FRH2O).
+! ISLTPK is soil type
+ bx = bb(isltyp(iCell))
+ smcmax = maxsmc(isltyp(iCell))
+ psisat = satpsi(isltyp(iCell))
+ if(bx > blim) bx = blim
+ fk = (((hlice/(grav*(-psisat))) * &
+ ((tslb(iCell,ns)-t0)/tslb(iCell,ns)) )**(-1/bx) )*smcmax
+ if (fk < 0.02) fk = 0.02
+ sh2o(iCell,ns) = min(fk,smois(iCell,ns))
+! ----------------------------------------------------------------------
+! now use iterative solution for liquid soil water content using
+! FUNCTION FRH2O with the initial guess for SH2O from above explicit
+! first guess.
+ call frh2o(free,tslb(iCell,ns),smois(iCell,ns),sh2o(iCell,ns), &
+ smcmax,bx,psisat)
+ sh2o(iCell,ns) = free
+ else ! of if (tslb(i,ns,j)
+! ----------------------------------------------------------------------
+! SH2O = SMOIS ( for T => 273.149K (-0.001C)
+ sh2o(iCell,ns)=smois(iCell,ns)
+! ----------------------------------------------------------------------
+ endif ! of if (tslb(i,ns,j)
+ enddo
+ else ! of if ((bx > 0.0)
+ do ns = 1, nSoilLevels
+ sh2o(iCell,ns)=smois(iCell,ns)
+ enddo
+ endif ! of if ((bx > 0.0)
+ enddo ! do iCell
+
+!initialize physical snow height SNOWH:
+ if(.not.fndsnowh)then
+!if no snowh do the following:
+ call physics_message( 'SNOW HEIGHT NOT FOUND - VALUE DEFINED IN LSMINIT' )
+ do iCell = 1, nCells
+ snowh(iCell)=snow(iCell)*0.005 ! snow in mm and snowh in m
+ enddo
+ endif
+ 110 continue
+
+ endif
+
+ end subroutine lsminit
+
+!=============================================================================================
+ subroutine soil_veg_gen_parm(dminfo,mminlu,mminsl)
+!=============================================================================================
+
+!input arguments:
+ type(dm_info),intent(in):: dminfo
+ character(len=15),intent(in):: mminlu, mminsl
+
+!local variables:
+ character*128:: mess , message
+
+ integer,parameter:: open_ok = 0
+ integer:: lumatch,iindex,lc,num_slope
+ integer:: istat
+
+!-----SPECIFY VEGETATION RELATED CHARACTERISTICS :
+!ALBBCK: SFC albedo (in percentage)
+!Z0 : Roughness length (m)
+!SHDFAC: Green vegetation fraction (in percentage)
+
+!Note: The ALBEDO, Z0, and SHDFAC values read from the following table
+!ALBEDO,and Z0 are specified in LAND-USE TABLE. SHDFAC is the monthly green vegetation data.
+!CMXTBL: MAX CNPY Capacity (m)
+!NROTBL: Rooting depth (layer)
+!RSMIN : Mimimum stomatal resistance (s m-1)
+!RSMAX : Max. stomatal resistance (s m-1)
+!RGL : Parameters used in radiation stress function
+!HS : Parameter used in vapor pressure deficit function
+!TOPT : Optimum transpiration air temperature. (K)
+!CMCMAX: Maximum canopy water capacity
+!CFACTR: Parameter used in the canopy inteception calculation
+!SNUP : Threshold snow depth (in water equivalent m) that implies 100% snow cover.
+!LAI : Leaf area index (dimensionless)
+!MAXALB: Upper bound on maximum albedo over deep snow
+!
+!-----READ IN VEGETATION PROPERTIES FROM VEGPARM.TBL
+
+!---------------------------------------------------------------------------------------------
+
+ write(0,*) ' enter subroutine soil_veg_gen_parm:'
+
+!read in the vegetation properties from vegparm.tbl:
+
+ if(dminfo % my_proc_id == IO_NODE) then
+ open(16,file='VEGPARM.TBL',form='FORMATTED',status='OLD',iostat=istat)
+ if(istat /= open_ok) &
+ call physics_error_fatal(istat,'subroutine soil_veg_gen_arm: ' // &
+ 'failure opening VEGPARM.TBL')
+
+ lumatch=0
+ find_lutype : do while (lumatch == 0)
+ read(16,*,end=2002)
+ read(16,*,end=2002) lutype
+ read(16,*) lucats,iindex
+
+ if(lutype.eq.trim(mminlu))then
+ write(mess,*) 'landuse type = ' // trim ( lutype ) // ' found', lucats,' categories'
+ call physics_message(mess)
+ lumatch=1
+ else
+ call physics_message("skipping over lutype = " // trim ( lutype ))
+ do lc = 1, lucats+12
+ read(16,*)
+ enddo
+ endif
+ enddo find_lutype
+
+!prevent possible array overwrite, Bill Bovermann, IBM, May 6, 2008:
+ if(size(shdtbl) < lucats .or. &
+ size(nrotbl) < lucats .or. &
+ size(rstbl) < lucats .or. &
+ size(rgltbl) < lucats .or. &
+ size(hstbl) < lucats .or. &
+ size(snuptbl) < lucats .or. &
+ size(maxalb) < lucats .or. &
+ size(laimintbl) < lucats .or. &
+ size(laimaxtbl) < lucats .or. &
+ size(z0mintbl) < lucats .or. &
+ size(z0maxtbl) < lucats .or. &
+ size(albedomintbl) < lucats .or. &
+ size(albedomaxtbl) < lucats .or. &
+ size(emissmintbl ) < lucats .or. &
+ size(emissmaxtbl ) < lucats) then
+! call wrf_error_fatal('table sizes too small for value of lucats in module_sf_noahdrv.f')
+ endif
+
+ if(lutype.eq.mminlu)then
+ do lc = 1, lucats
+ read(16,*) iindex,shdtbl(lc),nrotbl(lc),rstbl(lc),rgltbl(lc),hstbl(lc),snuptbl(lc), &
+ maxalb(lc),laimintbl(lc),laimaxtbl(lc),emissmintbl(lc),emissmaxtbl(lc), &
+ albedomintbl(lc),albedomaxtbl(lc),z0mintbl(lc),z0maxtbl(lc)
+ enddo
+ read (16,*)
+ read (16,*)topt_data
+ read (16,*)
+ read (16,*)cmcmax_data
+ read (16,*)
+ read (16,*)cfactr_data
+ read (16,*)
+ read (16,*)rsmax_data
+ read (16,*)
+ read (16,*)bare
+ read (16,*)
+ read (16,*)natural
+ endif
+
+ 2002 continue
+ close (16)
+ if(lumatch == 0) then
+! call wrf_error_fatal ("land use dataset '"//mminlu//"' not found in vegparm.tbl.")
+ endif
+
+ endif ! end dminfo
+
+!distribute data to all processors:
+ DM_BCAST_CHAR(lutype)
+ DM_BCAST_INTEGER(lucats)
+ DM_BCAST_INTEGER(iindex)
+ DM_BCAST_INTEGER(lumatch)
+ DM_BCAST_REALS(shdtbl)
+ DM_BCAST_INTEGERS(nrotbl)
+ DM_BCAST_REALS(rstbl)
+ DM_BCAST_REALS(rgltbl)
+ DM_BCAST_REALS(hstbl)
+ DM_BCAST_REALS(snuptbl)
+ DM_BCAST_REALS(laimintbl)
+ DM_BCAST_REALS(laimaxtbl)
+ DM_BCAST_REALS(z0mintbl)
+ DM_BCAST_REALS(z0maxtbl)
+ DM_BCAST_REALS(emissmintbl)
+ DM_BCAST_REALS(emissmaxtbl)
+ DM_BCAST_REALS(albedomintbl)
+ DM_BCAST_REALS(albedomaxtbl)
+ DM_BCAST_REALS(maxalb)
+ DM_BCAST_REAL(topt_data)
+ DM_BCAST_REAL(cmcmax_data)
+ DM_BCAST_REAL(cfactr_data)
+ DM_BCAST_REAL(rsmax_data)
+ DM_BCAST_INTEGER(bare)
+ DM_BCAST_INTEGER(natural)
+
+ write(0,*) ' LUTYPE = ',trim(lutype)
+ write(0,*) ' LUCATS = ',lucats
+ write(0,*) ' IINDEX = ',iindex
+ write(0,*) ' LUMATCH = ',lumatch
+
+ write(0,*) ' TOPT_DATA = ',topt_data
+ write(0,*) ' CMCMAX_DATA = ',cmcmax_data
+ write(0,*) ' CFACTR_DATA = ',cfactr_data
+ write(0,*) ' RSMAX_DATA = ',rsmax_data
+ write(0,*) ' BARE = ',bare
+ write(0,*) ' NATURAL = ',natural
+
+ write(0,*)
+ do lc = 1, lucats
+ write(0,101) lc,shdtbl(lc),float(nrotbl(lc)),rstbl(lc),rgltbl(lc),hstbl(lc),snuptbl(lc), &
+ maxalb(lc),laimintbl(lc),laimaxtbl(lc),emissmintbl(lc),emissmaxtbl(lc), &
+ albedomintbl(lc),albedomaxtbl(lc),z0mintbl(lc),z0maxtbl(lc)
+ enddo
+ 101 format(i4,15(3x,f9.5))
+
+!read in soil properties from soilparm.tbl:
+
+ if(dminfo % my_proc_id == IO_NODE) then
+ open(16,file='SOILPARM.TBL',form='FORMATTED',status='OLD',iostat=istat)
+ if(istat /= open_ok) &
+ call physics_error_fatal(istat,'module_sf_noahlsm.F: soil_veg_gen_parm: ' // &
+ 'failure opening SOILPARM.TBL' )
+
+ write(0,*)
+ write(mess,*) 'input soil texture classification = ', trim (mminsl)
+ call physics_message(mess)
+
+ lumatch=0
+ read(16,*)
+ read(16,2000,end=2003) sltype
+ 2000 format(a4)
+ read(16,*)slcats,iindex
+ if(sltype.eq.mminsl)then
+ write(mess,*) 'soil texture classification = ', trim ( sltype ) , ' found', &
+ slcats,' categories'
+! call wrf_message ( mess )
+ lumatch=1
+ endif
+
+!prevent possible array overwrite, Bill Bovermann, IBM, May 6, 2008:
+ if(size(bb) < slcats .or. &
+ size(drysmc) < slcats .or. &
+ size(f11 ) < slcats .or. &
+ size(maxsmc) < slcats .or. &
+ size(refsmc) < slcats .or. &
+ size(satpsi) < slcats .or. &
+ size(satdk ) < slcats .or. &
+ size(satdw ) < slcats .or. &
+ size(wltsmc) < slcats .or. &
+ size(qtz ) < slcats) then
+! call wrf_error_fatal('table sizes too small for value of slcats in module_sf_noahdrv.f')
+ endif
+ if(sltype.eq.mminsl) then
+ do lc = 1, slcats
+ read(16,*) iindex,bb(lc),drysmc(lc),f11(lc),maxsmc(lc),refsmc(lc),satpsi(lc), &
+ satdk(lc),satdw(lc),wltsmc(lc),qtz(lc)
+ enddo
+ endif
+
+ 2003 continue
+ close(16)
+ if(lumatch.eq.0)then
+ call physics_message( 'soil texture in input file does not ' )
+ call physics_message( 'match soilparm table' )
+! call wrf_error_fatal ( 'inconsistent or missing soilparm file' )
+ endif
+
+ endif
+
+!distribute data to all processors:
+ DM_BCAST_INTEGER(lumatch)
+ DM_BCAST_CHAR(sltype)
+!DM_BCAST_CHAR(trim(mminsl))
+ DM_BCAST_INTEGER(slcats)
+ DM_BCAST_INTEGER(iindex)
+ DM_BCAST_REALS(bb)
+ DM_BCAST_REALS(drysmc)
+ DM_BCAST_REALS(f11)
+ DM_BCAST_REALS(maxsmc)
+ DM_BCAST_REALS(refsmc)
+ DM_BCAST_REALS(satpsi)
+ DM_BCAST_REALS(satdk)
+ DM_BCAST_REALS(satdw)
+ DM_BCAST_REALS(wltsmc)
+ DM_BCAST_REALS(qtz)
+
+ write(0,*) ' LUMATCH=',lumatch
+ write(0,*) ' SLTYPE =',trim(sltype)
+ write(0,*) ' MMINSL =',mminsl
+ write(0,*) ' SLCATS =',slcats
+ write(0,*) ' IINDEX =',iindex
+ write(0,*)
+ do lc = 1, slcats
+ write(0,101) lc,bb(lc),drysmc(lc),f11(lc),maxsmc(lc),refsmc(lc),satpsi(lc), &
+ satdk(lc),satdw(lc),wltsmc(lc),qtz(lc)
+ enddo
+
+!read in general parameters from genparm.tbl:
+
+!if(dminfo % my_proc_id == IO_NODE) then
+! open(16,file='GENPARM.TBL',form='FORMATTED',status='OLD',iostat=istat)
+! if(istat /= open_ok) &
+! call physics_error_fatal(istat,'module_sf_noahlsm.F: soil_veg_gen_parm: ' // &
+! 'failure opening GENPARM.TBL' )
+! read(16,*)
+! read(16,*)
+! read(16,*) num_slope
+
+! slpcats=num_slope
+!prevent possible array overwrite, Bill Bovermann, IBM, May 6, 2008:
+! if(size(slope_data) < num_slope) &
+! call wrf_error_fatal('NUM_SLOPE too large for slope_data array in module_sf_noahdrv')
+
+! do lc = 1, slpcats
+! read(16,*)slope_data(lc)
+! enddo
+! read(16,*)
+! read(16,*)sbeta_data
+! read(16,*)
+! read(16,*)fxexp_data
+! read(16,*)
+! read(16,*)csoil_data
+! read(16,*)
+! read(16,*)salp_data
+! read(16,*)
+! read(16,*)refdk_data
+! read(16,*)
+! read(16,*)refkdt_data
+! read(16,*)
+! read(16,*)frzk_data
+! read(16,*)
+! read(16,*)zbot_data
+! read(16,*)
+! read(16,*)czil_data
+! read(16,*)
+! read(16,*)smlow_data
+! read(16,*)
+! read(16,*)smhigh_data
+! read(16,*)
+! read(16,*)lvcoef_data
+! close(16)
+!endif
+
+!DM_BCAST_INTEGER(num_slope)
+!DM_BCAST_INTEGER(slpcats)
+!DM_BCAST_REALS(slope_data)
+!DM_BCAST_REAL(sbeta_data)
+!DM_BCAST_REAL(fxexp_data)
+!DM_BCAST_REAL(csoil_data)
+!DM_BCAST_REAL(salp_data)
+!DM_BCAST_REAL(refdk_data)
+!DM_BCAST_REAL(refkdt_data)
+!DM_BCAST_REAL(frzk_data)
+!DM_BCAST_REAL(zbot_data)
+!DM_BCAST_REAL(czil_data)
+!DM_BCAST_REAL(smlow_data)
+!DM_BCAST_REAL(smhigh_data)
+!DM_BCAST_REAL(lvcoef_data)
+
+ write(0,*) ' end subroutine soil_veg_gen_parm:'
+
+ end subroutine soil_veg_gen_parm
+
+!=============================================================================================
+ end module module_physics_lsm_noahinit
+!=============================================================================================
</font>
</pre>