From 213adc69497543b90eb6cfff4008858722c2a389 Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Sun, 22 May 2022 11:08:44 -0600 Subject: [PATCH 01/14] Update .gitmodules and submodule pointer for ccpp-physics for code review and testing --- .gitmodules | 6 ++++-- physics/rte-rrtmgp | 2 +- 2 files changed, 5 insertions(+), 3 deletions(-) diff --git a/.gitmodules b/.gitmodules index 75e5ea836..dc5a31dd7 100644 --- a/.gitmodules +++ b/.gitmodules @@ -1,4 +1,6 @@ [submodule "physics/rte-rrtmgp"] path = physics/rte-rrtmgp - url = https://github.com/earth-system-radiation/rte-rrtmgp - branch = dtc/ccpp + #url = https://github.com/earth-system-radiation/rte-rrtmgp + #branch = dtc/ccpp + url = https://github.com/climbfuji/rte-rrtmgp + branch = feature/single_precision diff --git a/physics/rte-rrtmgp b/physics/rte-rrtmgp index cec1e8e12..c2a8f5751 160000 --- a/physics/rte-rrtmgp +++ b/physics/rte-rrtmgp @@ -1 +1 @@ -Subproject commit cec1e8e12d969c3c8c76574dbe4f40b366419cc7 +Subproject commit c2a8f57518c2af0789e6d04c37c5415a7c045dfa From a769f64115cd9f5282a6985196fa3c5e2d7858ae Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Sun, 22 May 2022 11:10:32 -0600 Subject: [PATCH 02/14] Interface changes in physics/aer_cloud.F to compile in single precision --- physics/aer_cloud.F | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/physics/aer_cloud.F b/physics/aer_cloud.F index 60df592b6..a334428d1 100644 --- a/physics/aer_cloud.F +++ b/physics/aer_cloud.F @@ -3477,7 +3477,8 @@ SUBROUTINE EMPIRICAL_PARAM_PHILLIPS(SI, SIW, SW, D_grid_dust, & D_grid_bio, n_grid_bio, ijstop_bio, A_solo, n_iw, DSH, & Nhet_dep,Nhet_dhf,fdust_dep,P_ice, T_ice) implicit none - real, intent(IN):: SI, SIW, SW, A_solo,P_ice, T_ice + real, intent(IN):: SI, SIW, SW, A_solo + real*8, intent(IN):: P_ice, T_ice real, dimension(:), intent(IN):: D_grid_dust, n_grid_dust, & D_grid_soot, n_grid_soot, D_grid_bio, n_grid_bio integer, intent(IN):: ijstop_dust, ijstop_soot, ijstop_bio @@ -3488,7 +3489,7 @@ SUBROUTINE EMPIRICAL_PARAM_PHILLIPS(SI, SIW, SW, D_grid_dust, & num_ic_solo_imm real, intent (inout) :: DSH, n_iw - real, intent (out) :: Nhet_dep,Nhet_dhf,fdust_dep + real*8, intent (out) :: Nhet_dep,Nhet_dhf,fdust_dep real :: dn_in_dust, dn_in_soot, dn_in_bio, dn_in_solo, dNall, & dNaux, naux, SS_w, dH_frac_dust, dH_frac_soot, dH_frac_solo, aux, From 35fdb06370e9ac61ec270a7e3357215aa0be6e2d Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Sun, 22 May 2022 11:11:27 -0600 Subject: [PATCH 03/14] Change working precision in physics/flake.F90 to support compiling physics in single precision --- physics/flake.F90 | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/physics/flake.F90 b/physics/flake.F90 index 2c2e7218c..78eb82632 100644 --- a/physics/flake.F90 +++ b/physics/flake.F90 @@ -36,6 +36,7 @@ MODULE data_parameters ! Description: ! Global parameters for the program are defined. ! + USE, INTRINSIC :: iso_fortran_env IMPLICIT NONE @@ -44,11 +45,11 @@ MODULE data_parameters ! Parameters for the Program: INTEGER, PARAMETER :: & - ireals = SELECTED_REAL_KIND (12,200), & - ! number of desired significant digits for - ! real variables - ! corresponds to 8 byte real variables - +#ifdef SINGLE_PREC + ireals = REAL32, & +#else + ireals = REAL64, & +#endif iintegers = KIND (1) ! kind-type parameter of the integer values ! corresponds to the default integers From a35dfda1762ac4dcc37957a15b076104599a4d87 Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Sun, 22 May 2022 11:11:54 -0600 Subject: [PATCH 04/14] Interface changes in physics/rrtmgp_lw_cloud_sampling.F90 and physics/rrtmgp_sw_cloud_sampling.F90 to support single precision builds --- physics/rrtmgp_lw_cloud_sampling.F90 | 20 ++++++++++---------- physics/rrtmgp_sw_cloud_sampling.F90 | 18 +++++++++--------- 2 files changed, 19 insertions(+), 19 deletions(-) diff --git a/physics/rrtmgp_lw_cloud_sampling.F90 b/physics/rrtmgp_lw_cloud_sampling.F90 index cb11607dc..cf43bb184 100644 --- a/physics/rrtmgp_lw_cloud_sampling.F90 +++ b/physics/rrtmgp_lw_cloud_sampling.F90 @@ -1,9 +1,9 @@ module rrtmgp_lw_cloud_sampling - use machine, only: kind_phys + use machine, only: kind_phys, kind_dbl_prec use mo_gas_optics_rrtmgp, only: ty_gas_optics_rrtmgp use mo_optical_props, only: ty_optical_props_2str use rrtmgp_sampling, only: sampled_mask, draw_samples - use mersenne_twister, only: random_setseed, random_number, random_stat + use mersenne_twister, only: random_setseed, random_number, random_stat use radiation_tools, only: check_error_msg use rrtmgp_lw_gas_optics, only: lw_gas_props use netcdf @@ -75,9 +75,9 @@ subroutine rrtmgp_lw_cloud_sampling_run(doLWrad, nCol, nLev, icseed_lw, iovr,iov integer :: iCol, iLay, iBand integer,dimension(ncol) :: ipseed_lw type(random_stat) :: rng_stat - real(kind_phys), dimension(lw_gas_props%get_ngpt(),nLev,ncol) :: rng3D,rng3D2 - real(kind_phys), dimension(lw_gas_props%get_ngpt()*nLev) :: rng2D - real(kind_phys), dimension(lw_gas_props%get_ngpt()) :: rng1D + real(kind_dbl_prec), dimension(lw_gas_props%get_ngpt(),nLev,ncol) :: rng3D,rng3D2 + real(kind_dbl_prec), dimension(lw_gas_props%get_ngpt()*nLev) :: rng2D + real(kind_dbl_prec), dimension(lw_gas_props%get_ngpt()) :: rng1D logical, dimension(ncol,nLev,lw_gas_props%get_ngpt()) :: maskMCICA ! Initialize CCPP error handling variables @@ -127,7 +127,7 @@ subroutine rrtmgp_lw_cloud_sampling_run(doLWrad, nCol, nLev, icseed_lw, iovr,iov ! Cloud-overlap. ! Maximum-random, random or maximum. if (iovr == iovr_maxrand .or. iovr == iovr_rand .or. iovr == iovr_max) then - call sampled_mask(rng3D, cld_frac, maskMCICA) + call sampled_mask(real(rng3D, kind=kind_phys), cld_frac, maskMCICA) endif ! Exponential decorrelation length overlap if (iovr == iovr_dcorr) then @@ -137,14 +137,14 @@ subroutine rrtmgp_lw_cloud_sampling_run(doLWrad, nCol, nLev, icseed_lw, iovr,iov call random_number(rng2D,rng_stat) rng3D2(:,:,iCol) = reshape(source = rng2D,shape=[lw_gas_props%get_ngpt(),nLev]) enddo - call sampled_mask(rng3D, cld_frac, maskMCICA, & + call sampled_mask(real(rng3D, kind=kind_phys), cld_frac, maskMCICA, & overlap_param = cloud_overlap_param(:,1:nLev-1), & - randoms2 = rng3D2) + randoms2 = real(rng3D2, kind=kind_phys)) endif ! Exponential or Exponential-random if (iovr == iovr_exp .or. iovr == iovr_exprand) then - call sampled_mask(rng3D, cld_frac, maskMCICA, & - overlap_param = cloud_overlap_param(:,1:nLev-1)) + call sampled_mask(real(rng3D, kind=kind_phys), cld_frac, maskMCICA, & + overlap_param = cloud_overlap_param(:,1:nLev-1)) endif ! diff --git a/physics/rrtmgp_sw_cloud_sampling.F90 b/physics/rrtmgp_sw_cloud_sampling.F90 index c4a5de4c8..a5b090149 100644 --- a/physics/rrtmgp_sw_cloud_sampling.F90 +++ b/physics/rrtmgp_sw_cloud_sampling.F90 @@ -1,5 +1,5 @@ module rrtmgp_sw_cloud_sampling - use machine, only: kind_phys + use machine, only: kind_phys, kind_dbl_prec use mo_gas_optics_rrtmgp, only: ty_gas_optics_rrtmgp use mo_optical_props, only: ty_optical_props_2str use rrtmgp_sampling, only: sampled_mask, draw_samples @@ -80,9 +80,9 @@ subroutine rrtmgp_sw_cloud_sampling_run(doSWrad, nCol, nDay, nLev, idxday, iovr, integer,dimension(nday) :: ipseed_sw type(random_stat) :: rng_stat real(kind_phys) :: tauloc,asyloc,ssaloc - real(kind_phys), dimension(sw_gas_props%get_ngpt(),nLev,nday) :: rng3D,rng3D2 - real(kind_phys), dimension(sw_gas_props%get_ngpt()*nLev) :: rng2D - real(kind_phys), dimension(sw_gas_props%get_ngpt()) :: rng1D + real(kind_dbl_prec), dimension(sw_gas_props%get_ngpt(),nLev,nday) :: rng3D,rng3D2 + real(kind_dbl_prec), dimension(sw_gas_props%get_ngpt()*nLev) :: rng2D + real(kind_dbl_prec), dimension(sw_gas_props%get_ngpt()) :: rng1D logical, dimension(nday,nLev,sw_gas_props%get_ngpt()) :: maskMCICA ! Initialize CCPP error handling variables @@ -131,7 +131,7 @@ subroutine rrtmgp_sw_cloud_sampling_run(doSWrad, nCol, nDay, nLev, idxday, iovr, ! Cloud overlap. ! Maximum-random, random, or maximum cloud overlap if (iovr == iovr_maxrand .or. iovr == iovr_max .or. iovr == iovr_rand) then - call sampled_mask(rng3D, cld_frac(idxday(1:nDay),:), maskMCICA) + call sampled_mask(real(rng3D, kind=kind_phys), cld_frac(idxday(1:nDay),:), maskMCICA) endif ! Decorrelation-length overlap if (iovr == iovr_dcorr) then @@ -140,13 +140,13 @@ subroutine rrtmgp_sw_cloud_sampling_run(doSWrad, nCol, nDay, nLev, idxday, iovr, call random_number(rng2D,rng_stat) rng3D2(:,:,iday) = reshape(source = rng2D,shape=[sw_gas_props%get_ngpt(),nLev]) enddo - call sampled_mask(rng3D, cld_frac(idxday(1:nDay),:), maskMCICA, & - overlap_param = cloud_overlap_param(idxday(1:nDay),1:nLev-1),& - randoms2 = rng3D2) + call sampled_mask(real(rng3D, kind=kind_phys), cld_frac(idxday(1:nDay),:), maskMCICA, & + overlap_param = cloud_overlap_param(idxday(1:nDay),1:nLev-1), & + randoms2 = real(rng3D2, kind=kind_phys)) endif ! Exponential or exponential-random cloud overlap if (iovr == iovr_exp .or. iovr == iovr_exprand) then - call sampled_mask(rng3D, cld_frac(idxday(1:nDay),:), maskMCICA, & + call sampled_mask(real(rng3D, kind=kind_phys), cld_frac(idxday(1:nDay),:), maskMCICA, & overlap_param = cloud_overlap_param(idxday(1:nDay),1:nLev-1)) endif From 82250de857dab51b8de4bcc40ec147f9e229a4e6 Mon Sep 17 00:00:00 2001 From: Dusan Jovic Date: Thu, 26 May 2022 18:47:34 +0000 Subject: [PATCH 05/14] Fix argument mismatch in w3nco subroutine calls --- physics/GFS_time_vary_pre.fv3.F90 | 22 ++++++++++++---- physics/GFS_time_vary_pre.scm.F90 | 24 +++++++++++++----- physics/aerinterp.F90 | 6 +++-- physics/cires_tauamf_data.F90 | 3 ++- physics/h2ointerp.f90 | 3 ++- physics/iccninterp.F90 | 3 ++- physics/machine.F | 2 +- physics/ozinterp.f90 | 3 ++- physics/sfcsub.F | 42 ++++++++++++++++--------------- 9 files changed, 70 insertions(+), 38 deletions(-) diff --git a/physics/GFS_time_vary_pre.fv3.F90 b/physics/GFS_time_vary_pre.fv3.F90 index 4b7648c38..86dfe3f40 100644 --- a/physics/GFS_time_vary_pre.fv3.F90 +++ b/physics/GFS_time_vary_pre.fv3.F90 @@ -70,7 +70,7 @@ subroutine GFS_time_vary_pre_timestep_init (jdat, idat, dtp, nsswr, nslwr, nhfrad, idate, debug, me, master, nscyc, sec, phour, zhour, fhour, & kdt, julian, yearlen, ipt, lprnt, lssav, lsswr, lslwr, solhr, errmsg, errflg) - use machine, only: kind_phys + use machine, only: kind_phys, kind_dbl_prec, kind_sngl_prec implicit none @@ -92,8 +92,10 @@ subroutine GFS_time_vary_pre_timestep_init (jdat, idat, dtp, nsswr, real(kind=kind_phys), parameter :: con_24 = 24.0_kind_phys real(kind=kind_phys), parameter :: con_hr = 3600.0_kind_phys - real(kind=kind_phys) :: rinc(5) + real(kind=kind_sngl_prec) :: rinc4(5) + real(kind=kind_dbl_prec) :: rinc8(5) + integer :: w3kindreal,w3kindint integer :: iw3jdn integer :: jd0, jd1 real :: fjd @@ -111,9 +113,19 @@ subroutine GFS_time_vary_pre_timestep_init (jdat, idat, dtp, nsswr, !--- jdat is being updated directly inside of FV3GFS_cap.F90 !--- update calendars and triggers - rinc(1:5) = 0 - call w3difdat(jdat,idat,4,rinc) - sec = rinc(4) + call w3kind(w3kindreal,w3kindint) + if (w3kindreal == 8) then + rinc8(1:5) = 0 + call w3difdat(jdat,idat,4,rinc8) + sec = rinc8(4) + else if (w3kindreal == 4) then + rinc4(1:5) = 0 + call w3difdat(jdat,idat,4,rinc4) + sec = rinc4(4) + else + write(0,*)' FATAL ERROR: Invalid w3kindreal' + call abort + endif phour = sec/con_hr !--- set current bucket hour zhour = phour diff --git a/physics/GFS_time_vary_pre.scm.F90 b/physics/GFS_time_vary_pre.scm.F90 index 0c34ca735..2bb6b3ceb 100644 --- a/physics/GFS_time_vary_pre.scm.F90 +++ b/physics/GFS_time_vary_pre.scm.F90 @@ -69,7 +69,7 @@ subroutine GFS_time_vary_pre_timestep_init (jdat, idat, dtp, nsswr, & nslwr, idate, debug, me, master, nscyc, sec, phour, zhour, fhour, kdt, & julian, yearlen, ipt, lprnt, lssav, lsswr, lslwr, solhr, errmsg, errflg) - use machine, only: kind_phys + use machine, only: kind_phys, kind_dbl_prec, kind_sngl_prec implicit none @@ -91,8 +91,10 @@ subroutine GFS_time_vary_pre_timestep_init (jdat, idat, dtp, nsswr, & real(kind=kind_phys), parameter :: con_24 = 24.0_kind_phys real(kind=kind_phys), parameter :: con_hr = 3600.0_kind_phys - real(kind=kind_phys) :: rinc(5) - + real(kind=kind_sngl_prec) :: rinc4(5) + real(kind=kind_dbl_prec) :: rinc8(5) + + integer :: w3kindreal,w3kindint integer :: iw3jdn integer :: jd0, jd1 real :: fjd @@ -112,9 +114,19 @@ subroutine GFS_time_vary_pre_timestep_init (jdat, idat, dtp, nsswr, & !--- jdat is being updated directly inside of the time integration !--- loop of scm.F90 !--- update calendars and triggers - rinc(1:5) = 0 - call w3difdat(jdat,idat,4,rinc) - sec = rinc(4) + call w3kind(w3kindreal,w3kindint) + if (w3kindreal == 8) then + rinc8(1:5) = 0 + call w3difdat(jdat,idat,4,rinc8) + sec = rinc8(4) + else if (w3kindreal == 4) then + rinc4(1:5) = 0 + call w3difdat(jdat,idat,4,rinc4) + sec = rina4c(4) + else + write(0,*)' FATAL ERROR: Invalid w3kindreal' + call abort + endif phour = sec/con_hr !--- set current bucket hour zhour = phour diff --git a/physics/aerinterp.F90 b/physics/aerinterp.F90 index 30ff97dff..8ad446f30 100644 --- a/physics/aerinterp.F90 +++ b/physics/aerinterp.F90 @@ -110,7 +110,8 @@ SUBROUTINE read_aerdataf ( me, master, iflip, idate, FHOUR, errmsg, errflg) integer :: i, j, k, n, ii, imon, klev, n1, n2 logical :: file_exist integer IDAT(8),JDAT(8) - real(kind=kind_phys) RINC(5), rjday + real(kind=kind_phys) rjday + real(8) RINC(5) integer jdow, jdoy, jday real(4) rinc4(5) integer w3kindreal,w3kindint @@ -244,8 +245,9 @@ SUBROUTINE aerinterpol( me,master,nthrds,npts,IDATE,FHOUR,iflip, jindx1,jindx2, real(kind=kind_phys) aerout(npts,lev,ntrcaer) real(kind=kind_phys) aerpm(npts,levsaer,ntrcaer) real(kind=kind_phys) prsl(npts,lev), aerpres(npts,levsaer) - real(kind=kind_phys) RINC(5), rjday + real(kind=kind_phys) rjday integer jdow, jdoy, jday + real(8) RINC(5) real(4) rinc4(5) integer w3kindreal,w3kindint diff --git a/physics/cires_tauamf_data.F90 b/physics/cires_tauamf_data.F90 index e0d43e74e..77c74b1a9 100644 --- a/physics/cires_tauamf_data.F90 +++ b/physics/cires_tauamf_data.F90 @@ -176,8 +176,9 @@ subroutine gfs_idate_calendar(idate, fhour, ddd, fddd) ! !locals ! - real(kind=kind_phys) :: rinc(5), rjday + real(kind=kind_phys) :: rjday integer :: jdow, jdoy, jday + real(8) :: rinc(5) real(4) :: rinc4(5) integer :: w3kindreal, w3kindint diff --git a/physics/h2ointerp.f90 b/physics/h2ointerp.f90 index f26ae6c0c..c4fb355fc 100644 --- a/physics/h2ointerp.f90 +++ b/physics/h2ointerp.f90 @@ -144,8 +144,9 @@ subroutine h2ointerpol(me,npts,idate,fhour,jindx1,jindx2,h2oplout,ddy) ! real(kind=kind_phys) ddy(npts) real(kind=kind_phys) h2oplout(npts,levh2o,h2o_coeff) - real(kind=kind_phys) rinc(5), rjday + real(kind=kind_phys) rjday integer jdow, jdoy, jday + real(8) rinc(5) real(4) rinc4(5) integer w3kindreal, w3kindint ! diff --git a/physics/iccninterp.F90 b/physics/iccninterp.F90 index a3a08dee8..c08b5c5e5 100644 --- a/physics/iccninterp.F90 +++ b/physics/iccninterp.F90 @@ -143,8 +143,9 @@ SUBROUTINE ciinterpol(me,npts,IDATE,FHOUR,jindx1,jindx2,ddy, & real(kind=kind_phys) ciplout(npts,lev),cipm(npts,kcipl) real(kind=kind_phys) ccnout(npts,lev),ccnpm(npts,kcipl) real(kind=kind_phys) cipres(npts,kcipl), prsl(npts,lev) - real(kind=kind_phys) RINC(5), rjday + real(kind=kind_phys) rjday integer jdow, jdoy, jday + real(8) RINC(5) real(4) rinc4(5) integer w3kindreal,w3kindint ! diff --git a/physics/machine.F b/physics/machine.F index 9b09d235c..ac881bc62 100644 --- a/physics/machine.F +++ b/physics/machine.F @@ -23,7 +23,7 @@ module machine &, kind_INTEGER = 4 ! -,,- #else - integer, parameter :: kind_io4 = 4, kind_io8 = 8 , kind_ior = 8 & + integer, parameter :: kind_io4 = 4, kind_io8 = 4 , kind_ior = 8 & &, kind_evod = 4, kind_dbl_prec = 8 & &, kind_sngl_prec = 4 # ifdef __PGI diff --git a/physics/ozinterp.f90 b/physics/ozinterp.f90 index 6fe86c8e1..5b3149d61 100644 --- a/physics/ozinterp.f90 +++ b/physics/ozinterp.f90 @@ -147,8 +147,9 @@ SUBROUTINE ozinterpol(me,npts,IDATE,FHOUR,jindx1,jindx2,ozplout,ddy) ! real(kind=kind_phys) DDY(npts) real(kind=kind_phys) ozplout(npts,levozp,oz_coeff) - real(kind=kind_phys) RINC(5), rjday + real(kind=kind_phys) rjday integer jdow, jdoy, jday + real(8) rinc(5) real(4) rinc4(5) integer w3kindreal,w3kindint ! diff --git a/physics/sfcsub.F b/physics/sfcsub.F index 78e5201be..dae710760 100644 --- a/physics/sfcsub.F +++ b/physics/sfcsub.F @@ -2736,7 +2736,7 @@ subroutine hmskrd(lugb,imsk,jmsk,fnmskh, & !>\ingroup mod_sfcsub subroutine fixrdg(lugb,idim,jdim,fngrib, & & kpds5,gdata,gaus,blno,blto,me) - use machine , only : kind_io8,kind_io4 + use machine , only : kind_io8,kind_dbl_prec,kind_sngl_prec use sfccyc_module, only : mdata implicit none integer lgrib,n,lskip,jret,j,ndata,lugi,jdim,idim,lugb, @@ -2747,8 +2747,8 @@ subroutine fixrdg(lugb,idim,jdim,fngrib, & real (kind=kind_io8) gdata(idim*jdim) logical gaus real (kind=kind_io8) blno,blto - real (kind=kind_io8), allocatable :: data8(:) - real (kind=kind_io4), allocatable :: data4(:) + real (kind=kind_dbl_prec), allocatable :: data8(:) + real (kind=kind_sngl_prec), allocatable :: data4(:) ! logical*1, allocatable :: lbms(:) ! @@ -2815,7 +2815,7 @@ subroutine fixrdg(lugb,idim,jdim,fngrib, & allocate(data4(1:idim*jdim)) call getgb(lugb,lugi,kdata,lskip,jpds,jgds,ndata,lskip, & kpds,kgds,lbms,data4,jret) - data8 = real(data4, kind=kind_io8) + data8(1:ndata) = real(data4(1:ndata), kind=kind_dbl_prec) deallocate(data4) else write(0,*)' FATAL ERROR: Invalid w3kindreal' @@ -6165,7 +6165,7 @@ subroutine qcbyfc(tsffcs,snofcs,qctsfs,qcsnos,qctsfi, & subroutine setrmsk(kpds5,slmask,igaul,jgaul,wlon,rnlat, & & data,imax,jmax,rlnout,rltout,lmask,rslmsk & &, gaus,blno, blto, kgds1, kpds4, lbms) - use machine , only : kind_io8,kind_io4 + use machine , only : kind_io8,kind_io4,kind_dbl_prec use sfccyc_module implicit none real (kind=kind_io8) blno,blto,wlon,rnlat,crit,data_max @@ -6177,7 +6177,8 @@ subroutine setrmsk(kpds5,slmask,igaul,jgaul,wlon,rnlat, & real (kind=kind_io8) slmask(igaul,jgaul) real (kind=kind_io8) data(imax,jmax),rslmsk(imax,jmax) &, rlnout(imax), rltout(jmax) - real (kind=kind_io8) a(jmax), w(jmax), radi, dlat, dlon + real (kind=kind_io8) radi, dlat, dlon + real (kind=kind_dbl_prec) a(jmax), w(jmax) logical lmask, gaus ! ! set the longitude and latitudes for the grib file @@ -6650,7 +6651,7 @@ subroutine setrmsk(kpds5,slmask,igaul,jgaul,wlon,rnlat, & !! This subroutine interpolates from lat/lon grid to other lat/lon grid. subroutine ga2la(gauin,imxin,jmxin,regout,imxout,jmxout, & & wlon,rnlat,rlnout,rltout,gaus,blno, blto) - use machine , only : kind_io8,kind_io4 + use machine , only : kind_io8,kind_io4,kind_dbl_prec use sfccyc_module , only : num_threads implicit none integer i1,i2,j2,ishft,i,jj,j1,jtem,jmxout,imxin,jmxin,imxout, & @@ -6673,7 +6674,7 @@ subroutine ga2la(gauin,imxin,jmxin,regout,imxout,jmxout, & data jmxsav/0/ save jmxsav, gaul, dlati real (kind=kind_io8) radi - real (kind=kind_io8) a(jmxin), w(jmxin) + real (kind=kind_dbl_prec) a(jmxin), w(jmxin) ! ! logical first @@ -6917,11 +6918,12 @@ end subroutine landtyp !>\ingroup mod_sfcsub subroutine gaulat(gaul,k) ! - use machine , only : kind_io8,kind_io4 + use machine , only : kind_io8,kind_io4,kind_dbl_prec implicit none integer n,k real (kind=kind_io8) radi - real (kind=kind_io8) a(k), w(k), gaul(k) + real (kind=kind_io8) gaul(k) + real (kind=kind_dbl_prec) a(k), w(k) ! call splat(4, k, a, w) ! @@ -7040,7 +7042,7 @@ subroutine clima(lugb,iy,im,id,ih,fh,len,lsoil,slmskl,slmskw, & data dayhf/ 15.5, 45.0, 74.5,105.0,135.5,166.0, & 196.5,227.5,258.0,288.5,319.0,349.5,380.5/ ! - real (kind=kind_io8) fha(5) + real(8) fha(5) real(4) fha4(5) integer w3kindreal,w3kindint integer ida(8),jda(8),ivtyp, kpd7 @@ -8305,7 +8307,7 @@ subroutine fixrdc(lugb,fngrib,kpds5,kpds7,mon,slmask, & & gdata,len,iret & &, imsk, jmsk, slmskh, gaus,blno, blto & &, outlat, outlon, me) - use machine , only : kind_io8,kind_io4 + use machine , only : kind_io8,kind_dbl_prec,kind_sngl_prec use sfccyc_module, only : mdata implicit none integer imax,jmax,ijmax,i,j,n,jret,inttyp,iret,imsk, & @@ -8321,8 +8323,8 @@ subroutine fixrdc(lugb,fngrib,kpds5,kpds7,mon,slmask, & ! real (kind=kind_io8) gdata(len), slmask(len) real (kind=kind_io8), allocatable :: data(:,:), rslmsk(:,:) - real (kind=kind_io8), allocatable :: data8(:) - real (kind=kind_io4), allocatable :: data4(:) + real (kind=kind_dbl_prec), allocatable :: data8(:) + real (kind=kind_sngl_prec), allocatable :: data4(:) real (kind=kind_io8), allocatable :: rlngrb(:), rltgrb(:) ! logical lmask, yr2kc, gaus, ijordr @@ -8398,7 +8400,7 @@ subroutine fixrdc(lugb,fngrib,kpds5,kpds7,mon,slmask, & allocate(data4(1:mdata)) call getgb(lugb,lugi,mdata,lskip,jpds,jgds,ndata,lskip, & kpds,kgds,lbms,data4,jret) - data8 = real(data4, kind=kind_io8) + data8(1:ndata) = real(data4(1:ndata), kind=kind_dbl_prec) deallocate(data4) endif if (me .eq. 0) write(6,*) ' input grib file dates=', @@ -8479,7 +8481,7 @@ subroutine fixrda(lugb,fngrib,kpds5,slmask, & & iy,im,id,ih,fh,gdata,len,iret & &, imsk, jmsk, slmskh, gaus,blno, blto & &, outlat, outlon, me) - use machine , only : kind_io8,kind_io4 + use machine , only : kind_io8,kind_dbl_prec,kind_sngl_prec use sfccyc_module, only : mdata implicit none integer nrepmx,nvalid,imo,iyr,idy,jret,ihr,nrept,lskip,lugi, & @@ -8506,8 +8508,8 @@ subroutine fixrda(lugb,fngrib,kpds5,slmask, & ! real (kind=kind_io8) gdata(len), slmask(len) real (kind=kind_io8), allocatable :: data(:,:),rslmsk(:,:) - real (kind=kind_io8), allocatable :: data8(:) - real (kind=kind_io4), allocatable :: data4(:) + real (kind=kind_dbl_prec), allocatable :: data8(:) + real (kind=kind_sngl_prec), allocatable :: data4(:) real (kind=kind_io8), allocatable :: rlngrb(:), rltgrb(:) ! logical lmask, yr2kc, gaus, ijordr @@ -8528,7 +8530,7 @@ subroutine fixrda(lugb,fngrib,kpds5,slmask, & integer mjday(12) data mjday/31,28,31,30,31,30,31,31,30,31,30,31/ ! - real (kind=kind_io8) fha(5) + real(8) fha(5) real(4) fha4(5) integer ida(8),jda(8) ! @@ -8645,7 +8647,7 @@ subroutine fixrda(lugb,fngrib,kpds5,slmask, & allocate (data4(1:mdata)) call getgb(lugb,lugi,mdata,lskip,jpds,jgds,ndata,lskip, & kpds,kgds,lbms,data4,jret) - data8 = real(data4, kind=kind_io8) + data8(1:ndata) = real(data4(1:ndata), kind=kind_dbl_prec) deallocate(data4) endif if (me .eq. 0) write(6,*) ' input grib file dates=', From d6978a2f05c41aa001607997c671bccff699f894 Mon Sep 17 00:00:00 2001 From: "samuel.trahan" Date: Wed, 22 Jun 2022 16:36:22 +0000 Subject: [PATCH 06/14] Use max(...,1e-7) in 1-exner comparision for single precision --- physics/GFS_suite_interstitial_3.F90 | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/physics/GFS_suite_interstitial_3.F90 b/physics/GFS_suite_interstitial_3.F90 index 79ab481ec..2fb3190b9 100644 --- a/physics/GFS_suite_interstitial_3.F90 +++ b/physics/GFS_suite_interstitial_3.F90 @@ -121,11 +121,19 @@ subroutine GFS_suite_interstitial_3_run (otsptflag, & do k=1,levs do i=1,im kk = max(10,kpbl(i)) +#ifdef SINGLE_PREC + if (k < kk) then + tem = rhcbot - (rhcbot-rhcpbl) * (one-prslk(i,k)) / max(one-prslk(i,kk),1e-7) + else + tem = rhcpbl - (rhcpbl-rhctop) * (prslk(i,kk)-prslk(i,k)) / max(prslk(i,kk),1e-7) + endif +#else if (k < kk) then tem = rhcbot - (rhcbot-rhcpbl) * (one-prslk(i,k)) / (one-prslk(i,kk)) else tem = rhcpbl - (rhcpbl-rhctop) * (prslk(i,kk)-prslk(i,k)) / prslk(i,kk) endif +#endif tem = rhcmax * work1(i) + tem * work2(i) rhc(i,k) = max(zero, min(one,tem)) enddo @@ -192,4 +200,4 @@ subroutine GFS_suite_interstitial_3_run (otsptflag, & end subroutine GFS_suite_interstitial_3_run - end module GFS_suite_interstitial_3 \ No newline at end of file + end module GFS_suite_interstitial_3 From a0f46cbabf917dc9fff107f60fd737bcc4828f9f Mon Sep 17 00:00:00 2001 From: "samuel.trahan" Date: Wed, 22 Jun 2022 16:36:37 +0000 Subject: [PATCH 07/14] flake gets its types from machine.F --- physics/flake.F90 | 18 ++---------------- 1 file changed, 2 insertions(+), 16 deletions(-) diff --git a/physics/flake.F90 b/physics/flake.F90 index 2c2e7218c..557e22949 100644 --- a/physics/flake.F90 +++ b/physics/flake.F90 @@ -35,23 +35,9 @@ MODULE data_parameters ! ! Description: ! Global parameters for the program are defined. +! Actually, scratch that. We'll import them from machine.F instead. ! - -IMPLICIT NONE - -!======================================================================= -! Global (i.e. public) Declarations: -! Parameters for the Program: - - INTEGER, PARAMETER :: & - ireals = SELECTED_REAL_KIND (12,200), & - ! number of desired significant digits for - ! real variables - ! corresponds to 8 byte real variables - - iintegers = KIND (1) - ! kind-type parameter of the integer values - ! corresponds to the default integers + use machine, only: ireals=>kind_phys, iintegers=>kind_INTEGER !======================================================================= From 2c720da4241739c534b46b61de984f7e2bf14c4f Mon Sep 17 00:00:00 2001 From: "samuel.trahan" Date: Wed, 22 Jun 2022 16:36:55 +0000 Subject: [PATCH 08/14] module_mp_fer_hires gets its types from machine.F --- physics/module_MP_FER_HIRES.F90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/physics/module_MP_FER_HIRES.F90 b/physics/module_MP_FER_HIRES.F90 index 092a2f941..380ba9de5 100644 --- a/physics/module_MP_FER_HIRES.F90 +++ b/physics/module_MP_FER_HIRES.F90 @@ -2249,9 +2249,10 @@ REAL FUNCTION DEPOSIT (PP,Tdum,WVdum,RHgrd,I,J,L) !-- Debug 20120111 !--- Also uses the Asai (1965) algorithm, but uses a different target ! vapor pressure for the adjustment ! + use machine, only: HIGH_PRES => kind_dbl_prec IMPLICIT NONE ! - INTEGER, PARAMETER :: HIGH_PRES=Selected_Real_Kind(15) + !INTEGER, PARAMETER :: HIGH_PRES=Selected_Real_Kind(15) REAL (KIND=HIGH_PRES), PARAMETER :: RHLIMIT=.001, & & RHLIMIT1=-RHLIMIT REAL (KIND=HIGH_PRES) :: DEP, SSAT From a9c97bba61311deef45974a2b0a60f8092fea263 Mon Sep 17 00:00:00 2001 From: "samuel.trahan" Date: Wed, 22 Jun 2022 18:34:04 +0000 Subject: [PATCH 09/14] Clean up machine.F and get rid of kind_evod & kind_rad --- physics/date_def.f | 10 ++++---- physics/machine.F | 42 ++++++++++++++++++------------- physics/module_nst_parameters.f90 | 3 +-- 3 files changed, 30 insertions(+), 25 deletions(-) diff --git a/physics/date_def.f b/physics/date_def.f index 2907d7416..fceb4334f 100644 --- a/physics/date_def.f +++ b/physics/date_def.f @@ -1,13 +1,13 @@ module date_def - use machine, ONLY: kind_evod + use machine, ONLY: kind_phys implicit none !jw integer idate(4) -!jw real(kind=kind_evod) fhour,shour,thour,z00 - real(kind=kind_evod) shour,thour,z00 - real(kind=kind_evod),target :: fhour, zhour +!jw real(kind=kind_phys) fhour,shour,thour,z00 + real(kind=kind_phys) shour,thour,z00 + real(kind=kind_phys),target :: fhour, zhour integer,target :: idate(4),idate7(7) ! - REAL(KIND=KIND_EVOD) ,ALLOCATABLE :: spdmax(:) + REAL(KIND=KIND_PHYS) ,ALLOCATABLE :: spdmax(:) end module date_def diff --git a/physics/machine.F b/physics/machine.F index eb1dcd257..0269db53c 100644 --- a/physics/machine.F +++ b/physics/machine.F @@ -5,32 +5,38 @@ module machine !! implicit none - - integer, parameter :: kind_io4 = 4, kind_io8 = 8 , kind_ior = 8 & - &, kind_evod = 8, kind_dbl_prec = 8 & - &, kind_sngl_prec = 4, kind_INTEGER = 4 & - &, kind_LOGICAL = 4 -#ifdef SINGLE_PREC - integer, parameter :: kind_rad = kind_sngl_prec & - &, kind_phys = kind_sngl_prec & - &, kind_grid = kind_dbl_prec &! atmos_cubed_sphere requres kind_grid=8 - &, kind_REAL = kind_sngl_prec ! used in cmp_comm + integer, parameter :: kind_sngl_prec = 4 & + &, kind_dbl_prec = 8 & +# ifdef __PGI + &, kind_qdt_prec = 8 & +# else + &, kind_qdt_prec = 16 & +# endif + &, kind_integer = 4 & + &, kind_logical = 4 + &, kind_io4 = kind_sngl_prec & + &, kind_ior = kind_dbl_prec & + &, kind_grid = kind_dbl_prec + +! Physics single precision flag +#ifndef SINGLE_PREC + integer, parameter :: kind_phys = kind_dbl_prec #else - integer, parameter :: kind_rad = kind_dbl_prec & - &, kind_phys = kind_dbl_prec & - &, kind_grid = kind_dbl_prec &! atmos_cubed_sphere requres kind_grid=8 - &, kind_REAL = kind_dbl_prec ! used in cmp_comm + integer, parameter :: kind_phys = kind_sngl_prec #endif + integer, parameter :: kind_io8 = kind_phys &! Note kind_io8 is not always 8 bytes + +! Dynamics single precision flag #ifdef OVERLOAD_R4 - integer, parameter :: kind_dyn = 4 + integer, parameter :: kind_dyn = kind_sngl_prec #else - integer, parameter :: kind_dyn = 8 + integer, parameter :: kind_dyn = kind_dbl_prec #endif ! - real(kind=kind_evod), parameter :: mprec = 1.e-12 ! machine precision to restrict dep - real(kind=kind_evod), parameter :: grib_undef = 9.99e20 ! grib undefine value + real(kind=kind_phys), parameter :: mprec = 1.e-12 ! machine precision to restrict dep + real(kind=kind_phys), parameter :: grib_undef = 9.99e20 ! grib undefine value ! end module machine diff --git a/physics/module_nst_parameters.f90 b/physics/module_nst_parameters.f90 index 6d8c8794b..f877bd520 100644 --- a/physics/module_nst_parameters.f90 +++ b/physics/module_nst_parameters.f90 @@ -9,8 +9,7 @@ !! history: !! 20210305: X.Li, reduce z_w_max from 30 m to 20 m module module_nst_parameters - use machine, only : kind_phys & - ,kind_rad ! for astronomy (date) calculations + use machine, only : kind_phys ! ! air constants and coefficients from the atmospehric model use physcons, only: & From 78caef910104df9c4ee44ab010d67e19f2322a2d Mon Sep 17 00:00:00 2001 From: "samuel.trahan" Date: Thu, 23 Jun 2022 14:47:14 +0000 Subject: [PATCH 10/14] point to sam's rte-rrtmgp --- .gitmodules | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.gitmodules b/.gitmodules index dc5a31dd7..835fd5c17 100644 --- a/.gitmodules +++ b/.gitmodules @@ -2,5 +2,5 @@ path = physics/rte-rrtmgp #url = https://github.com/earth-system-radiation/rte-rrtmgp #branch = dtc/ccpp - url = https://github.com/climbfuji/rte-rrtmgp - branch = feature/single_precision + url = https://github.com/SamuelTrahanNOAA/rte-rrtmgp + branch = sing_prec_from_main From 855bf8686642361bf1f2308916297e165aa3c4b0 Mon Sep 17 00:00:00 2001 From: "samuel.trahan" Date: Thu, 23 Jun 2022 19:58:18 +0000 Subject: [PATCH 11/14] revert physics/aer_cloud.F to working_32bit version --- physics/aer_cloud.F | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/physics/aer_cloud.F b/physics/aer_cloud.F index e8218cbb7..a334428d1 100644 --- a/physics/aer_cloud.F +++ b/physics/aer_cloud.F @@ -3476,10 +3476,9 @@ SUBROUTINE EMPIRICAL_PARAM_PHILLIPS(SI, SIW, SW, D_grid_dust, & n_grid_dust, ijstop_dust, D_grid_soot, n_grid_soot, ijstop_soot, & D_grid_bio, n_grid_bio, ijstop_bio, A_solo, n_iw, DSH, & Nhet_dep,Nhet_dhf,fdust_dep,P_ice, T_ice) - use machine, only: kind_dbl_prec implicit none real, intent(IN):: SI, SIW, SW, A_solo - real(kind_dbl_prec), intent(IN):: P_ice, T_ice + real*8, intent(IN):: P_ice, T_ice real, dimension(:), intent(IN):: D_grid_dust, n_grid_dust, & D_grid_soot, n_grid_soot, D_grid_bio, n_grid_bio integer, intent(IN):: ijstop_dust, ijstop_soot, ijstop_bio @@ -3490,7 +3489,7 @@ SUBROUTINE EMPIRICAL_PARAM_PHILLIPS(SI, SIW, SW, D_grid_dust, & num_ic_solo_imm real, intent (inout) :: DSH, n_iw - real(kind_dbl_prec), intent (out) :: Nhet_dep,Nhet_dhf,fdust_dep + real*8, intent (out) :: Nhet_dep,Nhet_dhf,fdust_dep real :: dn_in_dust, dn_in_soot, dn_in_bio, dn_in_solo, dNall, & dNaux, naux, SS_w, dH_frac_dust, dH_frac_soot, dH_frac_solo, aux, From d7a244fccc62afced1d9a8cf761c44e025058f24 Mon Sep 17 00:00:00 2001 From: "samuel.trahan" Date: Mon, 27 Jun 2022 17:55:14 +0000 Subject: [PATCH 12/14] point to authoritative repository for rte-rrtmgp --- .gitmodules | 6 ++---- physics/rte-rrtmgp | 2 +- 2 files changed, 3 insertions(+), 5 deletions(-) diff --git a/.gitmodules b/.gitmodules index 835fd5c17..75e5ea836 100644 --- a/.gitmodules +++ b/.gitmodules @@ -1,6 +1,4 @@ [submodule "physics/rte-rrtmgp"] path = physics/rte-rrtmgp - #url = https://github.com/earth-system-radiation/rte-rrtmgp - #branch = dtc/ccpp - url = https://github.com/SamuelTrahanNOAA/rte-rrtmgp - branch = sing_prec_from_main + url = https://github.com/earth-system-radiation/rte-rrtmgp + branch = dtc/ccpp diff --git a/physics/rte-rrtmgp b/physics/rte-rrtmgp index c848637d2..7f01618c9 160000 --- a/physics/rte-rrtmgp +++ b/physics/rte-rrtmgp @@ -1 +1 @@ -Subproject commit c848637d2158b7a4947ce3f97d8eda4b28b4c556 +Subproject commit 7f01618c92409658bddd3afa9acb004c608f6a0d From 70cdc31c664c5749b1088d9bff835ab332ab7bd8 Mon Sep 17 00:00:00 2001 From: "samuel.trahan" Date: Fri, 1 Jul 2022 14:33:35 +0000 Subject: [PATCH 13/14] correct a type mismatch in a call in module_sf_noahmplsm --- physics/module_sf_noahmplsm.f90 | 13 ++++++++++--- 1 file changed, 10 insertions(+), 3 deletions(-) diff --git a/physics/module_sf_noahmplsm.f90 b/physics/module_sf_noahmplsm.f90 index 3ea8203ab..19f0be523 100644 --- a/physics/module_sf_noahmplsm.f90 +++ b/physics/module_sf_noahmplsm.f90 @@ -10823,6 +10823,7 @@ real*8 function zolri(ri,za,z0,zt,zol1,psi_opt) real (kind=kind_phys) :: x1,x2,fx1,fx2 integer :: n integer, parameter :: nmax = 20 + real(kind=kind_phys) rphys_temp !real, dimension(nmax):: zlhux ! real :: zolri2 @@ -10855,7 +10856,9 @@ real*8 function zolri(ri,za,z0,zt,zol1,psi_opt) if (n==nmax .and. abs(x1 - x2) >= 0.01) then !if convergence fails, use approximate values: - call li_etal_2010(zolri, ri, za/z0, z0/zt) + rphys_temp = zolri + call li_etal_2010(rphys_temp, ri, za/z0, z0/zt) + zolri = rphys_temp !zlhux(n)=zolri !print*,"iter fail, n=",n," ri=",ri," z0=",z0 else @@ -10921,7 +10924,7 @@ real*8 function zolrib(ri,za,z0,zt,logz0,logzt,zol1,psi_opt) integer :: n integer, parameter :: nmax = 20 real (kind=kind_phys), dimension(nmax):: zlhux - real (kind=kind_phys) :: psit2,psix2 + real (kind=kind_phys) :: psit2,psix2,phys_temp ! real :: psim_unstable, psim_stable ! real :: psih_unstable, psih_stable @@ -10972,14 +10975,18 @@ real*8 function zolrib(ri,za,z0,zt,logz0,logzt,zol1,psi_opt) if (n==nmax .and. abs(zolold - zolrib) > 0.01 ) then !print*,"iter fail, n=",n," ri=",ri," z/l=",zolri !if convergence fails, use approximate values: - call li_etal_2010(zolrib, ri, za/z0, z0/zt) + phys_temp = zolrib + call li_etal_2010(phys_temp, ri, za/z0, z0/zt) + zolrib = phys_temp zlhux(n)=zolrib !print*,"failed, n=",n," ri=",ri," z0=",z0 !print*,"z/l=",zlhux(1:nmax) else !if(zolrib*ri .lt. 0.) then ! !print*,"end: wrong quadrants: z/l=",zolrib," ri=",ri + ! !phys_temp = zolrib ! !call li_etal_2010(zolrib, ri, za/z0, z0/zt) + ! !zolrib = phys_temp !endif !print*,"success,n=",n," ri=",ri," z0=",z0 endif From 18e35c6bda63d292ed441c62dbc8f8831120dd16 Mon Sep 17 00:00:00 2001 From: "samuel.trahan" Date: Tue, 5 Jul 2022 16:48:11 +0000 Subject: [PATCH 14/14] better names for new variables in physics/module_sf_noahmplsm.f90 --- physics/module_sf_noahmplsm.f90 | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/physics/module_sf_noahmplsm.f90 b/physics/module_sf_noahmplsm.f90 index 19f0be523..652db602d 100644 --- a/physics/module_sf_noahmplsm.f90 +++ b/physics/module_sf_noahmplsm.f90 @@ -10823,7 +10823,7 @@ real*8 function zolri(ri,za,z0,zt,zol1,psi_opt) real (kind=kind_phys) :: x1,x2,fx1,fx2 integer :: n integer, parameter :: nmax = 20 - real(kind=kind_phys) rphys_temp + real(kind=kind_phys) zolri_iteration !real, dimension(nmax):: zlhux ! real :: zolri2 @@ -10856,9 +10856,9 @@ real*8 function zolri(ri,za,z0,zt,zol1,psi_opt) if (n==nmax .and. abs(x1 - x2) >= 0.01) then !if convergence fails, use approximate values: - rphys_temp = zolri - call li_etal_2010(rphys_temp, ri, za/z0, z0/zt) - zolri = rphys_temp + zolri_iteration= zolri + call li_etal_2010(zolri_iteration, ri, za/z0, z0/zt) + zolri = zolri_iteration !zlhux(n)=zolri !print*,"iter fail, n=",n," ri=",ri," z0=",z0 else @@ -10924,7 +10924,7 @@ real*8 function zolrib(ri,za,z0,zt,logz0,logzt,zol1,psi_opt) integer :: n integer, parameter :: nmax = 20 real (kind=kind_phys), dimension(nmax):: zlhux - real (kind=kind_phys) :: psit2,psix2,phys_temp + real (kind=kind_phys) :: psit2,psix2,zolrib_iteration ! real :: psim_unstable, psim_stable ! real :: psih_unstable, psih_stable @@ -10975,9 +10975,9 @@ real*8 function zolrib(ri,za,z0,zt,logz0,logzt,zol1,psi_opt) if (n==nmax .and. abs(zolold - zolrib) > 0.01 ) then !print*,"iter fail, n=",n," ri=",ri," z/l=",zolri !if convergence fails, use approximate values: - phys_temp = zolrib - call li_etal_2010(phys_temp, ri, za/z0, z0/zt) - zolrib = phys_temp + zolrib_iteration = zolrib + call li_etal_2010(zolrib_iteration, ri, za/z0, z0/zt) + zolrib = zolrib_iteration zlhux(n)=zolrib !print*,"failed, n=",n," ri=",ri," z0=",z0 !print*,"z/l=",zlhux(1:nmax)