From 84fb687b2a9e423b6f0ad3dbf5ffacbd2769da49 Mon Sep 17 00:00:00 2001 From: MinsukJi-NOAA Date: Fri, 8 Oct 2021 15:56:30 -0400 Subject: [PATCH 1/8] Modify set_tracer_profile in tracer_manager.F90 --- tracer_manager/tracer_manager.F90 | 50 ++++++++++++++++++++++++------- 1 file changed, 39 insertions(+), 11 deletions(-) diff --git a/tracer_manager/tracer_manager.F90 b/tracer_manager/tracer_manager.F90 index 79ea8ac623..afb510cef5 100644 --- a/tracer_manager/tracer_manager.F90 +++ b/tracer_manager/tracer_manager.F90 @@ -76,6 +76,8 @@ module tracer_manager_mod fm_exists, & MODEL_NAMES +use platform_mod + implicit none private @@ -1037,7 +1039,7 @@ subroutine set_tracer_profile(model, n, tracer, err_msg) integer, intent(in) :: model !< Parameter representing component model in use integer, intent(in) :: n !< Tracer number -real, intent(inout), dimension(:,:,:) :: tracer !< Initialized tracer array +class(*), intent(inout), dimension(:,:,:) :: tracer !< Initialized tracer array character(len=*), intent(out), optional :: err_msg real :: surf_value, multiplier @@ -1063,7 +1065,12 @@ subroutine set_tracer_profile(model, n, tracer, err_msg) bottom_value = surf_value multiplier = 1.0 -tracer = surf_value +select type (tracer) +type is (real(kind=r4_kind)) + tracer = surf_value +type is (real(kind=r8_kind)) + tracer = surf_value +end select if ( query_method ( 'profile_type',model,n,scheme,control)) then !Change the tracer_number to the tracer_manager version @@ -1072,7 +1079,12 @@ subroutine set_tracer_profile(model, n, tracer, err_msg) profile_type = 'Fixed' flag =parse(control,'surface_value',surf_value) multiplier = 1.0 - tracer = surf_value + select type (tracer) + type is (real(kind=r4_kind)) + tracer = surf_value + type is (real(kind=r8_kind)) + tracer = surf_value + end select endif if(lowercase(trim(scheme(1:7))).eq.'profile') then @@ -1105,16 +1117,32 @@ subroutine set_tracer_profile(model, n, tracer, err_msg) select case (tracers(n1)%model) case (MODEL_ATMOS) multiplier = exp( log (top_value/surf_value) /numlevels) - tracer(:,:,1) = surf_value - do k = 2, size(tracer,3) - tracer(:,:,k) = tracer(:,:,k-1) * multiplier - enddo + select type (tracer) + type is (real(kind=r4_kind)) + tracer(:,:,1) = surf_value + do k = 2, size(tracer,3) + tracer(:,:,k) = tracer(:,:,k-1) * multiplier + enddo + type is (real(kind=r8_kind)) + tracer(:,:,1) = surf_value + do k = 2, size(tracer,3) + tracer(:,:,k) = tracer(:,:,k-1) * multiplier + enddo + end select case (MODEL_OCEAN) multiplier = exp( log (bottom_value/surf_value) /numlevels) - tracer(:,:,size(tracer,3)) = surf_value - do k = size(tracer,3) - 1, 1, -1 - tracer(:,:,k) = tracer(:,:,k+1) * multiplier - enddo + select type (tracer) + type is (real(kind=r4_kind)) + tracer(:,:,size(tracer,3)) = surf_value + do k = size(tracer,3) - 1, 1, -1 + tracer(:,:,k) = tracer(:,:,k+1) * multiplier + enddo + type is (real(kind=r8_kind)) + tracer(:,:,size(tracer,3)) = surf_value + do k = size(tracer,3) - 1, 1, -1 + tracer(:,:,k) = tracer(:,:,k+1) * multiplier + enddo + end select case default end select endif !scheme.eq.profile From ba92db8f5cadab7ab97740fc02f069b89bffd229 Mon Sep 17 00:00:00 2001 From: MinsukJi-NOAA Date: Fri, 8 Oct 2021 16:08:27 -0400 Subject: [PATCH 2/8] Specify platform_mod variables used --- tracer_manager/tracer_manager.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tracer_manager/tracer_manager.F90 b/tracer_manager/tracer_manager.F90 index afb510cef5..a8a30a2a39 100644 --- a/tracer_manager/tracer_manager.F90 +++ b/tracer_manager/tracer_manager.F90 @@ -76,7 +76,7 @@ module tracer_manager_mod fm_exists, & MODEL_NAMES -use platform_mod +use platform_mod, only: r4_kind, r8_kind implicit none private From 4b63ab91029716ca9e22a8b2cdb0ef40613c71ff Mon Sep 17 00:00:00 2001 From: MinsukJi-NOAA Date: Fri, 8 Oct 2021 21:41:27 -0400 Subject: [PATCH 3/8] Modify real_to_time_type in time_manager.F90 --- time_manager/time_manager.F90 | 22 ++++++++++++++++++---- 1 file changed, 18 insertions(+), 4 deletions(-) diff --git a/time_manager/time_manager.F90 b/time_manager/time_manager.F90 index 2a9599260a..53752b98a3 100644 --- a/time_manager/time_manager.F90 +++ b/time_manager/time_manager.F90 @@ -87,7 +87,7 @@ module time_manager_mod -use platform_mod, only: r8_kind +use platform_mod, only: r4_kind, r8_kind use constants_mod, only: rseconds_per_day=>seconds_per_day use fms_mod, only: error_mesg, FATAL, WARNING, write_version_number, stdout @@ -1206,7 +1206,7 @@ end function time_type_to_real !! @return A filled time type variable, and an error message if an !! error occurs. function real_to_time_type(x,err_msg) result(t) - real,intent(in) :: x !< Number of seconds. + class(*),intent(in) :: x !< Number of seconds. character(len=*),intent(out),optional :: err_msg !< Error message. type(time_type) :: t integer :: days @@ -1217,9 +1217,23 @@ function real_to_time_type(x,err_msg) result(t) real :: tps real :: a tps = real(ticks_per_second) - a = x/spd + + select type (x) + type is (real(kind=r4_kind)) + a = x/spd + type is (real(kind=r8_kind)) + a = x/spd + end select + days = safe_rtoi(a,do_floor) - a = x - real(days)*spd + + select type (x) + type is (real(kind=r4_kind)) + a = x - real(days)*spd + type is (real(kind=r8_kind)) + a = x - real(days)*spd + end select + seconds = safe_rtoi(a,do_floor) a = (a - real(seconds))*tps ticks = safe_rtoi(a,do_nearest) From a9d16385895d38a576f4c87b58f488d55e9f6318 Mon Sep 17 00:00:00 2001 From: MinsukJi-NOAA Date: Sun, 10 Oct 2021 11:15:34 -0400 Subject: [PATCH 4/8] Modify files in the sat_vapor_pres directory --- sat_vapor_pres/sat_vapor_pres.F90 | 237 ++- sat_vapor_pres/sat_vapor_pres_k.F90 | 2394 +++++++++++++++++++-------- 2 files changed, 1904 insertions(+), 727 deletions(-) diff --git a/sat_vapor_pres/sat_vapor_pres.F90 b/sat_vapor_pres/sat_vapor_pres.F90 index c92e134a94..901dfd1122 100644 --- a/sat_vapor_pres/sat_vapor_pres.F90 +++ b/sat_vapor_pres/sat_vapor_pres.F90 @@ -194,6 +194,8 @@ module sat_vapor_pres_mod lookup_des3_k, lookup_es3_des3_k, & compute_qs_k, compute_mrs_k + use platform_mod, only: r4_kind, r8_kind + implicit none private @@ -739,8 +741,8 @@ module sat_vapor_pres_mod ! subroutine lookup_es_0d ( temp, esat, err_msg ) - real, intent(in) :: temp - real, intent(out) :: esat + class(*), intent(in) :: temp + class(*), intent(out) :: esat character(len=*), intent(out), optional :: err_msg integer :: nbad @@ -771,8 +773,8 @@ end subroutine lookup_es_0d ! subroutine lookup_es_1d ( temp, esat, err_msg ) - real, intent(in) :: temp(:) - real, intent(out) :: esat(:) + class(*), intent(in) :: temp(:) + class(*), intent(out) :: esat(:) character(len=*), intent(out), optional :: err_msg character(len=54) :: err_msg_local @@ -807,8 +809,8 @@ end subroutine lookup_es_1d ! subroutine lookup_es_2d ( temp, esat, err_msg ) - real, intent(in) :: temp(:,:) - real, intent(out) :: esat(:,:) + class(*), intent(in) :: temp(:,:) + class(*), intent(out) :: esat(:,:) character(len=*), intent(out), optional :: err_msg character(len=54) :: err_msg_local @@ -843,8 +845,8 @@ end subroutine lookup_es_2d ! subroutine lookup_es_3d ( temp, esat, err_msg ) - real, intent(in) :: temp(:,:,:) - real, intent(out) :: esat(:,:,:) + class(*), intent(in) :: temp(:,:,:) + class(*), intent(out) :: esat(:,:,:) character(len=*), intent(out), optional :: err_msg integer :: nbad @@ -1975,10 +1977,10 @@ end subroutine lookup_es3_des3_3d subroutine compute_qs_0d ( temp, press, qsat, q, hc, dqsdT, esat, & err_msg, es_over_liq, es_over_liq_and_ice ) - real, intent(in) :: temp, press - real, intent(out) :: qsat - real, intent(in), optional :: q, hc - real, intent(out), optional :: dqsdT, esat + class(*), intent(in) :: temp, press + class(*), intent(out) :: qsat + class(*), intent(in), optional :: q, hc + class(*), intent(out), optional :: dqsdT, esat character(len=*), intent(out), optional :: err_msg logical,intent(in), optional :: es_over_liq logical,intent(in), optional :: es_over_liq_and_ice @@ -2033,11 +2035,11 @@ end subroutine compute_qs_0d subroutine compute_qs_1d ( temp, press, qsat, q, hc, dqsdT, esat, & err_msg, es_over_liq, es_over_liq_and_ice ) - real, intent(in) :: temp(:), press(:) - real, intent(out) :: qsat(:) - real, intent(in), optional :: q(:) -real, intent(in), optional :: hc - real, intent(out), optional :: dqsdT(:), esat(:) + class(*), intent(in) :: temp(:), press(:) + class(*), intent(out) :: qsat(:) + class(*), intent(in), optional :: q(:) + class(*), intent(in), optional :: hc + class(*), intent(out), optional :: dqsdT(:), esat(:) character(len=*), intent(out), optional :: err_msg logical,intent(in), optional :: es_over_liq logical,intent(in), optional :: es_over_liq_and_ice @@ -2095,11 +2097,11 @@ end subroutine compute_qs_1d subroutine compute_qs_2d ( temp, press, qsat, q, hc, dqsdT, esat, & err_msg, es_over_liq, es_over_liq_and_ice ) - real, intent(in) :: temp(:,:), press(:,:) - real, intent(out) :: qsat(:,:) - real, intent(in), optional :: q(:,:) - real, intent(in), optional :: hc - real, intent(out), optional :: dqsdT(:,:), esat(:,:) + class(*), intent(in) :: temp(:,:), press(:,:) + class(*), intent(out) :: qsat(:,:) + class(*), intent(in), optional :: q(:,:) + class(*), intent(in), optional :: hc + class(*), intent(out), optional :: dqsdT(:,:), esat(:,:) character(len=*), intent(out), optional :: err_msg logical,intent(in), optional :: es_over_liq logical,intent(in), optional :: es_over_liq_and_ice @@ -2156,11 +2158,11 @@ end subroutine compute_qs_2d subroutine compute_qs_3d ( temp, press, qsat, q, hc, dqsdT, esat, & err_msg, es_over_liq, es_over_liq_and_ice ) - real, intent(in) :: temp(:,:,:), press(:,:,:) - real, intent(out) :: qsat(:,:,:) - real, intent(in), optional :: q(:,:,:) - real, intent(in), optional :: hc - real, intent(out), optional :: dqsdT(:,:,:), esat(:,:,:) + class(*), intent(in) :: temp(:,:,:), press(:,:,:) + class(*), intent(out) :: qsat(:,:,:) + class(*), intent(in), optional :: q(:,:,:) + class(*), intent(in), optional :: hc + class(*), intent(out), optional :: dqsdT(:,:,:), esat(:,:,:) character(len=*), intent(out), optional :: err_msg logical,intent(in), optional :: es_over_liq logical,intent(in), optional :: es_over_liq_and_ice @@ -2608,74 +2610,120 @@ end subroutine sat_vapor_pres_init !####################################################################### function check_1d ( temp ) result ( nbad ) - real , intent(in) :: temp(:) + class(*), intent(in) :: temp(:) integer :: nbad, ind, i nbad = 0 - do i = 1, size(temp,1) - ind = int(dtinv*(temp(i)-tmin+teps)) - if (ind < 0 .or. ind > nlim) nbad = nbad+1 - enddo + + select type (temp) + type is (real(kind=r4_kind)) + do i = 1, size(temp,1) + ind = int(dtinv*(temp(i)-tmin+teps)) + if (ind < 0 .or. ind > nlim) nbad = nbad+1 + enddo + type is (real(kind=r8_kind)) + do i = 1, size(temp,1) + ind = int(dtinv*(temp(i)-tmin+teps)) + if (ind < 0 .or. ind > nlim) nbad = nbad+1 + enddo + end select end function check_1d !------------------------------------------------ function check_2d ( temp ) result ( nbad ) - real , intent(in) :: temp(:,:) + class(*), intent(in) :: temp(:,:) integer :: nbad integer :: j - nbad = 0 - do j = 1, size(temp,2) - nbad = nbad + check_1d ( temp(:,j) ) - enddo + nbad = 0 + + select type (temp) + type is (real(kind=r4_kind)) + do j = 1, size(temp,2) + nbad = nbad + check_1d ( temp(:,j) ) + enddo + type is (real(kind=r8_kind)) + do j = 1, size(temp,2) + nbad = nbad + check_1d ( temp(:,j) ) + enddo + end select + end function check_2d !####################################################################### subroutine temp_check_1d ( temp ) - real , intent(in) :: temp(:) + class(*), intent(in) :: temp(:) integer :: i, unit unit = stdoutunit - write(unit,*) 'Bad temperatures (dimension 1): ', (check_1d(temp(i:i)),i=1,size(temp,1)) + + select type (temp) + type is (real(kind=r4_kind)) + write(unit,*) 'Bad temperatures (dimension 1): ', (check_1d(temp(i:i)),i=1,size(temp,1)) + type is (real(kind=r8_kind)) + write(unit,*) 'Bad temperatures (dimension 1): ', (check_1d(temp(i:i)),i=1,size(temp,1)) + end select end subroutine temp_check_1d !-------------------------------------------------------------- subroutine temp_check_2d ( temp ) - real , intent(in) :: temp(:,:) + class(*), intent(in) :: temp(:,:) integer :: i, j, unit unit = stdoutunit - write(unit,*) 'Bad temperatures (dimension 1): ', (check_1d(temp(i,:)),i=1,size(temp,1)) - write(unit,*) 'Bad temperatures (dimension 2): ', (check_1d(temp(:,j)),j=1,size(temp,2)) + + select type (temp) + type is (real(kind=r4_kind)) + write(unit,*) 'Bad temperatures (dimension 1): ', (check_1d(temp(i,:)),i=1,size(temp,1)) + write(unit,*) 'Bad temperatures (dimension 2): ', (check_1d(temp(:,j)),j=1,size(temp,2)) + type is (real(kind=r8_kind)) + write(unit,*) 'Bad temperatures (dimension 1): ', (check_1d(temp(i,:)),i=1,size(temp,1)) + write(unit,*) 'Bad temperatures (dimension 2): ', (check_1d(temp(:,j)),j=1,size(temp,2)) + end select end subroutine temp_check_2d !-------------------------------------------------------------- subroutine temp_check_3d ( temp ) - real, intent(in) :: temp(:,:,:) + class(*), intent(in) :: temp(:,:,:) integer :: i, j, k, unit unit = stdoutunit - write(unit,*) 'Bad temperatures (dimension 1): ', (check_2d(temp(i,:,:)),i=1,size(temp,1)) - write(unit,*) 'Bad temperatures (dimension 2): ', (check_2d(temp(:,j,:)),j=1,size(temp,2)) - write(unit,*) 'Bad temperatures (dimension 3): ', (check_2d(temp(:,:,k)),k=1,size(temp,3)) + + select type (temp) + type is (real(kind=r4_kind)) + write(unit,*) 'Bad temperatures (dimension 1): ', (check_2d(temp(i,:,:)),i=1,size(temp,1)) + write(unit,*) 'Bad temperatures (dimension 2): ', (check_2d(temp(:,j,:)),j=1,size(temp,2)) + write(unit,*) 'Bad temperatures (dimension 3): ', (check_2d(temp(:,:,k)),k=1,size(temp,3)) + type is (real(kind=r8_kind)) + write(unit,*) 'Bad temperatures (dimension 1): ', (check_2d(temp(i,:,:)),i=1,size(temp,1)) + write(unit,*) 'Bad temperatures (dimension 2): ', (check_2d(temp(:,j,:)),j=1,size(temp,2)) + write(unit,*) 'Bad temperatures (dimension 3): ', (check_2d(temp(:,:,k)),k=1,size(temp,3)) + end select end subroutine temp_check_3d !####################################################################### subroutine show_all_bad_0d ( temp ) - real , intent(in) :: temp + class(*), intent(in) :: temp integer :: ind, unit unit = stdoutunit - ind = int(dtinv*(temp-tmin+teps)) + + select type (temp) + type is (real(kind=r4_kind)) + ind = int(dtinv*(temp-tmin+teps)) + type is (real(kind=r8_kind)) + ind = int(dtinv*(temp-tmin+teps)) + end select + if (ind < 0 .or. ind > nlim) then write(unit,'(a,e10.3,a,i6)') 'Bad temperature=',temp,' pe=',mpp_pe() endif @@ -2685,54 +2733,93 @@ end subroutine show_all_bad_0d !-------------------------------------------------------------- subroutine show_all_bad_1d ( temp ) - real , intent(in) :: temp(:) + class(*), intent(in) :: temp(:) integer :: i, ind, unit unit = stdoutunit - do i=1,size(temp) - ind = int(dtinv*(temp(i)-tmin+teps)) - if (ind < 0 .or. ind > nlim) then - write(unit,'(a,e10.3,a,i4,a,i6)') 'Bad temperature=',temp(i),' at i=',i,' pe=',mpp_pe() - endif - enddo + + select type (temp) + type is (real(kind=r4_kind)) + do i=1,size(temp) + ind = int(dtinv*(temp(i)-tmin+teps)) + if (ind < 0 .or. ind > nlim) then + write(unit,'(a,e10.3,a,i4,a,i6)') 'Bad temperature=',temp(i),' at i=',i,' pe=',mpp_pe() + endif + enddo + type is (real(kind=r8_kind)) + do i=1,size(temp) + ind = int(dtinv*(temp(i)-tmin+teps)) + if (ind < 0 .or. ind > nlim) then + write(unit,'(a,e10.3,a,i4,a,i6)') 'Bad temperature=',temp(i),' at i=',i,' pe=',mpp_pe() + endif + enddo + end select end subroutine show_all_bad_1d !-------------------------------------------------------------- subroutine show_all_bad_2d ( temp ) - real , intent(in) :: temp(:,:) + class(*), intent(in) :: temp(:,:) integer :: i, j, ind, unit unit = stdoutunit - do j=1,size(temp,2) - do i=1,size(temp,1) - ind = int(dtinv*(temp(i,j)-tmin+teps)) - if (ind < 0 .or. ind > nlim) then - write(unit,'(a,e10.3,a,i4,a,i4,a,i6)') 'Bad temperature=',temp(i,j),' at i=',i,' j=',j,' pe=',mpp_pe() - endif - enddo - enddo + + select type (temp) + type is (real(kind=r4_kind)) + do j=1,size(temp,2) + do i=1,size(temp,1) + ind = int(dtinv*(temp(i,j)-tmin+teps)) + if (ind < 0 .or. ind > nlim) then + write(unit,'(a,e10.3,a,i4,a,i4,a,i6)') 'Bad temperature=',temp(i,j),' at i=',i,' j=',j,' pe=',mpp_pe() + endif + enddo + enddo + type is (real(kind=r8_kind)) + do j=1,size(temp,2) + do i=1,size(temp,1) + ind = int(dtinv*(temp(i,j)-tmin+teps)) + if (ind < 0 .or. ind > nlim) then + write(unit,'(a,e10.3,a,i4,a,i4,a,i6)') 'Bad temperature=',temp(i,j),' at i=',i,' j=',j,' pe=',mpp_pe() + endif + enddo + enddo + end select end subroutine show_all_bad_2d !-------------------------------------------------------------- subroutine show_all_bad_3d ( temp ) - real, intent(in) :: temp(:,:,:) + class(*), intent(in) :: temp(:,:,:) integer :: i, j, k, ind, unit unit = stdoutunit - do k=1,size(temp,3) - do j=1,size(temp,2) - do i=1,size(temp,1) - ind = int(dtinv*(temp(i,j,k)-tmin+teps)) - if (ind < 0 .or. ind > nlim) then - write(unit,'(a,e10.3,a,i4,a,i4,a,i4,a,i6)') 'Bad temperature=',temp(i,j,k),' at i=',i,' j=',j,' k=',k,' pe=',mpp_pe() - endif - enddo - enddo - enddo + + select type (temp) + type is (real(kind=r4_kind)) + do k=1,size(temp,3) + do j=1,size(temp,2) + do i=1,size(temp,1) + ind = int(dtinv*(temp(i,j,k)-tmin+teps)) + if (ind < 0 .or. ind > nlim) then + write(unit,'(a,e10.3,a,i4,a,i4,a,i4,a,i6)') 'Bad temperature=',temp(i,j,k),' at i=',i,' j=',j,' k=',k,' pe=',mpp_pe() + endif + enddo + enddo + enddo + type is (real(kind=r8_kind)) + do k=1,size(temp,3) + do j=1,size(temp,2) + do i=1,size(temp,1) + ind = int(dtinv*(temp(i,j,k)-tmin+teps)) + if (ind < 0 .or. ind > nlim) then + write(unit,'(a,e10.3,a,i4,a,i4,a,i4,a,i6)') 'Bad temperature=',temp(i,j,k),' at i=',i,' j=',j,' k=',k,' pe=',mpp_pe() + endif + enddo + enddo + enddo + end select end subroutine show_all_bad_3d diff --git a/sat_vapor_pres/sat_vapor_pres_k.F90 b/sat_vapor_pres/sat_vapor_pres_k.F90 index a9662a7d3b..ca93cc0a37 100644 --- a/sat_vapor_pres/sat_vapor_pres_k.F90 +++ b/sat_vapor_pres/sat_vapor_pres_k.F90 @@ -50,6 +50,8 @@ module sat_vapor_pres_k_mod ! not be a fortran module. This complicates things greatly for questionable ! benefit and could be done as a second step anyway, if necessary. + use platform_mod, only: r4_kind, r8_kind + implicit none private @@ -475,15 +477,15 @@ end function compute_es_liq_ice_k subroutine compute_qs_k_3d (temp, press, eps, zvir, qs, nbad, q, hc, & dqsdT, esat, es_over_liq, es_over_liq_and_ice) - real, intent(in), dimension(:,:,:) :: temp, press - real, intent(in) :: eps, zvir - real, intent(out), dimension(:,:,:) :: qs - integer, intent(out) :: nbad - real, intent(in), dimension(:,:,:), optional :: q - real, intent(in), optional :: hc - real, intent(out), dimension(:,:,:), optional :: dqsdT, esat - logical,intent(in), optional :: es_over_liq - logical,intent(in), optional :: es_over_liq_and_ice + class(*), intent(in), dimension(:,:,:) :: temp, press + real, intent(in) :: eps, zvir + class(*), intent(out), dimension(:,:,:) :: qs + integer, intent(out) :: nbad + class(*), intent(in), dimension(:,:,:), optional :: q + class(*), intent(in), optional :: hc + class(*), intent(out), dimension(:,:,:), optional :: dqsdT, esat + logical,intent(in), optional :: es_over_liq + logical,intent(in), optional :: es_over_liq_and_ice real, dimension(size(temp,1), size(temp,2), size(temp,3)) :: & esloc, desat, denom @@ -491,66 +493,144 @@ subroutine compute_qs_k_3d (temp, press, eps, zvir, qs, nbad, q, hc, & real :: hc_loc if (present(hc)) then - hc_loc = hc + select type (hc) + type is (real(kind=r4_kind)) + hc_loc = hc + type is (real(kind=r8_kind)) + hc_loc = hc + end select else hc_loc = 1.0 endif - if (present(es_over_liq)) then - if (present (dqsdT)) then - call lookup_es2_des2_k (temp, esloc, desat, nbad) - desat = desat*hc_loc - else - call lookup_es2_k (temp, esloc, nbad) - endif - else if (present(es_over_liq_and_ice)) then - if (present (dqsdT)) then - call lookup_es3_des3_k (temp, esloc, desat, nbad) - desat = desat*hc_loc - else - call lookup_es3_k (temp, esloc, nbad) - endif - else - if (present (dqsdT)) then - call lookup_es_des_k (temp, esloc, desat, nbad) - desat = desat*hc_loc + + if (present(es_over_liq)) then + if (present (dqsdT)) then + call lookup_es2_des2_k (temp, esloc, desat, nbad) + desat = desat*hc_loc + else + call lookup_es2_k (temp, esloc, nbad) + endif + else if (present(es_over_liq_and_ice)) then + if (present (dqsdT)) then + call lookup_es3_des3_k (temp, esloc, desat, nbad) + desat = desat*hc_loc + else + call lookup_es3_k (temp, esloc, nbad) + endif else - call lookup_es_k (temp, esloc, nbad) + if (present (dqsdT)) then + call lookup_es_des_k (temp, esloc, desat, nbad) + desat = desat*hc_loc + else + call lookup_es_k (temp, esloc, nbad) + endif endif - endif + esloc = esloc*hc_loc if (present (esat)) then - esat = esloc + select type (esat) + type is (real(kind=r4_kind)) + esat = esloc + type is (real(kind=r8_kind)) + esat = esloc + end select endif + if (nbad == 0) then - if (present (q) .and. use_exact_qs) then - qs = (1.0 + zvir*q)*eps*esloc/press - if (present (dqsdT)) then - dqsdT = (1.0 + zvir*q)*eps*desat/press - endif - else ! (present(q)) - denom = press - (1.0 - eps)*esloc - do k=1,size(qs,3) - do j=1,size(qs,2) - do i=1,size(qs,1) - if (denom(i,j,k) > 0.0) then - qs(i,j,k) = eps*esloc(i,j,k)/denom(i,j,k) - else - qs(i,j,k) = eps + select type (press) + type is (real(kind=r4_kind)) + select type (qs) + type is (real(kind=r4_kind)) + if (present (q) .and. use_exact_qs) then + select type (q) + type is (real(kind=r4_kind)) + qs = (1.0 + zvir*q)*eps*esloc/press + if (present (dqsdT)) then + select type (dqsdT) + type is (real(kind=r4_kind)) + dqsdT = (1.0 + zvir*q)*eps*desat/press + end select endif + end select + else ! (present(q)) + denom = press - (1.0 - eps)*esloc + do k=1,size(qs,3) + do j=1,size(qs,2) + do i=1,size(qs,1) + if (denom(i,j,k) > 0.0) then + qs(i,j,k) = eps*esloc(i,j,k)/denom(i,j,k) + else + qs(i,j,k) = eps + endif + end do + end do end do - end do - end do - if (present (dqsdT)) then - dqsdT = eps*press*desat/denom**2 - endif - endif ! (present(q)) + if (present (dqsdT)) then + select type (dqsdT) + type is (real(kind=r4_kind)) + dqsdT = eps*press*desat/denom**2 + end select + endif + endif ! (present(q)) + end select + type is (real(kind=r8_kind)) + select type (qs) + type is (real(kind=r8_kind)) + if (present (q) .and. use_exact_qs) then + select type (q) + type is (real(kind=r8_kind)) + qs = (1.0 + zvir*q)*eps*esloc/press + if (present (dqsdT)) then + select type (dqsdT) + type is (real(kind=r8_kind)) + dqsdT = (1.0 + zvir*q)*eps*desat/press + end select + endif + end select + else ! (present(q)) + denom = press - (1.0 - eps)*esloc + do k=1,size(qs,3) + do j=1,size(qs,2) + do i=1,size(qs,1) + if (denom(i,j,k) > 0.0) then + qs(i,j,k) = eps*esloc(i,j,k)/denom(i,j,k) + else + qs(i,j,k) = eps + endif + end do + end do + end do + if (present (dqsdT)) then + select type (dqsdT) + type is (real(kind=r8_kind)) + dqsdT = eps*press*desat/denom**2 + end select + endif + endif ! (present(q)) + end select + end select else ! (nbad = 0) - qs = -999. + select type (qs) + type is (real(kind=r4_kind)) + qs = -999. + type is (real(kind=r8_kind)) + qs = -999. + end select if (present (dqsdT)) then - dqsdT = -999. + select type (dqsdT) + type is (real(kind=r4_kind)) + dqsdT = -999. + type is (real(kind=r8_kind)) + dqsdT = -999. + end select endif if (present (esat)) then - esat = -999. + select type (esat) + type is (real(kind=r4_kind)) + esat = -999. + type is (real(kind=r8_kind)) + esat = -999. + end select endif endif ! (nbad = 0) @@ -562,80 +642,155 @@ end subroutine compute_qs_k_3d subroutine compute_qs_k_2d (temp, press, eps, zvir, qs, nbad, q, hc, & dqsdT, esat, es_over_liq, es_over_liq_and_ice) - real, intent(in), dimension(:,:) :: temp, press - real, intent(in) :: eps, zvir - real, intent(out), dimension(:,:) :: qs - integer, intent(out) :: nbad - real, intent(in), dimension(:,:), optional :: q - real, intent(in), optional :: hc - real, intent(out), dimension(:,:), optional :: dqsdT, esat - logical,intent(in), optional :: es_over_liq - logical,intent(in), optional :: es_over_liq_and_ice + class(*), intent(in), dimension(:,:) :: temp, press + real, intent(in) :: eps, zvir + class(*), intent(out), dimension(:,:) :: qs + integer, intent(out) :: nbad + class(*), intent(in), dimension(:,:), optional :: q + class(*), intent(in), optional :: hc + class(*), intent(out), dimension(:,:), optional :: dqsdT, esat + logical,intent(in), optional :: es_over_liq + logical,intent(in), optional :: es_over_liq_and_ice real, dimension(size(temp,1), size(temp,2)) :: esloc, desat, denom integer :: i, j real :: hc_loc if (present(hc)) then - hc_loc = hc + select type (hc) + type is (real(kind=r4_kind)) + hc_loc = hc + type is (real(kind=r8_kind)) + hc_loc = hc + end select else hc_loc = 1.0 endif - if (present(es_over_liq)) then - if (present (dqsdT)) then - call lookup_es2_des2_k (temp, esloc, desat, nbad) - desat = desat*hc_loc - else - call lookup_es2_k (temp, esloc, nbad) - endif - else if (present(es_over_liq_and_ice)) then - if (present (dqsdT)) then - call lookup_es3_des3_k (temp, esloc, desat, nbad) - desat = desat*hc_loc - else - call lookup_es3_k (temp, esloc, nbad) - endif - else - if (present (dqsdT)) then - call lookup_es_des_k (temp, esloc, desat, nbad) - desat = desat*hc_loc + if (present(es_over_liq)) then + if (present (dqsdT)) then + call lookup_es2_des2_k (temp, esloc, desat, nbad) + desat = desat*hc_loc + else + call lookup_es2_k (temp, esloc, nbad) + endif + else if (present(es_over_liq_and_ice)) then + if (present (dqsdT)) then + call lookup_es3_des3_k (temp, esloc, desat, nbad) + desat = desat*hc_loc + else + call lookup_es3_k (temp, esloc, nbad) + endif else - call lookup_es_k (temp, esloc, nbad) + if (present (dqsdT)) then + call lookup_es_des_k (temp, esloc, desat, nbad) + desat = desat*hc_loc + else + call lookup_es_k (temp, esloc, nbad) + endif endif - endif + esloc = esloc*hc_loc if (present (esat)) then - esat = esloc + select type (esat) + type is (real(kind=r4_kind)) + esat = esloc + type is (real(kind=r8_kind)) + esat = esloc + end select endif + if (nbad == 0) then - if (present (q) .and. use_exact_qs) then - qs = (1.0 + zvir*q)*eps*esloc/press - if (present (dqsdT)) then - dqsdT = (1.0 + zvir*q)*eps*desat/press - endif - else ! (present(q)) - denom = press - (1.0 - eps)*esloc - do j=1,size(qs,2) - do i=1,size(qs,1) - if (denom(i,j) > 0.0) then - qs(i,j) = eps*esloc(i,j)/denom(i,j) - else - qs(i,j) = eps - endif - end do - end do - if (present (dqsdT)) then - dqsdT = eps*press*desat/denom**2 - endif - endif ! (present(q)) + select type (press) + type is (real(kind=r4_kind)) + select type (qs) + type is (real(kind=r4_kind)) + if (present (q) .and. use_exact_qs) then + select type (q) + type is (real(kind=r4_kind)) + qs = (1.0 + zvir*q)*eps*esloc/press + if (present (dqsdT)) then + select type (dqsdT) + type is (real(kind=r4_kind)) + dqsdT = (1.0 + zvir*q)*eps*desat/press + end select + endif + end select + else ! (present(q)) + denom = press - (1.0 - eps)*esloc + do j=1,size(qs,2) + do i=1,size(qs,1) + if (denom(i,j) > 0.0) then + qs(i,j) = eps*esloc(i,j)/denom(i,j) + else + qs(i,j) = eps + endif + end do + end do + if (present (dqsdT)) then + select type (dqsdT) + type is (real(kind=r4_kind)) + dqsdT = eps*press*desat/denom**2 + end select + endif + endif ! (present(q)) + end select + type is (real(kind=r8_kind)) + select type (qs) + type is (real(kind=r8_kind)) + if (present (q) .and. use_exact_qs) then + select type (q) + type is (real(kind=r8_kind)) + qs = (1.0 + zvir*q)*eps*esloc/press + if (present (dqsdT)) then + select type (dqsdT) + type is (real(kind=r8_kind)) + dqsdT = (1.0 + zvir*q)*eps*desat/press + end select + endif + end select + else ! (present(q)) + denom = press - (1.0 - eps)*esloc + do j=1,size(qs,2) + do i=1,size(qs,1) + if (denom(i,j) > 0.0) then + qs(i,j) = eps*esloc(i,j)/denom(i,j) + else + qs(i,j) = eps + endif + end do + end do + if (present (dqsdT)) then + select type (dqsdT) + type is (real(kind=r8_kind)) + dqsdT = eps*press*desat/denom**2 + end select + endif + endif ! (present(q)) + end select + end select else ! (nbad = 0) - qs = -999. + select type (qs) + type is (real(kind=r4_kind)) + qs = -999. + type is (real(kind=r8_kind)) + qs = -999. + end select if (present (dqsdT)) then - dqsdT = -999. + select type (dqsdT) + type is (real(kind=r4_kind)) + dqsdT = -999. + type is (real(kind=r8_kind)) + dqsdT = -999. + end select endif if (present (esat)) then - esat = -999. + select type (esat) + type is (real(kind=r4_kind)) + esat = -999. + type is (real(kind=r8_kind)) + esat = -999. + end select endif endif ! (nbad = 0) @@ -647,78 +802,151 @@ end subroutine compute_qs_k_2d subroutine compute_qs_k_1d (temp, press, eps, zvir, qs, nbad, q, hc, & dqsdT, esat, es_over_liq, es_over_liq_and_ice) - real, intent(in), dimension(:) :: temp, press - real, intent(in) :: eps, zvir - real, intent(out), dimension(:) :: qs - integer, intent(out) :: nbad - real, intent(in), dimension(:), optional :: q - real, intent(in), optional :: hc - real, intent(out), dimension(:), optional :: dqsdT, esat - logical,intent(in), optional :: es_over_liq - logical,intent(in), optional :: es_over_liq_and_ice + class(*), intent(in), dimension(:) :: temp, press + real, intent(in) :: eps, zvir + class(*), intent(out),dimension(:) :: qs + integer, intent(out) :: nbad + class(*), intent(in), dimension(:), optional :: q + class(*), intent(in), optional :: hc + class(*), intent(out), dimension(:),optional :: dqsdT, esat + logical,intent(in), optional :: es_over_liq + logical,intent(in), optional :: es_over_liq_and_ice real, dimension(size(temp,1)) :: esloc, desat, denom integer :: i real :: hc_loc if (present(hc)) then - hc_loc = hc + select type (hc) + type is (real(kind=r4_kind)) + hc_loc = hc + type is (real(kind=r4_kind)) + hc_loc = hc + end select else hc_loc = 1.0 endif - if (present(es_over_liq)) then - if (present (dqsdT)) then - call lookup_es2_des2_k (temp, esloc, desat, nbad) - desat = desat*hc_loc - else - call lookup_es2_k (temp, esloc, nbad) - endif - else if (present(es_over_liq_and_ice)) then - if (present (dqsdT)) then - call lookup_es3_des3_k (temp, esloc, desat, nbad) - desat = desat*hc_loc - else - call lookup_es3_k (temp, esloc, nbad) - endif - else - if (present (dqsdT)) then - call lookup_es_des_k (temp, esloc, desat, nbad) - desat = desat*hc_loc + if (present(es_over_liq)) then + if (present (dqsdT)) then + call lookup_es2_des2_k (temp, esloc, desat, nbad) + desat = desat*hc_loc + else + call lookup_es2_k (temp, esloc, nbad) + endif + else if (present(es_over_liq_and_ice)) then + if (present (dqsdT)) then + call lookup_es3_des3_k (temp, esloc, desat, nbad) + desat = desat*hc_loc + else + call lookup_es3_k (temp, esloc, nbad) + endif else - call lookup_es_k (temp, esloc, nbad) + if (present (dqsdT)) then + call lookup_es_des_k (temp, esloc, desat, nbad) + desat = desat*hc_loc + else + call lookup_es_k (temp, esloc, nbad) + endif endif - endif + esloc = esloc*hc_loc if (present (esat)) then - esat = esloc + select type (esat) + type is (real(kind=r4_kind)) + esat = esloc + type is (real(kind=r8_kind)) + esat = esloc + end select endif + if (nbad == 0) then - if (present (q) .and. use_exact_qs) then - qs = (1.0 + zvir*q)*eps*esloc/press - if (present (dqsdT)) then - dqsdT = (1.0 + zvir*q)*eps*desat/press - endif - else ! (present(q)) - denom = press - (1.0 - eps)*esloc - do i=1,size(qs,1) - if (denom(i) > 0.0) then - qs(i) = eps*esloc(i)/denom(i) - else - qs(i) = eps - endif - end do - if (present (dqsdT)) then - dqsdT = eps*press*desat/denom**2 - endif - endif ! (present(q)) + select type (press) + type is (real(kind=r4_kind)) + select type (qs) + type is (real(kind=r4_kind)) + if (present (q) .and. use_exact_qs) then + select type (q) + type is (real(kind=r4_kind)) + qs = (1.0 + zvir*q)*eps*esloc/press + if (present (dqsdT)) then + select type (dqsdT) + type is (real(kind=r4_kind)) + dqsdT = (1.0 + zvir*q)*eps*desat/press + end select + endif + end select + else ! (present(q)) + denom = press - (1.0 - eps)*esloc + do i=1,size(qs,1) + if (denom(i) > 0.0) then + qs(i) = eps*esloc(i)/denom(i) + else + qs(i) = eps + endif + end do + if (present (dqsdT)) then + select type (dqsdT) + type is (real(kind=r4_kind)) + dqsdT = eps*press*desat/denom**2 + end select + endif + endif ! (present(q)) + end select + type is (real(kind=r8_kind)) + select type (qs) + type is (real(kind=r8_kind)) + if (present (q) .and. use_exact_qs) then + select type (q) + type is (real(kind=r8_kind)) + qs = (1.0 + zvir*q)*eps*esloc/press + if (present (dqsdT)) then + select type (dqsdT) + type is (real(kind=r8_kind)) + dqsdT = (1.0 + zvir*q)*eps*desat/press + end select + endif + end select + else ! (present(q)) + denom = press - (1.0 - eps)*esloc + do i=1,size(qs,1) + if (denom(i) > 0.0) then + qs(i) = eps*esloc(i)/denom(i) + else + qs(i) = eps + endif + end do + if (present (dqsdT)) then + select type (dqsdT) + type is (real(kind=r8_kind)) + dqsdT = eps*press*desat/denom**2 + end select + endif + endif ! (present(q)) + end select + end select else ! (nbad = 0) - qs = -999. + select type (qs) + type is (real(kind=r4_kind)) + qs = -999. + type is (real(kind=r8_kind)) + qs = -999. + end select if (present (dqsdT)) then - dqsdT = -999. + select type (dqsdT) + type is (real(kind=r4_kind)) + dqsdT = -999. + type is (real(kind=r8_kind)) + dqsdT = -999. + end select endif if (present (esat)) then - esat = -999. + select type (esat) + type is (real(kind=r4_kind)) + esat = -999. + type is (real(kind=r8_kind)) + esat = -999. + end select endif endif ! (nbad = 0) @@ -730,79 +958,149 @@ end subroutine compute_qs_k_1d subroutine compute_qs_k_0d (temp, press, eps, zvir, qs, nbad, q, hc, & dqsdT, esat, es_over_liq, es_over_liq_and_ice) - real, intent(in) :: temp, press + class(*), intent(in) :: temp, press real, intent(in) :: eps, zvir - real, intent(out) :: qs + class(*), intent(out) :: qs integer, intent(out) :: nbad - real, intent(in), optional :: q - real, intent(in), optional :: hc - real, intent(out), optional :: dqsdT, esat + class(*), intent(in), optional :: q + class(*), intent(in), optional :: hc + class(*), intent(out), optional :: dqsdT, esat logical,intent(in), optional :: es_over_liq - logical,intent(in), optional :: es_over_liq_and_ice + logical,intent(in), optional :: es_over_liq_and_ice real :: esloc, desat, denom real :: hc_loc if (present(hc)) then - hc_loc = hc + select type (hc) + type is (real(kind=r4_kind)) + hc_loc = hc + type is (real(kind=r8_kind)) + hc_loc = hc + end select else hc_loc = 1.0 endif - if (present(es_over_liq)) then - if (present (dqsdT)) then - call lookup_es2_des2_k (temp, esloc, desat, nbad) - desat = desat*hc_loc - else - call lookup_es2_k (temp, esloc, nbad) - endif - else if (present(es_over_liq_and_ice)) then - if (present (dqsdT)) then - call lookup_es3_des3_k (temp, esloc, desat, nbad) - desat = desat*hc_loc - else - call lookup_es3_k (temp, esloc, nbad) - endif - else - if (present (dqsdT)) then - call lookup_es_des_k (temp, esloc, desat, nbad) - desat = desat*hc_loc + if (present(es_over_liq)) then + if (present (dqsdT)) then + call lookup_es2_des2_k (temp, esloc, desat, nbad) + desat = desat*hc_loc + else + call lookup_es2_k (temp, esloc, nbad) + endif + else if (present(es_over_liq_and_ice)) then + if (present (dqsdT)) then + call lookup_es3_des3_k (temp, esloc, desat, nbad) + desat = desat*hc_loc + else + call lookup_es3_k (temp, esloc, nbad) + endif else - call lookup_es_k (temp, esloc, nbad) + if (present (dqsdT)) then + call lookup_es_des_k (temp, esloc, desat, nbad) + desat = desat*hc_loc + else + call lookup_es_k (temp, esloc, nbad) + endif endif - endif + esloc = esloc*hc_loc if (present (esat)) then - esat = esloc + select type (esat) + type is (real(kind=r4_kind)) + esat = esloc + type is (real(kind=r8_kind)) + esat = esloc + end select endif + if (nbad == 0) then - if (present (q) .and. use_exact_qs) then - qs = (1.0 + zvir*q)*eps*esloc/press - if (present (dqsdT)) then - dqsdT = (1.0 + zvir*q)*eps*desat/press - endif - else ! (present(q)) - denom = press - (1.0 - eps)*esloc - if (denom > 0.0) then - qs = eps*esloc/denom - else - qs = eps - endif - if (present (dqsdT)) then - dqsdT = eps*press*desat/denom**2 - endif - endif ! (present(q)) + select type (press) + type is (real(kind=r4_kind)) + select type (qs) + type is (real(kind=r4_kind)) + if (present (q) .and. use_exact_qs) then + select type (q) + type is (real(kind=r4_kind)) + qs = (1.0 + zvir*q)*eps*esloc/press + if (present (dqsdT)) then + select type (dqsdT) + type is (real(kind=r4_kind)) + dqsdT = (1.0 + zvir*q)*eps*desat/press + end select + endif + end select + else ! (present(q)) + denom = press - (1.0 - eps)*esloc + if (denom > 0.0) then + qs = eps*esloc/denom + else + qs = eps + endif + if (present (dqsdT)) then + select type (dqsdT) + type is (real(kind=r4_kind)) + dqsdT = eps*press*desat/denom**2 + end select + endif + endif ! (present(q)) + end select + type is (real(kind=r8_kind)) + select type (qs) + type is (real(kind=r8_kind)) + if (present (q) .and. use_exact_qs) then + select type (q) + type is (real(kind=r8_kind)) + qs = (1.0 + zvir*q)*eps*esloc/press + if (present (dqsdT)) then + select type (dqsdT) + type is (real(kind=r8_kind)) + dqsdT = (1.0 + zvir*q)*eps*desat/press + end select + endif + end select + else ! (present(q)) + denom = press - (1.0 - eps)*esloc + if (denom > 0.0) then + qs = eps*esloc/denom + else + qs = eps + endif + if (present (dqsdT)) then + select type (dqsdT) + type is (real(kind=r8_kind)) + dqsdT = eps*press*desat/denom**2 + end select + endif + endif ! (present(q)) + end select + end select else ! (nbad = 0) - qs = -999. + select type (qs) + type is (real(kind=r4_kind)) + qs = -999. + type is (real(kind=r8_kind)) + qs = -999. + end select if (present (dqsdT)) then - dqsdT = -999. + select type (dqsdT) + type is (real(kind=r4_kind)) + dqsdT = -999. + type is (real(kind=r8_kind)) + dqsdT = -999. + end select endif if (present (esat)) then - esat = -999. + select type (esat) + type is (real(kind=r4_kind)) + esat = -999. + type is (real(kind=r8_kind)) + esat = -999. + end select endif endif ! (nbad = 0) - end subroutine compute_qs_k_0d !####################################################################### @@ -1148,107 +1446,211 @@ end subroutine compute_mrs_k_0d !####################################################################### subroutine lookup_es_des_k_3d (temp, esat, desat, nbad) - real, intent(in), dimension(:,:,:) :: temp - real, intent(out), dimension(:,:,:) :: esat, desat + class(*), intent(in), dimension(:,:,:) :: temp + class(*), intent(out), dimension(:,:,:) :: esat, desat integer, intent(out) :: nbad real :: tmp, del integer :: ind, i, j, k nbad = 0 - do k = 1, size(temp,3) - do j = 1, size(temp,2) - do i = 1, size(temp,1) - tmp = temp(i,j,k)-tminl - ind = int(dtinvl*(tmp+tepsl)) - if (ind < 0 .or. ind >= table_siz) then - nbad = nbad+1 - else - del = tmp-dtres*real(ind) - esat(i,j,k) = TABLE(ind+1) + & - del*(DTABLE(ind+1) + del*D2TABLE(ind+1)) - desat(i,j,k) = DTABLE(ind+1) + 2.*del*D2TABLE(ind+1) - endif - enddo - enddo - enddo + + select type (temp) + type is (real(kind=r4_kind)) + select type (esat) + type is (real(kind=r4_kind)) + select type (desat) + type is (real(kind=r4_kind)) + do k = 1, size(temp,3) + do j = 1, size(temp,2) + do i = 1, size(temp,1) + tmp = temp(i,j,k)-tminl + ind = int(dtinvl*(tmp+tepsl)) + if (ind < 0 .or. ind >= table_siz) then + nbad = nbad+1 + else + del = tmp-dtres*real(ind) + esat(i,j,k) = TABLE(ind+1) + del*(DTABLE(ind+1)+del*D2TABLE(ind+1)) + desat(i,j,k) = DTABLE(ind+1) + 2.*del*D2TABLE(ind+1) + endif + enddo + enddo + enddo + end select + end select + type is (real(kind=r8_kind)) + select type (esat) + type is (real(kind=r8_kind)) + select type (desat) + type is (real(kind=r8_kind)) + do k = 1, size(temp,3) + do j = 1, size(temp,2) + do i = 1, size(temp,1) + tmp = temp(i,j,k)-tminl + ind = int(dtinvl*(tmp+tepsl)) + if (ind < 0 .or. ind >= table_siz) then + nbad = nbad+1 + else + del = tmp-dtres*real(ind) + esat(i,j,k) = TABLE(ind+1) + del*(DTABLE(ind+1)+del*D2TABLE(ind+1)) + desat(i,j,k) = DTABLE(ind+1) + 2.*del*D2TABLE(ind+1) + endif + enddo + enddo + enddo + end select + end select + end select end subroutine lookup_es_des_k_3d !####################################################################### subroutine lookup_es_des_k_2d (temp, esat, desat, nbad) - real, intent(in), dimension(:,:) :: temp - real, intent(out), dimension(:,:) :: esat, desat + class(*), intent(in), dimension(:,:) :: temp + class(*), intent(out), dimension(:,:) :: esat, desat integer, intent(out) :: nbad real :: tmp, del integer :: ind, i, j nbad = 0 - do j = 1, size(temp,2) - do i = 1, size(temp,1) - tmp = temp(i,j)-tminl - ind = int(dtinvl*(tmp+tepsl)) - if (ind < 0 .or. ind >= table_siz) then - nbad = nbad+1 - else - del = tmp-dtres*real(ind) - esat(i,j) = TABLE(ind+1) + & - del*(DTABLE(ind+1) + del*D2TABLE(ind+1)) - desat(i,j) = DTABLE(ind+1) + 2.*del*D2TABLE(ind+1) - endif - enddo - enddo + + select type (temp) + type is (real(kind=r4_kind)) + select type (esat) + type is (real(kind=r4_kind)) + select type (desat) + type is (real(kind=r4_kind)) + do j = 1, size(temp,2) + do i = 1, size(temp,1) + tmp = temp(i,j)-tminl + ind = int(dtinvl*(tmp+tepsl)) + if (ind < 0 .or. ind >= table_siz) then + nbad = nbad+1 + else + del = tmp-dtres*real(ind) + esat(i,j) = TABLE(ind+1) + del*(DTABLE(ind+1)+del*D2TABLE(ind+1)) + desat(i,j) = DTABLE(ind+1) + 2.*del*D2TABLE(ind+1) + endif + enddo + enddo + end select + end select + type is (real(kind=r8_kind)) + select type (esat) + type is (real(kind=r8_kind)) + select type (desat) + type is (real(kind=r8_kind)) + do j = 1, size(temp,2) + do i = 1, size(temp,1) + tmp = temp(i,j)-tminl + ind = int(dtinvl*(tmp+tepsl)) + if (ind < 0 .or. ind >= table_siz) then + nbad = nbad+1 + else + del = tmp-dtres*real(ind) + esat(i,j) = TABLE(ind+1) + del*(DTABLE(ind+1)+del*D2TABLE(ind+1)) + desat(i,j) = DTABLE(ind+1) + 2.*del*D2TABLE(ind+1) + endif + enddo + enddo + end select + end select + end select end subroutine lookup_es_des_k_2d !####################################################################### subroutine lookup_es_des_k_1d (temp, esat, desat, nbad) - real, intent(in), dimension(:) :: temp - real, intent(out), dimension(:) :: esat, desat + class(*), intent(in), dimension(:) :: temp + class(*), intent(out), dimension(:) :: esat, desat integer, intent(out) :: nbad real :: tmp, del integer :: ind, i nbad = 0 - do i = 1, size(temp,1) - tmp = temp(i)-tminl - ind = int(dtinvl*(tmp+tepsl)) - if (ind < 0 .or. ind >= table_siz) then - nbad = nbad+1 - else - del = tmp-dtres*real(ind) - esat(i) = TABLE(ind+1) + & - del*(DTABLE(ind+1) + del*D2TABLE(ind+1)) - desat(i) = DTABLE(ind+1) + 2.*del*D2TABLE(ind+1) - endif - enddo + + select type (temp) + type is (real(kind=r4_kind)) + select type (esat) + type is (real(kind=r4_kind)) + select type (desat) + type is (real(kind=r4_kind)) + do i = 1, size(temp,1) + tmp = temp(i)-tminl + ind = int(dtinvl*(tmp+tepsl)) + if (ind < 0 .or. ind >= table_siz) then + nbad = nbad+1 + else + del = tmp-dtres*real(ind) + esat(i) = TABLE(ind+1) + del*(DTABLE(ind+1)+del*D2TABLE(ind+1)) + desat(i) = DTABLE(ind+1) + 2.*del*D2TABLE(ind+1) + endif + enddo + end select + end select + type is (real(kind=r8_kind)) + select type (esat) + type is (real(kind=r8_kind)) + select type (desat) + type is (real(kind=r8_kind)) + do i = 1, size(temp,1) + tmp = temp(i)-tminl + ind = int(dtinvl*(tmp+tepsl)) + if (ind < 0 .or. ind >= table_siz) then + nbad = nbad+1 + else + del = tmp-dtres*real(ind) + esat(i) = TABLE(ind+1) + del*(DTABLE(ind+1)+del*D2TABLE(ind+1)) + desat(i) = DTABLE(ind+1) + 2.*del*D2TABLE(ind+1) + endif + enddo + end select + end select + end select end subroutine lookup_es_des_k_1d !####################################################################### subroutine lookup_es_des_k_0d (temp, esat, desat, nbad) - real, intent(in) :: temp - real, intent(out) :: esat, desat + class(*), intent(in) :: temp + class(*), intent(out) :: esat, desat integer, intent(out) :: nbad real :: tmp, del integer :: ind nbad = 0 - tmp = temp-tminl + + select type (temp) + type is (real(kind=r4_kind)) + tmp = temp-tminl + type is (real(kind=r8_kind)) + tmp = temp-tminl + end select + ind = int(dtinvl*(tmp+tepsl)) if (ind < 0 .or. ind >= table_siz) then nbad = nbad+1 else del = tmp-dtres*real(ind) - esat = TABLE(ind+1) + & - del*(DTABLE(ind+1) + del*D2TABLE(ind+1)) - desat = DTABLE(ind+1) + 2.*del*D2TABLE(ind+1) + select type (esat) + type is (real(kind=r4_kind)) + esat = TABLE(ind+1)+del*(DTABLE(ind+1)+del*D2TABLE(ind+1)) + type is (real(kind=r8_kind)) + esat = TABLE(ind+1)+del*(DTABLE(ind+1)+del*D2TABLE(ind+1)) + end select + + select type (desat) + type is (real(kind=r4_kind)) + desat = DTABLE(ind+1) + 2.*del*D2TABLE(ind+1) + type is (real(kind=r8_kind)) + desat = DTABLE(ind+1) + 2.*del*D2TABLE(ind+1) + end select endif end subroutine lookup_es_des_k_0d @@ -1256,289 +1658,553 @@ end subroutine lookup_es_des_k_0d !####################################################################### subroutine lookup_es_k_3d(temp, esat, nbad) - real, intent(in), dimension(:,:,:) :: temp - real, intent(out), dimension(:,:,:) :: esat + class(*), intent(in), dimension(:,:,:) :: temp + class(*), intent(out), dimension(:,:,:) :: esat integer, intent(out) :: nbad real :: tmp, del integer :: ind, i, j, k nbad = 0 - do k = 1, size(temp,3) - do j = 1, size(temp,2) - do i = 1, size(temp,1) - tmp = temp(i,j,k)-tminl - ind = int(dtinvl*(tmp+tepsl)) - if (ind < 0 .or. ind >= table_siz) then - nbad = nbad+1 - else - del = tmp-dtres*real(ind) - esat(i,j,k) = TABLE(ind+1) + & - del*(DTABLE(ind+1) + del*D2TABLE(ind+1)) - endif - enddo - enddo - enddo + + select type (temp) + type is (real(kind=r4_kind)) + select type (esat) + type is (real(kind=r4_kind)) + do k = 1, size(temp,3) + do j = 1, size(temp,2) + do i = 1, size(temp,1) + tmp = temp(i,j,k)-tminl + ind = int(dtinvl*(tmp+tepsl)) + if (ind < 0 .or. ind >= table_siz) then + nbad = nbad+1 + else + del = tmp-dtres*real(ind) + esat(i,j,k) = TABLE(ind+1)+del*(DTABLE(ind+1)+del*D2TABLE(ind+1)) + endif + enddo + enddo + enddo + end select + type is (real(kind=r8_kind)) + select type (esat) + type is (real(kind=r8_kind)) + do k = 1, size(temp,3) + do j = 1, size(temp,2) + do i = 1, size(temp,1) + tmp = temp(i,j,k)-tminl + ind = int(dtinvl*(tmp+tepsl)) + if (ind < 0 .or. ind >= table_siz) then + nbad = nbad+1 + else + del = tmp-dtres*real(ind) + esat(i,j,k) = TABLE(ind+1)+del*(DTABLE(ind+1)+del*D2TABLE(ind+1)) + endif + enddo + enddo + enddo + end select + end select end subroutine lookup_es_k_3d !####################################################################### subroutine lookup_des_k_3d(temp, desat, nbad) - real, intent(in), dimension(:,:,:) :: temp - real, intent(out), dimension(:,:,:) :: desat + class(*), intent(in), dimension(:,:,:) :: temp + class(*), intent(out), dimension(:,:,:) :: desat integer, intent(out) :: nbad real :: tmp, del integer :: ind, i, j, k nbad = 0 - do k = 1, size(temp,3) - do j = 1, size(temp,2) - do i = 1, size(temp,1) - tmp = temp(i,j,k)-tminl - ind = int(dtinvl*(tmp+tepsl)) - if (ind < 0 .or. ind >= table_siz) then - nbad = nbad+1 - else - del = tmp-dtres*real(ind) - desat(i,j,k) = DTABLE(ind+1) + 2.*del*D2TABLE(ind+1) - endif - enddo - enddo - enddo + + select type (temp) + type is (real(kind=r4_kind)) + select type (desat) + type is (real(kind=r4_kind)) + do k = 1, size(temp,3) + do j = 1, size(temp,2) + do i = 1, size(temp,1) + tmp = temp(i,j,k)-tminl + ind = int(dtinvl*(tmp+tepsl)) + if (ind < 0 .or. ind >= table_siz) then + nbad = nbad+1 + else + del = tmp-dtres*real(ind) + desat(i,j,k) = DTABLE(ind+1) + 2.*del*D2TABLE(ind+1) + endif + enddo + enddo + enddo + end select + type is (real(kind=r8_kind)) + select type (desat) + type is (real(kind=r8_kind)) + do k = 1, size(temp,3) + do j = 1, size(temp,2) + do i = 1, size(temp,1) + tmp = temp(i,j,k)-tminl + ind = int(dtinvl*(tmp+tepsl)) + if (ind < 0 .or. ind >= table_siz) then + nbad = nbad+1 + else + del = tmp-dtres*real(ind) + desat(i,j,k) = DTABLE(ind+1) + 2.*del*D2TABLE(ind+1) + endif + enddo + enddo + enddo + end select + end select end subroutine lookup_des_k_3d !####################################################################### subroutine lookup_des_k_2d(temp, desat, nbad) - real, intent(in), dimension(:,:) :: temp - real, intent(out), dimension(:,:) :: desat + class(*), intent(in), dimension(:,:) :: temp + class(*), intent(out), dimension(:,:) :: desat integer, intent(out) :: nbad real :: tmp, del integer :: ind, i, j nbad = 0 - do j = 1, size(temp,2) - do i = 1, size(temp,1) - tmp = temp(i,j)-tminl - ind = int(dtinvl*(tmp+tepsl)) - if (ind < 0 .or. ind >= table_siz) then - nbad = nbad+1 - else - del = tmp-dtres*real(ind) - desat(i,j) = DTABLE(ind+1) + 2.*del*D2TABLE(ind+1) - endif - enddo - enddo + + select type (temp) + type is (real(kind=r4_kind)) + select type (desat) + type is (real(kind=r4_kind)) + do j = 1, size(temp,2) + do i = 1, size(temp,1) + tmp = temp(i,j)-tminl + ind = int(dtinvl*(tmp+tepsl)) + if (ind < 0 .or. ind >= table_siz) then + nbad = nbad+1 + else + del = tmp-dtres*real(ind) + desat(i,j) = DTABLE(ind+1) + 2.*del*D2TABLE(ind+1) + endif + enddo + enddo + end select + type is (real(kind=r8_kind)) + select type (desat) + type is (real(kind=r8_kind)) + do j = 1, size(temp,2) + do i = 1, size(temp,1) + tmp = temp(i,j)-tminl + ind = int(dtinvl*(tmp+tepsl)) + if (ind < 0 .or. ind >= table_siz) then + nbad = nbad+1 + else + del = tmp-dtres*real(ind) + desat(i,j) = DTABLE(ind+1) + 2.*del*D2TABLE(ind+1) + endif + enddo + enddo + end select + end select end subroutine lookup_des_k_2d !####################################################################### subroutine lookup_es_k_2d(temp, esat, nbad) - real, intent(in), dimension(:,:) :: temp - real, intent(out), dimension(:,:) :: esat + class(*), intent(in), dimension(:,:) :: temp + class(*), intent(out), dimension(:,:) :: esat integer, intent(out) :: nbad real :: tmp, del integer :: ind, i, j nbad = 0 - do j = 1, size(temp,2) - do i = 1, size(temp,1) - tmp = temp(i,j)-tminl - ind = int(dtinvl*(tmp+tepsl)) - if (ind < 0 .or. ind >= table_siz) then - nbad = nbad+1 - else - del = tmp-dtres*real(ind) - esat(i,j) = TABLE(ind+1) + del*(DTABLE(ind+1) + & - del*D2TABLE(ind+1)) - endif - enddo - enddo + + select type (temp) + type is (real(kind=r4_kind)) + select type (esat) + type is (real(kind=r4_kind)) + do j = 1, size(temp,2) + do i = 1, size(temp,1) + tmp = temp(i,j)-tminl + ind = int(dtinvl*(tmp+tepsl)) + if (ind < 0 .or. ind >= table_siz) then + nbad = nbad+1 + else + del = tmp-dtres*real(ind) + esat(i,j) = TABLE(ind+1)+del*(DTABLE(ind+1)+del*D2TABLE(ind+1)) + endif + enddo + enddo + end select + type is (real(kind=r8_kind)) + select type (esat) + type is (real(kind=r8_kind)) + do j = 1, size(temp,2) + do i = 1, size(temp,1) + tmp = temp(i,j)-tminl + ind = int(dtinvl*(tmp+tepsl)) + if (ind < 0 .or. ind >= table_siz) then + nbad = nbad+1 + else + del = tmp-dtres*real(ind) + esat(i,j) = TABLE(ind+1)+del*(DTABLE(ind+1)+del*D2TABLE(ind+1)) + endif + enddo + enddo + end select + end select end subroutine lookup_es_k_2d !####################################################################### subroutine lookup_des_k_1d(temp, desat, nbad) - real, intent(in), dimension(:) :: temp - real, intent(out), dimension(:) :: desat + class(*), intent(in), dimension(:) :: temp + class(*), intent(out), dimension(:) :: desat integer, intent(out) :: nbad real :: tmp, del integer :: ind, i nbad = 0 - do i = 1, size(temp,1) - tmp = temp(i)-tminl - ind = int(dtinvl*(tmp+tepsl)) - if (ind < 0 .or. ind >= table_siz) then - nbad = nbad+1 - else - del = tmp-dtres*real(ind) - desat(i) = DTABLE(ind+1) + 2.*del*D2TABLE(ind+1) - endif - enddo + + select type (temp) + type is (real(kind=r4_kind)) + select type (desat) + type is (real(kind=r4_kind)) + do i = 1, size(temp,1) + tmp = temp(i)-tminl + ind = int(dtinvl*(tmp+tepsl)) + if (ind < 0 .or. ind >= table_siz) then + nbad = nbad+1 + else + del = tmp-dtres*real(ind) + desat(i) = DTABLE(ind+1) + 2.*del*D2TABLE(ind+1) + endif + enddo + end select + type is (real(kind=r8_kind)) + select type (desat) + type is (real(kind=r8_kind)) + do i = 1, size(temp,1) + tmp = temp(i)-tminl + ind = int(dtinvl*(tmp+tepsl)) + if (ind < 0 .or. ind >= table_siz) then + nbad = nbad+1 + else + del = tmp-dtres*real(ind) + desat(i) = DTABLE(ind+1) + 2.*del*D2TABLE(ind+1) + endif + enddo + end select + end select end subroutine lookup_des_k_1d !####################################################################### subroutine lookup_es_k_1d(temp, esat, nbad) - real, intent(in), dimension(:) :: temp - real, intent(out), dimension(:) :: esat + class(*), intent(in), dimension(:) :: temp + class(*), intent(out), dimension(:) :: esat integer, intent(out) :: nbad real :: tmp, del integer :: ind, i nbad = 0 - do i = 1, size(temp,1) - tmp = temp(i)-tminl - ind = int(dtinvl*(tmp+tepsl)) - if (ind < 0 .or. ind >= table_siz) then - nbad = nbad+1 - else - del = tmp-dtres*real(ind) - esat(i) = TABLE(ind+1) + del*(DTABLE(ind+1) + del*D2TABLE(ind+1)) - endif - enddo + + select type (temp) + type is (real(kind=r4_kind)) + select type (esat) + type is (real(kind=r4_kind)) + do i = 1, size(temp,1) + tmp = temp(i)-tminl + ind = int(dtinvl*(tmp+tepsl)) + if (ind < 0 .or. ind >= table_siz) then + nbad = nbad+1 + else + del = tmp-dtres*real(ind) + esat(i) = TABLE(ind+1) + del*(DTABLE(ind+1) + del*D2TABLE(ind+1)) + endif + enddo + end select + type is (real(kind=r8_kind)) + select type (esat) + type is (real(kind=r8_kind)) + do i = 1, size(temp,1) + tmp = temp(i)-tminl + ind = int(dtinvl*(tmp+tepsl)) + if (ind < 0 .or. ind >= table_siz) then + nbad = nbad+1 + else + del = tmp-dtres*real(ind) + esat(i) = TABLE(ind+1) + del*(DTABLE(ind+1) + del*D2TABLE(ind+1)) + endif + enddo + end select + end select end subroutine lookup_es_k_1d !####################################################################### subroutine lookup_des_k_0d(temp, desat, nbad) - real, intent(in) :: temp - real, intent(out) :: desat + class(*), intent(in) :: temp + class(*), intent(out) :: desat integer, intent(out) :: nbad real :: tmp, del integer :: ind nbad = 0 - tmp = temp-tminl + + select type (temp) + type is (real(kind=r4_kind)) + tmp = temp-tminl + type is (real(kind=r8_kind)) + tmp = temp-tminl + end select + ind = int(dtinvl*(tmp+tepsl)) if (ind < 0 .or. ind >= table_siz) then nbad = nbad+1 else del = tmp-dtres*real(ind) - desat = DTABLE(ind+1) + 2.*del*D2TABLE(ind+1) + select type (desat) + type is (real(kind=r4_kind)) + desat = DTABLE(ind+1) + 2.*del*D2TABLE(ind+1) + type is (real(kind=r8_kind)) + desat = DTABLE(ind+1) + 2.*del*D2TABLE(ind+1) + end select endif end subroutine lookup_des_k_0d !####################################################################### subroutine lookup_es_k_0d(temp, esat, nbad) - real, intent(in) :: temp - real, intent(out) :: esat + class(*), intent(in) :: temp + class(*), intent(out) :: esat integer, intent(out) :: nbad real :: tmp, del integer :: ind nbad = 0 - tmp = temp-tminl + + select type (temp) + type is (real(kind=r4_kind)) + tmp = temp-tminl + type is (real(kind=r8_kind)) + tmp = temp-tminl + end select + ind = int(dtinvl*(tmp+tepsl)) if (ind < 0 .or. ind >= table_siz) then nbad = nbad+1 else del = tmp-dtres*real(ind) - esat = TABLE(ind+1) + del*(DTABLE(ind+1) + del*D2TABLE(ind+1)) + select type (esat) + type is (real(kind=r4_kind)) + esat = TABLE(ind+1) + del*(DTABLE(ind+1) + del*D2TABLE(ind+1)) + type is (real(kind=r8_kind)) + esat = TABLE(ind+1) + del*(DTABLE(ind+1) + del*D2TABLE(ind+1)) + end select endif end subroutine lookup_es_k_0d !####################################################################### subroutine lookup_es2_des2_k_3d (temp, esat, desat, nbad) - real, intent(in), dimension(:,:,:) :: temp - real, intent(out), dimension(:,:,:) :: esat, desat + class(*), intent(in), dimension(:,:,:) :: temp + class(*), intent(out), dimension(:,:,:) :: esat, desat integer, intent(out) :: nbad real :: tmp, del integer :: ind, i, j, k nbad = 0 - do k = 1, size(temp,3) - do j = 1, size(temp,2) - do i = 1, size(temp,1) - tmp = temp(i,j,k)-tminl - ind = int(dtinvl*(tmp+tepsl)) - if (ind < 0 .or. ind >= table_siz) then - nbad = nbad+1 - else - del = tmp-dtres*real(ind) - esat(i,j,k) = TABLE2(ind+1) + & - del*(DTABLE2(ind+1) + del*D2TABLE2(ind+1)) - desat(i,j,k) = DTABLE2(ind+1) + 2.*del*D2TABLE2(ind+1) - endif - enddo - enddo - enddo + + select type (temp) + type is (real(kind=r4_kind)) + select type (esat) + type is (real(kind=r4_kind)) + select type (desat) + type is (real(kind=r4_kind)) + do k = 1, size(temp,3) + do j = 1, size(temp,2) + do i = 1, size(temp,1) + tmp = temp(i,j,k)-tminl + ind = int(dtinvl*(tmp+tepsl)) + if (ind < 0 .or. ind >= table_siz) then + nbad = nbad+1 + else + del = tmp-dtres*real(ind) + esat(i,j,k) = TABLE2(ind+1) + del*(DTABLE2(ind+1)+del*D2TABLE2(ind+1)) + desat(i,j,k) = DTABLE2(ind+1) + 2.*del*D2TABLE2(ind+1) + endif + enddo + enddo + enddo + end select + end select + type is (real(kind=r8_kind)) + select type (esat) + type is (real(kind=r8_kind)) + select type (desat) + type is (real(kind=r8_kind)) + do k = 1, size(temp,3) + do j = 1, size(temp,2) + do i = 1, size(temp,1) + tmp = temp(i,j,k)-tminl + ind = int(dtinvl*(tmp+tepsl)) + if (ind < 0 .or. ind >= table_siz) then + nbad = nbad+1 + else + del = tmp-dtres*real(ind) + esat(i,j,k) = TABLE2(ind+1) + del*(DTABLE2(ind+1)+del*D2TABLE2(ind+1)) + desat(i,j,k) = DTABLE2(ind+1) + 2.*del*D2TABLE2(ind+1) + endif + enddo + enddo + enddo + end select + end select + end select end subroutine lookup_es2_des2_k_3d !####################################################################### subroutine lookup_es2_des2_k_2d (temp, esat, desat, nbad) - real, intent(in), dimension(:,:) :: temp - real, intent(out), dimension(:,:) :: esat, desat + class(*), intent(in), dimension(:,:) :: temp + class(*), intent(out), dimension(:,:) :: esat, desat integer, intent(out) :: nbad real :: tmp, del integer :: ind, i, j nbad = 0 - do j = 1, size(temp,2) - do i = 1, size(temp,1) - tmp = temp(i,j)-tminl - ind = int(dtinvl*(tmp+tepsl)) - if (ind < 0 .or. ind >= table_siz) then - nbad = nbad+1 - else - del = tmp-dtres*real(ind) - esat(i,j) = TABLE2(ind+1) + & - del*(DTABLE2(ind+1) + del*D2TABLE2(ind+1)) - desat(i,j) = DTABLE2(ind+1) + 2.*del*D2TABLE2(ind+1) - endif - enddo - enddo + + select type (temp) + type is (real(kind=r4_kind)) + select type (esat) + type is (real(kind=r4_kind)) + select type (desat) + type is (real(kind=r4_kind)) + do j = 1, size(temp,2) + do i = 1, size(temp,1) + tmp = temp(i,j)-tminl + ind = int(dtinvl*(tmp+tepsl)) + if (ind < 0 .or. ind >= table_siz) then + nbad = nbad+1 + else + del = tmp-dtres*real(ind) + esat(i,j) = TABLE2(ind+1) + del*(DTABLE2(ind+1)+del*D2TABLE2(ind+1)) + desat(i,j) = DTABLE2(ind+1) + 2.*del*D2TABLE2(ind+1) + endif + enddo + enddo + end select + end select + type is (real(kind=r8_kind)) + select type (esat) + type is (real(kind=r8_kind)) + select type (desat) + type is (real(kind=r8_kind)) + do j = 1, size(temp,2) + do i = 1, size(temp,1) + tmp = temp(i,j)-tminl + ind = int(dtinvl*(tmp+tepsl)) + if (ind < 0 .or. ind >= table_siz) then + nbad = nbad+1 + else + del = tmp-dtres*real(ind) + esat(i,j) = TABLE2(ind+1) + del*(DTABLE2(ind+1)+del*D2TABLE2(ind+1)) + desat(i,j) = DTABLE2(ind+1) + 2.*del*D2TABLE2(ind+1) + endif + enddo + enddo + end select + end select + end select end subroutine lookup_es2_des2_k_2d !####################################################################### subroutine lookup_es2_des2_k_1d (temp, esat, desat, nbad) - real, intent(in), dimension(:) :: temp - real, intent(out), dimension(:) :: esat, desat + class(*), intent(in), dimension(:) :: temp + class(*), intent(out), dimension(:) :: esat, desat integer, intent(out) :: nbad real :: tmp, del integer :: ind, i nbad = 0 - do i = 1, size(temp,1) - tmp = temp(i)-tminl - ind = int(dtinvl*(tmp+tepsl)) - if (ind < 0 .or. ind >= table_siz) then - nbad = nbad+1 - else - del = tmp-dtres*real(ind) - esat(i) = TABLE2(ind+1) + & - del*(DTABLE2(ind+1) + del*D2TABLE2(ind+1)) - desat(i) = DTABLE2(ind+1) + 2.*del*D2TABLE2(ind+1) - endif - enddo + + select type (temp) + type is (real(kind=r4_kind)) + select type (esat) + type is (real(kind=r4_kind)) + select type (desat) + type is (real(kind=r4_kind)) + do i = 1, size(temp,1) + tmp = temp(i)-tminl + ind = int(dtinvl*(tmp+tepsl)) + if (ind < 0 .or. ind >= table_siz) then + nbad = nbad+1 + else + del = tmp-dtres*real(ind) + esat(i) = TABLE2(ind+1) + del*(DTABLE2(ind+1)+del*D2TABLE2(ind+1)) + desat(i) = DTABLE2(ind+1) + 2.*del*D2TABLE2(ind+1) + endif + enddo + end select + end select + type is (real(kind=r8_kind)) + select type (esat) + type is (real(kind=r8_kind)) + select type (desat) + type is (real(kind=r8_kind)) + do i = 1, size(temp,1) + tmp = temp(i)-tminl + ind = int(dtinvl*(tmp+tepsl)) + if (ind < 0 .or. ind >= table_siz) then + nbad = nbad+1 + else + del = tmp-dtres*real(ind) + esat(i) = TABLE2(ind+1) + del*(DTABLE2(ind+1)+del*D2TABLE2(ind+1)) + desat(i) = DTABLE2(ind+1) + 2.*del*D2TABLE2(ind+1) + endif + enddo + end select + end select + end select end subroutine lookup_es2_des2_k_1d !####################################################################### subroutine lookup_es2_des2_k_0d (temp, esat, desat, nbad) - real, intent(in) :: temp - real, intent(out) :: esat, desat + class(*), intent(in) :: temp + class(*), intent(out) :: esat, desat integer, intent(out) :: nbad real :: tmp, del integer :: ind nbad = 0 - tmp = temp-tminl + + select type (temp) + type is (real(kind=r4_kind)) + tmp = temp-tminl + type is (real(kind=r8_kind)) + tmp = temp-tminl + end select + ind = int(dtinvl*(tmp+tepsl)) if (ind < 0 .or. ind >= table_siz) then nbad = nbad+1 else del = tmp-dtres*real(ind) - esat = TABLE2(ind+1) + & - del*(DTABLE2(ind+1) + del*D2TABLE2(ind+1)) - desat = DTABLE2(ind+1) + 2.*del*D2TABLE2(ind+1) + select type (esat) + type is (real(kind=r4_kind)) + esat = TABLE2(ind+1) + del*(DTABLE2(ind+1)+del*D2TABLE2(ind+1)) + type is (real(kind=r8_kind)) + esat = TABLE2(ind+1) + del*(DTABLE2(ind+1)+del*D2TABLE2(ind+1)) + end select + + select type (desat) + type is (real(kind=r4_kind)) + desat = DTABLE2(ind+1) + 2.*del*D2TABLE2(ind+1) + type is (real(kind=r8_kind)) + desat = DTABLE2(ind+1) + 2.*del*D2TABLE2(ind+1) + end select endif end subroutine lookup_es2_des2_k_0d @@ -1546,182 +2212,342 @@ end subroutine lookup_es2_des2_k_0d !####################################################################### subroutine lookup_es2_k_3d(temp, esat, nbad) - real, intent(in), dimension(:,:,:) :: temp - real, intent(out), dimension(:,:,:) :: esat + class(*), intent(in), dimension(:,:,:) :: temp + class(*), intent(out), dimension(:,:,:) :: esat integer, intent(out) :: nbad real :: tmp, del integer :: ind, i, j, k nbad = 0 - do k = 1, size(temp,3) - do j = 1, size(temp,2) - do i = 1, size(temp,1) - tmp = temp(i,j,k)-tminl - ind = int(dtinvl*(tmp+tepsl)) - if (ind < 0 .or. ind >= table_siz) then - nbad = nbad+1 - else - del = tmp-dtres*real(ind) - esat(i,j,k) = TABLE2(ind+1) + & - del*(DTABLE2(ind+1) + del*D2TABLE2(ind+1)) - endif - enddo - enddo - enddo + + select type (temp) + type is (real(kind=r4_kind)) + select type (esat) + type is (real(kind=r4_kind)) + do k = 1, size(temp,3) + do j = 1, size(temp,2) + do i = 1, size(temp,1) + tmp = temp(i,j,k)-tminl + ind = int(dtinvl*(tmp+tepsl)) + if (ind < 0 .or. ind >= table_siz) then + nbad = nbad+1 + else + del = tmp-dtres*real(ind) + esat(i,j,k) = TABLE2(ind+1)+del*(DTABLE2(ind+1)+del*D2TABLE2(ind+1)) + endif + enddo + enddo + enddo + end select + type is (real(kind=r8_kind)) + select type (esat) + type is (real(kind=r8_kind)) + do k = 1, size(temp,3) + do j = 1, size(temp,2) + do i = 1, size(temp,1) + tmp = temp(i,j,k)-tminl + ind = int(dtinvl*(tmp+tepsl)) + if (ind < 0 .or. ind >= table_siz) then + nbad = nbad+1 + else + del = tmp-dtres*real(ind) + esat(i,j,k) = TABLE2(ind+1)+del*(DTABLE2(ind+1)+del*D2TABLE2(ind+1)) + endif + enddo + enddo + enddo + end select + end select end subroutine lookup_es2_k_3d !####################################################################### subroutine lookup_des2_k_3d(temp, desat, nbad) - real, intent(in), dimension(:,:,:) :: temp - real, intent(out), dimension(:,:,:) :: desat + class(*), intent(in), dimension(:,:,:) :: temp + class(*), intent(out), dimension(:,:,:) :: desat integer, intent(out) :: nbad real :: tmp, del integer :: ind, i, j, k nbad = 0 - do k = 1, size(temp,3) - do j = 1, size(temp,2) - do i = 1, size(temp,1) - tmp = temp(i,j,k)-tminl - ind = int(dtinvl*(tmp+tepsl)) - if (ind < 0 .or. ind >= table_siz) then - nbad = nbad+1 - else - del = tmp-dtres*real(ind) - desat(i,j,k) = DTABLE2(ind+1) + 2.*del*D2TABLE2(ind+1) - endif - enddo - enddo - enddo + + select type (temp) + type is (real(kind=r4_kind)) + select type (desat) + type is (real(kind=r4_kind)) + do k = 1, size(temp,3) + do j = 1, size(temp,2) + do i = 1, size(temp,1) + tmp = temp(i,j,k)-tminl + ind = int(dtinvl*(tmp+tepsl)) + if (ind < 0 .or. ind >= table_siz) then + nbad = nbad+1 + else + del = tmp-dtres*real(ind) + desat(i,j,k) = DTABLE2(ind+1) + 2.*del*D2TABLE2(ind+1) + endif + enddo + enddo + enddo + end select + type is (real(kind=r8_kind)) + select type (desat) + type is (real(kind=r8_kind)) + do k = 1, size(temp,3) + do j = 1, size(temp,2) + do i = 1, size(temp,1) + tmp = temp(i,j,k)-tminl + ind = int(dtinvl*(tmp+tepsl)) + if (ind < 0 .or. ind >= table_siz) then + nbad = nbad+1 + else + del = tmp-dtres*real(ind) + desat(i,j,k) = DTABLE2(ind+1) + 2.*del*D2TABLE2(ind+1) + endif + enddo + enddo + enddo + end select + end select end subroutine lookup_des2_k_3d !####################################################################### subroutine lookup_des2_k_2d(temp, desat, nbad) - real, intent(in), dimension(:,:) :: temp - real, intent(out), dimension(:,:) :: desat + class(*), intent(in), dimension(:,:) :: temp + class(*), intent(out), dimension(:,:) :: desat integer, intent(out) :: nbad real :: tmp, del integer :: ind, i, j nbad = 0 - do j = 1, size(temp,2) - do i = 1, size(temp,1) - tmp = temp(i,j)-tminl - ind = int(dtinvl*(tmp+tepsl)) - if (ind < 0 .or. ind >= table_siz) then - nbad = nbad+1 - else - del = tmp-dtres*real(ind) - desat(i,j) = DTABLE2(ind+1) + 2.*del*D2TABLE2(ind+1) - endif - enddo - enddo + + select type (temp) + type is (real(kind=r4_kind)) + select type (desat) + type is (real(kind=r4_kind)) + do j = 1, size(temp,2) + do i = 1, size(temp,1) + tmp = temp(i,j)-tminl + ind = int(dtinvl*(tmp+tepsl)) + if (ind < 0 .or. ind >= table_siz) then + nbad = nbad+1 + else + del = tmp-dtres*real(ind) + desat(i,j) = DTABLE2(ind+1) + 2.*del*D2TABLE2(ind+1) + endif + enddo + enddo + end select + type is (real(kind=r8_kind)) + select type (desat) + type is (real(kind=r8_kind)) + do j = 1, size(temp,2) + do i = 1, size(temp,1) + tmp = temp(i,j)-tminl + ind = int(dtinvl*(tmp+tepsl)) + if (ind < 0 .or. ind >= table_siz) then + nbad = nbad+1 + else + del = tmp-dtres*real(ind) + desat(i,j) = DTABLE2(ind+1) + 2.*del*D2TABLE2(ind+1) + endif + enddo + enddo + end select + end select end subroutine lookup_des2_k_2d !####################################################################### subroutine lookup_es2_k_2d(temp, esat, nbad) - real, intent(in), dimension(:,:) :: temp - real, intent(out), dimension(:,:) :: esat + class(*), intent(in), dimension(:,:) :: temp + class(*), intent(out), dimension(:,:) :: esat integer, intent(out) :: nbad real :: tmp, del integer :: ind, i, j nbad = 0 - do j = 1, size(temp,2) - do i = 1, size(temp,1) - tmp = temp(i,j)-tminl - ind = int(dtinvl*(tmp+tepsl)) - if (ind < 0 .or. ind >= table_siz) then - nbad = nbad+1 - else - del = tmp-dtres*real(ind) - esat(i,j) = TABLE2(ind+1) + del*(DTABLE2(ind+1) + & - del*D2TABLE2(ind+1)) - endif - enddo - enddo + + select type (temp) + type is (real(kind=r4_kind)) + select type (esat) + type is (real(kind=r4_kind)) + do j = 1, size(temp,2) + do i = 1, size(temp,1) + tmp = temp(i,j)-tminl + ind = int(dtinvl*(tmp+tepsl)) + if (ind < 0 .or. ind >= table_siz) then + nbad = nbad+1 + else + del = tmp-dtres*real(ind) + esat(i,j) = TABLE2(ind+1)+del*(DTABLE2(ind+1)+del*D2TABLE2(ind+1)) + endif + enddo + enddo + end select + type is (real(kind=r8_kind)) + select type (esat) + type is (real(kind=r8_kind)) + do j = 1, size(temp,2) + do i = 1, size(temp,1) + tmp = temp(i,j)-tminl + ind = int(dtinvl*(tmp+tepsl)) + if (ind < 0 .or. ind >= table_siz) then + nbad = nbad+1 + else + del = tmp-dtres*real(ind) + esat(i,j) = TABLE2(ind+1)+del*(DTABLE2(ind+1)+del*D2TABLE2(ind+1)) + endif + enddo + enddo + end select + end select end subroutine lookup_es2_k_2d !####################################################################### subroutine lookup_des2_k_1d(temp, desat, nbad) - real, intent(in), dimension(:) :: temp - real, intent(out), dimension(:) :: desat + class(*), intent(in), dimension(:) :: temp + class(*), intent(out), dimension(:) :: desat integer, intent(out) :: nbad real :: tmp, del integer :: ind, i nbad = 0 - do i = 1, size(temp,1) - tmp = temp(i)-tminl - ind = int(dtinvl*(tmp+tepsl)) - if (ind < 0 .or. ind >= table_siz) then - nbad = nbad+1 - else - del = tmp-dtres*real(ind) - desat(i) = DTABLE2(ind+1) + 2.*del*D2TABLE2(ind+1) - endif - enddo + + select type (temp) + type is (real(kind=r4_kind)) + select type (desat) + type is (real(kind=r4_kind)) + do i = 1, size(temp,1) + tmp = temp(i)-tminl + ind = int(dtinvl*(tmp+tepsl)) + if (ind < 0 .or. ind >= table_siz) then + nbad = nbad+1 + else + del = tmp-dtres*real(ind) + desat(i) = DTABLE2(ind+1) + 2.*del*D2TABLE2(ind+1) + endif + enddo + end select + type is (real(kind=r8_kind)) + select type (desat) + type is (real(kind=r8_kind)) + do i = 1, size(temp,1) + tmp = temp(i)-tminl + ind = int(dtinvl*(tmp+tepsl)) + if (ind < 0 .or. ind >= table_siz) then + nbad = nbad+1 + else + del = tmp-dtres*real(ind) + desat(i) = DTABLE2(ind+1) + 2.*del*D2TABLE2(ind+1) + endif + enddo + end select + end select end subroutine lookup_des2_k_1d !####################################################################### subroutine lookup_es2_k_1d(temp, esat, nbad) - real, intent(in), dimension(:) :: temp - real, intent(out), dimension(:) :: esat + class(*), intent(in), dimension(:) :: temp + class(*), intent(out), dimension(:) :: esat integer, intent(out) :: nbad real :: tmp, del integer :: ind, i nbad = 0 - do i = 1, size(temp,1) - tmp = temp(i)-tminl - ind = int(dtinvl*(tmp+tepsl)) - if (ind < 0 .or. ind >= table_siz) then - nbad = nbad+1 - else - del = tmp-dtres*real(ind) - esat(i) = TABLE2(ind+1) + del*(DTABLE2(ind+1) + del*D2TABLE2(ind+1)) - endif - enddo + select type (temp) + type is (real(kind=r4_kind)) + select type (esat) + type is (real(kind=r4_kind)) + do i = 1, size(temp,1) + tmp = temp(i)-tminl + ind = int(dtinvl*(tmp+tepsl)) + if (ind < 0 .or. ind >= table_siz) then + nbad = nbad+1 + else + del = tmp-dtres*real(ind) + esat(i) = TABLE2(ind+1) + del*(DTABLE2(ind+1) + del*D2TABLE2(ind+1)) + endif + enddo + end select + type is (real(kind=r8_kind)) + select type (esat) + type is (real(kind=r8_kind)) + do i = 1, size(temp,1) + tmp = temp(i)-tminl + ind = int(dtinvl*(tmp+tepsl)) + if (ind < 0 .or. ind >= table_siz) then + nbad = nbad+1 + else + del = tmp-dtres*real(ind) + esat(i) = TABLE2(ind+1) + del*(DTABLE2(ind+1) + del*D2TABLE2(ind+1)) + endif + enddo + end select + end select + end subroutine lookup_es2_k_1d !####################################################################### subroutine lookup_des2_k_0d(temp, desat, nbad) - real, intent(in) :: temp - real, intent(out) :: desat + class(*), intent(in) :: temp + class(*), intent(out) :: desat integer, intent(out) :: nbad real :: tmp, del integer :: ind nbad = 0 - tmp = temp-tminl + + select type (temp) + type is (real(kind=r4_kind)) + tmp = temp-tminl + type is (real(kind=r8_kind)) + tmp = temp-tminl + end select + ind = int(dtinvl*(tmp+tepsl)) if (ind < 0 .or. ind >= table_siz) then nbad = nbad+1 else del = tmp-dtres*real(ind) - desat = DTABLE2(ind+1) + 2.*del*D2TABLE2(ind+1) + select type (desat) + type is (real(kind=r4_kind)) + desat = DTABLE2(ind+1) + 2.*del*D2TABLE2(ind+1) + type is (real(kind=r8_kind)) + desat = DTABLE2(ind+1) + 2.*del*D2TABLE2(ind+1) + end select endif end subroutine lookup_des2_k_0d !####################################################################### subroutine lookup_es2_k_0d(temp, esat, nbad) - real, intent(in) :: temp - real, intent(out) :: esat + class(*), intent(in) :: temp + class(*), intent(out) :: esat integer, intent(out) :: nbad real :: tmp, del integer :: ind nbad = 0 - tmp = temp-tminl + + select type (temp) + type is (real(kind=r4_kind)) + tmp = temp-tminl + type is (real(kind=r8_kind)) + tmp = temp-tminl + end select + ind = int(dtinvl*(tmp+tepsl)) if (ind < 0 .or. ind >= table_siz) then nbad = nbad+1 else del = tmp-dtres*real(ind) - esat = TABLE2(ind+1) + del*(DTABLE2(ind+1) + del*D2TABLE2(ind+1)) + select type (esat) + type is (real(kind=r4_kind)) + esat = TABLE2(ind+1) + del*(DTABLE2(ind+1) + del*D2TABLE2(ind+1)) + type is (real(kind=r8_kind)) + esat = TABLE2(ind+1) + del*(DTABLE2(ind+1) + del*D2TABLE2(ind+1)) + end select endif end subroutine lookup_es2_k_0d @@ -1730,107 +2556,211 @@ end subroutine lookup_es2_k_0d !####################################################################### subroutine lookup_es3_des3_k_3d (temp, esat, desat, nbad) - real, intent(in), dimension(:,:,:) :: temp - real, intent(out), dimension(:,:,:) :: esat, desat + class(*), intent(in), dimension(:,:,:) :: temp + class(*), intent(out), dimension(:,:,:) :: esat, desat integer, intent(out) :: nbad real :: tmp, del integer :: ind, i, j, k nbad = 0 - do k = 1, size(temp,3) - do j = 1, size(temp,2) - do i = 1, size(temp,1) - tmp = temp(i,j,k)-tminl - ind = int(dtinvl*(tmp+tepsl)) - if (ind < 0 .or. ind >= table_siz) then - nbad = nbad+1 - else - del = tmp-dtres*real(ind) - esat(i,j,k) = TABLE3(ind+1) + & - del*(DTABLE3(ind+1) + del*D2TABLE3(ind+1)) - desat(i,j,k) = DTABLE3(ind+1) + 2.*del*D2TABLE3(ind+1) - endif - enddo - enddo - enddo + + select type (temp) + type is (real(kind=r4_kind)) + select type (esat) + type is (real(kind=r4_kind)) + select type (desat) + type is (real(kind=r4_kind)) + do k = 1, size(temp,3) + do j = 1, size(temp,2) + do i = 1, size(temp,1) + tmp = temp(i,j,k)-tminl + ind = int(dtinvl*(tmp+tepsl)) + if (ind < 0 .or. ind >= table_siz) then + nbad = nbad+1 + else + del = tmp-dtres*real(ind) + esat(i,j,k) = TABLE3(ind+1) + del*(DTABLE3(ind+1)+del*D2TABLE3(ind+1)) + desat(i,j,k) = DTABLE3(ind+1) + 2.*del*D2TABLE3(ind+1) + endif + enddo + enddo + enddo + end select + end select + type is (real(kind=r8_kind)) + select type (esat) + type is (real(kind=r8_kind)) + select type (desat) + type is (real(kind=r8_kind)) + do k = 1, size(temp,3) + do j = 1, size(temp,2) + do i = 1, size(temp,1) + tmp = temp(i,j,k)-tminl + ind = int(dtinvl*(tmp+tepsl)) + if (ind < 0 .or. ind >= table_siz) then + nbad = nbad+1 + else + del = tmp-dtres*real(ind) + esat(i,j,k) = TABLE3(ind+1) + del*(DTABLE3(ind+1)+del*D2TABLE3(ind+1)) + desat(i,j,k) = DTABLE3(ind+1) + 2.*del*D2TABLE3(ind+1) + endif + enddo + enddo + enddo + end select + end select + end select end subroutine lookup_es3_des3_k_3d !####################################################################### subroutine lookup_es3_des3_k_2d (temp, esat, desat, nbad) - real, intent(in), dimension(:,:) :: temp - real, intent(out), dimension(:,:) :: esat, desat + class(*), intent(in), dimension(:,:) :: temp + class(*), intent(out), dimension(:,:) :: esat, desat integer, intent(out) :: nbad real :: tmp, del integer :: ind, i, j nbad = 0 - do j = 1, size(temp,2) - do i = 1, size(temp,1) - tmp = temp(i,j)-tminl - ind = int(dtinvl*(tmp+tepsl)) - if (ind < 0 .or. ind >= table_siz) then - nbad = nbad+1 - else - del = tmp-dtres*real(ind) - esat(i,j) = TABLE3(ind+1) + & - del*(DTABLE3(ind+1) + del*D2TABLE3(ind+1)) - desat(i,j) = DTABLE3(ind+1) + 2.*del*D2TABLE3(ind+1) - endif - enddo - enddo + + select type (temp) + type is (real(kind=r4_kind)) + select type (esat) + type is (real(kind=r4_kind)) + select type (desat) + type is (real(kind=r4_kind)) + do j = 1, size(temp,2) + do i = 1, size(temp,1) + tmp = temp(i,j)-tminl + ind = int(dtinvl*(tmp+tepsl)) + if (ind < 0 .or. ind >= table_siz) then + nbad = nbad+1 + else + del = tmp-dtres*real(ind) + esat(i,j) = TABLE3(ind+1) + del*(DTABLE3(ind+1)+del*D2TABLE3(ind+1)) + desat(i,j) = DTABLE3(ind+1) + 2.*del*D2TABLE3(ind+1) + endif + enddo + enddo + end select + end select + type is (real(kind=r8_kind)) + select type (esat) + type is (real(kind=r8_kind)) + select type (desat) + type is (real(kind=r8_kind)) + do j = 1, size(temp,2) + do i = 1, size(temp,1) + tmp = temp(i,j)-tminl + ind = int(dtinvl*(tmp+tepsl)) + if (ind < 0 .or. ind >= table_siz) then + nbad = nbad+1 + else + del = tmp-dtres*real(ind) + esat(i,j) = TABLE3(ind+1) + del*(DTABLE3(ind+1)+del*D2TABLE3(ind+1)) + desat(i,j) = DTABLE3(ind+1) + 2.*del*D2TABLE3(ind+1) + endif + enddo + enddo + end select + end select + end select end subroutine lookup_es3_des3_k_2d !####################################################################### subroutine lookup_es3_des3_k_1d (temp, esat, desat, nbad) - real, intent(in), dimension(:) :: temp - real, intent(out), dimension(:) :: esat, desat + class(*), intent(in), dimension(:) :: temp + class(*), intent(out), dimension(:) :: esat, desat integer, intent(out) :: nbad real :: tmp, del integer :: ind, i nbad = 0 - do i = 1, size(temp,1) - tmp = temp(i)-tminl - ind = int(dtinvl*(tmp+tepsl)) - if (ind < 0 .or. ind >= table_siz) then - nbad = nbad+1 - else - del = tmp-dtres*real(ind) - esat(i) = TABLE3(ind+1) + & - del*(DTABLE3(ind+1) + del*D2TABLE3(ind+1)) - desat(i) = DTABLE3(ind+1) + 2.*del*D2TABLE3(ind+1) - endif - enddo + + select type (temp) + type is (real(kind=r4_kind)) + select type (esat) + type is (real(kind=r4_kind)) + select type (desat) + type is (real(kind=r4_kind)) + do i = 1, size(temp,1) + tmp = temp(i)-tminl + ind = int(dtinvl*(tmp+tepsl)) + if (ind < 0 .or. ind >= table_siz) then + nbad = nbad+1 + else + del = tmp-dtres*real(ind) + esat(i) = TABLE3(ind+1) + del*(DTABLE3(ind+1)+del*D2TABLE3(ind+1)) + desat(i) = DTABLE3(ind+1) + 2.*del*D2TABLE3(ind+1) + endif + enddo + end select + end select + type is (real(kind=r8_kind)) + select type (esat) + type is (real(kind=r8_kind)) + select type (desat) + type is (real(kind=r8_kind)) + do i = 1, size(temp,1) + tmp = temp(i)-tminl + ind = int(dtinvl*(tmp+tepsl)) + if (ind < 0 .or. ind >= table_siz) then + nbad = nbad+1 + else + del = tmp-dtres*real(ind) + esat(i) = TABLE3(ind+1) + del*(DTABLE3(ind+1)+del*D2TABLE3(ind+1)) + desat(i) = DTABLE3(ind+1) + 2.*del*D2TABLE3(ind+1) + endif + enddo + end select + end select + end select end subroutine lookup_es3_des3_k_1d !####################################################################### subroutine lookup_es3_des3_k_0d (temp, esat, desat, nbad) - real, intent(in) :: temp - real, intent(out) :: esat, desat + class(*), intent(in) :: temp + class(*), intent(out) :: esat, desat integer, intent(out) :: nbad real :: tmp, del integer :: ind nbad = 0 - tmp = temp-tminl + + select type (temp) + type is (real(kind=r4_kind)) + tmp = temp-tminl + type is (real(kind=r8_kind)) + tmp = temp-tminl + end select + ind = int(dtinvl*(tmp+tepsl)) if (ind < 0 .or. ind >= table_siz) then nbad = nbad+1 else del = tmp-dtres*real(ind) - esat = TABLE3(ind+1) + & - del*(DTABLE3(ind+1) + del*D2TABLE3(ind+1)) - desat = DTABLE3(ind+1) + 2.*del*D2TABLE3(ind+1) + select type (esat) + type is (real(kind=r4_kind)) + esat = TABLE3(ind+1) + del*(DTABLE3(ind+1)+del*D2TABLE3(ind+1)) + type is (real(kind=r8_kind)) + esat = TABLE3(ind+1) + del*(DTABLE3(ind+1)+del*D2TABLE3(ind+1)) + end select + + select type (desat) + type is (real(kind=r4_kind)) + desat = DTABLE3(ind+1) + 2.*del*D2TABLE3(ind+1) + type is (real(kind=r8_kind)) + desat = DTABLE3(ind+1) + 2.*del*D2TABLE3(ind+1) + end select endif end subroutine lookup_es3_des3_k_0d @@ -1838,182 +2768,342 @@ end subroutine lookup_es3_des3_k_0d !####################################################################### subroutine lookup_es3_k_3d(temp, esat, nbad) - real, intent(in), dimension(:,:,:) :: temp - real, intent(out), dimension(:,:,:) :: esat + class(*), intent(in), dimension(:,:,:) :: temp + class(*), intent(out), dimension(:,:,:) :: esat integer, intent(out) :: nbad real :: tmp, del integer :: ind, i, j, k nbad = 0 - do k = 1, size(temp,3) - do j = 1, size(temp,2) - do i = 1, size(temp,1) - tmp = temp(i,j,k)-tminl - ind = int(dtinvl*(tmp+tepsl)) - if (ind < 0 .or. ind >= table_siz) then - nbad = nbad+1 - else - del = tmp-dtres*real(ind) - esat(i,j,k) = TABLE3(ind+1) + & - del*(DTABLE3(ind+1) + del*D2TABLE3(ind+1)) - endif - enddo - enddo - enddo + + select type (temp) + type is (real(kind=r4_kind)) + select type (esat) + type is (real(kind=r4_kind)) + do k = 1, size(temp,3) + do j = 1, size(temp,2) + do i = 1, size(temp,1) + tmp = temp(i,j,k)-tminl + ind = int(dtinvl*(tmp+tepsl)) + if (ind < 0 .or. ind >= table_siz) then + nbad = nbad+1 + else + del = tmp-dtres*real(ind) + esat(i,j,k) = TABLE3(ind+1)+del*(DTABLE3(ind+1)+del*D2TABLE3(ind+1)) + endif + enddo + enddo + enddo + end select + type is (real(kind=r8_kind)) + select type (esat) + type is (real(kind=r8_kind)) + do k = 1, size(temp,3) + do j = 1, size(temp,2) + do i = 1, size(temp,1) + tmp = temp(i,j,k)-tminl + ind = int(dtinvl*(tmp+tepsl)) + if (ind < 0 .or. ind >= table_siz) then + nbad = nbad+1 + else + del = tmp-dtres*real(ind) + esat(i,j,k) = TABLE3(ind+1)+del*(DTABLE3(ind+1)+del*D2TABLE3(ind+1)) + endif + enddo + enddo + enddo + end select + end select end subroutine lookup_es3_k_3d !####################################################################### subroutine lookup_des3_k_3d(temp, desat, nbad) - real, intent(in), dimension(:,:,:) :: temp - real, intent(out), dimension(:,:,:) :: desat + class(*), intent(in), dimension(:,:,:) :: temp + class(*), intent(out), dimension(:,:,:) :: desat integer, intent(out) :: nbad real :: tmp, del integer :: ind, i, j, k nbad = 0 - do k = 1, size(temp,3) - do j = 1, size(temp,2) - do i = 1, size(temp,1) - tmp = temp(i,j,k)-tminl - ind = int(dtinvl*(tmp+tepsl)) - if (ind < 0 .or. ind >= table_siz) then - nbad = nbad+1 - else - del = tmp-dtres*real(ind) - desat(i,j,k) = DTABLE3(ind+1) + 2.*del*D2TABLE3(ind+1) - endif - enddo - enddo - enddo + + select type (temp) + type is (real(kind=r4_kind)) + select type (desat) + type is (real(kind=r4_kind)) + do k = 1, size(temp,3) + do j = 1, size(temp,2) + do i = 1, size(temp,1) + tmp = temp(i,j,k)-tminl + ind = int(dtinvl*(tmp+tepsl)) + if (ind < 0 .or. ind >= table_siz) then + nbad = nbad+1 + else + del = tmp-dtres*real(ind) + desat(i,j,k) = DTABLE3(ind+1) + 2.*del*D2TABLE3(ind+1) + endif + enddo + enddo + enddo + end select + type is (real(kind=r8_kind)) + select type (desat) + type is (real(kind=r8_kind)) + do k = 1, size(temp,3) + do j = 1, size(temp,2) + do i = 1, size(temp,1) + tmp = temp(i,j,k)-tminl + ind = int(dtinvl*(tmp+tepsl)) + if (ind < 0 .or. ind >= table_siz) then + nbad = nbad+1 + else + del = tmp-dtres*real(ind) + desat(i,j,k) = DTABLE3(ind+1) + 2.*del*D2TABLE3(ind+1) + endif + enddo + enddo + enddo + end select + end select end subroutine lookup_des3_k_3d !####################################################################### subroutine lookup_des3_k_2d(temp, desat, nbad) - real, intent(in), dimension(:,:) :: temp - real, intent(out), dimension(:,:) :: desat + class(*), intent(in), dimension(:,:) :: temp + class(*), intent(out), dimension(:,:) :: desat integer, intent(out) :: nbad real :: tmp, del integer :: ind, i, j nbad = 0 - do j = 1, size(temp,2) - do i = 1, size(temp,1) - tmp = temp(i,j)-tminl - ind = int(dtinvl*(tmp+tepsl)) - if (ind < 0 .or. ind >= table_siz) then - nbad = nbad+1 - else - del = tmp-dtres*real(ind) - desat(i,j) = DTABLE3(ind+1) + 2.*del*D2TABLE3(ind+1) - endif - enddo - enddo + + select type (temp) + type is (real(kind=r4_kind)) + select type (desat) + type is (real(kind=r4_kind)) + do j = 1, size(temp,2) + do i = 1, size(temp,1) + tmp = temp(i,j)-tminl + ind = int(dtinvl*(tmp+tepsl)) + if (ind < 0 .or. ind >= table_siz) then + nbad = nbad+1 + else + del = tmp-dtres*real(ind) + desat(i,j) = DTABLE3(ind+1) + 2.*del*D2TABLE3(ind+1) + endif + enddo + enddo + end select + type is (real(kind=r8_kind)) + select type (desat) + type is (real(kind=r8_kind)) + do j = 1, size(temp,2) + do i = 1, size(temp,1) + tmp = temp(i,j)-tminl + ind = int(dtinvl*(tmp+tepsl)) + if (ind < 0 .or. ind >= table_siz) then + nbad = nbad+1 + else + del = tmp-dtres*real(ind) + desat(i,j) = DTABLE3(ind+1) + 2.*del*D2TABLE3(ind+1) + endif + enddo + enddo + end select + end select end subroutine lookup_des3_k_2d !####################################################################### subroutine lookup_es3_k_2d(temp, esat, nbad) - real, intent(in), dimension(:,:) :: temp - real, intent(out), dimension(:,:) :: esat + class(*), intent(in), dimension(:,:) :: temp + class(*), intent(out), dimension(:,:) :: esat integer, intent(out) :: nbad real :: tmp, del integer :: ind, i, j nbad = 0 - do j = 1, size(temp,2) - do i = 1, size(temp,1) - tmp = temp(i,j)-tminl - ind = int(dtinvl*(tmp+tepsl)) - if (ind < 0 .or. ind >= table_siz) then - nbad = nbad+1 - else - del = tmp-dtres*real(ind) - esat(i,j) = TABLE3(ind+1) + del*(DTABLE3(ind+1) + & - del*D2TABLE3(ind+1)) - endif - enddo - enddo + + select type (temp) + type is (real(kind=r4_kind)) + select type (esat) + type is (real(kind=r4_kind)) + do j = 1, size(temp,2) + do i = 1, size(temp,1) + tmp = temp(i,j)-tminl + ind = int(dtinvl*(tmp+tepsl)) + if (ind < 0 .or. ind >= table_siz) then + nbad = nbad+1 + else + del = tmp-dtres*real(ind) + esat(i,j) = TABLE3(ind+1)+del*(DTABLE3(ind+1)+del*D2TABLE3(ind+1)) + endif + enddo + enddo + end select + type is (real(kind=r8_kind)) + select type (esat) + type is (real(kind=r8_kind)) + do j = 1, size(temp,2) + do i = 1, size(temp,1) + tmp = temp(i,j)-tminl + ind = int(dtinvl*(tmp+tepsl)) + if (ind < 0 .or. ind >= table_siz) then + nbad = nbad+1 + else + del = tmp-dtres*real(ind) + esat(i,j) = TABLE3(ind+1)+del*(DTABLE3(ind+1)+del*D2TABLE3(ind+1)) + endif + enddo + enddo + end select + end select end subroutine lookup_es3_k_2d !####################################################################### subroutine lookup_des3_k_1d(temp, desat, nbad) - real, intent(in), dimension(:) :: temp - real, intent(out), dimension(:) :: desat + class(*), intent(in), dimension(:) :: temp + class(*), intent(out), dimension(:) :: desat integer, intent(out) :: nbad real :: tmp, del integer :: ind, i nbad = 0 - do i = 1, size(temp,1) - tmp = temp(i)-tminl - ind = int(dtinvl*(tmp+tepsl)) - if (ind < 0 .or. ind >= table_siz) then - nbad = nbad+1 - else - del = tmp-dtres*real(ind) - desat(i) = DTABLE3(ind+1) + 2.*del*D2TABLE3(ind+1) - endif - enddo + + select type (temp) + type is (real(kind=r4_kind)) + select type (desat) + type is (real(kind=r4_kind)) + do i = 1, size(temp,1) + tmp = temp(i)-tminl + ind = int(dtinvl*(tmp+tepsl)) + if (ind < 0 .or. ind >= table_siz) then + nbad = nbad+1 + else + del = tmp-dtres*real(ind) + desat(i) = DTABLE3(ind+1) + 2.*del*D2TABLE3(ind+1) + endif + enddo + end select + type is (real(kind=r8_kind)) + select type (desat) + type is (real(kind=r8_kind)) + do i = 1, size(temp,1) + tmp = temp(i)-tminl + ind = int(dtinvl*(tmp+tepsl)) + if (ind < 0 .or. ind >= table_siz) then + nbad = nbad+1 + else + del = tmp-dtres*real(ind) + desat(i) = DTABLE3(ind+1) + 2.*del*D2TABLE3(ind+1) + endif + enddo + end select + end select end subroutine lookup_des3_k_1d !####################################################################### subroutine lookup_es3_k_1d(temp, esat, nbad) - real, intent(in), dimension(:) :: temp - real, intent(out), dimension(:) :: esat + class(*), intent(in), dimension(:) :: temp + class(*), intent(out), dimension(:) :: esat integer, intent(out) :: nbad real :: tmp, del integer :: ind, i nbad = 0 - do i = 1, size(temp,1) - tmp = temp(i)-tminl - ind = int(dtinvl*(tmp+tepsl)) - if (ind < 0 .or. ind >= table_siz) then - nbad = nbad+1 - else - del = tmp-dtres*real(ind) - esat(i) = TABLE3(ind+1) + del*(DTABLE3(ind+1) + del*D2TABLE3(ind+1)) - endif - enddo + + select type (temp) + type is (real(kind=r4_kind)) + select type (esat) + type is (real(kind=r4_kind)) + do i = 1, size(temp,1) + tmp = temp(i)-tminl + ind = int(dtinvl*(tmp+tepsl)) + if (ind < 0 .or. ind >= table_siz) then + nbad = nbad+1 + else + del = tmp-dtres*real(ind) + esat(i) = TABLE3(ind+1) + del*(DTABLE3(ind+1) + del*D2TABLE3(ind+1)) + endif + enddo + end select + type is (real(kind=r8_kind)) + select type (esat) + type is (real(kind=r8_kind)) + do i = 1, size(temp,1) + tmp = temp(i)-tminl + ind = int(dtinvl*(tmp+tepsl)) + if (ind < 0 .or. ind >= table_siz) then + nbad = nbad+1 + else + del = tmp-dtres*real(ind) + esat(i) = TABLE3(ind+1) + del*(DTABLE3(ind+1) + del*D2TABLE3(ind+1)) + endif + enddo + end select + end select end subroutine lookup_es3_k_1d !####################################################################### subroutine lookup_des3_k_0d(temp, desat, nbad) - real, intent(in) :: temp - real, intent(out) :: desat + class(*), intent(in) :: temp + class(*), intent(out) :: desat integer, intent(out) :: nbad real :: tmp, del integer :: ind nbad = 0 - tmp = temp-tminl + + select type (temp) + type is (real(kind=r4_kind)) + tmp = temp-tminl + type is (real(kind=r8_kind)) + tmp = temp-tminl + end select + ind = int(dtinvl*(tmp+tepsl)) if (ind < 0 .or. ind >= table_siz) then nbad = nbad+1 else del = tmp-dtres*real(ind) - desat = DTABLE3(ind+1) + 2.*del*D2TABLE3(ind+1) + select type (desat) + type is (real(kind=r4_kind)) + desat = DTABLE3(ind+1) + 2.*del*D2TABLE3(ind+1) + type is (real(kind=r8_kind)) + desat = DTABLE3(ind+1) + 2.*del*D2TABLE3(ind+1) + end select endif end subroutine lookup_des3_k_0d !####################################################################### subroutine lookup_es3_k_0d(temp, esat, nbad) - real, intent(in) :: temp - real, intent(out) :: esat + class(*), intent(in) :: temp + class(*), intent(out) :: esat integer, intent(out) :: nbad real :: tmp, del integer :: ind nbad = 0 - tmp = temp-tminl + + select type (temp) + type is (real(kind=r4_kind)) + tmp = temp-tminl + type is (real(kind=r8_kind)) + tmp = temp-tminl + end select + ind = int(dtinvl*(tmp+tepsl)) if (ind < 0 .or. ind >= table_siz) then nbad = nbad+1 else del = tmp-dtres*real(ind) - esat = TABLE3(ind+1) + del*(DTABLE3(ind+1) + del*D2TABLE3(ind+1)) + select type (esat) + type is (real(kind=r4_kind)) + esat = TABLE3(ind+1) + del*(DTABLE3(ind+1) + del*D2TABLE3(ind+1)) + type is (real(kind=r8_kind)) + esat = TABLE3(ind+1) + del*(DTABLE3(ind+1) + del*D2TABLE3(ind+1)) + end select endif end subroutine lookup_es3_k_0d From 07ef6564038df32d29f04a46c74a6f9bcc028888 Mon Sep 17 00:00:00 2001 From: MinsukJi-NOAA Date: Sun, 10 Oct 2021 11:19:52 -0400 Subject: [PATCH 5/8] Fix select type statement in sat_vapor_pres_k.F90 --- sat_vapor_pres/sat_vapor_pres_k.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sat_vapor_pres/sat_vapor_pres_k.F90 b/sat_vapor_pres/sat_vapor_pres_k.F90 index ca93cc0a37..95d281e082 100644 --- a/sat_vapor_pres/sat_vapor_pres_k.F90 +++ b/sat_vapor_pres/sat_vapor_pres_k.F90 @@ -820,7 +820,7 @@ subroutine compute_qs_k_1d (temp, press, eps, zvir, qs, nbad, q, hc, & select type (hc) type is (real(kind=r4_kind)) hc_loc = hc - type is (real(kind=r4_kind)) + type is (real(kind=r8_kind)) hc_loc = hc end select else From d68a06c39bd9702f7db027c2791f26a11311be65 Mon Sep 17 00:00:00 2001 From: MinsukJi-NOAA Date: Sun, 10 Oct 2021 11:24:43 -0400 Subject: [PATCH 6/8] Move write temp statement inside select type --- sat_vapor_pres/sat_vapor_pres.F90 | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/sat_vapor_pres/sat_vapor_pres.F90 b/sat_vapor_pres/sat_vapor_pres.F90 index 901dfd1122..ce25eb7048 100644 --- a/sat_vapor_pres/sat_vapor_pres.F90 +++ b/sat_vapor_pres/sat_vapor_pres.F90 @@ -2720,14 +2720,16 @@ subroutine show_all_bad_0d ( temp ) select type (temp) type is (real(kind=r4_kind)) ind = int(dtinv*(temp-tmin+teps)) + if (ind < 0 .or. ind > nlim) then + write(unit,'(a,e10.3,a,i6)') 'Bad temperature=',temp,' pe=',mpp_pe() + endif type is (real(kind=r8_kind)) ind = int(dtinv*(temp-tmin+teps)) + if (ind < 0 .or. ind > nlim) then + write(unit,'(a,e10.3,a,i6)') 'Bad temperature=',temp,' pe=',mpp_pe() + endif end select - if (ind < 0 .or. ind > nlim) then - write(unit,'(a,e10.3,a,i6)') 'Bad temperature=',temp,' pe=',mpp_pe() - endif - end subroutine show_all_bad_0d !-------------------------------------------------------------- From 906b3a65ae7d0cf2d26fb4fc7efc98e90ad73bb1 Mon Sep 17 00:00:00 2001 From: MinsukJi-NOAA Date: Thu, 21 Oct 2021 21:16:03 +0000 Subject: [PATCH 7/8] Fix range size check if statements --- diag_manager/diag_manager.F90 | 24 +++++++++++++++--------- 1 file changed, 15 insertions(+), 9 deletions(-) diff --git a/diag_manager/diag_manager.F90 b/diag_manager/diag_manager.F90 index 131e401f5c..9e98c90768 100644 --- a/diag_manager/diag_manager.F90 +++ b/diag_manager/diag_manager.F90 @@ -383,9 +383,11 @@ INTEGER FUNCTION register_diag_field_scalar(module_name, field_name, init_time, IF ( PRESENT(err_msg) ) err_msg = '' ! Fatal error if range is present and its extent is not 2. - IF ( PRESENT(range) .AND. (SIZE(range) .NE. 2) ) THEN - ! extent of range should be 2 - CALL error_mesg ('diag_manager_mod::register_diag_field', 'extent of range should be 2', FATAL) + IF ( PRESENT(range) ) THEN + IF ( SIZE(range) .NE. 2 ) THEN + ! extent of range should be 2 + CALL error_mesg ('diag_manager_mod::register_diag_field', 'extent of range should be 2', FATAL) + END IF END IF IF ( PRESENT(init_time) ) THEN @@ -449,9 +451,11 @@ INTEGER FUNCTION register_diag_field_array(module_name, field_name, axes, init_t IF ( PRESENT(err_msg) ) err_msg = '' ! Fatal error if range is present and its extent is not 2. - IF ( PRESENT(range) .AND. (SIZE(range) .NE. 2) ) THEN - ! extent of range should be 2 - CALL error_mesg ('diag_manager_mod::register_diag_field', 'extent of range should be 2', FATAL) + IF ( PRESENT(range) ) THEN + IF ( SIZE(range) .NE. 2 ) THEN + ! extent of range should be 2 + CALL error_mesg ('diag_manager_mod::register_diag_field', 'extent of range should be 2', FATAL) + END IF END IF ! Call register static, then set static back to false @@ -671,9 +675,11 @@ INTEGER FUNCTION register_static_field(module_name, field_name, axes, long_name, END IF ! Fatal error if range is present and its extent is not 2. - IF ( PRESENT(range) .AND. (SIZE(range) .NE. 2) ) THEN - ! extent of range should be 2 - CALL error_mesg ('diag_manager_mod::register_static_field', 'extent of range should be 2', FATAL) + IF ( PRESENT(range) ) THEN + IF ( SIZE(range) .NE. 2 ) THEN + ! extent of range should be 2 + CALL error_mesg ('diag_manager_mod::register_static_field', 'extent of range should be 2', FATAL) + END IF END IF ! Namelist do_diag_field_log is by default false. Thus to log the From 6029a7dcba7836833240a144263528cd6410ea02 Mon Sep 17 00:00:00 2001 From: MinsukJi-NOAA Date: Sun, 31 Oct 2021 15:43:24 +0000 Subject: [PATCH 8/8] Modify sat_vapor_pres_k.F90 to fix FV3 32 bit atmos_4xdaily* file NaN values --- sat_vapor_pres/sat_vapor_pres_k.F90 | 412 +++++++++++++++++++++------- 1 file changed, 311 insertions(+), 101 deletions(-) diff --git a/sat_vapor_pres/sat_vapor_pres_k.F90 b/sat_vapor_pres/sat_vapor_pres_k.F90 index 95d281e082..8aa8f6a710 100644 --- a/sat_vapor_pres/sat_vapor_pres_k.F90 +++ b/sat_vapor_pres/sat_vapor_pres_k.F90 @@ -487,11 +487,22 @@ subroutine compute_qs_k_3d (temp, press, eps, zvir, qs, nbad, q, hc, & logical,intent(in), optional :: es_over_liq logical,intent(in), optional :: es_over_liq_and_ice - real, dimension(size(temp,1), size(temp,2), size(temp,3)) :: & - esloc, desat, denom + real(kind=r4_kind), allocatable, dimension(:,:,:) :: esloc_r4, desat_r4, denom_r4 + real(kind=r8_kind), allocatable, dimension(:,:,:) :: esloc_r8, desat_r8, denom_r8 integer :: i, j, k real :: hc_loc + select type (temp) + type is (real(kind=r4_kind)) + allocate(esloc_r4(size(temp,1), size(temp,2), size(temp,3))) + allocate(desat_r4(size(temp,1), size(temp,2), size(temp,3))) + allocate(denom_r4(size(temp,1), size(temp,2), size(temp,3))) + type is (real(kind=r8_kind)) + allocate(esloc_r8(size(temp,1), size(temp,2), size(temp,3))) + allocate(desat_r8(size(temp,1), size(temp,2), size(temp,3))) + allocate(denom_r8(size(temp,1), size(temp,2), size(temp,3))) + end select + if (present(hc)) then select type (hc) type is (real(kind=r4_kind)) @@ -505,34 +516,73 @@ subroutine compute_qs_k_3d (temp, press, eps, zvir, qs, nbad, q, hc, & if (present(es_over_liq)) then if (present (dqsdT)) then - call lookup_es2_des2_k (temp, esloc, desat, nbad) - desat = desat*hc_loc + select type (temp) + type is (real(kind=r4_kind)) + call lookup_es2_des2_k (temp, esloc_r4, desat_r4, nbad) + desat_r4 = desat_r4*hc_loc + type is (real(kind=r8_kind)) + call lookup_es2_des2_k (temp, esloc_r8, desat_r8, nbad) + desat_r8 = desat_r8*hc_loc + end select else - call lookup_es2_k (temp, esloc, nbad) + select type (temp) + type is (real(kind=r4_kind)) + call lookup_es2_k (temp, esloc_r4, nbad) + type is (real(kind=r8_kind)) + call lookup_es2_k (temp, esloc_r8, nbad) + end select endif else if (present(es_over_liq_and_ice)) then if (present (dqsdT)) then - call lookup_es3_des3_k (temp, esloc, desat, nbad) - desat = desat*hc_loc + select type (temp) + type is (real(kind=r4_kind)) + call lookup_es3_des3_k (temp, esloc_r4, desat_r4, nbad) + desat_r4 = desat_r4*hc_loc + type is (real(kind=r8_kind)) + call lookup_es3_des3_k (temp, esloc_r8, desat_r8, nbad) + desat_r8 = desat_r8*hc_loc + end select else - call lookup_es3_k (temp, esloc, nbad) + select type (temp) + type is (real(kind=r4_kind)) + call lookup_es3_k (temp, esloc_r4, nbad) + type is (real(kind=r8_kind)) + call lookup_es3_k (temp, esloc_r8, nbad) + end select endif else if (present (dqsdT)) then - call lookup_es_des_k (temp, esloc, desat, nbad) - desat = desat*hc_loc + select type (temp) + type is (real(kind=r4_kind)) + call lookup_es_des_k (temp, esloc_r4, desat_r4, nbad) + desat_r4 = desat_r4*hc_loc + type is (real(kind=r8_kind)) + call lookup_es_des_k (temp, esloc_r8, desat_r8, nbad) + desat_r8 = desat_r8*hc_loc + end select else - call lookup_es_k (temp, esloc, nbad) + select type (temp) + type is (real(kind=r4_kind)) + call lookup_es_k (temp, esloc_r4, nbad) + type is (real(kind=r8_kind)) + call lookup_es_k (temp, esloc_r8, nbad) + end select endif endif - esloc = esloc*hc_loc + select type (temp) + type is (real(kind=r4_kind)) + esloc_r4 = esloc_r4*hc_loc + type is (real(kind=r8_kind)) + esloc_r8 = esloc_r8*hc_loc + end select + if (present (esat)) then select type (esat) type is (real(kind=r4_kind)) - esat = esloc + esat = esloc_r4 type is (real(kind=r8_kind)) - esat = esloc + esat = esloc_r8 end select endif @@ -544,21 +594,21 @@ subroutine compute_qs_k_3d (temp, press, eps, zvir, qs, nbad, q, hc, & if (present (q) .and. use_exact_qs) then select type (q) type is (real(kind=r4_kind)) - qs = (1.0 + zvir*q)*eps*esloc/press + qs = (1.0 + zvir*q)*eps*esloc_r4/press if (present (dqsdT)) then select type (dqsdT) type is (real(kind=r4_kind)) - dqsdT = (1.0 + zvir*q)*eps*desat/press + dqsdT = (1.0 + zvir*q)*eps*desat_r4/press end select endif end select else ! (present(q)) - denom = press - (1.0 - eps)*esloc + denom_r4 = press - (1.0 - eps)*esloc_r4 do k=1,size(qs,3) do j=1,size(qs,2) do i=1,size(qs,1) - if (denom(i,j,k) > 0.0) then - qs(i,j,k) = eps*esloc(i,j,k)/denom(i,j,k) + if (denom_r4(i,j,k) > 0.0) then + qs(i,j,k) = eps*esloc_r4(i,j,k)/denom_r4(i,j,k) else qs(i,j,k) = eps endif @@ -568,7 +618,7 @@ subroutine compute_qs_k_3d (temp, press, eps, zvir, qs, nbad, q, hc, & if (present (dqsdT)) then select type (dqsdT) type is (real(kind=r4_kind)) - dqsdT = eps*press*desat/denom**2 + dqsdT = eps*press*desat_r4/denom_r4**2 end select endif endif ! (present(q)) @@ -579,21 +629,21 @@ subroutine compute_qs_k_3d (temp, press, eps, zvir, qs, nbad, q, hc, & if (present (q) .and. use_exact_qs) then select type (q) type is (real(kind=r8_kind)) - qs = (1.0 + zvir*q)*eps*esloc/press + qs = (1.0 + zvir*q)*eps*esloc_r8/press if (present (dqsdT)) then select type (dqsdT) type is (real(kind=r8_kind)) - dqsdT = (1.0 + zvir*q)*eps*desat/press + dqsdT = (1.0 + zvir*q)*eps*desat_r8/press end select endif end select else ! (present(q)) - denom = press - (1.0 - eps)*esloc + denom_r8 = press - (1.0 - eps)*esloc_r8 do k=1,size(qs,3) do j=1,size(qs,2) do i=1,size(qs,1) - if (denom(i,j,k) > 0.0) then - qs(i,j,k) = eps*esloc(i,j,k)/denom(i,j,k) + if (denom_r8(i,j,k) > 0.0) then + qs(i,j,k) = eps*esloc_r8(i,j,k)/denom_r8(i,j,k) else qs(i,j,k) = eps endif @@ -603,7 +653,7 @@ subroutine compute_qs_k_3d (temp, press, eps, zvir, qs, nbad, q, hc, & if (present (dqsdT)) then select type (dqsdT) type is (real(kind=r8_kind)) - dqsdT = eps*press*desat/denom**2 + dqsdT = eps*press*desat_r8/denom_r8**2 end select endif endif ! (present(q)) @@ -634,6 +684,12 @@ subroutine compute_qs_k_3d (temp, press, eps, zvir, qs, nbad, q, hc, & endif endif ! (nbad = 0) + select type (temp) + type is (real(kind=r4_kind)) + deallocate(esloc_r4, desat_r4, denom_r4) + type is (real(kind=r8_kind)) + deallocate(esloc_r8, desat_r8, denom_r8) + end select end subroutine compute_qs_k_3d @@ -652,10 +708,22 @@ subroutine compute_qs_k_2d (temp, press, eps, zvir, qs, nbad, q, hc, & logical,intent(in), optional :: es_over_liq logical,intent(in), optional :: es_over_liq_and_ice - real, dimension(size(temp,1), size(temp,2)) :: esloc, desat, denom + real(kind=r4_kind), allocatable, dimension(:,:) :: esloc_r4, desat_r4, denom_r4 + real(kind=r8_kind), allocatable, dimension(:,:) :: esloc_r8, desat_r8, denom_r8 integer :: i, j real :: hc_loc + select type (temp) + type is (real(kind=r4_kind)) + allocate(esloc_r4(size(temp,1), size(temp,2))) + allocate(desat_r4(size(temp,1), size(temp,2))) + allocate(denom_r4(size(temp,1), size(temp,2))) + type is (real(kind=r8_kind)) + allocate(esloc_r8(size(temp,1), size(temp,2))) + allocate(desat_r8(size(temp,1), size(temp,2))) + allocate(denom_r8(size(temp,1), size(temp,2))) + end select + if (present(hc)) then select type (hc) type is (real(kind=r4_kind)) @@ -669,34 +737,73 @@ subroutine compute_qs_k_2d (temp, press, eps, zvir, qs, nbad, q, hc, & if (present(es_over_liq)) then if (present (dqsdT)) then - call lookup_es2_des2_k (temp, esloc, desat, nbad) - desat = desat*hc_loc + select type (temp) + type is (real(kind=r4_kind)) + call lookup_es2_des2_k (temp, esloc_r4, desat_r4, nbad) + desat_r4 = desat_r4*hc_loc + type is (real(kind=r8_kind)) + call lookup_es2_des2_k (temp, esloc_r8, desat_r8, nbad) + desat_r8 = desat_r8*hc_loc + end select else - call lookup_es2_k (temp, esloc, nbad) + select type (temp) + type is (real(kind=r4_kind)) + call lookup_es2_k (temp, esloc_r4, nbad) + type is (real(kind=r8_kind)) + call lookup_es2_k (temp, esloc_r8, nbad) + end select endif else if (present(es_over_liq_and_ice)) then if (present (dqsdT)) then - call lookup_es3_des3_k (temp, esloc, desat, nbad) - desat = desat*hc_loc + select type (temp) + type is (real(kind=r4_kind)) + call lookup_es3_des3_k (temp, esloc_r4, desat_r4, nbad) + desat_r4 = desat_r4*hc_loc + type is (real(kind=r8_kind)) + call lookup_es3_des3_k (temp, esloc_r8, desat_r8, nbad) + desat_r8 = desat_r8*hc_loc + end select else - call lookup_es3_k (temp, esloc, nbad) + select type (temp) + type is (real(kind=r4_kind)) + call lookup_es3_k (temp, esloc_r4, nbad) + type is (real(kind=r8_kind)) + call lookup_es3_k (temp, esloc_r8, nbad) + end select endif else if (present (dqsdT)) then - call lookup_es_des_k (temp, esloc, desat, nbad) - desat = desat*hc_loc + select type (temp) + type is (real(kind=r4_kind)) + call lookup_es_des_k (temp, esloc_r4, desat_r4, nbad) + desat_r4 = desat_r4*hc_loc + type is (real(kind=r8_kind)) + call lookup_es_des_k (temp, esloc_r8, desat_r8, nbad) + desat_r8 = desat_r8*hc_loc + end select else - call lookup_es_k (temp, esloc, nbad) + select type (temp) + type is (real(kind=r4_kind)) + call lookup_es_k (temp, esloc_r4, nbad) + type is (real(kind=r8_kind)) + call lookup_es_k (temp, esloc_r8, nbad) + end select endif endif - esloc = esloc*hc_loc + select type (temp) + type is (real(kind=r4_kind)) + esloc_r4 = esloc_r4*hc_loc + type is (real(kind=r8_kind)) + esloc_r8 = esloc_r8*hc_loc + end select + if (present (esat)) then select type (esat) type is (real(kind=r4_kind)) - esat = esloc + esat = esloc_r4 type is (real(kind=r8_kind)) - esat = esloc + esat = esloc_r8 end select endif @@ -708,20 +815,20 @@ subroutine compute_qs_k_2d (temp, press, eps, zvir, qs, nbad, q, hc, & if (present (q) .and. use_exact_qs) then select type (q) type is (real(kind=r4_kind)) - qs = (1.0 + zvir*q)*eps*esloc/press + qs = (1.0 + zvir*q)*eps*esloc_r4/press if (present (dqsdT)) then select type (dqsdT) type is (real(kind=r4_kind)) - dqsdT = (1.0 + zvir*q)*eps*desat/press + dqsdT = (1.0 + zvir*q)*eps*desat_r4/press end select endif end select else ! (present(q)) - denom = press - (1.0 - eps)*esloc + denom_r4 = press - (1.0 - eps)*esloc_r4 do j=1,size(qs,2) do i=1,size(qs,1) - if (denom(i,j) > 0.0) then - qs(i,j) = eps*esloc(i,j)/denom(i,j) + if (denom_r4(i,j) > 0.0) then + qs(i,j) = eps*esloc_r4(i,j)/denom_r4(i,j) else qs(i,j) = eps endif @@ -730,7 +837,7 @@ subroutine compute_qs_k_2d (temp, press, eps, zvir, qs, nbad, q, hc, & if (present (dqsdT)) then select type (dqsdT) type is (real(kind=r4_kind)) - dqsdT = eps*press*desat/denom**2 + dqsdT = eps*press*desat_r4/denom_r4**2 end select endif endif ! (present(q)) @@ -741,20 +848,20 @@ subroutine compute_qs_k_2d (temp, press, eps, zvir, qs, nbad, q, hc, & if (present (q) .and. use_exact_qs) then select type (q) type is (real(kind=r8_kind)) - qs = (1.0 + zvir*q)*eps*esloc/press + qs = (1.0 + zvir*q)*eps*esloc_r8/press if (present (dqsdT)) then select type (dqsdT) type is (real(kind=r8_kind)) - dqsdT = (1.0 + zvir*q)*eps*desat/press + dqsdT = (1.0 + zvir*q)*eps*desat_r8/press end select endif end select else ! (present(q)) - denom = press - (1.0 - eps)*esloc + denom_r8 = press - (1.0 - eps)*esloc_r8 do j=1,size(qs,2) do i=1,size(qs,1) - if (denom(i,j) > 0.0) then - qs(i,j) = eps*esloc(i,j)/denom(i,j) + if (denom_r8(i,j) > 0.0) then + qs(i,j) = eps*esloc_r8(i,j)/denom_r8(i,j) else qs(i,j) = eps endif @@ -763,7 +870,7 @@ subroutine compute_qs_k_2d (temp, press, eps, zvir, qs, nbad, q, hc, & if (present (dqsdT)) then select type (dqsdT) type is (real(kind=r8_kind)) - dqsdT = eps*press*desat/denom**2 + dqsdT = eps*press*desat_r8/denom_r8**2 end select endif endif ! (present(q)) @@ -794,6 +901,12 @@ subroutine compute_qs_k_2d (temp, press, eps, zvir, qs, nbad, q, hc, & endif endif ! (nbad = 0) + select type (temp) + type is (real(kind=r4_kind)) + deallocate(esloc_r4, desat_r4, denom_r4) + type is (real(kind=r8_kind)) + deallocate(esloc_r8, desat_r8, denom_r8) + end select end subroutine compute_qs_k_2d @@ -812,10 +925,22 @@ subroutine compute_qs_k_1d (temp, press, eps, zvir, qs, nbad, q, hc, & logical,intent(in), optional :: es_over_liq logical,intent(in), optional :: es_over_liq_and_ice - real, dimension(size(temp,1)) :: esloc, desat, denom + real(kind=r4_kind), allocatable, dimension(:) :: esloc_r4, desat_r4, denom_r4 + real(kind=r8_kind), allocatable, dimension(:) :: esloc_r8, desat_r8, denom_r8 integer :: i real :: hc_loc + select type (temp) + type is (real(kind=r4_kind)) + allocate(esloc_r4(size(temp,1))) + allocate(desat_r4(size(temp,1))) + allocate(denom_r4(size(temp,1))) + type is (real(kind=r8_kind)) + allocate(esloc_r8(size(temp,1))) + allocate(desat_r8(size(temp,1))) + allocate(denom_r8(size(temp,1))) + end select + if (present(hc)) then select type (hc) type is (real(kind=r4_kind)) @@ -829,34 +954,73 @@ subroutine compute_qs_k_1d (temp, press, eps, zvir, qs, nbad, q, hc, & if (present(es_over_liq)) then if (present (dqsdT)) then - call lookup_es2_des2_k (temp, esloc, desat, nbad) - desat = desat*hc_loc + select type (temp) + type is (real(kind=r4_kind)) + call lookup_es2_des2_k (temp, esloc_r4, desat_r4, nbad) + desat_r4 = desat_r4*hc_loc + type is (real(kind=r8_kind)) + call lookup_es2_des2_k (temp, esloc_r8, desat_r8, nbad) + desat_r8 = desat_r8*hc_loc + end select else - call lookup_es2_k (temp, esloc, nbad) + select type (temp) + type is (real(kind=r4_kind)) + call lookup_es2_k (temp, esloc_r4, nbad) + type is (real(kind=r8_kind)) + call lookup_es2_k (temp, esloc_r8, nbad) + end select endif else if (present(es_over_liq_and_ice)) then if (present (dqsdT)) then - call lookup_es3_des3_k (temp, esloc, desat, nbad) - desat = desat*hc_loc + select type (temp) + type is (real(kind=r4_kind)) + call lookup_es3_des3_k (temp, esloc_r4, desat_r4, nbad) + desat_r4 = desat_r4*hc_loc + type is (real(kind=r8_kind)) + call lookup_es3_des3_k (temp, esloc_r8, desat_r8, nbad) + desat_r8 = desat_r8*hc_loc + end select else - call lookup_es3_k (temp, esloc, nbad) + select type (temp) + type is (real(kind=r4_kind)) + call lookup_es3_k (temp, esloc_r4, nbad) + type is (real(kind=r8_kind)) + call lookup_es3_k (temp, esloc_r8, nbad) + end select endif else if (present (dqsdT)) then - call lookup_es_des_k (temp, esloc, desat, nbad) - desat = desat*hc_loc + select type (temp) + type is (real(kind=r4_kind)) + call lookup_es_des_k (temp, esloc_r4, desat_r4, nbad) + desat_r4 = desat_r4*hc_loc + type is (real(kind=r8_kind)) + call lookup_es_des_k (temp, esloc_r8, desat_r8, nbad) + desat_r8 = desat_r8*hc_loc + end select else - call lookup_es_k (temp, esloc, nbad) + select type (temp) + type is (real(kind=r4_kind)) + call lookup_es_k (temp, esloc_r4, nbad) + type is (real(kind=r8_kind)) + call lookup_es_k (temp, esloc_r8, nbad) + end select endif endif - esloc = esloc*hc_loc + select type (temp) + type is (real(kind=r4_kind)) + esloc_r4 = esloc_r4*hc_loc + type is (real(kind=r8_kind)) + esloc_r8 = esloc_r8*hc_loc + end select + if (present (esat)) then select type (esat) type is (real(kind=r4_kind)) - esat = esloc + esat = esloc_r4 type is (real(kind=r8_kind)) - esat = esloc + esat = esloc_r8 end select endif @@ -868,19 +1032,19 @@ subroutine compute_qs_k_1d (temp, press, eps, zvir, qs, nbad, q, hc, & if (present (q) .and. use_exact_qs) then select type (q) type is (real(kind=r4_kind)) - qs = (1.0 + zvir*q)*eps*esloc/press + qs = (1.0 + zvir*q)*eps*esloc_r4/press if (present (dqsdT)) then select type (dqsdT) type is (real(kind=r4_kind)) - dqsdT = (1.0 + zvir*q)*eps*desat/press + dqsdT = (1.0 + zvir*q)*eps*desat_r4/press end select endif end select else ! (present(q)) - denom = press - (1.0 - eps)*esloc + denom_r4 = press - (1.0 - eps)*esloc_r4 do i=1,size(qs,1) - if (denom(i) > 0.0) then - qs(i) = eps*esloc(i)/denom(i) + if (denom_r4(i) > 0.0) then + qs(i) = eps*esloc_r4(i)/denom_r4(i) else qs(i) = eps endif @@ -888,7 +1052,7 @@ subroutine compute_qs_k_1d (temp, press, eps, zvir, qs, nbad, q, hc, & if (present (dqsdT)) then select type (dqsdT) type is (real(kind=r4_kind)) - dqsdT = eps*press*desat/denom**2 + dqsdT = eps*press*desat_r4/denom_r4**2 end select endif endif ! (present(q)) @@ -899,19 +1063,19 @@ subroutine compute_qs_k_1d (temp, press, eps, zvir, qs, nbad, q, hc, & if (present (q) .and. use_exact_qs) then select type (q) type is (real(kind=r8_kind)) - qs = (1.0 + zvir*q)*eps*esloc/press + qs = (1.0 + zvir*q)*eps*esloc_r8/press if (present (dqsdT)) then select type (dqsdT) type is (real(kind=r8_kind)) - dqsdT = (1.0 + zvir*q)*eps*desat/press + dqsdT = (1.0 + zvir*q)*eps*desat_r8/press end select endif end select else ! (present(q)) - denom = press - (1.0 - eps)*esloc + denom_r8 = press - (1.0 - eps)*esloc_r8 do i=1,size(qs,1) - if (denom(i) > 0.0) then - qs(i) = eps*esloc(i)/denom(i) + if (denom_r8(i) > 0.0) then + qs(i) = eps*esloc_r8(i)/denom_r8(i) else qs(i) = eps endif @@ -919,7 +1083,7 @@ subroutine compute_qs_k_1d (temp, press, eps, zvir, qs, nbad, q, hc, & if (present (dqsdT)) then select type (dqsdT) type is (real(kind=r8_kind)) - dqsdT = eps*press*desat/denom**2 + dqsdT = eps*press*desat_r8/denom_r8**2 end select endif endif ! (present(q)) @@ -950,6 +1114,12 @@ subroutine compute_qs_k_1d (temp, press, eps, zvir, qs, nbad, q, hc, & endif endif ! (nbad = 0) + select type (temp) + type is (real(kind=r4_kind)) + deallocate(esloc_r4, desat_r4, denom_r4) + type is (real(kind=r8_kind)) + deallocate(esloc_r8, desat_r8, denom_r8) + end select end subroutine compute_qs_k_1d @@ -968,7 +1138,8 @@ subroutine compute_qs_k_0d (temp, press, eps, zvir, qs, nbad, q, hc, & logical,intent(in), optional :: es_over_liq logical,intent(in), optional :: es_over_liq_and_ice - real :: esloc, desat, denom + real(kind=r4_kind) :: esloc_r4, desat_r4, denom_r4 + real(kind=r8_kind) :: esloc_r8, desat_r8, denom_r8 real :: hc_loc if (present(hc)) then @@ -984,34 +1155,73 @@ subroutine compute_qs_k_0d (temp, press, eps, zvir, qs, nbad, q, hc, & if (present(es_over_liq)) then if (present (dqsdT)) then - call lookup_es2_des2_k (temp, esloc, desat, nbad) - desat = desat*hc_loc + select type (temp) + type is (real(kind=r4_kind)) + call lookup_es2_des2_k (temp, esloc_r4, desat_r4, nbad) + desat_r4 = desat_r4*hc_loc + type is (real(kind=r8_kind)) + call lookup_es2_des2_k (temp, esloc_r8, desat_r8, nbad) + desat_r8 = desat_r8*hc_loc + end select else - call lookup_es2_k (temp, esloc, nbad) + select type (temp) + type is (real(kind=r4_kind)) + call lookup_es2_k (temp, esloc_r4, nbad) + type is (real(kind=r8_kind)) + call lookup_es2_k (temp, esloc_r8, nbad) + end select endif else if (present(es_over_liq_and_ice)) then if (present (dqsdT)) then - call lookup_es3_des3_k (temp, esloc, desat, nbad) - desat = desat*hc_loc + select type (temp) + type is (real(kind=r4_kind)) + call lookup_es3_des3_k (temp, esloc_r4, desat_r4, nbad) + desat_r4 = desat_r4*hc_loc + type is (real(kind=r8_kind)) + call lookup_es3_des3_k (temp, esloc_r8, desat_r8, nbad) + desat_r8 = desat_r8*hc_loc + end select else - call lookup_es3_k (temp, esloc, nbad) + select type (temp) + type is (real(kind=r4_kind)) + call lookup_es3_k (temp, esloc_r4, nbad) + type is (real(kind=r8_kind)) + call lookup_es3_k (temp, esloc_r8, nbad) + end select endif else if (present (dqsdT)) then - call lookup_es_des_k (temp, esloc, desat, nbad) - desat = desat*hc_loc + select type (temp) + type is (real(kind=r4_kind)) + call lookup_es_des_k (temp, esloc_r4, desat_r4, nbad) + desat_r4 = desat_r4*hc_loc + type is (real(kind=r8_kind)) + call lookup_es_des_k (temp, esloc_r8, desat_r8, nbad) + desat_r8 = desat_r8*hc_loc + end select else - call lookup_es_k (temp, esloc, nbad) + select type (temp) + type is (real(kind=r4_kind)) + call lookup_es_k (temp, esloc_r4, nbad) + type is (real(kind=r8_kind)) + call lookup_es_k (temp, esloc_r8, nbad) + end select endif endif - esloc = esloc*hc_loc + select type (temp) + type is (real(kind=r4_kind)) + esloc_r4 = esloc_r4*hc_loc + type is (real(kind=r8_kind)) + esloc_r8 = esloc_r8*hc_loc + end select + if (present (esat)) then select type (esat) type is (real(kind=r4_kind)) - esat = esloc + esat = esloc_r4 type is (real(kind=r8_kind)) - esat = esloc + esat = esloc_r8 end select endif @@ -1023,25 +1233,25 @@ subroutine compute_qs_k_0d (temp, press, eps, zvir, qs, nbad, q, hc, & if (present (q) .and. use_exact_qs) then select type (q) type is (real(kind=r4_kind)) - qs = (1.0 + zvir*q)*eps*esloc/press + qs = (1.0 + zvir*q)*eps*esloc_r4/press if (present (dqsdT)) then select type (dqsdT) type is (real(kind=r4_kind)) - dqsdT = (1.0 + zvir*q)*eps*desat/press + dqsdT = (1.0 + zvir*q)*eps*desat_r4/press end select endif end select else ! (present(q)) - denom = press - (1.0 - eps)*esloc - if (denom > 0.0) then - qs = eps*esloc/denom + denom_r4 = press - (1.0 - eps)*esloc_r4 + if (denom_r4 > 0.0) then + qs = eps*esloc_r4/denom_r4 else qs = eps endif if (present (dqsdT)) then select type (dqsdT) type is (real(kind=r4_kind)) - dqsdT = eps*press*desat/denom**2 + dqsdT = eps*press*desat_r4/denom_r4**2 end select endif endif ! (present(q)) @@ -1052,25 +1262,25 @@ subroutine compute_qs_k_0d (temp, press, eps, zvir, qs, nbad, q, hc, & if (present (q) .and. use_exact_qs) then select type (q) type is (real(kind=r8_kind)) - qs = (1.0 + zvir*q)*eps*esloc/press + qs = (1.0 + zvir*q)*eps*esloc_r8/press if (present (dqsdT)) then select type (dqsdT) type is (real(kind=r8_kind)) - dqsdT = (1.0 + zvir*q)*eps*desat/press + dqsdT = (1.0 + zvir*q)*eps*desat_r8/press end select endif end select else ! (present(q)) - denom = press - (1.0 - eps)*esloc - if (denom > 0.0) then - qs = eps*esloc/denom + denom_r8 = press - (1.0 - eps)*esloc_r8 + if (denom_r8 > 0.0) then + qs = eps*esloc_r8/denom_r8 else qs = eps endif if (present (dqsdT)) then select type (dqsdT) type is (real(kind=r8_kind)) - dqsdT = eps*press*desat/denom**2 + dqsdT = eps*press*desat_r8/denom_r8**2 end select endif endif ! (present(q))