diff --git a/ccpp/data/GFS_typedefs.F90 b/ccpp/data/GFS_typedefs.F90 index 2666d59b0..526f12b3a 100644 --- a/ccpp/data/GFS_typedefs.F90 +++ b/ccpp/data/GFS_typedefs.F90 @@ -915,6 +915,7 @@ module GFS_typedefs integer :: iopt_snf !rainfall & snowfall (1-jordan91; 2->bats; 3->noah) integer :: iopt_tbot !lower boundary of soil temperature (1->zero-flux; 2->noah) integer :: iopt_stc !snow/soil temperature time scheme (only layer 1) + integer :: iopt_trs !thermal roughness scheme (1-z0h=z0m; 2-czil; 3-ec;4-kb inversed) logical :: use_ufo !< flag for gcycle surface option @@ -3333,6 +3334,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & integer :: iopt_snf = 1 !rainfall & snowfall (1-jordan91; 2->bats; 3->noah) integer :: iopt_tbot = 2 !lower boundary of soil temperature (1->zero-flux; 2->noah) integer :: iopt_stc = 1 !snow/soil temperature time scheme (only layer 1) + integer :: iopt_trs = 2 !thermal roughness scheme (1-z0h=z0m; 2-czil; 3-ec;4-kb reversed) logical :: use_ufo = .false. !< flag for gcycle surface option @@ -3681,6 +3683,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & ! Noah MP options iopt_dveg,iopt_crs,iopt_btr,iopt_run,iopt_sfc, iopt_frz, & iopt_inf, iopt_rad,iopt_alb,iopt_snf,iopt_tbot,iopt_stc, & + iopt_trs, & ! GFDL surface layer options lcurr_sf, pert_cd, ntsflg, sfenth, & !--- lake model control @@ -4297,6 +4300,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & Model%iopt_snf = iopt_snf Model%iopt_tbot = iopt_tbot Model%iopt_stc = iopt_stc + Model%iopt_trs = iopt_trs !--- tuning parameters for physical parameterizations Model%ras = ras @@ -5133,6 +5137,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & print *,'iopt_snf = ', Model%iopt_snf print *,'iopt_tbot = ',Model%iopt_tbot print *,'iopt_stc = ', Model%iopt_stc + print *,'iopt_trs = ', Model%iopt_trs elseif (Model%lsm == Model%lsm_ruc) then print *,' RUC Land Surface Model used' else @@ -5957,6 +5962,7 @@ subroutine control_print(Model) print *, ' iopt_snf : ', Model%iopt_snf print *, ' iopt_tbot : ', Model%iopt_tbot print *, ' iopt_stc : ', Model%iopt_stc + print *, ' iopt_trs : ', Model%iopt_trs endif print *, ' use_ufo : ', Model%use_ufo print *, ' lcurr_sf : ', Model%lcurr_sf diff --git a/ccpp/data/GFS_typedefs.meta b/ccpp/data/GFS_typedefs.meta index f8966b817..e6012a92e 100644 --- a/ccpp/data/GFS_typedefs.meta +++ b/ccpp/data/GFS_typedefs.meta @@ -4027,6 +4027,12 @@ units = index dimensions = () type = integer +[iopt_trs] + standard_name = control_for_land_surface_scheme_surface_thermal_roughness + long_name = choice for surface thermal roughness option (see noahmp module for definition) + units = index + dimensions = () + type = integer [use_ufo] standard_name = flag_for_gcycle_surface_option long_name = flag for gcycle surface option diff --git a/ccpp/physics b/ccpp/physics index e8dc72338..3ee7678bb 160000 --- a/ccpp/physics +++ b/ccpp/physics @@ -1 +1 @@ -Subproject commit e8dc7233804baa920a3c044a39d29a56d9b18930 +Subproject commit 3ee7678bb3c9ab889c584fc975715f5760917802 diff --git a/io/module_wrt_grid_comp.F90 b/io/module_wrt_grid_comp.F90 index 23d3a4047..d687d803d 100644 --- a/io/module_wrt_grid_comp.F90 +++ b/io/module_wrt_grid_comp.F90 @@ -3540,45 +3540,55 @@ end subroutine rtll subroutine lambert(stlat1,stlat2,c_lat,c_lon,glon,glat,x,y,inv) !------------------------------------------------------------------------------- - real(ESMF_KIND_R8), intent(in) :: stlat1,stlat2,c_lat,c_lon - real(ESMF_KIND_R8), intent(inout) :: glon, glat + real(ESMF_KIND_R8), intent(in) :: stlat1,stlat2,c_lat,c_lon + real(ESMF_KIND_R8), intent(inout) :: glon, glat real(ESMF_KIND_R8), intent(inout) :: x, y - integer, intent(in) :: inv + integer, intent(in) :: inv !------------------------------------------------------------------------------- -! real(ESMF_KIND_R8), parameter :: pi=3.14159265358979323846 - real(ESMF_KIND_R8), parameter :: dtor=pi/180.0 - real(ESMF_KIND_R8), parameter :: rtod=180.0/pi +! real(ESMF_KIND_R8), parameter :: pi = 3.14159265358979323846 + real(ESMF_KIND_R8), parameter :: dtor = pi/180.0 + real(ESMF_KIND_R8), parameter :: rtod = 180.0/pi real(ESMF_KIND_R8), parameter :: a = 6371200.0 !------------------------------------------------------------------------------- -! inv == 1 (glon,glat) ---> (x,y) lat/lon to grid -! inv == -1 (x,y) ---> (glon,glat) grid to lat/lon +! inv == 1 (glon,glat) ---> (x,y) +! inv == -1 (x,y) ---> (glon,glat) - real(ESMF_KIND_R8) :: en,f,rho,rho0, dlon, theta + real(ESMF_KIND_R8) :: xp, yp, en, de, rho, rho0, rho2, dlon, theta, dr2 + real(ESMF_KIND_R8) :: h = 1.0 - IF (stlat1 == stlat2) THEN - en=sin(stlat1*dtor) - ELSE - en=log(cos(stlat1*dtor)/cos(stlat2*dtor))/ & - log(tan((45+0.5*stlat2)*dtor)/tan((45+0.5*stlat1)*dtor)) - ENDIF + ! For reference see: + ! John P. Snyder (1987), Map projections: A working manual (pp. 104-110) + ! https://doi.org/10.3133/pp1395 - f=(cos(stlat1*dtor)*tan((45+0.5*stlat1)*dtor)**en)/en - rho0=a*f/(tan((45+0.5*c_lat)*dtor)**en) + if (stlat1 == stlat2) then + en = sin(stlat1*dtor) + else + en = log(cos(stlat1*dtor)/cos(stlat2*dtor)) / & + log(tan((45+0.5*stlat2)*dtor)/tan((45+0.5*stlat1)*dtor)) ! (15-3) + endif + h = sign(1.0_ESMF_KIND_R8,en) + + de = a*(cos(stlat1*dtor)*tan((45+0.5*stlat1)*dtor)**en)/en ! (15-2) + rho0 = de/(tan((45+0.5*c_lat)*dtor)**en) ! (15-1a) if (inv == 1) then ! FORWARD TRANSFORMATION - rho=a*f/(tan((45+0.5*glat)*dtor)**en) - dlon=modulo(glon-c_lon+180+3600,360.)-180.D0 - theta=en*dlon*dtor - x=rho*sin(theta) - y=rho0-rho*cos(theta) + rho = de/(tan((45+0.5*glat)*dtor)**en) ! (15-1) + dlon = modulo(glon-c_lon+180.0+3600.0,360.0)-180.0 + theta = en*dlon*dtor ! (14-4) + x = rho*sin(theta) ! (14-1) + y = rho0-rho*cos(theta) ! (14-2) else if (inv == -1) then ! INVERSE TRANSFORMATION - y=rho0-y - rho = sqrt(x*x+y*y) - theta=atan2(x,y) - glon=c_lon+(theta/en)*rtod - glon=modulo(glon+180+3600,360.)-180.D0 -! glat=(2.0*atan((a*f/rho)**(1.0/en))-0.5*pi)*rtod - glat=(0.5*pi-2.0*atan((rho/(a*f))**(1.0/en)))*rtod + xp = h*x; + yp = h*(rho0-y) + theta = atan2(xp,yp) ! (14-11) + glon = c_lon+(theta/en)*rtod ! (14-9) + glon = modulo(glon+180.0+3600.0,360.0)-180.0 + rho2 = xp*xp+yp*yp ! (14-10) + if (rho2 == 0.0) then + glat = h*90.0 + else + glat = 2.0*atan((de*de/rho2)**(1.0/(2.0*en)))*rtod-90.0 ! (15-5) + endif else write (unit=*,fmt=*) " lambert: unknown inv argument" return