From 801352fff30aaa1c16dd6449a4a296d3093de4c6 Mon Sep 17 00:00:00 2001 From: ntlewis Date: Mon, 16 Jan 2023 09:47:42 +0000 Subject: [PATCH 01/19] thermodynamic sea ice --- .../two_stream_gray_rad.F90 | 9 +- .../driver/solo/idealized_moist_phys.F90 | 14 +- .../driver/solo/mixed_layer.F90 | 373 ++++++++++++++---- src/shared/constants/constants.F90 | 1 + 4 files changed, 310 insertions(+), 87 deletions(-) diff --git a/src/atmos_param/two_stream_gray_rad/two_stream_gray_rad.F90 b/src/atmos_param/two_stream_gray_rad/two_stream_gray_rad.F90 index 9849c30c0..da2cf4c60 100644 --- a/src/atmos_param/two_stream_gray_rad/two_stream_gray_rad.F90 +++ b/src/atmos_param/two_stream_gray_rad/two_stream_gray_rad.F90 @@ -138,6 +138,7 @@ module two_stream_gray_rad_mod character(len=256) :: co2_file='co2' ! file name of co2 file to read character(len=256) :: co2_variable_name='co2' ! file name of co2 file to read +logical :: do_toa_albedo = .false. namelist/two_stream_gray_rad_nml/ solar_constant, del_sol, & ir_tau_eq, ir_tau_pole, odp, atm_abs, sw_diff, & @@ -148,7 +149,8 @@ module two_stream_gray_rad_mod window, carbon_conc, rad_scheme, & do_read_co2, co2_file, co2_variable_name, solday, equinox_day, bog_a, bog_b, bog_mu, & use_time_average_coszen, dt_rad_avg,& - diabatic_acce !Schneider Liu values + diabatic_acce, & !Schneider Liu values + do_toa_albedo !================================================================================== !-------------------- diagnostics fields ------------------------------- @@ -453,6 +455,9 @@ subroutine two_stream_gray_rad_down (is, js, Time_diag, lat, lon, p_half, t, ! Default: Averaged Earth insolation at all longitudes p2 = (1. - 3.*sin(lat)**2)/4. insolation = 0.25 * solar_constant * (1.0 + del_sol * p2 + del_sw * sin(lat)) + if (do_toa_albedo) then + insolation = insolation * (0.77 + 0.125 * 2. * p2) + endif end if select case(sw_scheme) @@ -487,10 +492,12 @@ subroutine two_stream_gray_rad_down (is, js, Time_diag, lat, lon, p_half, t, end do ! compute downward shortwave flux + do k = 1, n+1 sw_down(:,:,k) = insolation(:,:) * exp(-sw_tau(:,:,k)) end do + case(B_SCHNEIDER_LIU) ! Schneider & Liu 2009 Giant planet scheme ! SW optical thickness diff --git a/src/atmos_spectral/driver/solo/idealized_moist_phys.F90 b/src/atmos_spectral/driver/solo/idealized_moist_phys.F90 index 1b86aff46..490c7a1e0 100644 --- a/src/atmos_spectral/driver/solo/idealized_moist_phys.F90 +++ b/src/atmos_spectral/driver/solo/idealized_moist_phys.F90 @@ -167,6 +167,7 @@ module idealized_moist_phys_mod real, allocatable, dimension(:,:) :: & z_surf, & ! surface height t_surf, & ! surface temperature + t_ml, h_thermo_ice, & ! NTL 01/23 thermodynamic sea ice q_surf, & ! surface moisture u_surf, & ! surface U wind v_surf, & ! surface V wind @@ -448,6 +449,10 @@ subroutine idealized_moist_phys_init(Time, Time_step_in, nhum, rad_lon_2d, rad_l allocate(depth_change_conv(is:ie, js:je)) allocate(z_surf (is:ie, js:je)) allocate(t_surf (is:ie, js:je)) +! NTL 01/23 thermodynamic sea ice +allocate(h_thermo_ice(is:ie, js:je)) +allocate(t_ml (is:ie, js:je)) +! allocate(q_surf (is:ie, js:je)); q_surf = 0.0 allocate(u_surf (is:ie, js:je)); u_surf = 0.0 allocate(v_surf (is:ie, js:je)); v_surf = 0.0 @@ -611,7 +616,9 @@ subroutine idealized_moist_phys_init(Time, Time_step_in, nhum, rad_lon_2d, rad_l ! to quickly enter the atmosphere avoiding problems with the convection scheme t_surf = t_surf_init + 1.0 - call mixed_layer_init(is, ie, js, je, num_levels, t_surf, bucket_depth, get_axis_id(), Time, albedo, rad_lonb_2d(:,:), rad_latb_2d(:,:), land, bucket) ! t_surf is intent(inout) !s albedo distribution set here. + call mixed_layer_init(is, ie, js, je, num_levels, t_surf, & + h_thermo_ice, t_ml, & ! NTL 01/23 thermodynamic sea ice + bucket_depth, get_axis_id(), Time, albedo, rad_lonb_2d(:,:), rad_latb_2d(:,:), land, bucket) ! t_surf is intent(inout) !s albedo distribution set here. elseif(gp_surface) then albedo=0.0 @@ -1248,6 +1255,7 @@ subroutine idealized_moist_phys(Time, p_half, p_full, z_half, z_full, ug, vg, tg js, & je, & t_surf(:,:), & ! t_surf is intent(inout) + h_thermo_ice(:,:), t_ml(:,:), & ! NTL 01/23 thermodynamic sea ice flux_t(:,:), & flux_q(:,:), & flux_r(:,:), & @@ -1347,7 +1355,9 @@ subroutine idealized_moist_phys_end call vert_turb_driver_end endif call lscale_cond_end -if(mixed_layer_bc) call mixed_layer_end(t_surf, bucket_depth, bucket) +if(mixed_layer_bc) call mixed_layer_end(t_surf, & + h_thermo_ice, t_ml, albedo, & ! NTL 01/23 thermodynamic sea ice + bucket_depth, bucket) if(do_damping) call damping_driver_end #ifdef SOC_NO_COMPILE diff --git a/src/atmos_spectral/driver/solo/mixed_layer.F90 b/src/atmos_spectral/driver/solo/mixed_layer.F90 index 4eee4ff1d..ca424c8d2 100644 --- a/src/atmos_spectral/driver/solo/mixed_layer.F90 +++ b/src/atmos_spectral/driver/solo/mixed_layer.F90 @@ -37,7 +37,8 @@ module mixed_layer_mod use tracer_manager_mod, only: get_tracer_names, get_number_tracers -use constants_mod, only: HLV, PI, RHO_CP, CP_AIR, KELVIN +use constants_mod, only: HLV, PI, RHO_CP, CP_AIR, KELVIN, & + HLF, DENS_ICE, TFREEZE ! NTL 01/23 for sea ice code use diag_manager_mod, only: register_diag_field, register_static_field, send_data @@ -128,6 +129,20 @@ module mixed_layer_mod character(len=256) :: flux_lhe_anom_file_name = 'INPUT/flux_lhe_anom.nc' character(len=256) :: flux_lhe_anom_field_name = 'flux_lhe_anom' +! THERMO ICE OPTIONS - NTL 01/23 +logical :: do_explicit = .false. ! explicit or implicit timestepping for ice +logical :: do_thermo_ice = .false. ! do thermodynamic sea ice? +real :: thermo_ice_albedo = 0.4 ! from Feldl & Merlis - original XZ is 0.8, M. England 0.6 +real :: L_thermo_ice = DENS_ICE * HLF ! latent heat of fusion (J/m^3) +real, parameter :: k_thermo_ice = 2 ! conductivity of ice (W/m/K) +real, parameter :: thermo_ice_basal_flux_const = 120 ! linear coefficient for heat flux from ocean to ice base (W/m^2/K) +real :: t_thermo_ice_base = TFREEZE ! temperature at base of ice, taken to be freshwater freezing point +real :: t_surf_freeze = TFREEZE + + + + + namelist/mixed_layer_nml/ evaporation, depth, qflux_amp, qflux_width, tconst,& delta_T, prescribe_initial_dist,albedo_value, & land_depth,trop_depth, & !mj @@ -145,7 +160,8 @@ module mixed_layer_mod ice_albedo_value, specify_sst_over_ocean_only, & ice_concentration_threshold, & add_latent_heat_flux_anom,flux_lhe_anom_file_name,& - flux_lhe_anom_field_name, do_ape_sst, qflux_field_name + flux_lhe_anom_field_name, do_ape_sst, qflux_field_name,& + do_explicit, do_thermo_ice, thermo_ice_albedo, t_surf_freeze ! ICE, NTL 01/23 !================================================================================================================================= @@ -163,7 +179,11 @@ module mixed_layer_mod id_heat_cap, & ! heat capacity id_albedo, & ! mj albedo id_ice_conc, & ! st ice concentration - id_delta_t_surf + id_delta_t_surf, & + ! NTL 01/23 thermodynamic sea ice + id_h_thermo_ice, & ! sea ice thickness + id_t_ml, & ! mixed layer temperature + id_flux_thermo_ice ! conductive heat flux through ice real, allocatable, dimension(:,:) :: & ocean_qflux, & ! Q-flux @@ -193,7 +213,14 @@ module mixed_layer_mod zsurf, & ! mj know about topography land_sea_heat_capacity,& sst_new, & ! mj input SST - albedo_initial + albedo_initial, & + ! NTL 01/23 thermodynamic sea ice + dFdt_surf, & ! d(corrected_flux)/d(t_surf), for calculation of ice sfc temperature + delta_t_ml, & ! Increment in mixed layer temperature + delta_h_thermo_ice, & ! Increment in sea ice thickness + delta_t_thermo_ice, & ! Increment in ice (steady-state) surface temperature + flux_thermo_ice, & ! Conductive heat flux through ice + corrected_flux_noq logical, allocatable, dimension(:,:) :: land_mask @@ -209,10 +236,13 @@ module mixed_layer_mod contains !================================================================================================================================= -subroutine mixed_layer_init(is, ie, js, je, num_levels, t_surf, bucket_depth, axes, Time, albedo, rad_lonb_2d,rad_latb_2d, land, restart_file_bucket_depth) +subroutine mixed_layer_init(is, ie, js, je, num_levels, t_surf, & + h_thermo_ice, t_ml, & ! NTL 01/23 thermodynamic sea ice + bucket_depth, axes, Time, albedo, rad_lonb_2d,rad_latb_2d, land, restart_file_bucket_depth) type(time_type), intent(in) :: Time -real, intent(out), dimension(:,:) :: t_surf, albedo +real, intent(out), dimension(:,:) :: t_surf, albedo +real, intent(out), dimension(:,:) :: h_thermo_ice, t_ml ! NTL 01/23 thermodynamic sea ice real, intent(out), dimension(:,:,:) :: bucket_depth integer, intent(in), dimension(4) :: axes real, intent(in), dimension(:,:) :: rad_lonb_2d, rad_latb_2d @@ -276,12 +306,21 @@ subroutine mixed_layer_init(is, ie, js, je, num_levels, t_surf, bucket_depth, ax allocate(delta_t_surf (is:ie, js:je)) allocate(eff_heat_capacity (is:ie, js:je)) allocate(corrected_flux (is:ie, js:je)) +allocate(corrected_flux_noq (is:ie, js:je)) ! NTL 01/23 thermodynamic sea ice allocate(t_surf_dependence (is:ie, js:je)) allocate (albedo_initial (is:ie, js:je)) allocate(land_sea_heat_capacity (is:ie, js:je)) allocate(zsurf (is:ie, js:je)) allocate(sst_new (is:ie, js:je)) allocate(land_mask (is:ie, js:je)); land_mask=land + +! NTL 01/23 thermodynamic sea ice +allocate(dFdt_surf (is:ie, js:je)) +allocate(delta_t_ml (is:ie, js:je)) +allocate(delta_h_thermo_ice (is:ie, js:je)) +allocate(delta_t_thermo_ice (is:ie, js:je)) +allocate(flux_thermo_ice (is:ie, js:je)) + ! !see if restart file exists for the surface temperature ! @@ -290,6 +329,12 @@ subroutine mixed_layer_init(is, ie, js, je, num_levels, t_surf, bucket_depth, ax ! latitude will be needed for oceanic q flux !s Moved up slightly so that rad_lat_2d can be used in initial temperature distribution if necessary. +if (do_thermo_ice .and. (.not. do_explicit)) then + call error_mesg('mixed_layer','setting do_explicit .true. for consistency with sea ice', WARNING) + do_explicit = .true. +endif + + call get_deg_lat(deg_lat) do j=js,je rad_lat_2d(:,j) = deg_lat(j)*PI/180. @@ -317,6 +362,11 @@ subroutine mixed_layer_init(is, ie, js, je, num_levels, t_surf, bucket_depth, ax call nullify_domain() call read_data(trim('INPUT/mixed_layer.res'), 't_surf', t_surf, grid_domain) + if (do_thermo_ice) then + call read_data(trim('INPUT/mixed_layer.res'), 'h_thermo_ice', h_thermo_ice, grid_domain) + call read_data(trim('INPUT/mixed_layer.res'), 't_ml', t_ml, grid_domain) + call read_data(trim('INPUT/mixed_layer.res'), 'albedo', albedo, grid_domain) + endif if (restart_file_bucket_depth) then call read_data(trim('INPUT/mixed_layer.res'), 'bucket_depth', bucket_depth, grid_domain) @@ -336,14 +386,26 @@ subroutine mixed_layer_init(is, ie, js, je, num_levels, t_surf, bucket_depth, ax elseif (prescribe_initial_dist) then ! call error_mesg('mixed_layer','mixed_layer restart file not found - initializing from prescribed distribution', WARNING) + t_surf(:,:) = tconst - delta_T*((3.*sin(rad_lat_2d)**2.)-1.)/3. - + if (do_thermo_ice) then ! NTL 01/23 thermodynamic sea ice + h_thermo_ice(:,:) = 0.0 + t_ml(:,:) = t_surf(:,:) + albedo(:,:) = albedo_value + where(land) albedo(:,:) = land_albedo_prefactor * albedo_value + albedo_initial(:,:) = albedo (:,:) + !where(.not. land) + ! where (t_surf .lt. TFREEZE) albedo(:,:) = thermo_ice_albedo + !endwhere + endif else call error_mesg('mixed_layer','mixed_layer restart file not found - initializing from lowest model level temp', WARNING) endif + + id_t_surf = register_diag_field(mod_name, 't_surf', & axes(1:2), Time, 'surface temperature','K') id_flux_t = register_diag_field(mod_name, 'flux_t', & @@ -356,11 +418,21 @@ subroutine mixed_layer_init(is, ie, js, je, num_levels, t_surf, bucket_depth, ax axes(1:2), 'mixed layer heat capacity','joules/m^2/deg C') id_delta_t_surf = register_diag_field(mod_name, 'delta_t_surf', & axes(1:2), Time, 'change in sst','K') -if (update_albedo_from_ice) then +if (update_albedo_from_ice .or. do_thermo_ice) then id_albedo = register_diag_field(mod_name, 'albedo', & axes(1:2), Time, 'surface albedo', 'none') - id_ice_conc = register_diag_field(mod_name, 'ice_conc', & - axes(1:2), Time, 'ice_concentration', 'none') + if (do_thermo_ice) then + ! NTL 01/23 thermodynamic sea ice + id_h_thermo_ice = register_diag_field(mod_name, 'h_thermo_ice', & + axes(1:2), Time, 'sea ice thickness','m') + id_t_ml = register_diag_field(mod_name, 't_ml', & + axes(1:2), Time, 'mixed layer tempeature','K') + id_flux_thermo_ice = register_diag_field(mod_name, 'flux_thermo_ice', & + axes(1:2), Time, 'conductive heat flux through sea ice','watts/m2') + else + id_ice_conc = register_diag_field(mod_name, 'ice_conc', & + axes(1:2), Time, 'ice_concentration', 'none') + endif else id_albedo = register_static_field(mod_name, 'albedo', & axes(1:2), 'surface albedo', 'none') @@ -418,65 +490,70 @@ subroutine mixed_layer_init(is, ie, js, je, num_levels, t_surf, bucket_depth, ax !s Prescribe albedo distribution here so that it will be the same in both two_stream_gray later and rrtmg radiation. -albedo(:,:) = albedo_value +if (.not. do_thermo_ice) then ! NTL 01/23 thermodynamic sea ice -- this is set / loaded earlier -if(trim(land_option) .eq. 'input') then + albedo(:,:) = albedo_value -where(land) albedo = land_albedo_prefactor * albedo + if(trim(land_option) .eq. 'input') then -endif + where(land) albedo = land_albedo_prefactor * albedo -!mj MiMA albedo choices. -select case (albedo_choice) - case (2) ! higher_albedo northward (lat_glacier>0) or southward (lat_glacier <0 ) of lat_glacier - do j = 1, size(t_surf,2) - lat = deg_lat(js+j-1) - ! mj SH or NH only - if ( lat_glacier .ge. 0. ) then - if ( lat > lat_glacier ) then - albedo(:,j) = higher_albedo - else - albedo(:,j) = albedo_value - endif - else - if ( lat < lat_glacier ) then + endif + + !mj MiMA albedo choices. + select case (albedo_choice) + case (2) ! higher_albedo northward (lat_glacier>0) or southward (lat_glacier <0 ) of lat_glacier + do j = 1, size(t_surf,2) + lat = deg_lat(js+j-1) + ! mj SH or NH only + if ( lat_glacier .ge. 0. ) then + if ( lat > lat_glacier ) then + albedo(:,j) = higher_albedo + else + albedo(:,j) = albedo_value + endif + else + if ( lat < lat_glacier ) then + albedo(:,j) = higher_albedo + else + albedo(:,j) = albedo_value + endif + endif + enddo + case (3) ! higher_albedo poleward of lat_glacier + do j = 1, size(t_surf,2) + lat = deg_lat(js+j-1) + if ( abs(lat) > lat_glacier ) then albedo(:,j) = higher_albedo - else + else albedo(:,j) = albedo_value - endif - endif - enddo - case (3) ! higher_albedo poleward of lat_glacier - do j = 1, size(t_surf,2) - lat = deg_lat(js+j-1) - if ( abs(lat) > lat_glacier ) then - albedo(:,j) = higher_albedo - else - albedo(:,j) = albedo_value - endif - enddo - case (4) ! exponential increase with albedo_exp - do j = 1, size(t_surf,2) - lat = abs(deg_lat(js+j-1)) - albedo(:,j) = albedo_value + (higher_albedo-albedo_value)*(lat/90.)**albedo_exp - enddo - case (5) ! tanh increase around albedo_cntr with albedo_wdth - do j = 1, size(t_surf,2) - lat = abs(deg_lat(js+j-1)) - albedo(:,j) = albedo_value + (higher_albedo-albedo_value)*& - 0.5*(1+tanh((lat-albedo_cntr)/albedo_wdth)) - enddo -end select - -albedo_initial=albedo - -if (update_albedo_from_ice) then - call interpolator_init( ice_interp, trim(ice_file_name)//'.nc', rad_lonb_2d, rad_latb_2d, data_out_of_bounds=(/CONSTANT/) ) - call read_ice_conc(Time) - call albedo_calc(albedo,Time) -else - if ( id_albedo > 0 ) used = send_data ( id_albedo, albedo ) -endif + endif + enddo + case (4) ! exponential increase with albedo_exp + do j = 1, size(t_surf,2) + lat = abs(deg_lat(js+j-1)) + albedo(:,j) = albedo_value + (higher_albedo-albedo_value)*(lat/90.)**albedo_exp + enddo + case (5) ! tanh increase around albedo_cntr with albedo_wdth + do j = 1, size(t_surf,2) + lat = abs(deg_lat(js+j-1)) + albedo(:,j) = albedo_value + (higher_albedo-albedo_value)*& + 0.5*(1+tanh((lat-albedo_cntr)/albedo_wdth)) + enddo + end select + + albedo_initial=albedo + + if (update_albedo_from_ice) then + call interpolator_init( ice_interp, trim(ice_file_name)//'.nc', rad_lonb_2d, rad_latb_2d, data_out_of_bounds=(/CONSTANT/) ) + call read_ice_conc(Time) + call albedo_calc(albedo,Time) + else + if ( id_albedo > 0 ) used = send_data ( id_albedo, albedo ) + endif +else ! NTL 01/23 thermodynamic sea ice + if ( id_albedo > 0 ) used = send_data ( id_albedo, albedo, Time) +endif ! if(do_thermo_ice) ! Note: do_calc_eff_heat_cap is true by default and control when the surface ! heat capacity is calculated (land and ocean). @@ -558,6 +635,10 @@ subroutine mixed_layer ( & Time_next, & js, je, & t_surf, & + ! NTL 01/23 thermodynamic sea ice + h_thermo_ice, & + t_ml, & + ! end ice addition flux_t, & flux_q, & flux_r, & @@ -579,6 +660,7 @@ subroutine mixed_layer ( & real, intent(in), dimension(:,:) :: net_surf_sw_down, surf_lw_down real, intent(in), dimension(:,:) :: flux_t, flux_q, flux_r real, intent(inout), dimension(:,:) :: t_surf +real, intent(inout), dimension(:,:) :: h_thermo_ice, t_ml ! NTL 01/23 thermodynamic sea ice real, intent(in), dimension(:,:) :: dhdt_surf, dedt_surf, dedq_surf, & drdt_surf, dhdt_atm, dedq_atm real, intent(in) :: dt @@ -604,7 +686,9 @@ subroutine mixed_layer ( & land_ice_mask=land_mask endif -call albedo_calc(albedo_out,Time_next) +if (.not. do_thermo_ice) then ! NTL 01/23 thermodynamic sea ice + call albedo_calc(albedo_out,Time_next) +endif !s Add latent heat flux anomalies before any of the calculations take place @@ -664,6 +748,14 @@ subroutine mixed_layer ( & t_surf_dependence = t_surf_dependence + beta_q * HLV endif +if (do_thermo_ice) then + ! NTL 01/23 thermodynamic sea ice + corrected_flux_noq = corrected_flux + ocean_qflux + flux_thermo_ice = 0 ! Initialize conductive heat flux through sea ice + ! for calculation of ice sfc temperature + dFdt_surf = dhdt_surf + drdt_surf + HLV * (dedt_surf) ! d(corrected_flux)/d(T_surf) +endif + !s Surface heat_capacity calculation based on that in MiMA by mj if(do_sc_sst) then !mj sst read from input file @@ -705,33 +797,134 @@ subroutine mixed_layer ( & endif -if (do_calc_eff_heat_cap) then - !s use the land_sea_heat_capacity calculated in mixed_layer_init +if (.not. do_thermo_ice) then ! NTL 01/23 thermodynamic sea ice switch - ! Now update the mixed layer surface temperature using an implicit step - ! - eff_heat_capacity = land_sea_heat_capacity + t_surf_dependence * dt !s need to investigate how this works + if (do_calc_eff_heat_cap) then + !s use the land_sea_heat_capacity calculated in mixed_layer_init - if (any(eff_heat_capacity .eq. 0.0)) then - write(*,*) 'mixed_layer: error', eff_heat_capacity - call error_mesg('mixed_layer', 'Avoiding division by zero',fatal) - end if + ! Now update the mixed layer surface temperature using an implicit step + ! + if (do_explicit) then ! NTL 01/23 thermodynamic sea ice -- add explicit option for consistency with thermo ice + eff_heat_capacity = land_sea_heat_capacity + else + eff_heat_capacity = land_sea_heat_capacity + t_surf_dependence * dt !s need to investigate how this works + endif - if(do_sc_sst.and.specify_sst_over_ocean_only) then - where (land_ice_mask) delta_t_surf = - corrected_flux * dt / eff_heat_capacity - where (land_ice_mask) t_surf = t_surf + delta_t_surf - else - delta_t_surf = - corrected_flux * dt / eff_heat_capacity - t_surf = t_surf + delta_t_surf - endif + if (any(eff_heat_capacity .eq. 0.0)) then + write(*,*) 'mixed_layer: error', eff_heat_capacity + call error_mesg('mixed_layer', 'Avoiding division by zero',fatal) + end if + + if(do_sc_sst.and.specify_sst_over_ocean_only) then + where (land_ice_mask) delta_t_surf = - corrected_flux * dt / eff_heat_capacity + where (land_ice_mask) t_surf = t_surf + delta_t_surf + else + delta_t_surf = - corrected_flux * dt / eff_heat_capacity + t_surf = t_surf + delta_t_surf + endif + + endif !s end of if(do_sc_sst). + +else ! if (do_thermo_ice) + where (land_ice_mask) ! name misleading, for thermo sea ice, land_ice_mask is just land mask + delta_t_ml = - corrected_flux * dt / land_sea_heat_capacity + delta_h_thermo_ice = 0 + ! t_surf and t_ml are equal (they differ only when mixed layer is ice-covered) + t_surf = t_ml + delta_t_ml + ! = update state = + t_ml = t_ml + delta_t_ml + h_thermo_ice = h_thermo_ice + delta_h_thermo_ice + albedo_out(:,:) = albedo_value * land_albedo_prefactor + elsewhere + + ! NTL 01/23 thermodynamic sea ice + ! COPIED FROM Sea ice by Ian Eisenman, XZ 02/2018 + ! === (3) surface temperature increment is calculated with + ! explicit forward time step using flux corrected for implicit atm; + ! next, (4) surface state is updated === + + ! ====================================================================== + ! + ! Ocean mixed layer and sea ice model equations + ! + ! ====================================================================== + ! + ! Mixed layer temperature is evolved is t_ml, and sea ice thickness is + ! h_ice. Atmosphere cares about surface temparature (t_surf). Where + ! h_ice=0, t_surf=t_ml, but where h_ice>0, t_surf=t_ice which is the + ! ice surface temperature (assumed to be in steady-state: "zero layer + ! model"). So h_ice, t_ml, and t_surf are all saved at each time step + ! (t_ice does not need to be saved). + ! + ! This model is equivalent to the Semtner (1976) "zero-layer model", + ! with several differences (no snow, ice latent heat is same at base + ! as surface, ice surface temperature calculated by including sensible + ! and latent heat derivatives in dF/dT_surf, inclusion of frazil + ! growth). + + ! calculate values of delta_h_ice, delta_t_ml, and delta_t_ice + + where ( h_thermo_ice .le. 0 ) ! = ice-free = [ should be equivalent to use .eq. instead of .le. ] + delta_t_ml = - ( corrected_flux - ocean_qflux) * dt / land_sea_heat_capacity + delta_h_thermo_ice = 0 + elsewhere ! = ice-covered = ( h>0 ) + delta_t_ml = - ( thermo_ice_basal_flux_const * ( t_ml - t_thermo_ice_base ) - ocean_qflux ) * dt / land_sea_heat_capacity + delta_h_thermo_ice = ( corrected_flux - thermo_ice_basal_flux_const * ( t_ml - t_thermo_ice_base ) ) * dt / L_thermo_ice + endwhere + where ( t_ml + delta_t_ml .lt. t_surf_freeze ) ! = frazil growth = + delta_h_thermo_ice = delta_h_thermo_ice - ( t_ml + delta_t_ml - t_surf_freeze ) * (land_sea_heat_capacity)/L_thermo_ice + delta_t_ml = t_surf_freeze - t_ml + endwhere + where ( ( h_thermo_ice .gt. 0 ) .and. ( h_thermo_ice + delta_h_thermo_ice .le. 0 ) ) ! = complete ablation = + delta_t_ml = delta_t_ml - ( h_thermo_ice + delta_h_thermo_ice ) * L_thermo_ice/land_sea_heat_capacity + delta_h_thermo_ice = - h_thermo_ice + endwhere + ! = update surface temperature = + ! compute surface melt + where ( h_thermo_ice + delta_h_thermo_ice .gt. 0 ) ! surface is ice-covered + ! calculate increment in steady-state ice surface temperature + flux_thermo_ice = k_thermo_ice / (h_thermo_ice + delta_h_thermo_ice) * ( t_thermo_ice_base - t_surf ) + delta_t_thermo_ice = ( - corrected_flux + flux_thermo_ice ) & + / ( k_thermo_ice / (h_thermo_ice + delta_h_thermo_ice) + dFdt_surf ) + where ( t_surf + delta_t_thermo_ice .gt. t_surf_freeze ) ! surface ablation + delta_t_thermo_ice = t_surf_freeze - t_surf + endwhere + ! surface is ice-covered, so update t_surf as ice surface temperature + t_surf = t_surf + delta_t_thermo_ice + elsewhere ! ice-free, so update t_surf as mixed layer temperature + t_surf = t_ml + delta_t_ml + endwhere + + + + + ! = update state = + t_ml = t_ml + delta_t_ml + h_thermo_ice = h_thermo_ice + delta_h_thermo_ice + + where ( h_thermo_ice .gt. 0 ) + albedo_out(:,:) = albedo_value + (thermo_ice_albedo - albedo_value) * tanh(h_thermo_ice/2.) + elsewhere + albedo_out(:,:) = albedo_value + endwhere + endwhere + + + if ( id_albedo > 0 ) used = send_data ( id_albedo, albedo_out, Time_next ) + +endif ! if (do_thermo_ice) -endif !s end of if(do_sc_sst). ! ! Finally calculate the increments for the lowest atmospheric layer ! -Tri_surf%delta_t = fn_t + en_t * delta_t_surf -if (evaporation) Tri_surf%delta_tr(:,:,nhum) = fn_q + en_q * delta_t_surf +if (do_explicit) then ! NTL 01/23 thermodynamic sea ice uses explicit step + Tri_surf%delta_t = fn_t + if (evaporation) Tri_surf%delta_tr(:,:,nhum) = fn_q +else + Tri_surf%delta_t = fn_t + en_t * delta_t_surf + if (evaporation) Tri_surf%delta_tr(:,:,nhum) = fn_q + en_q * delta_t_surf +endif ! ! Note: @@ -743,6 +936,10 @@ subroutine mixed_layer ( & if(id_flux_oceanq > 0) used = send_data(id_flux_oceanq, ocean_qflux, Time_next) if(id_delta_t_surf > 0) used = send_data(id_delta_t_surf, delta_t_surf, Time_next) +! NTL 01/23 thermodynamic sea ice +if(id_h_thermo_ice > 0) used = send_data(id_h_thermo_ice, h_thermo_ice, Time_next) +if(id_t_ml > 0) used = send_data(id_t_ml, t_ml, Time_next) +if(id_flux_thermo_ice > 0) used = send_data(id_flux_thermo_ice, flux_thermo_ice, Time_next) end subroutine mixed_layer @@ -780,9 +977,12 @@ subroutine read_ice_conc(Time) end subroutine read_ice_conc !================================================================================================================================= -subroutine mixed_layer_end(t_surf, bucket_depth, restart_file_bucket_depth) +subroutine mixed_layer_end(t_surf, & + h_thermo_ice, t_ml, albedo, & ! NTL 01/23 thermodynamic sea ice + bucket_depth, restart_file_bucket_depth) real, intent(inout), dimension(:,:) :: t_surf +real, intent(inout), dimension(:,:) :: h_thermo_ice, t_ml, albedo ! NTL 01/23 thermodynamic sea ice real, intent(inout), dimension(:,:,:) :: bucket_depth logical, intent(in) :: restart_file_bucket_depth integer:: unit @@ -792,6 +992,11 @@ subroutine mixed_layer_end(t_surf, bucket_depth, restart_file_bucket_depth) ! write a restart file for the surface temperature call nullify_domain() call write_data(trim('RESTART/mixed_layer.res'), 't_surf', t_surf, grid_domain) +if (do_thermo_ice) then ! NTL 01/23 thermodynamic sea ice + call write_data(trim('RESTART/mixed_layer.res'), 'h_thermo_ice', h_thermo_ice, grid_domain) + call write_data(trim('RESTART/mixed_layer.res'), 't_ml', t_ml, grid_domain) + call write_data(trim('RESTART/mixed_layer.res'), 'albedo', albedo, grid_domain) +endif if (restart_file_bucket_depth) then call write_data(trim('RESTART/mixed_layer.res'), 'bucket_depth', bucket_depth, grid_domain) endif diff --git a/src/shared/constants/constants.F90 b/src/shared/constants/constants.F90 index 42ee900a5..f7fa83c2d 100644 --- a/src/shared/constants/constants.F90 +++ b/src/shared/constants/constants.F90 @@ -120,6 +120,7 @@ module constants_mod real, public, parameter :: RVGAS = 461.50 real, public, parameter :: CP_VAPOR = 4.0*RVGAS real, public, parameter :: DENS_H2O = 1000. +real, public, parameter :: DENS_ICE = 917. ! NTL 01/23 for sea ice code real, public, parameter :: HLV = 2.500e6 real, public, parameter :: HLF = 3.34e5 real, public, parameter :: HLS = HLV + HLF From 3f76abb42f91a9af1bf80f52050ebb48df48d94d Mon Sep 17 00:00:00 2001 From: ntlewis Date: Mon, 16 Jan 2023 10:25:30 +0000 Subject: [PATCH 02/19] options for single precision (won't work with socrates) --- .../qe_moist_convection/qe_moist_convection.F90 | 5 +++-- src/extra/python/isca/codebase.py | 17 ++++++++++++++++- src/extra/python/isca/templates/compile.sh | 16 +++++++++------- .../python/isca/templates/mkmf.template.ia64 | 4 ++-- 4 files changed, 30 insertions(+), 12 deletions(-) diff --git a/src/atmos_param/qe_moist_convection/qe_moist_convection.F90 b/src/atmos_param/qe_moist_convection/qe_moist_convection.F90 index af299fcfa..aa55e0ec6 100644 --- a/src/atmos_param/qe_moist_convection/qe_moist_convection.F90 +++ b/src/atmos_param/qe_moist_convection/qe_moist_convection.F90 @@ -70,13 +70,14 @@ module qe_moist_convection_mod real :: val_inc = 0.01 real :: val_min = -1. ! calculated in get_lcl_temp_table_size real :: val_max = 1. ! calculated in get_lcl_temp_table_size + real :: precision = 1.e-7 real, parameter :: small = 1.e-10, & ! to avoid division by 0 in dry limit pref = 1.e5 real, allocatable, dimension(:) :: lcl_temp_table - namelist /qe_moist_convection_nml/ tau_bm, rhbm, Tmin, Tmax, val_inc + namelist /qe_moist_convection_nml/ tau_bm, rhbm, Tmin, Tmax, val_inc, precision !------------------------------------------------------------------ ! Description of namelist variables @@ -1084,7 +1085,7 @@ end subroutine get_lcl_temp real function lcl_temp(value, lcl_temp_guess) real, intent(in) :: value, lcl_temp_guess - real, parameter :: precision = 1.e-7 + ! real, parameter :: precision = 1.e-7 integer, parameter :: max_iter = 100 real :: T, dT integer :: iter diff --git a/src/extra/python/isca/codebase.py b/src/extra/python/isca/codebase.py index 09353e7ce..6c8b27370 100644 --- a/src/extra/python/isca/codebase.py +++ b/src/extra/python/isca/codebase.py @@ -120,6 +120,8 @@ def __init__(self, repo=None, commit=None, directory=None, storedir=P(GFDL_WORK, # read path names from the default file self.path_names = [] self.compile_flags = [] # users can append to this to add additional compiler options + self.precision_compile_flags = ['-DOVERLOAD_C8'] # default double precision + self.other_compile_flags = ['-r8', '-i4'] # default double precision @property def code_is_available(self): @@ -245,6 +247,7 @@ def compile(self, debug=False, optimisation=None): mkdir(self.builddir) compile_flags = [] + other_flags = [] # if debug: # compile_flags.append('-g') # compile_flags.append('-traceback') @@ -254,7 +257,11 @@ def compile(self, debug=False, optimisation=None): # compile_flags.append('-O%d' % optimisation) compile_flags.extend(self.compile_flags) + compile_flags.extend(self.precision_compile_flags) compile_flags_str = ' '.join(compile_flags) + + other_flags.extend(self.other_compile_flags) + other_flags_str = ' '.join(other_flags) # get path_names from the directory if not self.path_names: @@ -268,6 +275,7 @@ def compile(self, debug=False, optimisation=None): 'srcdir': self.srcdir, 'workdir': self.workdir, 'compile_flags': compile_flags_str, + 'extra_compile_flags': other_flags_str, 'env_source': env, 'path_names': path_names_str, 'executable_name': self.executable_name, @@ -280,6 +288,13 @@ def compile(self, debug=False, optimisation=None): self._log_line(line) self.log.info('Compilation complete.') + + def use_single_precision(self): + self.log.info('Using single precision') + self.precision_compile_flags = ['-DOVERLOAD_C4', '-DOVERLOAD_R4'] + self.other_compile_flags = ['-i4'] + self.executable_name = self.executable_name.strip('.x')+'_single.x' + self.builddir = P(self.workdir, 'build', self.executable_name.split('.')[0]) @@ -393,4 +408,4 @@ class DryCodeBase(GreyCodeBase): # """The Shallow Water Equations. # """ # name = 'shallow' -# executable_name = 'shallow.x' +# executable_name = 'shallow.x' \ No newline at end of file diff --git a/src/extra/python/isca/templates/compile.sh b/src/extra/python/isca/templates/compile.sh index d20377492..7811a643d 100755 --- a/src/extra/python/isca/templates/compile.sh +++ b/src/extra/python/isca/templates/compile.sh @@ -18,9 +18,9 @@ template_debug={{ template_dir }}/mkmf.template.debug execdir={{ execdir }} # where code is compiled and executable is created executable={{ executable_name }} -netcdf_flags=`nf-config --fflags --flibs` +#netcdf_flags=`nf-config --fflags --flibs` -ulimit -s unlimited # Set stack size to unlimited +#ulimit -s unlimited # Set stack size to unlimited export MALLOC_CHECK_=0 # 3. compile the mppncombine tool if it hasn't yet been done. @@ -53,14 +53,16 @@ if [ $debug == True ]; then echo "Compiling in debug mode" # execute mkmf to create makefile -cppDefs="-Duse_libMPI -Duse_netCDF -Duse_LARGEFILE -DINTERNAL_FILE_NML -DOVERLOAD_C8 {{compile_flags}}" -$mkmf -a $sourcedir -t $template_debug -p $executable -c "$cppDefs" $pathnames $sourcedir/shared/include $sourcedir/shared/mpp/include +cppDefs="-Duse_libMPI -Duse_netCDF -Duse_LARGEFILE -DINTERNAL_FILE_NML {{compile_flags}}" +othDefs="{{extra_compile_flags}}" +$mkmf -a $sourcedir -t $template_debug -p $executable -c "$cppDefs" -o "$othDefs" $pathnames $sourcedir/shared/include $sourcedir/shared/mpp/include else # execute mkmf to create makefile -cppDefs="-Duse_libMPI -Duse_netCDF -Duse_LARGEFILE -DINTERNAL_FILE_NML -DOVERLOAD_C8 {{compile_flags}}" -$mkmf -a $sourcedir -t $template -p $executable -c "$cppDefs" $pathnames $sourcedir/shared/include $sourcedir/shared/mpp/include +cppDefs="-Duse_libMPI -Duse_netCDF -Duse_LARGEFILE -DINTERNAL_FILE_NML ${CDEFS} {{compile_flags}}" +othDefs="{{extra_compile_flags}}" +$mkmf -a $sourcedir -t $template -p $executable -c "$cppDefs" -o "$othDefs" $pathnames $sourcedir/shared/include $sourcedir/shared/mpp/include fi @@ -80,4 +82,4 @@ make $executable if [ $? != 0 ]; then echo "ERROR: make failed for $executable" exit 1 -fi +fi \ No newline at end of file diff --git a/src/extra/python/isca/templates/mkmf.template.ia64 b/src/extra/python/isca/templates/mkmf.template.ia64 index dcb29e579..e3943e08e 100644 --- a/src/extra/python/isca/templates/mkmf.template.ia64 +++ b/src/extra/python/isca/templates/mkmf.template.ia64 @@ -18,11 +18,11 @@ NETCDF_LIBS = `nc-config --libs` # -diag-disable 6843: # This suppresses the warning: `warning #6843: A dummy argument with an explicit INTENT(OUT) declaration is not given an explicit value.` of which # there are a lot of instances in the GFDL codebase. -FFLAGS = $(CPPFLAGS) -fpp -stack_temps -safe_cray_ptr -ftz -assume byterecl -shared-intel -i4 -r8 -g -O2 -diag-disable 6843 -mcmodel large +FFLAGS = $(CPPFLAGS) -fpp -stack_temps -safe_cray_ptr -ftz -assume byterecl -shared-intel -g -O2 -diag-disable 6843 -mcmodel large #FFLAGS = $(CPPFLAGS) -fltconsistency -stack_temps -safe_cray_ptr -ftz -shared-intel -assume byterecl -g -O0 -i4 -r8 -check -warn -warn noerrors -debug variable_locations -inline_debug_info -traceback FC = $(F90) LD = $(F90) $(NETCDF_LIBS) #CC = mpicc LDFLAGS = -lnetcdff -lnetcdf -lmpi -shared-intel -CFLAGS = -D__IFC +CFLAGS = -D__IFC \ No newline at end of file From 420d99847855d95658fd960580bccfe999fcf199 Mon Sep 17 00:00:00 2001 From: ntlewis Date: Fri, 21 Jul 2023 10:06:38 +0100 Subject: [PATCH 03/19] changes to rrtm to have more control over coszen, and to allow for single precision --- .../rrtm_radiation/rrtm_radiation.f90 | 56 +++++++++++++++---- .../modules/{parkind.f90 => parkind.F90} | 12 +++- ...mg_lw_setcoef.f90 => rrtmg_lw_setcoef.F90} | 5 +- .../modules/{parkind.f90 => parkind.F90} | 13 ++++- .../rrtmg_sw/gcm_model/src/rrtmg_sw_rad.f90 | 8 ++- .../gcm_model/src/rrtmg_sw_rad.nomcica.f90 | 8 ++- .../gcm_model/src/rrtmg_sw_spcvmc.f90 | 5 +- .../gcm_model/src/rrtmg_sw_spcvrt.f90 | 5 +- 8 files changed, 89 insertions(+), 23 deletions(-) rename src/atmos_param/rrtm_radiation/rrtmg_lw/gcm_model/modules/{parkind.f90 => parkind.F90} (65%) rename src/atmos_param/rrtm_radiation/rrtmg_lw/gcm_model/src/{rrtmg_lw_setcoef.f90 => rrtmg_lw_setcoef.F90} (99%) rename src/atmos_param/rrtm_radiation/rrtmg_sw/gcm_model/modules/{parkind.f90 => parkind.F90} (65%) diff --git a/src/atmos_param/rrtm_radiation/rrtm_radiation.f90 b/src/atmos_param/rrtm_radiation/rrtm_radiation.f90 index 808d5239d..b547b9f05 100644 --- a/src/atmos_param/rrtm_radiation/rrtm_radiation.f90 +++ b/src/atmos_param/rrtm_radiation/rrtm_radiation.f90 @@ -79,7 +79,7 @@ module rrtm_vars ! cloud & aerosol optical depths, cloud and aerosol specific parameters. Set to zero real(kind=rb),allocatable,dimension(:,:,:) :: taucld,tauaer, sw_zro, zro_sw ! heating rates and fluxes, zenith angle when in-between radiation time steps - real(kind=rb),allocatable,dimension(:,:) :: sw_flux,lw_flux,zencos, olr, toa_sw! surface and TOA fluxes, cos(zenith angle) + real(kind=rb),allocatable,dimension(:,:) :: sw_flux,lw_flux,zencos, olr, toa_sw, toa_sw_down! surface and TOA fluxes, cos(zenith angle) ! dimension (lon x lat) real(kind=rb),allocatable,dimension(:,:,:) :: tdt_rad ! heating rate [K/s] ! dimension (lon x lat x pfull) @@ -162,7 +162,8 @@ module rrtm_vars ! for faster radiation calculation !!!!!! mp586 added for annual mean insolation!!!!! - + logical :: do_toa_albedo = .false. + real(kind=rb):: alb_scaler = 1. logical :: frierson_solar_rad =.false. real(kind=rb) :: del_sol = 0.95 ! frierson 2006 default = 1.4, but 0.95 gets the curve closer to the annual mean insolation real(kind=rb) :: del_sw = 0.0 !frierson 2006 default @@ -192,7 +193,7 @@ module rrtm_vars ! !-------------------- diagnostics fields ------------------------------- - integer :: id_tdt_rad, id_tdt_sw, id_tdt_lw, id_coszen, id_flux_sw, id_flux_lw, id_olr, id_toa_sw, id_albedo,id_ozone, id_co2, id_fracday, id_half_level_temp, id_full_level_temp + integer :: id_tdt_rad, id_tdt_sw, id_tdt_lw, id_coszen, id_flux_sw, id_flux_lw, id_olr, id_toa_sw, id_toa_sw_down, id_albedo,id_ozone, id_co2, id_fracday, id_half_level_temp, id_full_level_temp character(len=14), parameter :: mod_name_rad = 'rrtm_radiation' !s changed parameter name from mod_name to mod_name_rad as compiler objected, presumably because mod_name also defined in idealized_moist_physics.F90 after use rrtm_vars is included. real :: missing_value = -999. @@ -210,7 +211,7 @@ module rrtm_vars &lonstep, do_zm_tracers, do_zm_rad, & &do_precip_albedo, precip_albedo_mode, precip_albedo, precip_lat,& &do_read_co2, co2_file, co2_variable_name, use_dyofyr, solrad, & - &solday, equinox_day,solr_cnst + &solday, equinox_day,solr_cnst,do_toa_albedo,alb_scaler end module rrtm_vars !***************************************************************************************** @@ -305,6 +306,10 @@ subroutine rrtm_radiation_init(axes,Time,ncols,nlay,lonb,latb, Time_step) register_diag_field ( mod_name_rad, 'toa_sw', axes(1:2), Time, & 'Net TOA SW flux', & 'W/m2', missing_value=missing_value ) + id_toa_sw_down = & + register_diag_field ( mod_name_rad, 'toa_sw_down', axes(1:2), Time, & + 'Downward TOA SW flux', & + 'W/m2', missing_value=missing_value ) id_albedo = & register_diag_field ( mod_name_rad, 'rrtm_albedo', axes(1:2), Time, & 'Interactive albedo', & @@ -463,6 +468,8 @@ subroutine rrtm_radiation_init(axes,Time,ncols,nlay,lonb,latb, Time_step) allocate(olr(size(lonb,1)-1,size(latb,2)-1)) if(id_toa_sw > 0) & allocate(toa_sw(size(lonb,1)-1,size(latb,2)-1)) + if(id_toa_sw_down > 0) & + allocate(toa_sw_down(size(lonb,1)-1,size(latb,2)-1)) if(do_precip_albedo)allocate(rrtm_precip(size(lonb,1)-1,size(latb,2)-1)) if(store_intermediate_rad .or. id_tdt_rad > 0)& allocate(tdt_rad(size(lonb,1)-1,size(latb,2)-1,nlay)) @@ -591,10 +598,11 @@ subroutine run_rrtmg(is,js,Time,lat,lon,p_full,p_half,albedo,q,t,t_surf_rad,tdt, real(kind=rb),dimension(size(q,1)/lonstep,size(q,2),size(q,3) ) :: swijk,lwijk real(kind=rb),dimension(size(q,1)/lonstep,size(q,2)) :: swflxijk,lwflxijk real(kind=rb),dimension(ncols_rrt,nlay_rrt+1):: phalf,thalf - real(kind=rb),dimension(ncols_rrt) :: tsrf,cosz_rr,albedo_rr + real(kind=rb),dimension(ncols_rrt) :: tsrf,cosz_rr,albedo_rr, cloud_coalb_spatial_rr real(kind=rb) :: dlon,dlat,dj,di type(time_type) :: Time_loc real(kind=rb),dimension(size(q,1),size(q,2)) :: albedo_loc + real(kind=rb),dimension(size(q,1),size(q,2)) :: cloud_coalb_spatial real(kind=rb),dimension(size(q,1),size(q,2),size(q,3)) :: q_tmp, h2o_vmr real(kind=rb),dimension(size(q,1),size(q,2)) :: fracsun real(kind=rb),dimension(size(q,1),size(q,2)) :: p2 !mp586 addition for annual mean insolation @@ -648,8 +656,15 @@ subroutine run_rrtmg(is,js,Time,lat,lon,p_full,p_half,albedo,q,t,t_surf_rad,tdt, !!!!! mp586 addition for annual mean insolation !!!!! !!!! following https://github.com/sit23/Isca/blob/master/src/atmos_param/socrates/interface/socrates_interface.F90#L888 !!!! - if (frierson_solar_rad) then - p2 = (1. - 3.*sin(lat(:,:))**2)/4. + p2 = (1. - 3.*sin(lat(:,:))**2)/4. + if (do_toa_albedo) then + !cloud_coalb_spatial(:,:) = (0.75 + 0.15 * 2. * p2(:,:))*alb_scaler - (alb_scaler - 1.) + cloud_coalb_spatial(:,:) = (0.87 + 0.05 * 2. * p2(:,:))*alb_scaler - (alb_scaler - 1.) + else + cloud_coalb_spatial(:,:) = 1. + endif + + if (frierson_solar_rad) then coszen = 0.25 * (1.0 + del_sol * p2 + del_sw * sin(lat(:,:))) rrsun = 1 ! needs to be set, set to 1 so that stellar_radiation is unchanged in socrates_interface else @@ -808,6 +823,7 @@ subroutine run_rrtmg(is,js,Time,lat,lon,p_full,p_half,albedo,q,t,t_surf_rad,tdt, cosz_rr = reshape(coszen (1:si:lonstep,:),(/ si*sj/lonstep /)) albedo_rr = reshape(albedo_loc(1:si:lonstep,:),(/ si*sj/lonstep /)) tsrf = reshape(t_surf_rad(1:si:lonstep,:),(/ si*sj/lonstep /)) + cloud_coalb_spatial_rr = reshape(cloud_coalb_spatial (1:si:lonstep,:),(/ si*sj/lonstep /)) !--------------------------------------------------------------------------------------------------------------- ! now actually run RRTM @@ -834,7 +850,7 @@ subroutine run_rrtmg(is,js,Time,lat,lon,p_full,p_half,albedo,q,t,t_surf_rad,tdt, pfull , phalf , tfull , thalf , tsrf , & h2o , o3 , co2 , ch4_val*ones , n2o_val*ones , o2_val*ones , & albedo_rr , albedo_rr, albedo_rr, albedo_rr, & - cosz_rr , solrad , dyofyr , solr_cnst, & + cosz_rr , cloud_coalb_spatial_rr, solrad , dyofyr , solr_cnst, & inflglw , iceflglw , liqflglw , & ! cloud parameters zeros , taucld , sw_zro , sw_zro , sw_zro , & @@ -848,7 +864,7 @@ subroutine run_rrtmg(is,js,Time,lat,lon,p_full,p_half,albedo,q,t,t_surf_rad,tdt, pfull , phalf , tfull , thalf , tsrf , & h2o , o3 , co2 , zeros , zeros, zeros, & albedo_rr , albedo_rr, albedo_rr, albedo_rr, & - cosz_rr , solrad , dyofyr , solr_cnst, & + cosz_rr , cloud_coalb_spatial_rr, solrad , dyofyr , solr_cnst, & inflglw , iceflglw , liqflglw , & ! cloud parameters zeros , taucld , sw_zro , sw_zro , sw_zro , & @@ -1007,6 +1023,21 @@ subroutine run_rrtmg(is,js,Time,lat,lon,p_full,p_half,albedo,q,t,t_surf_rad,tdt, enddo enddo endif + if(id_toa_sw_down > 0)then + swflxijk = reshape(swdflx(:,sk+1),(/ si/lonstep,sj /)) ! net TOA SW flux, +ve down + dlon=1./lonstep + do i=1,size(swijk,1) + i1 = i+1 + ! close toroidally + if(i1 > size(swijk,1)) i1=1 + do ij=1,lonstep + di = (ij-1)*dlon + ij1 = (i-1)*lonstep + ij + toa_sw_down(ij1,:) = di*swflxijk(i1,:) + (1.-di)*swflxijk(i ,:) + enddo + enddo + endif + @@ -1029,7 +1060,8 @@ subroutine write_diag_rrtm(Time,is,js,ozone,cotwo,fracday,albedo_loc, t_full) use rrtm_vars,only: sw_flux,lw_flux,zencos,tdt_rad,tdt_sw_rad,tdt_lw_rad,t_half,& &id_tdt_rad,id_tdt_sw,id_tdt_lw,id_coszen,& &id_flux_sw,id_flux_lw,id_albedo,id_ozone, id_co2, id_fracday,& - &id_olr,id_toa_sw,olr,toa_sw, id_half_level_temp, id_full_level_temp + &id_olr,id_toa_sw,id_toa_sw_down,olr,toa_sw,toa_sw_down,& + id_half_level_temp, id_full_level_temp use diag_manager_mod, only: register_diag_field, send_data use time_manager_mod,only: time_type @@ -1081,6 +1113,10 @@ subroutine write_diag_rrtm(Time,is,js,ozone,cotwo,fracday,albedo_loc, t_full) if ( id_toa_sw > 0 ) then used = send_data ( id_toa_sw, toa_sw, Time) endif +!------- Downward SW toa flux ------------ + if ( id_toa_sw_down > 0 ) then + used = send_data ( id_toa_sw_down, toa_sw_down, Time) + endif !------- Interactive albedo ------------ if ( present(albedo_loc)) then ! used = send_data ( id_albedo, albedo_loc, Time, is, js ) diff --git a/src/atmos_param/rrtm_radiation/rrtmg_lw/gcm_model/modules/parkind.f90 b/src/atmos_param/rrtm_radiation/rrtmg_lw/gcm_model/modules/parkind.F90 similarity index 65% rename from src/atmos_param/rrtm_radiation/rrtmg_lw/gcm_model/modules/parkind.f90 rename to src/atmos_param/rrtm_radiation/rrtmg_lw/gcm_model/modules/parkind.F90 index 5667d460b..8719337f9 100644 --- a/src/atmos_param/rrtm_radiation/rrtmg_lw/gcm_model/modules/parkind.f90 +++ b/src/atmos_param/rrtm_radiation/rrtmg_lw/gcm_model/modules/parkind.F90 @@ -15,7 +15,11 @@ module parkind ! integer kinds ! ------------- ! - integer, parameter :: kind_ib = selected_int_kind(13) ! 8 byte integer +#ifdef OVERLOAD_R4 + integer, parameter :: kind_ib = selected_int_kind(6) ! 4 byte integer +#else + integer, parameter :: kind_ib = selected_int_kind(13) ! 8 byte integer +#endif integer, parameter :: kind_im = selected_int_kind(6) ! 4 byte integer integer, parameter :: kind_in = kind(1) ! native integer @@ -23,7 +27,11 @@ module parkind ! real kinds ! ---------- ! - integer, parameter :: kind_rb = selected_real_kind(12) ! 8 byte real +#ifdef OVERLOAD_R4 + integer, parameter :: kind_rb = selected_real_kind(6) ! 4 byte real +#else + integer, parameter :: kind_rb = selected_real_kind(13) ! 8 byte real +#endif integer, parameter :: kind_rm = selected_real_kind(6) ! 4 byte real integer, parameter :: kind_rn = kind(1.0) ! native real diff --git a/src/atmos_param/rrtm_radiation/rrtmg_lw/gcm_model/src/rrtmg_lw_setcoef.f90 b/src/atmos_param/rrtm_radiation/rrtmg_lw/gcm_model/src/rrtmg_lw_setcoef.F90 similarity index 99% rename from src/atmos_param/rrtm_radiation/rrtmg_lw/gcm_model/src/rrtmg_lw_setcoef.f90 rename to src/atmos_param/rrtm_radiation/rrtmg_lw/gcm_model/src/rrtmg_lw_setcoef.F90 index cc65c677a..3d1afcf8d 100644 --- a/src/atmos_param/rrtm_radiation/rrtmg_lw/gcm_model/src/rrtmg_lw_setcoef.f90 +++ b/src/atmos_param/rrtm_radiation/rrtmg_lw/gcm_model/src/rrtmg_lw_setcoef.F90 @@ -253,8 +253,11 @@ subroutine setcoef(nlayers, istart, pavel, tavel, tz, tbound, semiss, & ! layer pressure. Store them in JP and JP1. Store in FP the ! fraction of the difference (in ln(pressure)) between these ! two values that the layer pressure lies. -! plog = alog(pavel(lay)) +#ifdef OVERLOAD_R4 + plog = alog(pavel(lay)) +#else plog = dlog(pavel(lay)) +#endif jp(lay) = int(36._rb - 5*(plog+0.04_rb)) if (jp(lay) .lt. 1) then jp(lay) = 1 diff --git a/src/atmos_param/rrtm_radiation/rrtmg_sw/gcm_model/modules/parkind.f90 b/src/atmos_param/rrtm_radiation/rrtmg_sw/gcm_model/modules/parkind.F90 similarity index 65% rename from src/atmos_param/rrtm_radiation/rrtmg_sw/gcm_model/modules/parkind.f90 rename to src/atmos_param/rrtm_radiation/rrtmg_sw/gcm_model/modules/parkind.F90 index 5667d460b..50a2b1f26 100644 --- a/src/atmos_param/rrtm_radiation/rrtmg_sw/gcm_model/modules/parkind.f90 +++ b/src/atmos_param/rrtm_radiation/rrtmg_sw/gcm_model/modules/parkind.F90 @@ -15,7 +15,12 @@ module parkind ! integer kinds ! ------------- ! - integer, parameter :: kind_ib = selected_int_kind(13) ! 8 byte integer + +#ifdef OVERLOAD_R4 + integer, parameter :: kind_ib = selected_int_kind(6) ! 4 byte integer +#else + integer, parameter :: kind_ib = selected_int_kind(13) ! 8 byte integer +#endif integer, parameter :: kind_im = selected_int_kind(6) ! 4 byte integer integer, parameter :: kind_in = kind(1) ! native integer @@ -23,7 +28,11 @@ module parkind ! real kinds ! ---------- ! - integer, parameter :: kind_rb = selected_real_kind(12) ! 8 byte real +#ifdef OVERLOAD_R4 + integer, parameter :: kind_rb = selected_real_kind(6) ! 4 byte real +#else + integer, parameter :: kind_rb = selected_real_kind(13) ! 8 byte real +#endif integer, parameter :: kind_rm = selected_real_kind(6) ! 4 byte real integer, parameter :: kind_rn = kind(1.0) ! native real diff --git a/src/atmos_param/rrtm_radiation/rrtmg_sw/gcm_model/src/rrtmg_sw_rad.f90 b/src/atmos_param/rrtm_radiation/rrtmg_sw/gcm_model/src/rrtmg_sw_rad.f90 index 1ed25daad..a9071c78d 100644 --- a/src/atmos_param/rrtm_radiation/rrtmg_sw/gcm_model/src/rrtmg_sw_rad.f90 +++ b/src/atmos_param/rrtm_radiation/rrtmg_sw/gcm_model/src/rrtmg_sw_rad.f90 @@ -82,7 +82,8 @@ subroutine rrtmg_sw & play ,plev ,tlay ,tlev ,tsfc , & h2ovmr , o3vmr ,co2vmr ,ch4vmr ,n2ovmr ,o2vmr , & asdir ,asdif ,aldir ,aldif , & - coszen ,adjes ,dyofyr ,scon , & + coszen , cloud_coalb_spatial, & + adjes ,dyofyr ,scon , & inflgsw ,iceflgsw,liqflgsw,cldfmcl , & taucmcl ,ssacmcl ,asmcmcl ,fsfcmcl , & ciwpmcl ,clwpmcl ,reicmcl ,relqmcl , & @@ -242,6 +243,7 @@ subroutine rrtmg_sw & real(kind=rb), intent(in) :: adjes ! Flux adjustment for Earth/Sun distance real(kind=rb), intent(in) :: coszen(:) ! Cosine of solar zenith angle ! Dimensions: (ncol) + real(kind=rb), intent(in) :: cloud_coalb_spatial(:) ! prescribed coalbedo for clouds real(kind=rb), intent(in) :: scon ! Solar constant (W/m2) integer(kind=im), intent(in) :: inflgsw ! Flag for cloud optical properties @@ -340,6 +342,7 @@ subroutine rrtmg_sw & ! real(kind=rb) :: earth_sun ! function for Earth/Sun distance factor real(kind=rb) :: cossza ! Cosine of solar zenith angle + real(kind=rb) :: cloud_coalb ! co albedo for prescribed clouds real(kind=rb) :: adjflux(jpband) ! adjustment for current Earth/Sun distance real(kind=rb) :: solvar(jpband) ! solar constant scaling factor from rrtmg_sw ! default value of 1368.22 Wm-2 at 1 AU @@ -569,6 +572,7 @@ subroutine rrtmg_sw & ! is below horizon cossza = coszen(iplon) + cloud_coalb = cloud_coalb_spatial(iplon) if (cossza .lt. zepzen) cossza = zepzen @@ -695,7 +699,7 @@ subroutine rrtmg_sw & (nlayers, istart, iend, icpr, idelm, iout, & pavel, tavel, pz, tz, tbound, albdif, albdir, & zcldfmc, ztaucmc, zasycmc, zomgcmc, ztaormc, & - ztaua, zasya, zomga, cossza, coldry, wkl, adjflux, & + ztaua, zasya, zomga, cossza, cloud_coalb, coldry, wkl, adjflux, & laytrop, layswtch, laylow, jp, jt, jt1, & co2mult, colch4, colco2, colh2o, colmol, coln2o, colo2, colo3, & fac00, fac01, fac10, fac11, & diff --git a/src/atmos_param/rrtm_radiation/rrtmg_sw/gcm_model/src/rrtmg_sw_rad.nomcica.f90 b/src/atmos_param/rrtm_radiation/rrtmg_sw/gcm_model/src/rrtmg_sw_rad.nomcica.f90 index c2d1e68e5..8bf4ddabb 100644 --- a/src/atmos_param/rrtm_radiation/rrtmg_sw/gcm_model/src/rrtmg_sw_rad.nomcica.f90 +++ b/src/atmos_param/rrtm_radiation/rrtmg_sw/gcm_model/src/rrtmg_sw_rad.nomcica.f90 @@ -80,7 +80,8 @@ subroutine rrtmg_sw & play ,plev ,tlay ,tlev ,tsfc , & h2ovmr ,o3vmr ,co2vmr ,ch4vmr ,n2ovmr ,o2vmr, & asdir ,asdif ,aldir ,aldif , & - coszen ,adjes ,dyofyr ,scon , & + coszen ,cloud_coalb_spatial, & + adjes ,dyofyr ,scon , & inflgsw ,iceflgsw,liqflgsw,cldfr , & taucld ,ssacld ,asmcld ,fsfcld , & cicewp ,cliqwp ,reice ,reliq , & @@ -228,6 +229,7 @@ subroutine rrtmg_sw & real(kind=rb), intent(in) :: coszen(:) ! Cosine of solar zenith angle ! Dimensions: (ncol) + real(kind=rb), intent(in) :: cloud_coalb_spatial(:) ! prescribed coalbedo for clouds real(kind=rb), intent(in) :: adjes ! Flux adjustment for Earth/Sun distance integer(kind=im), intent(in) :: dyofyr ! Day of the year (used to get Earth/Sun ! distance if adjflx not provided) @@ -329,6 +331,7 @@ subroutine rrtmg_sw & ! real(kind=rb) :: earth_sun ! function for Earth/Sun distance factor real(kind=rb) :: cossza ! Cosine of solar zenith angle + real(kind=rb) :: cloud_coalb ! co albedo for prescribed clouds real(kind=rb) :: adjflux(jpband) ! adjustment for current Earth/Sun distance real(kind=rb) :: solvar(jpband) ! solar constant scaling factor from rrtmg_sw ! default value of 1368.22 Wm-2 at 1 AU @@ -557,6 +560,7 @@ subroutine rrtmg_sw & ! is below horizon cossza = coszen(iplon) + cloud_coalb = cloud_coalb_spatial(iplon) if (cossza .lt. zepzen) cossza = zepzen @@ -680,7 +684,7 @@ subroutine rrtmg_sw & (nlayers, istart, iend, icpr, idelm, iout, & pavel, tavel, pz, tz, tbound, albdif, albdir, & cldfrac, ztauc, zasyc, zomgc, ztaucorig, & - ztaua, zasya, zomga, cossza, coldry, wkl, adjflux, & + ztaua, zasya, zomga, cossza, cloud_coalb, coldry, wkl, adjflux, & laytrop, layswtch, laylow, jp, jt, jt1, & co2mult, colch4, colco2, colh2o, colmol, coln2o, colo2, colo3, & fac00, fac01, fac10, fac11, & diff --git a/src/atmos_param/rrtm_radiation/rrtmg_sw/gcm_model/src/rrtmg_sw_spcvmc.f90 b/src/atmos_param/rrtm_radiation/rrtmg_sw/gcm_model/src/rrtmg_sw_spcvmc.f90 index 749b117f5..d448648d5 100644 --- a/src/atmos_param/rrtm_radiation/rrtmg_sw/gcm_model/src/rrtmg_sw_spcvmc.f90 +++ b/src/atmos_param/rrtm_radiation/rrtmg_sw/gcm_model/src/rrtmg_sw_spcvmc.f90 @@ -35,7 +35,7 @@ subroutine spcvmc_sw & (nlayers, istart, iend, icpr, idelm, iout, & pavel, tavel, pz, tz, tbound, palbd, palbp, & pcldfmc, ptaucmc, pasycmc, pomgcmc, ptaormc, & - ptaua, pasya, pomga, prmu0, coldry, wkl, adjflux, & + ptaua, pasya, pomga, prmu0, cloud_coalb, coldry, wkl, adjflux, & laytrop, layswtch, laylow, jp, jt, jt1, & co2mult, colch4, colco2, colh2o, colmol, coln2o, colo2, colo3, & fac00, fac01, fac10, fac11, & @@ -124,6 +124,7 @@ subroutine spcvmc_sw & real(kind=rb), intent(in) :: palbp(:) ! surface albedo (direct) ! Dimensions: (nbndsw) real(kind=rb), intent(in) :: prmu0 ! cosine of solar zenith angle + real(kind=rb), intent(in) :: cloud_coalb ! prescribed cloud co-albedo real(kind=rb), intent(in) :: pcldfmc(:,:) ! cloud fraction [mcica] ! Dimensions: (nlayers,ngptsw) real(kind=rb), intent(in) :: ptaucmc(:,:) ! cloud optical depth [mcica] @@ -313,7 +314,7 @@ subroutine spcvmc_sw & iw = iw+1 ! Apply adjustment for correct Earth/Sun distance and zenith angle to incoming solar flux - zincflx(iw) = adjflux(jb) * zsflxzen(iw) * prmu0 + zincflx(iw) = adjflux(jb) * zsflxzen(iw) * prmu0 * cloud_coalb ! zincflux = zincflux + adjflux(jb) * zsflxzen(iw) * prmu0 ! inactive ! Compute layer reflectances and transmittances for direct and diffuse sources, diff --git a/src/atmos_param/rrtm_radiation/rrtmg_sw/gcm_model/src/rrtmg_sw_spcvrt.f90 b/src/atmos_param/rrtm_radiation/rrtmg_sw/gcm_model/src/rrtmg_sw_spcvrt.f90 index 06c67fb8d..b92030be0 100644 --- a/src/atmos_param/rrtm_radiation/rrtmg_sw/gcm_model/src/rrtmg_sw_spcvrt.f90 +++ b/src/atmos_param/rrtm_radiation/rrtmg_sw/gcm_model/src/rrtmg_sw_spcvrt.f90 @@ -35,7 +35,7 @@ subroutine spcvrt_sw & (nlayers, istart, iend, icpr, idelm, iout, & pavel, tavel, pz, tz, tbound, palbd, palbp, & pclfr, ptauc, pasyc, pomgc, ptaucorig, & - ptaua, pasya, pomga, prmu0, coldry, wkl, adjflux, & + ptaua, pasya, pomga, prmu0, cloud_coalb, coldry, wkl, adjflux, & laytrop, layswtch, laylow, jp, jt, jt1, & co2mult, colch4, colco2, colh2o, colmol, coln2o, colo2, colo3, & fac00, fac01, fac10, fac11, & @@ -121,6 +121,7 @@ subroutine spcvrt_sw & real(kind=rb), intent(in) :: palbp(:) ! surface albedo (direct) ! Dimensions: (nbndsw) real(kind=rb), intent(in) :: prmu0 ! cosine of solar zenith angle + real(kind=rb), intent(in) :: cloud_coalb ! prescribed cloud co-albedo real(kind=rb), intent(in) :: pclfr(:) ! cloud fraction ! Dimensions: (nlayers) real(kind=rb), intent(in) :: ptauc(:,:) ! cloud optical depth @@ -312,7 +313,7 @@ subroutine spcvrt_sw & iw = iw+1 ! Apply adjustments for correct Earth/Sun distance and zenith angle to incoming solar flux - zincflx(iw) = adjflux(jb) * zsflxzen(iw) * prmu0 + zincflx(iw) = adjflux(jb) * zsflxzen(iw) * prmu0 * cloud_coalb ! zincflux = zincflux + adjflux(jb) * zsflxzen(iw) * prmu0 ! inactive ! Compute layer reflectances and transmittances for direct and diffuse sources, From 12d0baed02e6f57b543892205038f549df69f79a Mon Sep 17 00:00:00 2001 From: ntlewis Date: Fri, 21 Jul 2023 10:07:11 +0100 Subject: [PATCH 04/19] let my25 read restart so it actually works --- src/atmos_param/my25_turb/my25_turb.F90 | 15 ++++++++++----- 1 file changed, 10 insertions(+), 5 deletions(-) diff --git a/src/atmos_param/my25_turb/my25_turb.F90 b/src/atmos_param/my25_turb/my25_turb.F90 index 989d342db..fd7f60358 100644 --- a/src/atmos_param/my25_turb/my25_turb.F90 +++ b/src/atmos_param/my25_turb/my25_turb.F90 @@ -27,14 +27,15 @@ MODULE MY25_TURB_MOD use mpp_mod, only : input_nml_file use fms_mod, only : file_exist, open_namelist_file, error_mesg, & - FATAL, close_file, note, read_data, & + FATAL, close_file, note, read_data, write_data, & check_nml_error, mpp_pe, mpp_root_pe, & - write_version_number, stdlog, open_restart_file + write_version_number, stdlog, open_restart_file, nullify_domain use fms_io_mod, only : register_restart_field, restart_file_type, & save_restart, restore_state use tridiagonal_mod, only : tri_invert, close_tridiagonal use constants_mod, only : grav, vonkarm use monin_obukhov_mod, only : mo_diff + use transforms_mod, only: grid_domain !--------------------------------------------------------------------- implicit none @@ -706,14 +707,16 @@ SUBROUTINE MY25_TURB_INIT( ix, jx, kx ) ! --- Input TKE !--------------------------------------------------------------------- - id_restart = register_restart_field(Tur_restart, 'my25_turb.res', 'TKE', TKE) + id_restart = register_restart_field(Tur_restart, 'my25_turb.res.nc', 'TKE', TKE) if (file_exist( 'INPUT/my25_turb.res.nc' )) then if (mpp_pe() == mpp_root_pe() ) then call error_mesg ('my25_turb_mod', 'MY25_TURB_INIT:& &Reading netCDF formatted restart file: & &INPUT/my25_turb.res.nc', NOTE) endif - call restore_state(Tur_restart) + !call restore_state(Tur_restart) + call nullify_domain() + call read_data(trim('INPUT/my25_turb.res'), 'TKE', TKE, grid_domain) else if( FILE_EXIST( 'INPUT/my25_turb.res' ) ) then @@ -781,7 +784,9 @@ subroutine my25_turb_restart(timestamp) ! write out the restart data which is always needed, regardless of ! when the first donner calculation step is after restart. !---------------------------------------------------------------------- - call save_restart(Tur_restart, timestamp) + !call save_restart(Tur_restart, timestamp) + call nullify_domain() + call write_data(trim('RESTART/my25_turb.res'), 'TKE', TKE, grid_domain) end subroutine my25_turb_restart ! NAME="my25_turb_restart" From 8021f3d40a6419245e4b683f90bd7e984ddc82f5 Mon Sep 17 00:00:00 2001 From: ntlewis Date: Fri, 21 Jul 2023 10:07:48 +0100 Subject: [PATCH 05/19] additions to two_stream_gray for tidally locked planets --- .../two_stream_gray_rad.F90 | 148 ++++++++++++++++-- 1 file changed, 134 insertions(+), 14 deletions(-) diff --git a/src/atmos_param/two_stream_gray_rad/two_stream_gray_rad.F90 b/src/atmos_param/two_stream_gray_rad/two_stream_gray_rad.F90 index da2cf4c60..0a3fd3865 100644 --- a/src/atmos_param/two_stream_gray_rad/two_stream_gray_rad.F90 +++ b/src/atmos_param/two_stream_gray_rad/two_stream_gray_rad.F90 @@ -28,7 +28,7 @@ module two_stream_gray_rad_mod mpp_pe, close_file, error_mesg, & NOTE, FATAL, uppercase - use constants_mod, only: stefan, cp_air, grav, pstd_mks, pstd_mks_earth, seconds_per_sol, orbital_period + use constants_mod, only: stefan, cp_air, grav, pstd_mks, pstd_mks_earth, seconds_per_sol, orbital_period, radius use diag_manager_mod, only: register_diag_field, send_data @@ -85,6 +85,11 @@ module two_stream_gray_rad_mod real :: equinox_day = 0.75 !s Fraction of year [0,1] where NH autumn equinox occurs (only really useful if calendar has defined months). logical :: use_time_average_coszen = .false. !s if .true., then time-averaging is done on coszen so that insolation doesn't depend on timestep real :: dt_rad_avg = -1 +logical :: do_tl = .false. +logical :: do_closein = .false. +logical :: do_normal_integration_method = .true. +real :: R_stellar = 476700000. +real :: d_stellar = 1117512000.0 character(len=32) :: rad_scheme = 'frierson' @@ -119,6 +124,7 @@ module two_stream_gray_rad_mod real :: bog_mu = 1.0 real, allocatable, dimension(:,:) :: insolation, p2, lw_tau_0, sw_tau_0 !s albedo now defined in mixed_layer_init +real, allocatable, dimension(:,:) :: lat_tl, lat_tl_insol real, allocatable, dimension(:,:) :: b_surf, b_surf_gp real, allocatable, dimension(:,:,:) :: b, tdt_rad, tdt_solar real, allocatable, dimension(:,:,:) :: lw_up, lw_down, lw_flux, sw_up, sw_down, sw_flux, rad_flux @@ -150,14 +156,14 @@ module two_stream_gray_rad_mod do_read_co2, co2_file, co2_variable_name, solday, equinox_day, bog_a, bog_b, bog_mu, & use_time_average_coszen, dt_rad_avg,& diabatic_acce, & !Schneider Liu values - do_toa_albedo + do_toa_albedo, do_tl, do_closein, R_stellar, d_stellar, do_normal_integration_method !================================================================================== !-------------------- diagnostics fields ------------------------------- integer :: id_olr, id_swdn_sfc, id_swdn_toa, id_net_lw_surf, id_lwdn_sfc, id_lwup_sfc, & id_tdt_rad, id_tdt_solar, id_flux_rad, id_flux_lw, id_flux_sw, id_coszen, id_fracsun, & - id_lw_dtrans, id_lw_dtrans_win, id_sw_dtrans, id_co2 + id_lw_dtrans, id_lw_dtrans_win, id_sw_dtrans, id_co2, id_nudge_flux, id_tau character(len=10), parameter :: mod_name = 'two_stream' @@ -268,6 +274,8 @@ subroutine two_stream_gray_rad_init(is, ie, js, je, num_levels, axes, Time, lonb allocate (olr (ie-is+1, je-js+1)) allocate (net_lw_surf (ie-is+1, je-js+1)) allocate (toa_sw_in (ie-is+1, je-js+1)) +allocate (lat_tl (ie-is+1, je-js+1)) +allocate (lat_tl_insol (ie-is+1, je-js+1)) allocate (insolation (ie-is+1, je-js+1)) allocate (p2 (ie-is+1, je-js+1)) @@ -354,6 +362,12 @@ subroutine two_stream_gray_rad_init(is, ie, js, je, num_levels, axes, Time, lonb register_diag_field ( mod_name, 'coszen', axes(1:2), Time, & 'cosine of zenith angle', & 'none', missing_value=missing_value ) + + id_tau = & + register_diag_field ( mod_name, 'tau', axes(1:2), Time, & + 'surface lw optical depth', & + 'none', missing_value=missing_value ) + id_fracsun = & register_diag_field ( mod_name, 'fracsun', axes(1:2), Time, & 'daylight fraction of time interval', & @@ -380,13 +394,17 @@ subroutine two_stream_gray_rad_init(is, ie, js, je, num_levels, axes, Time, lonb register_diag_field ( mod_name, 'lw_dtrans', axes(1:3), Time, & 'LW transmission (non window)', & 'none', missing_value=missing_value ) + id_nudge_flux = & + register_diag_field ( mod_name, 'nudge_flux', axes(1:2), Time, & + 'Sea ice nudging flux', & + 'W/m2', missing_value=missing_value ) return end subroutine two_stream_gray_rad_init ! ================================================================================== subroutine two_stream_gray_rad_down (is, js, Time_diag, lat, lon, p_half, t, & - net_surf_sw_down, surf_lw_down, albedo, q) + net_surf_sw_down, surf_lw_down, albedo, q, const_correct, nudge_out) ! NTL ICE CHANGE ! Begin the radiation calculation by computing downward fluxes. ! This part of the calculation does not depend on the surface temperature. @@ -397,12 +415,18 @@ subroutine two_stream_gray_rad_down (is, js, Time_diag, lat, lon, p_half, t, real, intent(out), dimension(:,:) :: net_surf_sw_down real, intent(out), dimension(:,:) :: surf_lw_down real, intent(in), dimension(:,:,:) :: t, q, p_half +real, intent(in), dimension(:,:) :: const_correct, nudge_out integer :: i, j, k, n, dyofyr integer :: seconds, year_in_s, days real :: r_seconds, frac_of_day, frac_of_year, gmt, time_since_ae, rrsun, day_in_s, r_solday, r_total_seconds, r_days, r_dt_rad_avg, dt_rad_radians logical :: used +! all for closein insolation +real :: ci_lat_f, ci_lat_d +real, dimension(size(q,1), size(q,2)) :: ci_rho_dist, ci_J_f, ci_alpha, ci_beta, ci_phi, & + ci_t, ci_l, ci_s, ci_delta1, ci_delta2, ci_J_p + real ,dimension(size(q,1),size(q,2),size(q,3)) :: co2f @@ -448,6 +472,63 @@ subroutine two_stream_gray_rad_down (is, js, Time_diag, lat, lon, p_half, t, end if insolation = solar_constant * coszen + + if (do_toa_albedo) then + p2 = (1. - 3.*sin(lat)**2)/4. + insolation = insolation * (0.75 + 0.15 * 2. * p2) + endif + + insolation = insolation + +else if (do_tl) then + + + if (do_closein) then + + lat_tl_insol = asin(cos(lat)*cos(lon)) + pi/2 + + ci_lat_f = acos((radius+R_stellar) / d_stellar) + ci_lat_d = acos((radius-R_stellar) / d_stellar) + + ci_rho_dist = sqrt(d_stellar**2 + radius**2 - 2*d_stellar*radius*cos(lat_tl_insol)) + + ! J_f + ci_J_f = (d_stellar * cos(lat_tl_insol) - radius)*R_stellar**2 / ci_rho_dist**3 + + ! J_p + ci_alpha = acos((d_stellar*cos(lat_tl_insol) - radius) / ci_rho_dist) + ci_beta = pi/2 - ci_alpha + ci_phi = asin(R_stellar / ci_rho_dist) + ci_t = ci_rho_dist * cos(ci_phi) + ci_l = ci_t * sin(ci_phi) + ci_s = ci_t * sin(ci_phi) * cos(ci_alpha) + ci_delta1 = acos(cos(ci_phi)/cos(ci_beta)) + ci_delta2 = acos(tan(ci_beta)/tan(ci_phi)) + ci_J_p = ci_l * ci_s / pi / ci_t**2 * (pi - ci_delta2 + sin(ci_delta2)*cos(ci_delta2)) + 1 / pi * (ci_delta1 - sin(ci_delta1)*cos(ci_delta1)) + + where (lat_tl_insol .le. ci_lat_f) + insolation = solar_constant * (d_stellar / R_stellar)**2 * ci_J_f + elsewhere ((lat_tl_insol .gt. ci_lat_f) .and. (lat_tl_insol .lt. ci_lat_d)) + insolation = solar_constant * (d_stellar / R_stellar)**2 * ci_J_p + elsewhere (lat_tl_insol .ge. ci_lat_d) + insolation = 0. + endwhere + + + + + else + + insolation = solar_constant * cos(lat) * cos(lon) + + where (insolation .le. 0. ) + insolation = 0. + endwhere + endif + + + coszen = insolation / solar_constant + else if (sw_scheme==B_SCHNEIDER_LIU) then insolation = (solar_constant/pi)*cos(lat) @@ -456,8 +537,9 @@ subroutine two_stream_gray_rad_down (is, js, Time_diag, lat, lon, p_half, t, p2 = (1. - 3.*sin(lat)**2)/4. insolation = 0.25 * solar_constant * (1.0 + del_sol * p2 + del_sw * sin(lat)) if (do_toa_albedo) then - insolation = insolation * (0.77 + 0.125 * 2. * p2) + insolation = insolation * (0.75 + 0.15 * 2. * p2) endif + insolation = insolation end if select case(sw_scheme) @@ -580,7 +662,15 @@ subroutine two_stream_gray_rad_down (is, js, Time_diag, lat, lon, p_half, t, case(B_FRIERSON) ! longwave optical thickness function of latitude and pressure - lw_tau_0 = ir_tau_eq + (ir_tau_pole - ir_tau_eq)*sin(lat)**2 + if (do_tl) then + lat_tl = pi/2 - asin(cos(lat)*cos(lon)) + lw_tau_0 = ir_tau_eq + (ir_tau_pole - ir_tau_eq)*sin(lat_tl)**2 + where (lat_tl .ge. (pi/2)) + lw_tau_0 = ir_tau_pole + endwhere + else + lw_tau_0 = ir_tau_eq + (ir_tau_pole - ir_tau_eq)*sin(lat)**2 + endif lw_tau_0 = lw_tau_0 * odp ! scale by optical depth parameter - default 1 ! compute optical depths for each model level @@ -596,9 +686,17 @@ subroutine two_stream_gray_rad_down (is, js, Time_diag, lat, lon, p_half, t, ! compute downward longwave flux by integrating downward lw_down(:,:,1) = 0. - do k = 1, n - lw_down(:,:,k+1) = lw_down(:,:,k)*lw_dtrans(:,:,k) + b(:,:,k)*(1. - lw_dtrans(:,:,k)) - end do + if (do_normal_integration_method) then + do k = 1, n + lw_down(:,:,k+1) = lw_down(:,:,k)*lw_dtrans(:,:,k) + b(:,:,k)*(1. - lw_dtrans(:,:,k)) + end do + + else + do k = 1, n + lw_down(:,:,k+1) = 2.*b(:,:,k) * (lw_tau(:,:,k+1) - lw_tau(:,:,k))/(2.+ (lw_tau(:,:,k+1) - lw_tau(:,:,k))) + & + lw_down(:,:,k) * (2.- (lw_tau(:,:,k+1) - lw_tau(:,:,k)))/(2.+ (lw_tau(:,:,k+1) - lw_tau(:,:,k))) + end do + endif case(B_SCHNEIDER_LIU) @@ -625,9 +723,9 @@ subroutine two_stream_gray_rad_down (is, js, Time_diag, lat, lon, p_half, t, end select ! ================================================================================= -surf_lw_down = lw_down(:, :, n+1) +surf_lw_down = lw_down(:, :, n+1) toa_sw_in = sw_down(:, :, 1) -net_surf_sw_down = sw_down(:, :, n+1) * (1. - albedo) +net_surf_sw_down = sw_down(:, :, n+1) * (1. - albedo) ! ================================================================================= if(lw_scheme.eq.B_SCHNEIDER_LIU) then @@ -639,10 +737,20 @@ subroutine two_stream_gray_rad_down (is, js, Time_diag, lat, lon, p_half, t, if ( id_lwdn_sfc > 0 ) then used = send_data ( id_lwdn_sfc, surf_lw_down, Time_diag) endif + +!surf_lw_down = surf_lw_down + const_correct + nudge_out ! NTL UPDATE NUDGING + !------- incoming sw flux toa ------- if ( id_swdn_toa > 0 ) then used = send_data ( id_swdn_toa, toa_sw_in, Time_diag) endif + +! NTL NUDGING +if( id_nudge_flux > 0) then + used = send_data ( id_nudge_flux, const_correct+nudge_out, Time_diag) +endif + + !------- downward sw flux surface ------- if ( id_swdn_sfc > 0 ) then used = send_data ( id_swdn_sfc, net_surf_sw_down, Time_diag) @@ -653,6 +761,10 @@ subroutine two_stream_gray_rad_down (is, js, Time_diag, lat, lon, p_half, t, used = send_data ( id_coszen, coszen, Time_diag) endif +if (id_tau > 0) then + used = send_data ( id_tau, lw_tau_0, Time_diag) +endif + !------- carbon dioxide concentration ------------ if ( id_co2 > 0 ) then used = send_data ( id_co2, carbon_conc, Time_diag) @@ -697,9 +809,17 @@ subroutine two_stream_gray_rad_up (is, js, Time_diag, lat, p_half, t_surf, t, td case(B_FRIERSON, B_BYRNE) ! compute upward longwave flux by integrating upward lw_up(:,:,n+1) = b_surf - do k = n, 1, -1 - lw_up(:,:,k) = lw_up(:,:,k+1)*lw_dtrans(:,:,k) + b(:,:,k)*(1.0 - lw_dtrans(:,:,k)) - end do + if (do_normal_integration_method) then + do k = n, 1, -1 + lw_up(:,:,k) = lw_up(:,:,k+1)*lw_dtrans(:,:,k) + b(:,:,k)*(1.0 - lw_dtrans(:,:,k)) + end do + + else + do k = n, 1, -1 + lw_up(:,:,k) = 2.*b(:,:,k) * -1.*(lw_tau(:,:,k) - lw_tau(:,:,k+1))/ (2.- (lw_tau(:,:,k) - lw_tau(:,:,k+1))) & + + lw_up(:,:,k+1)* (2. + (lw_tau(:,:,k) - lw_tau(:,:,k+1)))/(2. - (lw_tau(:,:,k) - lw_tau(:,:,k+1))) + end do + endif case(B_SCHNEIDER_LIU) ! compute upward longwave flux by integrating upward From 92f222dc6d616e0f915c198fd2255da9d44be8ea Mon Sep 17 00:00:00 2001 From: ntlewis Date: Fri, 21 Jul 2023 10:11:22 +0100 Subject: [PATCH 06/19] add additional dry convection scheme (+ other changes?) --- src/atmos_param/dry_adj/dry_adj.F90 | 327 ++++++++++++++++++ .../driver/solo/idealized_moist_phys.F90 | 76 +++- 2 files changed, 392 insertions(+), 11 deletions(-) create mode 100644 src/atmos_param/dry_adj/dry_adj.F90 diff --git a/src/atmos_param/dry_adj/dry_adj.F90 b/src/atmos_param/dry_adj/dry_adj.F90 new file mode 100644 index 000000000..a15c99cca --- /dev/null +++ b/src/atmos_param/dry_adj/dry_adj.F90 @@ -0,0 +1,327 @@ +MODULE DRY_ADJ_MOD + + !======================================================================= + ! DRY ADIABATIC ADJUSTMENT + !======================================================================= + + + + use fms_mod, ONLY: file_exist, error_mesg, open_file, & + check_nml_error, mpp_pe, FATAL, & + NOTE, close_file, uppercase + + use constants_mod, ONLY: grav, kappa + !--------------------------------------------------------------------- + implicit none + private + + public :: dry_adj, dry_adj_init, dry_adj_end, dry_adj_bdgt + + !--------------------------------------------------------------------- + + character(len=128) :: version = '$Id: dry_adj.F90,v 17.0.4.1 2010/08/30 20:33:34 wfc Exp $' + character(len=128) :: tag = '$Name: testing $' + logical :: module_is_initialized = .false. + + !--------------------------------------------------------------------- + + real, parameter :: p00 = 1000.0E2 + real :: kappa_use = 0.0 ! modified kappa + ! (kappa_use = kappa * alpha) + ! initialised in dry_adj_init + + !--------------------------------------------------------------------- + ! --- NAMELIST + !--------------------------------------------------------------------- + ! itermax - Max number of iterations + !--------------------------------------------------------------------- + + integer :: itermax = 5 + real :: small = 0.001 + logical :: do_mcm_dry_adj = .true. + logical :: do_relax = .false. + real :: tau_relax = 7200. + real :: alpha = 1.0 ! values less than one achieve stable column. + + namelist /dry_adj_nml/ itermax, small, do_mcm_dry_adj, do_relax, tau_relax!, alpha + + !--------------------------------------------------------------------- + + contains + + !####################################################################### + !####################################################################### + + SUBROUTINE DRY_ADJ ( temp0, pres, pres_int, dtemp, timestep, mask ) + + !======================================================================= + ! DRY ADIABATIC ADJUSTMENT + !======================================================================= + !--------------------------------------------------------------------- + ! Arguments (Intent in) + ! temp0 - Temperature + ! pres - Pressure + ! pres_int - Pressure at layer interface + ! mask - OPTIONAL; floating point mask (0. or 1.) designating + ! where data is present + !--------------------------------------------------------------------- + real, intent(in), dimension(:,:,:) :: temp0, pres, pres_int + real, intent(in) :: timestep + real, intent(in), OPTIONAL, dimension(:,:,:) :: mask + + !--------------------------------------------------------------------- + ! Arguments (Intent out) + ! dtemp - Change in temperature + !--------------------------------------------------------------------- + real, intent(out), dimension(:,:,:) :: dtemp + + !--------------------------------------------------------------------- + ! (Intent local) + !--------------------------------------------------------------------- + + real, dimension(size(temp0,1),size(temp0,2),size(temp0,3)) :: & + temp, pi, theta, pixdp, dpres, ppp + + real, dimension(size(temp0,1),size(temp0,2)) :: store, xxx + logical, dimension(size(temp0,1),size(temp0,2)) :: do_daa + + integer :: kmax, iter, k, i, j + logical :: do_any, did_adj + + !==================================================================== + + ! --- Check to see if dry_adj has been initialized + if(.not. module_is_initialized ) CALL error_mesg( 'DRY_ADJ', & + 'dry_adj_init has not been called', FATAL ) + + ! --- Set dimensions + kmax = size( temp0, 3 ) + + ! --- Compute pressure thickness of layers + dpres(:,:,1:kmax) = pres_int(:,:,2:kmax+1) - pres_int(:,:,1:kmax) + + ! --- Copy input temperature + temp = temp0 + + ! --- Compute exner function + !do i = 1,size(temp0,1) + ! do j = 1,size(temp0,2) + ! pi(i,j,:) = (pres(i,j,:) / pres(i,j,size(temp0,3))) ** kappa_use + ! end do + !end do + ! --- Compute exner function + pi = ( pres / p00 ) ** kappa_use + + ! --- Compute product of pi and dpres + pixdp = pi * dpres + + ! --- Compute potential temperature + theta = temp / pi + + if(do_mcm_dry_adj) then + do k = 2,kmax + xxx = 0.5*kappa_use*(pres(:,:,k) - pres(:,:,k-1))/pres_int(:,:,k) + ppp(:,:,k) = (1.0 + xxx)/(1.0 - xxx) + enddo + endif + + !----------------------------------------------------------------- + ! iteration loop starts + !----------------------------------------------------------------- + do iter = 1,itermax + !----------------------------------------------------------------- + + did_adj = .false. + + do k = 1,kmax - 1 + ! ---------------------------------------------- + + ! --- Flag layers needing adjustment + if(do_mcm_dry_adj) then + do_daa(:,:) = temp(:,:,k+1) > ( temp(:,:,k)*ppp(:,:,k+1) + small ) + else + do_daa(:,:) = ( theta(:,:,k+1) - theta(:,:,k) ) > small + endif + + if( PRESENT( mask ) ) then + do_daa(:,:) = do_daa(:,:) .and. ( mask(:,:,k+1) > 0.5 ) + endif + do_any = ANY( do_daa(:,:) ) + + ! --- Do adjustment + if ( do_any ) then + if(do_mcm_dry_adj) then + where ( do_daa ) + temp(:,:,k) = (temp(:,:,k) * dpres(:,:,k ) + & + temp(:,:,k+1)* dpres(:,:,k+1) ) & + /(dpres(:,:,k) + ppp(:,:,k+1)*dpres(:,:,k+1)) + temp(:,:,k+1) = temp(:,:,k)*ppp(:,:,k+1) + end where + did_adj = .true. + else + where ( do_daa ) + store(:,:) = ( theta(:,:,k ) * pixdp(:,:,k ) & + + theta(:,:,k+1) * pixdp(:,:,k+1) ) & + / ( pixdp(:,:,k ) + pixdp(:,:,k+1) ) + theta(:,:,k ) = store(:,:) + theta(:,:,k+1) = store(:,:) + temp(:,:,k ) = pi(:,:,k ) * theta(:,:,k ) + temp(:,:,k+1) = pi(:,:,k+1) * theta(:,:,k+1) + end where + did_adj = .true. + endif + end if + + ! ---------------------------------------------- + end do + + ! --- If no adjusment made this pass, exit iteration loop. + if ( .not. did_adj ) go to 900 + + !----------------------------------------------------------------- + end do + !----------------------------------------------------------------- + ! iteration loop ends + !----------------------------------------------------------------- + if(.not.do_mcm_dry_adj) then + call error_mesg ('DRY_ADJ', 'Non-convergence in dry_adj', NOTE) + endif + 900 continue + + ! --- Compute change in temperature + if (do_relax) then + dtemp = (temp - temp0) * (1 - exp(-timestep / tau_relax)) + else + dtemp = temp - temp0 + endif + + !======================================================================= + end SUBROUTINE DRY_ADJ + + !##################################################################### + !##################################################################### + + SUBROUTINE DRY_ADJ_INIT(alpha_in) + + !======================================================================= + ! ***** INITIALIZE RAS + !======================================================================= + + real, intent(in) :: alpha_in + !--------------------------------------------------------------------- + ! (Intent local) + !--------------------------------------------------------------------- + + integer :: unit, io, ierr + + !===================================================================== + + !--------------------------------------------------------------------- + ! --- READ NAMELIST + !--------------------------------------------------------------------- + + if (file_exist('input.nml')) then + unit = open_file (file='input.nml', action='read') + ierr = 1 + do while (ierr /= 0) + read (unit, nml=dry_adj_nml, iostat=io, end=10) + ierr = check_nml_error(io, 'dry_adj_nml') + end do + 10 call close_file(unit) + endif + + !------- write version number and namelist --------- + + unit = open_file (file='logfile.out', action='append') + if ( mpp_pe() == 0 ) then + write (unit,'(/,80("="),/(a))') trim(version), trim(tag) + write (unit,nml=dry_adj_nml) + endif + call close_file(unit) + + !------------------------------------------------------------------- + + alpha = alpha_in + + kappa_use = kappa * alpha + + module_is_initialized = .true. + + !===================================================================== + end SUBROUTINE DRY_ADJ_INIT + + + !####################################################################### + !####################################################################### + + SUBROUTINE DRY_ADJ_END() + + !------------------------------------------------------------------- + + module_is_initialized = .false. + + !===================================================================== + end SUBROUTINE DRY_ADJ_END + + + !####################################################################### + !####################################################################### + + SUBROUTINE DRY_ADJ_BDGT ( dtemp, pres_int ) + + !======================================================================= + ! Budget check for dry adiabatic adjustment - a debugging tool + !======================================================================= + + !--------------------------------------------------------------------- + ! Arguments (Intent in) + ! dtemp - Temperature change + ! pres_int - Pressure at layer interface + !--------------------------------------------------------------------- + real, intent(in), dimension(:,:,:) :: dtemp, pres_int + + !--------------------------------------------------------------------- + ! (Intent local) + !--------------------------------------------------------------------- + + real, dimension(size(dtemp,1),size(dtemp,2),size(dtemp,3)) :: dpres + real :: sum_dtemp + integer :: imax, jmax, kmax, i, j, k + + !======================================================================= + + imax = size ( dtemp, 1 ) + jmax = size ( dtemp, 2 ) + kmax = size ( dtemp, 3 ) + + ! --- Compute pressure thickness of layers + dpres(:,:,1:kmax) = pres_int(:,:,2:kmax+1) - pres_int(:,:,1:kmax) + + ! --- Check budget + + do j = 1,jmax + do i = 1,imax + + sum_dtemp = 0. + + do k = 1,kmax + sum_dtemp = sum_dtemp + dtemp(i,j,k)*dpres(i,j,k) / Grav + end do + + if ( abs( sum_dtemp ) > 1.0e-4 ) then + print * + print *, ' DRY ADIABATIC ADJUSTMENT BUDGET CHECK AT i,j = ', i,j + print *, ' sum_dtemp = ', sum_dtemp + print *, 'STOP' + ! STOP + endif + + end do + end do + + !======================================================================= + end SUBROUTINE DRY_ADJ_BDGT + + !####################################################################### + !####################################################################### + end MODULE DRY_ADJ_MOD \ No newline at end of file diff --git a/src/atmos_spectral/driver/solo/idealized_moist_phys.F90 b/src/atmos_spectral/driver/solo/idealized_moist_phys.F90 index 490c7a1e0..b1e8f1efc 100644 --- a/src/atmos_spectral/driver/solo/idealized_moist_phys.F90 +++ b/src/atmos_spectral/driver/solo/idealized_moist_phys.F90 @@ -33,6 +33,8 @@ module idealized_moist_phys_mod use dry_convection_mod, only: dry_convection_init, dry_convection +use dry_adj_mod, only: dry_adj_init, dry_adj, dry_adj_end + use diag_manager_mod, only: register_diag_field, send_data use transforms_mod, only: get_grid_domain @@ -102,13 +104,15 @@ module idealized_moist_phys_mod SIMPLE_BETTS_CONV = 1, & FULL_BETTS_MILLER_CONV = 2,& DRY_CONV = 3, & - RAS_CONV = 4 + RAS_CONV = 4, & + DRYADJ_CONV = 5 integer :: r_conv_scheme = UNSET ! the selected convection scheme logical :: lwet_convection = .false. logical :: do_bm = .false. logical :: do_ras = .false. +logical :: do_dryadj = .false. ! Cloud options logical :: do_cloud_simple = .false. ! by default the cloud scheme is off. @@ -147,7 +151,9 @@ module idealized_moist_phys_mod real :: raw_bucket = 0.53 ! default raw coefficient for bucket depth LJJ ! end RG Add bucket -namelist / idealized_moist_phys_nml / turb, lwet_convection, do_bm, do_ras, roughness_heat, & +real :: alpha = 1.0 + +namelist / idealized_moist_phys_nml / turb, lwet_convection, do_bm, do_ras, do_dryadj, roughness_heat, & do_cloud_simple, & two_stream_gray, do_rrtm_radiation, do_damping,& mixed_layer_bc, do_simple, & @@ -167,7 +173,7 @@ module idealized_moist_phys_mod real, allocatable, dimension(:,:) :: & z_surf, & ! surface height t_surf, & ! surface temperature - t_ml, h_thermo_ice, & ! NTL 01/23 thermodynamic sea ice + t_ml, h_thermo_ice, const_correct, nudge_out, & ! NTL 01/23 thermodynamic sea ice q_surf, & ! surface moisture u_surf, & ! surface U wind v_surf, & ! surface V wind @@ -372,6 +378,7 @@ subroutine idealized_moist_phys_init(Time, Time_step_in, nhum, rad_lon_2d, rad_l lwet_convection = .false. do_bm = .false. do_ras = .false. + do_dryadj = .false. call error_mesg('idealized_moist_phys','No convective adjustment scheme used.', NOTE) else if(uppercase(trim(convection_scheme)) == 'SIMPLE_BETTS_MILLER') then @@ -380,6 +387,7 @@ subroutine idealized_moist_phys_init(Time, Time_step_in, nhum, rad_lon_2d, rad_l lwet_convection = .true. do_bm = .false. do_ras = .false. + do_dryadj = .false. else if(uppercase(trim(convection_scheme)) == 'FULL_BETTS_MILLER') then r_conv_scheme = FULL_BETTS_MILLER_CONV @@ -387,6 +395,7 @@ subroutine idealized_moist_phys_init(Time, Time_step_in, nhum, rad_lon_2d, rad_l do_bm = .true. lwet_convection = .false. do_ras = .false. + do_dryadj = .false. else if(uppercase(trim(convection_scheme)) == 'RAS') then r_conv_scheme = RAS_CONV @@ -394,6 +403,7 @@ subroutine idealized_moist_phys_init(Time, Time_step_in, nhum, rad_lon_2d, rad_l do_ras = .true. do_bm = .false. lwet_convection = .false. + do_dryadj = .false. else if(uppercase(trim(convection_scheme)) == 'DRY') then r_conv_scheme = DRY_CONV @@ -401,6 +411,15 @@ subroutine idealized_moist_phys_init(Time, Time_step_in, nhum, rad_lon_2d, rad_l lwet_convection = .false. do_bm = .false. do_ras = .false. + do_dryadj = .false. + +else if(uppercase(trim(convection_scheme)) == 'DRYADJ') then + r_conv_scheme = DRYADJ_CONV + call error_mesg('idealized_moist_phys','Using AM3 dry adjustment convection scheme.', NOTE) + lwet_convection = .false. + do_bm = .false. + do_ras = .false. + do_dryadj = .true. else if(uppercase(trim(convection_scheme)) == 'UNSET') then call error_mesg('idealized_moist_phys','determining convection scheme from flags', NOTE) @@ -411,7 +430,11 @@ subroutine idealized_moist_phys_init(Time, Time_step_in, nhum, rad_lon_2d, rad_l if (do_bm) then r_conv_scheme = FULL_BETTS_MILLER_CONV call error_mesg('idealized_moist_phys','Using Betts-Miller convection scheme.', NOTE) - end if + end if + if (do_dryadj) then + r_conv_scheme = DRYADJ_CONV + call error_mesg('idealized_moist_phys','Using AM3 dry adjustment convection scheme.',NOTE) + end if if (do_ras) then r_conv_scheme = RAS_CONV call error_mesg('idealized_moist_phys','Using relaxed Arakawa Schubert convection scheme.', NOTE) @@ -424,11 +447,20 @@ subroutine idealized_moist_phys_init(Time, Time_step_in, nhum, rad_lon_2d, rad_l if(lwet_convection .and. do_bm) & call error_mesg('idealized_moist_phys','lwet_convection and do_bm cannot both be .true.',FATAL) +if(do_dryadj .and. do_bm) & + call error_mesg('idealized_moist_phys','do_dryadj and do_bm cannot both be .true.',FATAL) + if(lwet_convection .and. do_ras) & - call error_mesg('idealized_moist_phys','lwet_convection and do_ras cannot both be .true.',FATAL) + call error_mesg('idealized_moist_phys','lwet_convection and do_ras cannot both be .true.',FATAL) + +if(lwet_convection .and. do_dryadj) & + call error_mesg('idealized_moist_phys','lwet_convection and do_dryadj cannot both be .true.',FATAL) if(do_bm .and. do_ras) & call error_mesg('idealized_moist_phys','do_bm and do_ras cannot both be .true.',FATAL) + +if(do_dryadj .and. do_ras) & + call error_mesg('idealized_moist_phys', 'do_dryadj and do_ras cannot both be .true.',FATAL) nsphum = nhum Time_step = Time_step_in @@ -452,6 +484,8 @@ subroutine idealized_moist_phys_init(Time, Time_step_in, nhum, rad_lon_2d, rad_l ! NTL 01/23 thermodynamic sea ice allocate(h_thermo_ice(is:ie, js:je)) allocate(t_ml (is:ie, js:je)) +allocate(const_correct(is:ie, js:je)) +allocate(nudge_out(is:ie, js:je)) ! allocate(q_surf (is:ie, js:je)); q_surf = 0.0 allocate(u_surf (is:ie, js:je)); u_surf = 0.0 @@ -617,7 +651,7 @@ subroutine idealized_moist_phys_init(Time, Time_step_in, nhum, rad_lon_2d, rad_l t_surf = t_surf_init + 1.0 call mixed_layer_init(is, ie, js, je, num_levels, t_surf, & - h_thermo_ice, t_ml, & ! NTL 01/23 thermodynamic sea ice + h_thermo_ice, t_ml, const_correct,nudge_out, & ! NTL 01/23 thermodynamic sea ice bucket_depth, get_axis_id(), Time, albedo, rad_lonb_2d(:,:), rad_latb_2d(:,:), land, bucket) ! t_surf is intent(inout) !s albedo distribution set here. elseif(gp_surface) then @@ -731,6 +765,9 @@ subroutine idealized_moist_phys_init(Time, Time_step_in, nhum, rad_lon_2d, rad_l call ras_init (do_strat, axes,Time,tracers_in_ras) +case (DRYADJ_CONV) + call dry_adj_init(alpha) + end select !jp not sure why these diag_fields are fenced when condensation ones above are not... @@ -791,6 +828,8 @@ subroutine idealized_moist_phys_init(Time, Time_step_in, nhum, rad_lon_2d, rad_l id_rh = register_diag_field ( mod_name, 'rh', & axes(1:3), Time, 'relative humidity', 'percent') + + end subroutine idealized_moist_phys_init !================================================================================================================================= subroutine idealized_moist_phys(Time, p_half, p_full, z_half, z_full, ug, vg, tg, grid_tracers, & @@ -887,6 +926,17 @@ subroutine idealized_moist_phys(Time, p_half, p_full, z_half, z_full, ug, vg, tg if(id_cape > 0) used = send_data(id_cape, cape, Time) if(id_cin > 0) used = send_data(id_cin, cin, Time) +case(DRYADJ_CONV) + call dry_adj(tg(:, :, :, previous), & + p_full(:,:,:,previous), p_half(:,:,:,previous), & + conv_dt_tg, delta_t) + + tg_tmp = conv_dt_tg(:,:,:) + tg(:,:,:,previous) + conv_dt_tg = conv_dt_tg / delta_t + conv_dt_qg = 0.0 + qg_tmp = grid_tracers(:,:,:,previous,nsphum) + if(id_conv_dt_tg > 0) used = send_data(id_conv_dt_tg, conv_dt_tg, Time) + case(DRY_CONV) call dry_convection(Time, tg(:, :, :, previous), & p_full(:,:,:,previous), p_half(:,:,:,previous), & @@ -1006,7 +1056,7 @@ subroutine idealized_moist_phys(Time, p_half, p_full, z_half, z_full, ug, vg, tg tg(:,:,:,previous), & net_surf_sw_down(:,:), & surf_lw_down(:,:), albedo, & - grid_tracers(:,:,:,previous,nsphum)) + grid_tracers(:,:,:,previous,nsphum), const_correct(:,:), nudge_out(:,:)) !NTL_CHANGE end if if(.not.mixed_layer_bc) then @@ -1255,7 +1305,7 @@ subroutine idealized_moist_phys(Time, p_half, p_full, z_half, z_full, ug, vg, tg js, & je, & t_surf(:,:), & ! t_surf is intent(inout) - h_thermo_ice(:,:), t_ml(:,:), & ! NTL 01/23 thermodynamic sea ice + h_thermo_ice(:,:), t_ml(:,:), const_correct(:,:), nudge_out(:,:), & ! NTL 01/23 thermodynamic sea ice flux_t(:,:), & flux_q(:,:), & flux_r(:,:), & @@ -1282,8 +1332,11 @@ subroutine idealized_moist_phys(Time, p_half, p_full, z_half, z_full, ug, vg, tg endif ! if(turb) then !s Adding relative humidity calculation so as to allow comparison with Frierson's thesis. - call rh_calc (p_full(:,:,:,previous),tg_tmp,qg_tmp,RH) - if(id_rh >0) used = send_data(id_rh, RH*100., Time) + + if(id_rh >0) then + call rh_calc (p_full(:,:,:,previous),tg_tmp,qg_tmp,RH) + used = send_data(id_rh, RH*100., Time) + endif ! RG Add bucket @@ -1349,6 +1402,7 @@ subroutine idealized_moist_phys_end if(two_stream_gray) call two_stream_gray_rad_end if(lwet_convection) call qe_moist_convection_end if(do_ras) call ras_end +if(do_dryadj) call dry_adj_end if(turb) then call vert_diff_end @@ -1356,7 +1410,7 @@ subroutine idealized_moist_phys_end endif call lscale_cond_end if(mixed_layer_bc) call mixed_layer_end(t_surf, & - h_thermo_ice, t_ml, albedo, & ! NTL 01/23 thermodynamic sea ice + h_thermo_ice, t_ml, const_correct, nudge_out, albedo, & ! NTL 01/23 thermodynamic sea ice bucket_depth, bucket) if(do_damping) call damping_driver_end From f25a3f123ce4fcb4dfc42ad6a9af8b72a4e5b24b Mon Sep 17 00:00:00 2001 From: ntlewis Date: Fri, 21 Jul 2023 10:11:50 +0100 Subject: [PATCH 07/19] change valid range for wind for exoplanets --- src/atmos_spectral/model/spectral_dynamics.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/atmos_spectral/model/spectral_dynamics.F90 b/src/atmos_spectral/model/spectral_dynamics.F90 index ab2e0b9ae..e7983c477 100644 --- a/src/atmos_spectral/model/spectral_dynamics.F90 +++ b/src/atmos_spectral/model/spectral_dynamics.F90 @@ -1568,7 +1568,7 @@ subroutine spectral_diagnostics_init(Time) real,dimension(2) :: vrange character(len=128) :: tname, longname, units -vrange = (/ -400., 400. /) +vrange = (/ -10000., 10000. /) rad_to_deg = 180./pi call get_grid_boundaries(lonb,latb,global=.true.) From 990e90735a6070abb7293d857f79199aeecbf955 Mon Sep 17 00:00:00 2001 From: ntlewis Date: Fri, 21 Jul 2023 10:12:34 +0100 Subject: [PATCH 08/19] add path to new dry convection scheme --- src/extra/model/grey/path_names | 1 + src/extra/model/isca/path_names | 7 ++++--- src/extra/model/socrates/path_names | 1 + 3 files changed, 6 insertions(+), 3 deletions(-) diff --git a/src/extra/model/grey/path_names b/src/extra/model/grey/path_names index 51c284f52..6d629a670 100644 --- a/src/extra/model/grey/path_names +++ b/src/extra/model/grey/path_names @@ -10,6 +10,7 @@ atmos_param/sea_esf_rad/null/rad_utilities.F90 atmos_param/shallow_conv/shallow_conv.F90 atmos_param/stable_bl_turb/stable_bl_turb.F90 atmos_param/strat_cloud/null/strat_cloud.F90 +atmos_param/dry_adj/dry_adj.F90 atmos_param/cloud_simple/cloud_simple.F90 atmos_param/two_stream_gray_rad/two_stream_gray_rad.F90 atmos_param/qflux/qflux.f90 diff --git a/src/extra/model/isca/path_names b/src/extra/model/isca/path_names index eebdfd2b0..287cf0c0c 100644 --- a/src/extra/model/isca/path_names +++ b/src/extra/model/isca/path_names @@ -11,6 +11,7 @@ atmos_param/sea_esf_rad/null/rad_utilities.F90 atmos_param/shallow_conv/shallow_conv.F90 atmos_param/stable_bl_turb/stable_bl_turb.F90 atmos_param/strat_cloud/null/strat_cloud.F90 +atmos_param/dry_adj/dry_adj.F90 atmos_param/cloud_simple/cloud_simple.F90 atmos_param/two_stream_gray_rad/two_stream_gray_rad.F90 atmos_param/qflux/qflux.f90 @@ -19,7 +20,7 @@ atmos_param/monin_obukhov/monin_obukhov_kernel.F90 atmos_param/monin_obukhov/monin_obukhov.F90 atmos_param/dry_convection/dry_convection.f90 atmos_param/rayleigh_bottom_drag/rayleigh_bottom_drag.F90 -atmos_param/rrtm_radiation/rrtmg_lw/gcm_model/modules/parkind.f90 +atmos_param/rrtm_radiation/rrtmg_lw/gcm_model/modules/parkind.F90 atmos_param/rrtm_radiation/rrtmg_lw/gcm_model/modules/parrrtm.f90 atmos_param/rrtm_radiation/rrtmg_lw/gcm_model/modules/rrlw_cld.f90 atmos_param/rrtm_radiation/rrtmg_lw/gcm_model/modules/rrlw_con.f90 @@ -49,7 +50,7 @@ atmos_param/rrtm_radiation/rrtmg_lw/gcm_model/src/mcica_random_numbers.f90 atmos_param/rrtm_radiation/rrtmg_lw/gcm_model/src/rrtmg_lw_cldprop.f90 atmos_param/rrtm_radiation/rrtmg_lw/gcm_model/src/rrtmg_lw_rad.nomcica.f90 atmos_param/rrtm_radiation/rrtmg_lw/gcm_model/src/rrtmg_lw_rtrn.f90 -atmos_param/rrtm_radiation/rrtmg_lw/gcm_model/src/rrtmg_lw_setcoef.f90 +atmos_param/rrtm_radiation/rrtmg_lw/gcm_model/src/rrtmg_lw_setcoef.F90 atmos_param/rrtm_radiation/rrtmg_lw/gcm_model/src/mcica_subcol_gen_lw.f90 atmos_param/rrtm_radiation/rrtmg_lw/gcm_model/src/rrtmg_lw_init.f90 atmos_param/rrtm_radiation/rrtmg_lw/gcm_model/src/rrtmg_lw_rtrnmc.f90 @@ -57,7 +58,7 @@ atmos_param/rrtm_radiation/rrtmg_lw/gcm_model/src/rrtmg_lw_taumol.f90 atmos_param/rrtm_radiation/rrtmg_lw/gcm_model/src/rrtmg_lw_cldprmc.f90 atmos_param/rrtm_radiation/rrtmg_lw/gcm_model/src/rrtmg_lw_k_g.f90 atmos_param/rrtm_radiation/rrtmg_lw/gcm_model/src/rrtmg_lw_rtrnmr.f90 -atmos_param/rrtm_radiation/rrtmg_sw/gcm_model/modules/parkind.f90 +atmos_param/rrtm_radiation/rrtmg_sw/gcm_model/modules/parkind.F90 atmos_param/rrtm_radiation/rrtmg_sw/gcm_model/modules/parrrsw.f90 atmos_param/rrtm_radiation/rrtmg_sw/gcm_model/modules/rrsw_aer.f90 atmos_param/rrtm_radiation/rrtmg_sw/gcm_model/modules/rrsw_cld.f90 diff --git a/src/extra/model/socrates/path_names b/src/extra/model/socrates/path_names index c92ad053f..120c82696 100644 --- a/src/extra/model/socrates/path_names +++ b/src/extra/model/socrates/path_names @@ -11,6 +11,7 @@ atmos_param/sea_esf_rad/null/rad_utilities.F90 atmos_param/shallow_conv/shallow_conv.F90 atmos_param/stable_bl_turb/stable_bl_turb.F90 atmos_param/strat_cloud/null/strat_cloud.F90 +atmos_param/dry_adj/dry_adj.F90 atmos_param/cloud_simple/cloud_simple.F90 atmos_param/two_stream_gray_rad/two_stream_gray_rad.F90 atmos_param/qflux/qflux.f90 From abc824166031389e5e93e76167a5127419864c94 Mon Sep 17 00:00:00 2001 From: ntlewis Date: Fri, 21 Jul 2023 10:12:56 +0100 Subject: [PATCH 09/19] add nudging options to sea ice code --- .../driver/solo/mixed_layer.F90 | 413 ++++++++++++++++-- 1 file changed, 385 insertions(+), 28 deletions(-) diff --git a/src/atmos_spectral/driver/solo/mixed_layer.F90 b/src/atmos_spectral/driver/solo/mixed_layer.F90 index ca424c8d2..432d032db 100644 --- a/src/atmos_spectral/driver/solo/mixed_layer.F90 +++ b/src/atmos_spectral/driver/solo/mixed_layer.F90 @@ -44,7 +44,7 @@ module mixed_layer_mod use time_manager_mod, only: time_type -use transforms_mod, only: get_deg_lat, get_deg_lon, grid_domain +use transforms_mod, only: get_deg_lat, get_deg_lon, grid_domain, area_weighted_global_mean use vert_diff_mod, only: surf_diff_type @@ -138,10 +138,21 @@ module mixed_layer_mod real, parameter :: thermo_ice_basal_flux_const = 120 ! linear coefficient for heat flux from ocean to ice base (W/m^2/K) real :: t_thermo_ice_base = TFREEZE ! temperature at base of ice, taken to be freshwater freezing point real :: t_surf_freeze = TFREEZE - - - - +logical :: do_var_thermo_ice_albedo = .false. +real :: thermo_ice_nudge_timescale = 7. * 86400. +logical :: do_nudge_ice = .false. +logical :: do_nudging_correction = .false. +logical :: correct_nh_only = .false. +logical :: diff_correction = .false. +character(len=256) :: nudge_ice_file_name = 'nudge_ice' +logical :: time_varying_nudge_ice = .false. +logical :: do_prescribe_albedo = .false. +character(len=256) :: albedo_file_name = 'albedo' +logical :: time_varying_albedo = .false. +logical :: prescribe_nh_albedo_only = .false. +logical :: read_const_correct = .true. +logical :: read_nudge_out = .true. +logical :: albedo_from_nudge_file = .false. namelist/mixed_layer_nml/ evaporation, depth, qflux_amp, qflux_width, tconst,& delta_T, prescribe_initial_dist,albedo_value, & @@ -161,7 +172,11 @@ module mixed_layer_mod ice_concentration_threshold, & add_latent_heat_flux_anom,flux_lhe_anom_file_name,& flux_lhe_anom_field_name, do_ape_sst, qflux_field_name,& - do_explicit, do_thermo_ice, thermo_ice_albedo, t_surf_freeze ! ICE, NTL 01/23 + do_explicit, do_thermo_ice, thermo_ice_albedo, t_surf_freeze, t_thermo_ice_base, & + do_var_thermo_ice_albedo,& ! ICE, NTL 01/23 + do_nudge_ice, nudge_ice_file_name, thermo_ice_nudge_timescale, time_varying_nudge_ice, & + do_nudging_correction, correct_nh_only, do_prescribe_albedo, time_varying_albedo, albedo_file_name, & + prescribe_nh_albedo_only, diff_correction, read_const_correct, read_nudge_out, albedo_from_nudge_file !================================================================================================================================= @@ -183,7 +198,12 @@ module mixed_layer_mod ! NTL 01/23 thermodynamic sea ice id_h_thermo_ice, & ! sea ice thickness id_t_ml, & ! mixed layer temperature - id_flux_thermo_ice ! conductive heat flux through ice + id_flux_thermo_ice, & ! conductive heat flux through ice + id_ice_nudge_flux, & + id_ice_nudge_correction, & + id_ice_nudge_flux_avg, & + id_ice_nudge_correction_avg, & + id_ice_nudge_correction_chk !NTL_DELETE real, allocatable, dimension(:,:) :: & ocean_qflux, & ! Q-flux @@ -221,9 +241,16 @@ module mixed_layer_mod delta_t_thermo_ice, & ! Increment in ice (steady-state) surface temperature flux_thermo_ice, & ! Conductive heat flux through ice corrected_flux_noq +real, allocatable, dimension(:,:) :: h_thermo_ice_in, nudge_flux, nudge_correction, albedo_in +real, allocatable, dimension(:,:) :: delta_t_surf_correct, delta_t_surf_nocorrect !NTL_DELETE +real :: const_nudge_correction +type(interpolate_type),save :: nudge_ice_interp +type(interpolate_type),save :: albedo_interp logical, allocatable, dimension(:,:) :: land_mask + + !mj read sst from input file type(interpolate_type),save :: sst_interp type(interpolate_type),save :: qflux_interp @@ -237,12 +264,12 @@ module mixed_layer_mod !================================================================================================================================= subroutine mixed_layer_init(is, ie, js, je, num_levels, t_surf, & - h_thermo_ice, t_ml, & ! NTL 01/23 thermodynamic sea ice + h_thermo_ice, t_ml, const_correct, nudge_out, & ! NTL 01/23 thermodynamic sea ice bucket_depth, axes, Time, albedo, rad_lonb_2d,rad_latb_2d, land, restart_file_bucket_depth) type(time_type), intent(in) :: Time real, intent(out), dimension(:,:) :: t_surf, albedo -real, intent(out), dimension(:,:) :: h_thermo_ice, t_ml ! NTL 01/23 thermodynamic sea ice +real, intent(out), dimension(:,:) :: h_thermo_ice, t_ml, const_correct, nudge_out ! NTL 01/23 thermodynamic sea ice real, intent(out), dimension(:,:,:) :: bucket_depth integer, intent(in), dimension(4) :: axes real, intent(in), dimension(:,:) :: rad_lonb_2d, rad_latb_2d @@ -320,7 +347,12 @@ subroutine mixed_layer_init(is, ie, js, je, num_levels, t_surf, & allocate(delta_h_thermo_ice (is:ie, js:je)) allocate(delta_t_thermo_ice (is:ie, js:je)) allocate(flux_thermo_ice (is:ie, js:je)) - +allocate(h_thermo_ice_in (is:ie, js:je)) +allocate(albedo_in (is:ie, js:je)) +allocate(nudge_flux (is:ie, js:je)) +allocate(nudge_correction (is:ie, js:je)) +allocate(delta_t_surf_nocorrect (is:ie, js:je)) !NTL_DELETE +allocate(delta_t_surf_correct (is:ie, js:je)) !NTL_DELETE ! !see if restart file exists for the surface temperature ! @@ -365,6 +397,18 @@ subroutine mixed_layer_init(is, ie, js, je, num_levels, t_surf, & if (do_thermo_ice) then call read_data(trim('INPUT/mixed_layer.res'), 'h_thermo_ice', h_thermo_ice, grid_domain) call read_data(trim('INPUT/mixed_layer.res'), 't_ml', t_ml, grid_domain) + if (read_const_correct) then + call read_data(trim('INPUT/mixed_layer.res'), 'const_correct', const_correct, grid_domain) + else + const_correct(:,:) = 0.0 + endif + if (read_nudge_out) then + call read_data(trim('INPUT/mixed_layer.res'), 'nudge_out', nudge_out, grid_domain) + else + nudge_out(:,:) = 0.0 + endif + call read_data(trim('INPUT/mixed_layer.res'), 'albedo', albedo, grid_domain) + else if (do_prescribe_albedo) then call read_data(trim('INPUT/mixed_layer.res'), 'albedo', albedo, grid_domain) endif @@ -391,12 +435,18 @@ subroutine mixed_layer_init(is, ie, js, je, num_levels, t_surf, & if (do_thermo_ice) then ! NTL 01/23 thermodynamic sea ice h_thermo_ice(:,:) = 0.0 t_ml(:,:) = t_surf(:,:) + const_correct(:,:) = 0.0 + nudge_out(:,:) = 0.0 albedo(:,:) = albedo_value where(land) albedo(:,:) = land_albedo_prefactor * albedo_value albedo_initial(:,:) = albedo (:,:) !where(.not. land) ! where (t_surf .lt. TFREEZE) albedo(:,:) = thermo_ice_albedo !endwhere + else if (do_prescribe_albedo) then + albedo(:,:) = albedo_value + where(land) albedo(:,:) = land_albedo_prefactor * albedo_value + albedo_initial(:,:) = albedo (:,:) endif else @@ -418,7 +468,7 @@ subroutine mixed_layer_init(is, ie, js, je, num_levels, t_surf, & axes(1:2), 'mixed layer heat capacity','joules/m^2/deg C') id_delta_t_surf = register_diag_field(mod_name, 'delta_t_surf', & axes(1:2), Time, 'change in sst','K') -if (update_albedo_from_ice .or. do_thermo_ice) then +if (update_albedo_from_ice .or. do_thermo_ice .or. do_prescribe_albedo) then id_albedo = register_diag_field(mod_name, 'albedo', & axes(1:2), Time, 'surface albedo', 'none') if (do_thermo_ice) then @@ -429,6 +479,16 @@ subroutine mixed_layer_init(is, ie, js, je, num_levels, t_surf, & axes(1:2), Time, 'mixed layer tempeature','K') id_flux_thermo_ice = register_diag_field(mod_name, 'flux_thermo_ice', & axes(1:2), Time, 'conductive heat flux through sea ice','watts/m2') + id_ice_nudge_flux = register_diag_field(mod_name, 'ice_nudge_flux', & + axes(1:2), Time, 'nudging heat flux','watts/m2') + id_ice_nudge_correction = register_diag_field(mod_name, 'ice_nudge_correction', & + axes(1:2), Time, 'nudging heat flux correction','watts/m2') + id_ice_nudge_correction_chk = register_diag_field(mod_name, 'ice_nudge_correction_chk', & + axes(1:2), Time, 'nudging heat flux correction check','watts/m2') !NTL_DELETE + id_ice_nudge_flux_avg = register_diag_field(mod_name, 'ice_nudge_flux_avg', & + Time, 'nudging heat flux','watts/m2') + id_ice_nudge_correction_avg = register_diag_field(mod_name, 'ice_nudge_correction_avg', & + Time, 'nudging heat flux correction','watts/m2') else id_ice_conc = register_diag_field(mod_name, 'ice_conc', & axes(1:2), Time, 'ice_concentration', 'none') @@ -468,6 +528,63 @@ subroutine mixed_layer_init(is, ie, js, je, num_levels, t_surf, & endif endif + +h_thermo_ice_in = 999. ! large number stops this doing anything if do_nudge_ice is false +! load ice +if (do_nudge_ice .or. albedo_from_nudge_file) then + if (time_varying_nudge_ice .or. albedo_from_nudge_file) then + call interpolator_init( nudge_ice_interp, trim(nudge_ice_file_name)//'.nc', rad_lonb_2d, rad_latb_2d, data_out_of_bounds=(/CONSTANT/) ) + else + + if(file_exist(trim('INPUT/'//nudge_ice_file_name//'.nc'))) then + call mpp_get_global_domain(grid_domain, xsize=global_num_lon, ysize=global_num_lat) + call field_size(trim('INPUT/'//nudge_ice_file_name//'.nc'), trim(nudge_ice_file_name), siz) + if ( siz(1) == global_num_lon .or. siz(2) == global_num_lat ) then + call read_data(trim('INPUT/'//nudge_ice_file_name//'.nc'), trim(nudge_ice_file_name), h_thermo_ice_in, grid_domain) + else + write(ctmp1(1: 4),'(i4)') siz(1) + write(ctmp1(9:12),'(i4)') siz(2) + write(ctmp2(1: 4),'(i4)') global_num_lon + write(ctmp2(9:12),'(i4)') global_num_lat + call error_mesg ('get_ice','nudge_ice file contains data on a '// & + ctmp1//' grid, but atmos model grid is '//ctmp2, FATAL) + endif + else + call error_mesg('get_ice','do_nudge_ice="'//trim('True')//'"'// & + ' but '//trim(nudge_ice_file_name)//' does not exist', FATAL) + endif + + endif +endif + +albedo_in = 0.0 +if (do_prescribe_albedo) then + if (time_varying_albedo) then + call interpolator_init( albedo_interp, trim(albedo_file_name)//'.nc', rad_lonb_2d, rad_latb_2d, data_out_of_bounds=(/CONSTANT/) ) + else + + if(file_exist(trim('INPUT/'//albedo_file_name//'.nc'))) then + call mpp_get_global_domain(grid_domain, xsize=global_num_lon, ysize=global_num_lat) + call field_size(trim('INPUT/'//albedo_file_name//'.nc'), trim(albedo_file_name), siz) + if ( siz(1) == global_num_lon .or. siz(2) == global_num_lat ) then + call read_data(trim('INPUT/'//albedo_file_name//'.nc'), trim(albedo_file_name), albedo_in, grid_domain) + else + write(ctmp1(1: 4),'(i4)') siz(1) + write(ctmp1(9:12),'(i4)') siz(2) + write(ctmp2(1: 4),'(i4)') global_num_lon + write(ctmp2(9:12),'(i4)') global_num_lat + call error_mesg ('get_albedo','albedo file contains data on a '// & + ctmp1//' grid, but atmos model grid is '//ctmp2, FATAL) + endif + else + call error_mesg('get_albedo','do_prescribe_albedo="'//trim('True')//'"'// & + ' but '//trim(albedo_file_name)//' does not exist', FATAL) + endif + + endif +endif + + !s Adding MiMA options for qfluxes. if ( do_qflux .or. do_warmpool) then @@ -490,7 +607,7 @@ subroutine mixed_layer_init(is, ie, js, je, num_levels, t_surf, & !s Prescribe albedo distribution here so that it will be the same in both two_stream_gray later and rrtmg radiation. -if (.not. do_thermo_ice) then ! NTL 01/23 thermodynamic sea ice -- this is set / loaded earlier +if ((.not. do_thermo_ice) .and. (.not. do_prescribe_albedo)) then ! NTL 01/23 thermodynamic sea ice -- this is set / loaded earlier albedo(:,:) = albedo_value @@ -638,6 +755,8 @@ subroutine mixed_layer ( & ! NTL 01/23 thermodynamic sea ice h_thermo_ice, & t_ml, & + const_correct, & + nudge_out, & ! end ice addition flux_t, & flux_q, & @@ -652,7 +771,7 @@ subroutine mixed_layer ( & drdt_surf, & dhdt_atm, & dedq_atm, & - albedo_out) + albedo_out) ! ---- arguments ----------------------------------------------------------- type(time_type), intent(in) :: Time, Time_next @@ -660,7 +779,7 @@ subroutine mixed_layer ( & real, intent(in), dimension(:,:) :: net_surf_sw_down, surf_lw_down real, intent(in), dimension(:,:) :: flux_t, flux_q, flux_r real, intent(inout), dimension(:,:) :: t_surf -real, intent(inout), dimension(:,:) :: h_thermo_ice, t_ml ! NTL 01/23 thermodynamic sea ice +real, intent(inout), dimension(:,:) :: h_thermo_ice, t_ml, const_correct, nudge_out ! NTL 01/23 thermodynamic sea ice real, intent(in), dimension(:,:) :: dhdt_surf, dedt_surf, dedq_surf, & drdt_surf, dhdt_atm, dedq_atm real, intent(in) :: dt @@ -737,12 +856,79 @@ subroutine mixed_layer ( & endif +if((do_nudge_ice.and.time_varying_nudge_ice).or. albedo_from_nudge_file) then + call interpolator( nudge_ice_interp, Time_next, h_thermo_ice_in, trim(nudge_ice_file_name) ) +endif + +if(do_prescribe_albedo.and.time_varying_albedo) then + call interpolator( albedo_interp, Time, albedo_in, trim(albedo_file_name) ) +endif + + + + + + ! ! Implement mixed layer surface boundary condition ! corrected_flux = - net_surf_sw_down - surf_lw_down + alpha_t * CP_AIR + alpha_lw - ocean_qflux t_surf_dependence = beta_t * CP_AIR + beta_lw + +nudge_out = 0.0 +where ((h_thermo_ice .gt. 0) .and. (h_thermo_ice .le. 0.1)) +nudge_out = h_thermo_ice / (86400./4.) * L_thermo_ice +endwhere +corrected_flux = corrected_flux - nudge_out +const_nudge_correction = 0.0 +const_correct = 0.0 +const_nudge_correction = -area_weighted_global_mean(nudge_out) +if (.not. (const_nudge_correction .eq. 0.)) then + where (nudge_out .eq. 0.) + const_correct(:,:) = const_nudge_correction + endwhere + const_correct = const_correct * const_nudge_correction / area_weighted_global_mean(const_correct) + !const_correct(:,:) = -area_weighted_global_mean(nudge_out) !-area_weighted_global_mean(nudge_flux * L_thermo_ice/dt) + corrected_flux = corrected_flux - const_correct +endif + + + + +nudge_out = 0.0 +if (do_nudge_ice) then + where ((h_thermo_ice .gt. 0.1) .and. (h_thermo_ice_in .le. 0)) + nudge_out = h_thermo_ice / thermo_ice_nudge_timescale * L_thermo_ice + endwhere + corrected_flux = corrected_flux - nudge_out +endif + + + +! if (prescribe_nh_albedo_only) then +! where (rad_lat_2d .gt. 0) + +const_nudge_correction = 0.0 +const_correct = 0.0 +if (do_nudging_correction) then + const_nudge_correction = -area_weighted_global_mean(nudge_out) + if (.not. (const_nudge_correction .eq. 0.)) then + if (correct_nh_only) then + where ((nudge_out .eq. 0.) .and. (rad_lat_2d .gt. 0)) + const_correct(:,:) = const_nudge_correction + endwhere + else + where (nudge_out .eq. 0.) + const_correct(:,:) = const_nudge_correction + endwhere + endif + const_correct = const_correct * const_nudge_correction / area_weighted_global_mean(const_correct) + !const_correct(:,:) = -area_weighted_global_mean(nudge_out) !-area_weighted_global_mean(nudge_flux * L_thermo_ice/dt) + corrected_flux = corrected_flux - const_correct + endif +endif + if (evaporation) then corrected_flux = corrected_flux + alpha_q * HLV t_surf_dependence = t_surf_dependence + beta_q * HLV @@ -752,8 +938,13 @@ subroutine mixed_layer ( & ! NTL 01/23 thermodynamic sea ice corrected_flux_noq = corrected_flux + ocean_qflux flux_thermo_ice = 0 ! Initialize conductive heat flux through sea ice + nudge_flux = 0. ! Initialize nudging flux + nudge_correction = 0. ! Initializing nudging correction + const_nudge_correction = 0. ! for calculation of ice sfc temperature dFdt_surf = dhdt_surf + drdt_surf + HLV * (dedt_surf) ! d(corrected_flux)/d(T_surf) + delta_t_surf_correct = 0. !NTL_DELETE + delta_t_surf_nocorrect = 0. !NTL_DELETE endif !s Surface heat_capacity calculation based on that in MiMA by mj @@ -825,7 +1016,7 @@ subroutine mixed_layer ( & endif !s end of if(do_sc_sst). -else ! if (do_thermo_ice) +else if (do_thermo_ice) then where (land_ice_mask) ! name misleading, for thermo sea ice, land_ice_mask is just land mask delta_t_ml = - corrected_flux * dt / land_sea_heat_capacity delta_h_thermo_ice = 0 @@ -834,7 +1025,7 @@ subroutine mixed_layer ( & ! = update state = t_ml = t_ml + delta_t_ml h_thermo_ice = h_thermo_ice + delta_h_thermo_ice - albedo_out(:,:) = albedo_value * land_albedo_prefactor + !albedo_out(:,:) = albedo_value * land_albedo_prefactor elsewhere ! NTL 01/23 thermodynamic sea ice @@ -863,35 +1054,56 @@ subroutine mixed_layer ( & ! growth). ! calculate values of delta_h_ice, delta_t_ml, and delta_t_ice - + + + + where ( h_thermo_ice .le. 0 ) ! = ice-free = [ should be equivalent to use .eq. instead of .le. ] - delta_t_ml = - ( corrected_flux - ocean_qflux) * dt / land_sea_heat_capacity + delta_t_ml = - ( corrected_flux_noq - ocean_qflux) * dt / land_sea_heat_capacity delta_h_thermo_ice = 0 elsewhere ! = ice-covered = ( h>0 ) delta_t_ml = - ( thermo_ice_basal_flux_const * ( t_ml - t_thermo_ice_base ) - ocean_qflux ) * dt / land_sea_heat_capacity - delta_h_thermo_ice = ( corrected_flux - thermo_ice_basal_flux_const * ( t_ml - t_thermo_ice_base ) ) * dt / L_thermo_ice + delta_h_thermo_ice = ( corrected_flux_noq - thermo_ice_basal_flux_const * ( t_ml - t_thermo_ice_base ) ) * dt / L_thermo_ice + !where (h_thermo_ice_in .le. 0 ) ! by default h_thermo_ice_in is zero and this doesn't do anything + ! nudge_flux = - h_thermo_ice_in / thermo_ice_nudge_timescale * dt + ! delta_h_thermo_ice = delta_h_thermo_ice + nudge_flux + !endwhere endwhere + where ( t_ml + delta_t_ml .lt. t_surf_freeze ) ! = frazil growth = delta_h_thermo_ice = delta_h_thermo_ice - ( t_ml + delta_t_ml - t_surf_freeze ) * (land_sea_heat_capacity)/L_thermo_ice delta_t_ml = t_surf_freeze - t_ml endwhere + where ( ( h_thermo_ice .gt. 0 ) .and. ( h_thermo_ice + delta_h_thermo_ice .le. 0 ) ) ! = complete ablation = delta_t_ml = delta_t_ml - ( h_thermo_ice + delta_h_thermo_ice ) * L_thermo_ice/land_sea_heat_capacity delta_h_thermo_ice = - h_thermo_ice + + ! elsewhere ( ( h_thermo_ice .gt. 0 ) .and. ( h_thermo_ice + delta_h_thermo_ice - h_thermo_ice / thermo_ice_nudge_timescale * dt .le. 0 ) .and. (h_thermo_ice_in .le. 0 )) ! by default h_thermo_ice_in is 999 and this doesn't do anything + ! nudge_flux = (h_thermo_ice + delta_h_thermo_ice) + ! ! no need to update delta_t_ml as dT_ml = dT_ml - (h + dh + n_f), but this = dT_ml + ! delta_h_thermo_ice = - h_thermo_ice + ! elsewhere ( (h_thermo_ice .gt. 0) .and. ( h_thermo_ice + delta_h_thermo_ice - h_thermo_ice / thermo_ice_nudge_timescale * dt .gt. 0 ) .and. (h_thermo_ice_in .le. 0 )) ! by default h_thermo_ice_in is 999 and this doesn't do anything + ! nudge_flux = h_thermo_ice / thermo_ice_nudge_timescale * dt + ! delta_h_thermo_ice = delta_h_thermo_ice - nudge_flux + endwhere + ! = update surface temperature = ! compute surface melt where ( h_thermo_ice + delta_h_thermo_ice .gt. 0 ) ! surface is ice-covered ! calculate increment in steady-state ice surface temperature flux_thermo_ice = k_thermo_ice / (h_thermo_ice + delta_h_thermo_ice) * ( t_thermo_ice_base - t_surf ) - delta_t_thermo_ice = ( - corrected_flux + flux_thermo_ice ) & + delta_t_thermo_ice = ( - corrected_flux_noq + flux_thermo_ice ) & / ( k_thermo_ice / (h_thermo_ice + delta_h_thermo_ice) + dFdt_surf ) where ( t_surf + delta_t_thermo_ice .gt. t_surf_freeze ) ! surface ablation delta_t_thermo_ice = t_surf_freeze - t_surf endwhere ! surface is ice-covered, so update t_surf as ice surface temperature + delta_t_surf = delta_t_thermo_ice t_surf = t_surf + delta_t_thermo_ice elsewhere ! ice-free, so update t_surf as mixed layer temperature + delta_t_surf = delta_t_ml t_surf = t_ml + delta_t_ml endwhere @@ -901,19 +1113,155 @@ subroutine mixed_layer ( & ! = update state = t_ml = t_ml + delta_t_ml h_thermo_ice = h_thermo_ice + delta_h_thermo_ice - - where ( h_thermo_ice .gt. 0 ) - albedo_out(:,:) = albedo_value + (thermo_ice_albedo - albedo_value) * tanh(h_thermo_ice/2.) - elsewhere - albedo_out(:,:) = albedo_value - endwhere + + endwhere + ! delta_t_surf_nocorrect = delta_t_surf !NTL_DELETE + + ! ! do nudging correction to ice free area + ! if (do_nudging_correction) then + ! if(.not. diff_correction) then + ! const_nudge_correction = -area_weighted_global_mean(nudge_flux * L_thermo_ice/land_sea_heat_capacity) + + ! if (.not. (const_nudge_correction .eq. 0.)) then + ! where ( h_thermo_ice .le. 0. ) + ! nudge_correction = const_nudge_correction + ! endwhere + + ! nudge_correction = nudge_correction * const_nudge_correction / area_weighted_global_mean(nudge_correction) + + ! where ( h_thermo_ice .le. 0. ) + ! delta_t_surf = delta_t_surf + nudge_correction + ! t_surf = t_surf + nudge_correction + ! delta_t_ml = delta_t_ml + nudge_correction !NTL_NEWCODE + ! t_ml = t_ml + nudge_correction !NTL_NEWCODE + ! endwhere + ! endif + ! else if (diff_correction) then + + ! ! do nudging correction to all locs where there is no nudging + ! !if (do_nudging_correction) then + ! const_nudge_correction = -area_weighted_global_mean(nudge_flux * L_thermo_ice/land_sea_heat_capacity) + + ! if (.not. (const_nudge_correction .eq. 0.)) then + ! where ( nudge_flux .eq. 0. ) + ! nudge_correction = const_nudge_correction + ! endwhere + + ! nudge_correction = nudge_correction * const_nudge_correction / area_weighted_global_mean(nudge_correction) + + ! where ( nudge_flux .eq. 0. ) + ! where ( h_thermo_ice .le. 0. ) + ! delta_t_surf = delta_t_surf + nudge_correction + ! t_surf = t_surf + nudge_correction + ! delta_t_ml = delta_t_ml + nudge_correction !NTL_NEWCODE + ! t_ml = t_ml + nudge_correction !NTL_NEWCODE + ! elsewhere (h_thermo_ice .gt. 0.) + ! h_thermo_ice = h_thermo_ice - delta_h_thermo_ice + ! t_surf = t_surf - delta_t_thermo_ice + ! delta_h_thermo_ice = delta_h_thermo_ice - nudge_correction*depth*RHO_CP/L_thermo_ice + ! flux_thermo_ice = k_thermo_ice / (h_thermo_ice + delta_h_thermo_ice) * ( t_thermo_ice_base - t_surf ) + ! delta_t_thermo_ice = ( - corrected_flux_noq + flux_thermo_ice ) & + ! / ( k_thermo_ice / (h_thermo_ice + delta_h_thermo_ice) + dFdt_surf ) + ! where ( t_surf + delta_t_thermo_ice .gt. t_surf_freeze ) ! surface ablation + ! delta_t_thermo_ice = t_surf_freeze - t_surf + ! endwhere + ! ! surface is ice-covered, so update t_surf as ice surface temperature + ! delta_t_surf = delta_t_thermo_ice + ! t_surf = t_surf + delta_t_thermo_ice + ! h_thermo_ice = h_thermo_ice + delta_h_thermo_ice + ! endwhere + ! endwhere + ! endif + ! endif + ! endif + ! ! end of nudging correction + + ! elsewhere ( ( h_thermo_ice .gt. 0 ) .and. ( h_thermo_ice + delta_h_thermo_ice - h_thermo_ice / thermo_ice_nudge_timescale * dt .le. 0 ) .and. (h_thermo_ice_in .le. 0 )) ! by default h_thermo_ice_in is 999 and this doesn't do anything + ! nudge_flux = (h_thermo_ice + delta_h_thermo_ice) + ! ! no need to update delta_t_ml as dT_ml = dT_ml - (h + dh + n_f), but this = dT_ml + ! delta_h_thermo_ice = - h_thermo_ice + ! elsewhere ( (h_thermo_ice .gt. 0) .and. ( h_thermo_ice + delta_h_thermo_ice - h_thermo_ice / thermo_ice_nudge_timescale * dt .gt. 0 ) .and. (h_thermo_ice_in .le. 0 )) ! by default h_thermo_ice_in is 999 and this doesn't do anything + ! nudge_flux = h_thermo_ice / thermo_ice_nudge_timescale * dt + ! delta_h_thermo_ice = delta_h_thermo_ice - nudge_flux + + + + ! if (do_nudge_ice) then + ! where ((h_thermo_ice .gt. 0) .and. (h_thermo_ice_in .le. 0)) + ! nudge_out = h_thermo_ice / thermo_ice_nudge_timescale * L_thermo_ice + ! elsewhere + ! nudge_out = 0.0 + ! endwhere + ! endif + + + ! if (do_nudging_correction) then + ! const_nudge_correction = -area_weighted_global_mean(nudge_out) + ! where (nudge_out .eq. 0.) + ! const_correct(:,:) = const_nudge_correction + ! elsewhere + ! const_correct(:,:) = 0.0 + ! endwhere + ! const_correct = const_correct * const_nudge_correction / area_weighted_global_mean(const_correct) + ! !const_correct(:,:) = -area_weighted_global_mean(nudge_out) !-area_weighted_global_mean(nudge_flux * L_thermo_ice/dt) + ! else + ! const_correct(:,:) = 0.0 + ! endif + + + !delta_t_surf_correct = delta_t_surf - delta_t_surf_nocorrect !NTL_DELETE + + + ! compute albedo + if (do_var_thermo_ice_albedo) then + where (land_ice_mask) + albedo_out(:,:) = albedo_value * land_albedo_prefactor + elsewhere + where ( h_thermo_ice .gt. 0.1 ) + albedo_out(:,:) = albedo_value + (thermo_ice_albedo - albedo_value) * tanh((h_thermo_ice-0.1)/1.) + elsewhere + albedo_out(:,:) = albedo_value + endwhere + endwhere + else + where (land_ice_mask) + albedo_out(:,:) = albedo_value * land_albedo_prefactor + elsewhere + where ( h_thermo_ice .gt. 0.1 ) + albedo_out(:,:) = thermo_ice_albedo + elsewhere + albedo_out(:,:) = albedo_value + endwhere + endwhere + endif + + + if (do_prescribe_albedo) then + if (prescribe_nh_albedo_only) then + where (rad_lat_2d .gt. 0) + albedo_out(:,:) = albedo_in(:,:) + endwhere + else + albedo_out(:,:) = albedo_in(:,:) + endif + endif if ( id_albedo > 0 ) used = send_data ( id_albedo, albedo_out, Time_next ) +endif + +if ((do_prescribe_albedo).and.(.not.do_thermo_ice)) then + albedo_out(:,:) = albedo_in(:,:) + if ( id_albedo > 0 ) used = send_data ( id_albedo, albedo_out, Time_next ) endif ! if (do_thermo_ice) +if (albedo_from_nudge_file) then + where ((h_thermo_ice .gt. 0) .and. (h_thermo_ice_in .le. 0)) + albedo_out(:,:) = albedo_value + endwhere +endif ! ! Finally calculate the increments for the lowest atmospheric layer @@ -940,6 +1288,11 @@ subroutine mixed_layer ( & if(id_h_thermo_ice > 0) used = send_data(id_h_thermo_ice, h_thermo_ice, Time_next) if(id_t_ml > 0) used = send_data(id_t_ml, t_ml, Time_next) if(id_flux_thermo_ice > 0) used = send_data(id_flux_thermo_ice, flux_thermo_ice, Time_next) +if(id_ice_nudge_flux > 0) used = send_data(id_ice_nudge_flux, L_thermo_ice * nudge_flux / dt, Time_next) +if(id_ice_nudge_correction > 0) used = send_data(id_ice_nudge_correction, land_sea_heat_capacity * nudge_correction / dt, Time_next) +if(id_ice_nudge_correction_chk > 0) used = send_data(id_ice_nudge_correction_chk, land_sea_heat_capacity * delta_t_surf_correct / dt, Time_next) ! NTL_DELETE +if(id_ice_nudge_flux_avg > 0) used = send_data(id_ice_nudge_flux_avg, area_weighted_global_mean(L_thermo_ice * nudge_flux / dt), Time_next) +if(id_ice_nudge_correction_avg > 0) used = send_data(id_ice_nudge_correction_avg, area_weighted_global_mean(land_sea_heat_capacity * nudge_correction / dt), Time_next) end subroutine mixed_layer @@ -978,11 +1331,11 @@ end subroutine read_ice_conc !================================================================================================================================= subroutine mixed_layer_end(t_surf, & - h_thermo_ice, t_ml, albedo, & ! NTL 01/23 thermodynamic sea ice + h_thermo_ice, t_ml, const_correct, nudge_out, albedo, & ! NTL 01/23 thermodynamic sea ice bucket_depth, restart_file_bucket_depth) real, intent(inout), dimension(:,:) :: t_surf -real, intent(inout), dimension(:,:) :: h_thermo_ice, t_ml, albedo ! NTL 01/23 thermodynamic sea ice +real, intent(inout), dimension(:,:) :: h_thermo_ice, t_ml, const_correct, nudge_out, albedo ! NTL 01/23 thermodynamic sea ice real, intent(inout), dimension(:,:,:) :: bucket_depth logical, intent(in) :: restart_file_bucket_depth integer:: unit @@ -995,6 +1348,10 @@ subroutine mixed_layer_end(t_surf, & if (do_thermo_ice) then ! NTL 01/23 thermodynamic sea ice call write_data(trim('RESTART/mixed_layer.res'), 'h_thermo_ice', h_thermo_ice, grid_domain) call write_data(trim('RESTART/mixed_layer.res'), 't_ml', t_ml, grid_domain) + call write_data(trim('RESTART/mixed_layer.res'), 'const_correct', const_correct, grid_domain) + call write_data(trim('RESTART/mixed_layer.res'), 'nudge_out', nudge_out, grid_domain) + call write_data(trim('RESTART/mixed_layer.res'), 'albedo', albedo, grid_domain) +else if (do_prescribe_albedo) then call write_data(trim('RESTART/mixed_layer.res'), 'albedo', albedo, grid_domain) endif if (restart_file_bucket_depth) then From 3cbfbdc75fb09f1d71fb972881df2382be160420 Mon Sep 17 00:00:00 2001 From: ntlewis Date: Wed, 18 Oct 2023 12:26:42 +0100 Subject: [PATCH 10/19] compiler options in python interface --- ci/environment-py3.9.yml | 2 +- src/extra/python/isca/codebase.py | 26 +++++++++++++---- src/extra/python/isca/experiment.py | 5 +++- .../templates/mkmf.template.ubuntu_conda2 | 28 +++++++++++++++++++ 4 files changed, 54 insertions(+), 7 deletions(-) create mode 100644 src/extra/python/isca/templates/mkmf.template.ubuntu_conda2 diff --git a/ci/environment-py3.9.yml b/ci/environment-py3.9.yml index 169e3adef..345e74483 100644 --- a/ci/environment-py3.9.yml +++ b/ci/environment-py3.9.yml @@ -1,4 +1,4 @@ -name: isca_env +name: isca_env2 channels: - conda-forge dependencies: diff --git a/src/extra/python/isca/codebase.py b/src/extra/python/isca/codebase.py index 6c8b27370..ada94b55b 100644 --- a/src/extra/python/isca/codebase.py +++ b/src/extra/python/isca/codebase.py @@ -34,7 +34,8 @@ def from_repo(cls, repo, commit=None, **kwargs): def from_directory(cls, directory, **kwargs): return cls(directory=directory, **kwargs) - def __init__(self, repo=None, commit=None, directory=None, storedir=P(GFDL_WORK, 'codebase'), safe_mode=False): + def __init__(self, repo=None, commit=None, directory=None, storedir=P(GFDL_WORK, 'codebase'), safe_mode=False, + compiler='intel', env_name=None): """Create a new CodeBase object. A CodeBase can be created with either a git repository or a file directory as it's source. @@ -59,9 +60,13 @@ def __init__(self, repo=None, commit=None, directory=None, storedir=P(GFDL_WORK, self.log.error('Too many sources. Cannot create a CodeBase with both a source directory and a source repository.') raise AttributeError('Either repo= or directory= required to create CodeBase.') - + self.compiler = compiler + self.env_name = env_name self.safe_mode = safe_mode self.storedir = storedir + if compiler == 'gfortran': + self.executable_name = self.executable_name.strip('.x')+'_gfort.x' + #print('HERE!!', self.executable_name) if directory is not None: self.directory = directory @@ -121,7 +126,12 @@ def __init__(self, repo=None, commit=None, directory=None, storedir=P(GFDL_WORK, self.path_names = [] self.compile_flags = [] # users can append to this to add additional compiler options self.precision_compile_flags = ['-DOVERLOAD_C8'] # default double precision - self.other_compile_flags = ['-r8', '-i4'] # default double precision + print(compiler, self.executable_name) + if self.compiler == 'intel': + self.other_compile_flags = ['-r8', '-i4'] # default double precision + elif self.compiler == 'gfortran': + self.other_compile_flags = ['-fdefault-real-8', '-fdefault-double-8'] + print('here', compiler, self.executable_name) @property def code_is_available(self): @@ -243,7 +253,10 @@ def _log_line(self, line): @useworkdir @destructive def compile(self, debug=False, optimisation=None): - env = get_env_file() + if self.env_name is not None: + env = get_env_file(self.env_name) + else: + env = get_env_file() mkdir(self.builddir) compile_flags = [] @@ -292,7 +305,10 @@ def compile(self, debug=False, optimisation=None): def use_single_precision(self): self.log.info('Using single precision') self.precision_compile_flags = ['-DOVERLOAD_C4', '-DOVERLOAD_R4'] - self.other_compile_flags = ['-i4'] + if self.compiler == 'intel': + self.other_compile_flags = ['-i4'] + elif self.compiler == 'gfortran': + self.other_compile_flags = ['-freal-8-real-4'] self.executable_name = self.executable_name.strip('.x')+'_single.x' self.builddir = P(self.workdir, 'build', self.executable_name.split('.')[0]) diff --git a/src/extra/python/isca/experiment.py b/src/extra/python/isca/experiment.py index d176033d5..167078bda 100755 --- a/src/extra/python/isca/experiment.py +++ b/src/extra/python/isca/experiment.py @@ -74,7 +74,10 @@ def __init__(self, name, codebase, safe_mode=False, workbase=GFDL_WORK, database self.restartdir = P(self.datadir, 'restarts') # where restarts will be stored self.template_dir = P(_module_directory, 'templates') - self.env_source = get_env_file() + if self.codebase.env_name is not None: + self.env_source = get_env_file(self.codebase.env_name) + else: + self.env_source = get_env_file() self.templates = Environment(loader=FileSystemLoader(self.template_dir)) diff --git a/src/extra/python/isca/templates/mkmf.template.ubuntu_conda2 b/src/extra/python/isca/templates/mkmf.template.ubuntu_conda2 new file mode 100644 index 000000000..65d8eeb68 --- /dev/null +++ b/src/extra/python/isca/templates/mkmf.template.ubuntu_conda2 @@ -0,0 +1,28 @@ +# template for the gfortran compiler +# typical use with mkmf +# mkmf -t template.ifc -c"-Duse_libMPI -Duse_netCDF" path_names /usr/local/include +CPPFLAGS = `nc-config --cflags` +NC_INC=`nc-config --fflags` +NC_LIB=`nc-config --flibs` + +# FFLAGS: +# -cpp: Use the fortran preprocessor +# -ffree-line-length-none -fno-range-check: Allow arbitrarily long lines +# -fcray-pointer: Cray pointers don't alias other variables. +# -ftz: Denormal numbers are flushed to zero. +# -assume byterecl: Specifies the units for the OPEN statement as bytes. +# -shared-intel: Load intel libraries dynamically +# -i4: 4 byte integers +# -fdefault-real-8: 8 byte reals (compatability for some parts of GFDL code) +# -fdefault-double-8: 8 byte doubles (compat. with RRTM) +# -O2: Level 2 speed optimisations + +FFLAGS = $(CPPFLAGS) $(NC_LIB) -cpp -fcray-pointer \ + -O2 -ffree-line-length-none -fno-range-check \ + -fallow-invalid-boz -fallow-argument-mismatch + +FC = $(F90) +LD = $(F90) + +LDFLAGS = -lnetcdff -lnetcdf -lmpi +CFLAGS = -D__IFC From 272d573b635e8c9de40be0fc7bb46cbfdea6e6ae Mon Sep 17 00:00:00 2001 From: ntlewis Date: Wed, 18 Oct 2023 12:27:10 +0100 Subject: [PATCH 11/19] stuff to do with ozone in grey model --- .../two_stream_gray_rad.F90 | 98 +++++++++++++++++-- .../driver/solo/idealized_moist_phys.F90 | 2 +- 2 files changed, 93 insertions(+), 7 deletions(-) diff --git a/src/atmos_param/two_stream_gray_rad/two_stream_gray_rad.F90 b/src/atmos_param/two_stream_gray_rad/two_stream_gray_rad.F90 index 0a3fd3865..def71287e 100644 --- a/src/atmos_param/two_stream_gray_rad/two_stream_gray_rad.F90 +++ b/src/atmos_param/two_stream_gray_rad/two_stream_gray_rad.F90 @@ -94,7 +94,7 @@ module two_stream_gray_rad_mod character(len=32) :: rad_scheme = 'frierson' integer, parameter :: B_GEEN = 1, B_FRIERSON = 2, & - B_BYRNE = 3, B_SCHNEIDER_LIU=4 + B_BYRNE = 3, B_SCHNEIDER_LIU=4, B_DAVIS=5, B_NEIL=6 integer, private :: sw_scheme = B_FRIERSON integer, private :: lw_scheme = B_FRIERSON @@ -163,7 +163,7 @@ module two_stream_gray_rad_mod integer :: id_olr, id_swdn_sfc, id_swdn_toa, id_net_lw_surf, id_lwdn_sfc, id_lwup_sfc, & id_tdt_rad, id_tdt_solar, id_flux_rad, id_flux_lw, id_flux_sw, id_coszen, id_fracsun, & - id_lw_dtrans, id_lw_dtrans_win, id_sw_dtrans, id_co2, id_nudge_flux, id_tau + id_lw_dtrans, id_lw_dtrans_win, id_sw_dtrans, id_co2, id_nudge_flux, id_tau, id_sw_tau character(len=10), parameter :: mod_name = 'two_stream' @@ -229,6 +229,12 @@ subroutine two_stream_gray_rad_init(is, ie, js, je, num_levels, axes, Time, lonb else if(uppercase(trim(rad_scheme)) == 'BYRNE') then lw_scheme = B_BYRNE call error_mesg('two_stream_gray_rad','Using Byrne & OGorman (2013) radiation scheme.', NOTE) +else if(uppercase(trim(rad_scheme)) == 'DAVIS') then + lw_scheme = B_DAVIS + call error_mesg('two_stream_gray_rad','Using Davis and Birner (2019) radiation scheme.', NOTE) +else if(uppercase(trim(rad_scheme)) == 'NEIL') then + lw_scheme = B_NEIL + call error_mesg('two_stream_gray_rad','Using Neils radiation scheme.', NOTE) else if(uppercase(trim(rad_scheme)) == 'SCHNEIDER') then lw_scheme = B_SCHNEIDER_LIU call error_mesg('two_stream_gray_rad','Using Schneider & Liu (2009) radiation scheme for GIANT PLANETS.', NOTE) @@ -357,6 +363,12 @@ subroutine two_stream_gray_rad_init(is, ie, js, je, num_levels, axes, Time, lonb register_diag_field ( mod_name, 'flux_sw', axes(half), Time, & 'Net shortwave radiative flux (positive up)', & 'W/m^2', missing_value=missing_value ) + + + id_sw_tau = & + register_diag_field ( mod_name, 'sw_tau', axes(half), Time, & + 'shortwave optical depth', & + 'nondim', missing_value=missing_value ) id_coszen = & register_diag_field ( mod_name, 'coszen', axes(1:2), Time, & @@ -403,7 +415,7 @@ end subroutine two_stream_gray_rad_init ! ================================================================================== -subroutine two_stream_gray_rad_down (is, js, Time_diag, lat, lon, p_half, t, & +subroutine two_stream_gray_rad_down (is, js, Time_diag, lat, lon, p_full, p_half, t, & net_surf_sw_down, surf_lw_down, albedo, q, const_correct, nudge_out) ! NTL ICE CHANGE ! Begin the radiation calculation by computing downward fluxes. @@ -414,7 +426,7 @@ subroutine two_stream_gray_rad_down (is, js, Time_diag, lat, lon, p_half, t, real, intent(in), dimension(:,:) :: lat, lon, albedo real, intent(out), dimension(:,:) :: net_surf_sw_down real, intent(out), dimension(:,:) :: surf_lw_down -real, intent(in), dimension(:,:,:) :: t, q, p_half +real, intent(in), dimension(:,:,:) :: t, q, p_half, p_full real, intent(in), dimension(:,:) :: const_correct, nudge_out integer :: i, j, k, n, dyofyr @@ -427,9 +439,15 @@ subroutine two_stream_gray_rad_down (is, js, Time_diag, lat, lon, p_half, t, real, dimension(size(q,1), size(q,2)) :: ci_rho_dist, ci_J_f, ci_alpha, ci_beta, ci_phi, & ci_t, ci_l, ci_s, ci_delta1, ci_delta2, ci_J_p +! for davis ozone +real, dimension(size(q,1), size(q,2)) :: p_center, p_width real ,dimension(size(q,1),size(q,2),size(q,3)) :: co2f +! for neil ozone +real, dimension(size(q,1),size(q,2),size(q,3)) :: dist, Z +real :: H + n = size(t,3) @@ -563,12 +581,63 @@ subroutine two_stream_gray_rad_down (is, js, Time_diag, lat, lon, p_half, t, sw_down(:,:,k+1) = sw_down(:,:,k) * sw_dtrans(:,:,k) end do +case(B_DAVIS) + ! Default: Davis and Birner handling of SW radiation + ! SW optical thickness + sw_tau_0 = atm_abs * sqrt(pi) !/ 2. !(1.0 - sw_diff*sin(lat)**2)*atm_abs + p_center = (130. - 90. * cos(0.9*lat)**2.) * 100. + p_width = (20. + 40. * sin(lat)**2.) * 100. + + ! compute optical depths for each model level + ! solar_tau_0(j)*frac_oz*(.25*sqrt(pi)*(erf(((p_half(1,j,k)-140)/60))+1)+& + ! (cos(0.3*lat(1,j))**2)*0.3333*sqrt(pi)*(erf((((p_half(1,j,k)-p_center(j))/p_width(j))))+1)) + + do k = 1, n+1 + sw_tau(:,:,k) = sw_tau_0 * (1./4. * (erf_func((p_half(:,:,k) - 14000.)/6000.) + 1.) + & + 1./3. * cos(0.3*lat(:,:))**2. * (erf_func((p_half(:,:,k) - p_center(:,:))/p_width(:,:)) + 1.)) + end do + + ! compute downward shortwave flux + + do k = 1, n+1 + sw_down(:,:,k) = insolation(:,:) * exp(-sw_tau(:,:,k)) + end do + +case(B_NEIL) + + sw_tau_0 = atm_abs + + H = 287.*300./9.81 + do k = 1, n + Z(:,:,k) = -1 * H *log(p_full(:,:,k) / p_half(:,:,n+1)) + end do + + do k = 1, n + dist(:,:,k) = cos(-(Z(:,:,k) - 3.e4)*pi/3.e4) * (0.4 + 0.6*cos(lat(:,:))**2.) + end do + + where (Z .ge. (3.e4 + 3.e4/2.)) + dist(:,:,:) = 0. + elsewhere (Z .le. (3.e4 - 3.e4/2.)) + dist(:,:,:) = 0. + endwhere + + sw_tau(:,:,1) = 0. + do k = 1, n + sw_tau(:,:,k+1) = sw_tau_0*dist(:,:,k)*(p_half(:,:,k+1) - p_half(:,:,k))/grav + sw_tau(:,:,k) + end do + + do k = 1, n+1 + sw_down(:,:,k) = insolation(:,:) * exp(-sw_tau(:,:,k)) + end do + case(B_FRIERSON, B_BYRNE) ! Default: Frierson handling of SW radiation ! SW optical thickness sw_tau_0 = (1.0 - sw_diff*sin(lat)**2)*atm_abs ! compute optical depths for each model level + do k = 1, n+1 sw_tau(:,:,k) = sw_tau_0 * (p_half(:,:,k)/pstd_mks)**solar_exponent end do @@ -660,7 +729,7 @@ subroutine two_stream_gray_rad_down (is, js, Time_diag, lat, lon, p_half, t, lw_down(:,:,k+1) = lw_down(:,:,k)*lw_dtrans(:,:,k) + b(:,:,k)*(1. - lw_dtrans(:,:,k)) end do -case(B_FRIERSON) +case(B_FRIERSON,B_DAVIS,B_NEIL) ! longwave optical thickness function of latitude and pressure if (do_tl) then lat_tl = pi/2 - asin(cos(lat)*cos(lon)) @@ -765,6 +834,10 @@ subroutine two_stream_gray_rad_down (is, js, Time_diag, lat, lon, p_half, t, used = send_data ( id_tau, lw_tau_0, Time_diag) endif +if (id_sw_tau > 0) then + used = send_data ( id_sw_tau, sw_tau, Time_diag) +endif + !------- carbon dioxide concentration ------------ if ( id_co2 > 0 ) then used = send_data ( id_co2, carbon_conc, Time_diag) @@ -806,7 +879,7 @@ subroutine two_stream_gray_rad_up (is, js, Time_diag, lat, p_half, t_surf, t, td end do lw_up = lw_up + lw_up_win -case(B_FRIERSON, B_BYRNE) +case(B_FRIERSON, B_BYRNE, B_DAVIS, B_NEIL) ! compute upward longwave flux by integrating upward lw_up(:,:,n+1) = b_surf if (do_normal_integration_method) then @@ -931,6 +1004,19 @@ subroutine two_stream_gray_rad_end end subroutine two_stream_gray_rad_end +! =================================================================================== + +function erf_func(x) + + implicit none + + real, intent(in), dimension(:,:) :: x + real, dimension(size(x,1), size(x,2)) :: erf_func + + erf_func=(2/sqrt(pi))*(x/abs(x))*sqrt(1-exp(-(x**2)))*((sqrt(pi)/2)+ (31/200)*exp(-(x**2)) - (341/8000)*(exp(-2*(x**2)))) + +end function erf_func + ! ================================================================================== end module two_stream_gray_rad_mod diff --git a/src/atmos_spectral/driver/solo/idealized_moist_phys.F90 b/src/atmos_spectral/driver/solo/idealized_moist_phys.F90 index b1e8f1efc..acb50da81 100644 --- a/src/atmos_spectral/driver/solo/idealized_moist_phys.F90 +++ b/src/atmos_spectral/driver/solo/idealized_moist_phys.F90 @@ -1052,7 +1052,7 @@ subroutine idealized_moist_phys(Time, p_half, p_full, z_half, z_full, ug, vg, tg call two_stream_gray_rad_down(is, js, Time, & rad_lat(:,:), & rad_lon(:,:), & - p_half(:,:,:,current), & + p_full(:,:,:,current), p_half(:,:,:,current), & tg(:,:,:,previous), & net_surf_sw_down(:,:), & surf_lw_down(:,:), albedo, & From e218620115f0f105d5ae199dd1144bb4fbf34e80 Mon Sep 17 00:00:00 2001 From: ntlewis Date: Wed, 18 Oct 2023 12:28:35 +0100 Subject: [PATCH 12/19] remove modifications that only nudge ice for thickness greater than 0.1 --- .../driver/solo/mixed_layer.F90 | 40 +++++++++---------- 1 file changed, 20 insertions(+), 20 deletions(-) diff --git a/src/atmos_spectral/driver/solo/mixed_layer.F90 b/src/atmos_spectral/driver/solo/mixed_layer.F90 index 432d032db..ae7d96ab7 100644 --- a/src/atmos_spectral/driver/solo/mixed_layer.F90 +++ b/src/atmos_spectral/driver/solo/mixed_layer.F90 @@ -876,29 +876,29 @@ subroutine mixed_layer ( & t_surf_dependence = beta_t * CP_AIR + beta_lw -nudge_out = 0.0 -where ((h_thermo_ice .gt. 0) .and. (h_thermo_ice .le. 0.1)) -nudge_out = h_thermo_ice / (86400./4.) * L_thermo_ice -endwhere -corrected_flux = corrected_flux - nudge_out -const_nudge_correction = 0.0 -const_correct = 0.0 -const_nudge_correction = -area_weighted_global_mean(nudge_out) -if (.not. (const_nudge_correction .eq. 0.)) then - where (nudge_out .eq. 0.) - const_correct(:,:) = const_nudge_correction - endwhere - const_correct = const_correct * const_nudge_correction / area_weighted_global_mean(const_correct) - !const_correct(:,:) = -area_weighted_global_mean(nudge_out) !-area_weighted_global_mean(nudge_flux * L_thermo_ice/dt) - corrected_flux = corrected_flux - const_correct -endif +!nudge_out = 0.0 +!where ((h_thermo_ice .gt. 0) .and. (h_thermo_ice .le. 0.1)) +!nudge_out = h_thermo_ice / (86400./4.) * L_thermo_ice +!endwhere +!corrected_flux = corrected_flux - nudge_out +!const_nudge_correction = 0.0 +!const_correct = 0.0 +!const_nudge_correction = -area_weighted_global_mean(nudge_out) +!if (.not. (const_nudge_correction .eq. 0.)) then +! where (nudge_out .eq. 0.) +! const_correct(:,:) = const_nudge_correction +! endwhere +! const_correct = const_correct * const_nudge_correction / area_weighted_global_mean(const_correct) +! !const_correct(:,:) = -area_weighted_global_mean(nudge_out) !-area_weighted_global_mean(nudge_flux * L_thermo_ice/dt) +! corrected_flux = corrected_flux - const_correct +!endif nudge_out = 0.0 if (do_nudge_ice) then - where ((h_thermo_ice .gt. 0.1) .and. (h_thermo_ice_in .le. 0)) + where ((h_thermo_ice .gt. 0.0) .and. (h_thermo_ice_in .le. 0)) nudge_out = h_thermo_ice / thermo_ice_nudge_timescale * L_thermo_ice endwhere corrected_flux = corrected_flux - nudge_out @@ -1219,8 +1219,8 @@ subroutine mixed_layer ( & where (land_ice_mask) albedo_out(:,:) = albedo_value * land_albedo_prefactor elsewhere - where ( h_thermo_ice .gt. 0.1 ) - albedo_out(:,:) = albedo_value + (thermo_ice_albedo - albedo_value) * tanh((h_thermo_ice-0.1)/1.) + where ( h_thermo_ice .gt. 0.0 ) + albedo_out(:,:) = albedo_value + (thermo_ice_albedo - albedo_value) * tanh((h_thermo_ice)/1.) elsewhere albedo_out(:,:) = albedo_value endwhere @@ -1229,7 +1229,7 @@ subroutine mixed_layer ( & where (land_ice_mask) albedo_out(:,:) = albedo_value * land_albedo_prefactor elsewhere - where ( h_thermo_ice .gt. 0.1 ) + where ( h_thermo_ice .gt. 0.0 ) albedo_out(:,:) = thermo_ice_albedo elsewhere albedo_out(:,:) = albedo_value From 2630ea65ff15ac94d5769207d442a127c1ce548a Mon Sep 17 00:00:00 2001 From: ntlewis Date: Fri, 31 May 2024 11:13:27 +0100 Subject: [PATCH 13/19] tidy sea ice branch, so only code necessary for thermodynamic sea ice is included --- ci/environment-py3.9.yml | 2 +- src/atmos_param/dry_adj/dry_adj.F90 | 327 ------------------ src/atmos_param/my25_turb/my25_turb.F90 | 15 +- .../qe_moist_convection.F90 | 5 +- .../rrtm_radiation/rrtm_radiation.f90 | 56 +-- .../modules/{parkind.F90 => parkind.f90} | 12 +- ...mg_lw_setcoef.F90 => rrtmg_lw_setcoef.f90} | 5 +- .../modules/{parkind.F90 => parkind.f90} | 13 +- .../rrtmg_sw/gcm_model/src/rrtmg_sw_rad.f90 | 8 +- .../gcm_model/src/rrtmg_sw_rad.nomcica.f90 | 8 +- .../gcm_model/src/rrtmg_sw_spcvmc.f90 | 5 +- .../gcm_model/src/rrtmg_sw_spcvrt.f90 | 5 +- .../two_stream_gray_rad.F90 | 218 +----------- .../driver/solo/idealized_moist_phys.F90 | 61 +--- .../driver/solo/mixed_layer.F90 | 213 ++---------- .../model/spectral_dynamics.F90 | 2 +- src/extra/model/grey/path_names | 1 - src/extra/model/isca/path_names | 7 +- src/extra/model/socrates/path_names | 1 - src/extra/python/isca/codebase.py | 39 +-- src/extra/python/isca/experiment.py | 5 +- src/extra/python/isca/templates/compile.sh | 16 +- .../python/isca/templates/mkmf.template.ia64 | 4 +- 23 files changed, 91 insertions(+), 937 deletions(-) delete mode 100644 src/atmos_param/dry_adj/dry_adj.F90 rename src/atmos_param/rrtm_radiation/rrtmg_lw/gcm_model/modules/{parkind.F90 => parkind.f90} (65%) rename src/atmos_param/rrtm_radiation/rrtmg_lw/gcm_model/src/{rrtmg_lw_setcoef.F90 => rrtmg_lw_setcoef.f90} (99%) rename src/atmos_param/rrtm_radiation/rrtmg_sw/gcm_model/modules/{parkind.F90 => parkind.f90} (65%) diff --git a/ci/environment-py3.9.yml b/ci/environment-py3.9.yml index 345e74483..169e3adef 100644 --- a/ci/environment-py3.9.yml +++ b/ci/environment-py3.9.yml @@ -1,4 +1,4 @@ -name: isca_env2 +name: isca_env channels: - conda-forge dependencies: diff --git a/src/atmos_param/dry_adj/dry_adj.F90 b/src/atmos_param/dry_adj/dry_adj.F90 deleted file mode 100644 index a15c99cca..000000000 --- a/src/atmos_param/dry_adj/dry_adj.F90 +++ /dev/null @@ -1,327 +0,0 @@ -MODULE DRY_ADJ_MOD - - !======================================================================= - ! DRY ADIABATIC ADJUSTMENT - !======================================================================= - - - - use fms_mod, ONLY: file_exist, error_mesg, open_file, & - check_nml_error, mpp_pe, FATAL, & - NOTE, close_file, uppercase - - use constants_mod, ONLY: grav, kappa - !--------------------------------------------------------------------- - implicit none - private - - public :: dry_adj, dry_adj_init, dry_adj_end, dry_adj_bdgt - - !--------------------------------------------------------------------- - - character(len=128) :: version = '$Id: dry_adj.F90,v 17.0.4.1 2010/08/30 20:33:34 wfc Exp $' - character(len=128) :: tag = '$Name: testing $' - logical :: module_is_initialized = .false. - - !--------------------------------------------------------------------- - - real, parameter :: p00 = 1000.0E2 - real :: kappa_use = 0.0 ! modified kappa - ! (kappa_use = kappa * alpha) - ! initialised in dry_adj_init - - !--------------------------------------------------------------------- - ! --- NAMELIST - !--------------------------------------------------------------------- - ! itermax - Max number of iterations - !--------------------------------------------------------------------- - - integer :: itermax = 5 - real :: small = 0.001 - logical :: do_mcm_dry_adj = .true. - logical :: do_relax = .false. - real :: tau_relax = 7200. - real :: alpha = 1.0 ! values less than one achieve stable column. - - namelist /dry_adj_nml/ itermax, small, do_mcm_dry_adj, do_relax, tau_relax!, alpha - - !--------------------------------------------------------------------- - - contains - - !####################################################################### - !####################################################################### - - SUBROUTINE DRY_ADJ ( temp0, pres, pres_int, dtemp, timestep, mask ) - - !======================================================================= - ! DRY ADIABATIC ADJUSTMENT - !======================================================================= - !--------------------------------------------------------------------- - ! Arguments (Intent in) - ! temp0 - Temperature - ! pres - Pressure - ! pres_int - Pressure at layer interface - ! mask - OPTIONAL; floating point mask (0. or 1.) designating - ! where data is present - !--------------------------------------------------------------------- - real, intent(in), dimension(:,:,:) :: temp0, pres, pres_int - real, intent(in) :: timestep - real, intent(in), OPTIONAL, dimension(:,:,:) :: mask - - !--------------------------------------------------------------------- - ! Arguments (Intent out) - ! dtemp - Change in temperature - !--------------------------------------------------------------------- - real, intent(out), dimension(:,:,:) :: dtemp - - !--------------------------------------------------------------------- - ! (Intent local) - !--------------------------------------------------------------------- - - real, dimension(size(temp0,1),size(temp0,2),size(temp0,3)) :: & - temp, pi, theta, pixdp, dpres, ppp - - real, dimension(size(temp0,1),size(temp0,2)) :: store, xxx - logical, dimension(size(temp0,1),size(temp0,2)) :: do_daa - - integer :: kmax, iter, k, i, j - logical :: do_any, did_adj - - !==================================================================== - - ! --- Check to see if dry_adj has been initialized - if(.not. module_is_initialized ) CALL error_mesg( 'DRY_ADJ', & - 'dry_adj_init has not been called', FATAL ) - - ! --- Set dimensions - kmax = size( temp0, 3 ) - - ! --- Compute pressure thickness of layers - dpres(:,:,1:kmax) = pres_int(:,:,2:kmax+1) - pres_int(:,:,1:kmax) - - ! --- Copy input temperature - temp = temp0 - - ! --- Compute exner function - !do i = 1,size(temp0,1) - ! do j = 1,size(temp0,2) - ! pi(i,j,:) = (pres(i,j,:) / pres(i,j,size(temp0,3))) ** kappa_use - ! end do - !end do - ! --- Compute exner function - pi = ( pres / p00 ) ** kappa_use - - ! --- Compute product of pi and dpres - pixdp = pi * dpres - - ! --- Compute potential temperature - theta = temp / pi - - if(do_mcm_dry_adj) then - do k = 2,kmax - xxx = 0.5*kappa_use*(pres(:,:,k) - pres(:,:,k-1))/pres_int(:,:,k) - ppp(:,:,k) = (1.0 + xxx)/(1.0 - xxx) - enddo - endif - - !----------------------------------------------------------------- - ! iteration loop starts - !----------------------------------------------------------------- - do iter = 1,itermax - !----------------------------------------------------------------- - - did_adj = .false. - - do k = 1,kmax - 1 - ! ---------------------------------------------- - - ! --- Flag layers needing adjustment - if(do_mcm_dry_adj) then - do_daa(:,:) = temp(:,:,k+1) > ( temp(:,:,k)*ppp(:,:,k+1) + small ) - else - do_daa(:,:) = ( theta(:,:,k+1) - theta(:,:,k) ) > small - endif - - if( PRESENT( mask ) ) then - do_daa(:,:) = do_daa(:,:) .and. ( mask(:,:,k+1) > 0.5 ) - endif - do_any = ANY( do_daa(:,:) ) - - ! --- Do adjustment - if ( do_any ) then - if(do_mcm_dry_adj) then - where ( do_daa ) - temp(:,:,k) = (temp(:,:,k) * dpres(:,:,k ) + & - temp(:,:,k+1)* dpres(:,:,k+1) ) & - /(dpres(:,:,k) + ppp(:,:,k+1)*dpres(:,:,k+1)) - temp(:,:,k+1) = temp(:,:,k)*ppp(:,:,k+1) - end where - did_adj = .true. - else - where ( do_daa ) - store(:,:) = ( theta(:,:,k ) * pixdp(:,:,k ) & - + theta(:,:,k+1) * pixdp(:,:,k+1) ) & - / ( pixdp(:,:,k ) + pixdp(:,:,k+1) ) - theta(:,:,k ) = store(:,:) - theta(:,:,k+1) = store(:,:) - temp(:,:,k ) = pi(:,:,k ) * theta(:,:,k ) - temp(:,:,k+1) = pi(:,:,k+1) * theta(:,:,k+1) - end where - did_adj = .true. - endif - end if - - ! ---------------------------------------------- - end do - - ! --- If no adjusment made this pass, exit iteration loop. - if ( .not. did_adj ) go to 900 - - !----------------------------------------------------------------- - end do - !----------------------------------------------------------------- - ! iteration loop ends - !----------------------------------------------------------------- - if(.not.do_mcm_dry_adj) then - call error_mesg ('DRY_ADJ', 'Non-convergence in dry_adj', NOTE) - endif - 900 continue - - ! --- Compute change in temperature - if (do_relax) then - dtemp = (temp - temp0) * (1 - exp(-timestep / tau_relax)) - else - dtemp = temp - temp0 - endif - - !======================================================================= - end SUBROUTINE DRY_ADJ - - !##################################################################### - !##################################################################### - - SUBROUTINE DRY_ADJ_INIT(alpha_in) - - !======================================================================= - ! ***** INITIALIZE RAS - !======================================================================= - - real, intent(in) :: alpha_in - !--------------------------------------------------------------------- - ! (Intent local) - !--------------------------------------------------------------------- - - integer :: unit, io, ierr - - !===================================================================== - - !--------------------------------------------------------------------- - ! --- READ NAMELIST - !--------------------------------------------------------------------- - - if (file_exist('input.nml')) then - unit = open_file (file='input.nml', action='read') - ierr = 1 - do while (ierr /= 0) - read (unit, nml=dry_adj_nml, iostat=io, end=10) - ierr = check_nml_error(io, 'dry_adj_nml') - end do - 10 call close_file(unit) - endif - - !------- write version number and namelist --------- - - unit = open_file (file='logfile.out', action='append') - if ( mpp_pe() == 0 ) then - write (unit,'(/,80("="),/(a))') trim(version), trim(tag) - write (unit,nml=dry_adj_nml) - endif - call close_file(unit) - - !------------------------------------------------------------------- - - alpha = alpha_in - - kappa_use = kappa * alpha - - module_is_initialized = .true. - - !===================================================================== - end SUBROUTINE DRY_ADJ_INIT - - - !####################################################################### - !####################################################################### - - SUBROUTINE DRY_ADJ_END() - - !------------------------------------------------------------------- - - module_is_initialized = .false. - - !===================================================================== - end SUBROUTINE DRY_ADJ_END - - - !####################################################################### - !####################################################################### - - SUBROUTINE DRY_ADJ_BDGT ( dtemp, pres_int ) - - !======================================================================= - ! Budget check for dry adiabatic adjustment - a debugging tool - !======================================================================= - - !--------------------------------------------------------------------- - ! Arguments (Intent in) - ! dtemp - Temperature change - ! pres_int - Pressure at layer interface - !--------------------------------------------------------------------- - real, intent(in), dimension(:,:,:) :: dtemp, pres_int - - !--------------------------------------------------------------------- - ! (Intent local) - !--------------------------------------------------------------------- - - real, dimension(size(dtemp,1),size(dtemp,2),size(dtemp,3)) :: dpres - real :: sum_dtemp - integer :: imax, jmax, kmax, i, j, k - - !======================================================================= - - imax = size ( dtemp, 1 ) - jmax = size ( dtemp, 2 ) - kmax = size ( dtemp, 3 ) - - ! --- Compute pressure thickness of layers - dpres(:,:,1:kmax) = pres_int(:,:,2:kmax+1) - pres_int(:,:,1:kmax) - - ! --- Check budget - - do j = 1,jmax - do i = 1,imax - - sum_dtemp = 0. - - do k = 1,kmax - sum_dtemp = sum_dtemp + dtemp(i,j,k)*dpres(i,j,k) / Grav - end do - - if ( abs( sum_dtemp ) > 1.0e-4 ) then - print * - print *, ' DRY ADIABATIC ADJUSTMENT BUDGET CHECK AT i,j = ', i,j - print *, ' sum_dtemp = ', sum_dtemp - print *, 'STOP' - ! STOP - endif - - end do - end do - - !======================================================================= - end SUBROUTINE DRY_ADJ_BDGT - - !####################################################################### - !####################################################################### - end MODULE DRY_ADJ_MOD \ No newline at end of file diff --git a/src/atmos_param/my25_turb/my25_turb.F90 b/src/atmos_param/my25_turb/my25_turb.F90 index fd7f60358..989d342db 100644 --- a/src/atmos_param/my25_turb/my25_turb.F90 +++ b/src/atmos_param/my25_turb/my25_turb.F90 @@ -27,15 +27,14 @@ MODULE MY25_TURB_MOD use mpp_mod, only : input_nml_file use fms_mod, only : file_exist, open_namelist_file, error_mesg, & - FATAL, close_file, note, read_data, write_data, & + FATAL, close_file, note, read_data, & check_nml_error, mpp_pe, mpp_root_pe, & - write_version_number, stdlog, open_restart_file, nullify_domain + write_version_number, stdlog, open_restart_file use fms_io_mod, only : register_restart_field, restart_file_type, & save_restart, restore_state use tridiagonal_mod, only : tri_invert, close_tridiagonal use constants_mod, only : grav, vonkarm use monin_obukhov_mod, only : mo_diff - use transforms_mod, only: grid_domain !--------------------------------------------------------------------- implicit none @@ -707,16 +706,14 @@ SUBROUTINE MY25_TURB_INIT( ix, jx, kx ) ! --- Input TKE !--------------------------------------------------------------------- - id_restart = register_restart_field(Tur_restart, 'my25_turb.res.nc', 'TKE', TKE) + id_restart = register_restart_field(Tur_restart, 'my25_turb.res', 'TKE', TKE) if (file_exist( 'INPUT/my25_turb.res.nc' )) then if (mpp_pe() == mpp_root_pe() ) then call error_mesg ('my25_turb_mod', 'MY25_TURB_INIT:& &Reading netCDF formatted restart file: & &INPUT/my25_turb.res.nc', NOTE) endif - !call restore_state(Tur_restart) - call nullify_domain() - call read_data(trim('INPUT/my25_turb.res'), 'TKE', TKE, grid_domain) + call restore_state(Tur_restart) else if( FILE_EXIST( 'INPUT/my25_turb.res' ) ) then @@ -784,9 +781,7 @@ subroutine my25_turb_restart(timestamp) ! write out the restart data which is always needed, regardless of ! when the first donner calculation step is after restart. !---------------------------------------------------------------------- - !call save_restart(Tur_restart, timestamp) - call nullify_domain() - call write_data(trim('RESTART/my25_turb.res'), 'TKE', TKE, grid_domain) + call save_restart(Tur_restart, timestamp) end subroutine my25_turb_restart ! NAME="my25_turb_restart" diff --git a/src/atmos_param/qe_moist_convection/qe_moist_convection.F90 b/src/atmos_param/qe_moist_convection/qe_moist_convection.F90 index aa55e0ec6..af299fcfa 100644 --- a/src/atmos_param/qe_moist_convection/qe_moist_convection.F90 +++ b/src/atmos_param/qe_moist_convection/qe_moist_convection.F90 @@ -70,14 +70,13 @@ module qe_moist_convection_mod real :: val_inc = 0.01 real :: val_min = -1. ! calculated in get_lcl_temp_table_size real :: val_max = 1. ! calculated in get_lcl_temp_table_size - real :: precision = 1.e-7 real, parameter :: small = 1.e-10, & ! to avoid division by 0 in dry limit pref = 1.e5 real, allocatable, dimension(:) :: lcl_temp_table - namelist /qe_moist_convection_nml/ tau_bm, rhbm, Tmin, Tmax, val_inc, precision + namelist /qe_moist_convection_nml/ tau_bm, rhbm, Tmin, Tmax, val_inc !------------------------------------------------------------------ ! Description of namelist variables @@ -1085,7 +1084,7 @@ end subroutine get_lcl_temp real function lcl_temp(value, lcl_temp_guess) real, intent(in) :: value, lcl_temp_guess - ! real, parameter :: precision = 1.e-7 + real, parameter :: precision = 1.e-7 integer, parameter :: max_iter = 100 real :: T, dT integer :: iter diff --git a/src/atmos_param/rrtm_radiation/rrtm_radiation.f90 b/src/atmos_param/rrtm_radiation/rrtm_radiation.f90 index b547b9f05..808d5239d 100644 --- a/src/atmos_param/rrtm_radiation/rrtm_radiation.f90 +++ b/src/atmos_param/rrtm_radiation/rrtm_radiation.f90 @@ -79,7 +79,7 @@ module rrtm_vars ! cloud & aerosol optical depths, cloud and aerosol specific parameters. Set to zero real(kind=rb),allocatable,dimension(:,:,:) :: taucld,tauaer, sw_zro, zro_sw ! heating rates and fluxes, zenith angle when in-between radiation time steps - real(kind=rb),allocatable,dimension(:,:) :: sw_flux,lw_flux,zencos, olr, toa_sw, toa_sw_down! surface and TOA fluxes, cos(zenith angle) + real(kind=rb),allocatable,dimension(:,:) :: sw_flux,lw_flux,zencos, olr, toa_sw! surface and TOA fluxes, cos(zenith angle) ! dimension (lon x lat) real(kind=rb),allocatable,dimension(:,:,:) :: tdt_rad ! heating rate [K/s] ! dimension (lon x lat x pfull) @@ -162,8 +162,7 @@ module rrtm_vars ! for faster radiation calculation !!!!!! mp586 added for annual mean insolation!!!!! - logical :: do_toa_albedo = .false. - real(kind=rb):: alb_scaler = 1. + logical :: frierson_solar_rad =.false. real(kind=rb) :: del_sol = 0.95 ! frierson 2006 default = 1.4, but 0.95 gets the curve closer to the annual mean insolation real(kind=rb) :: del_sw = 0.0 !frierson 2006 default @@ -193,7 +192,7 @@ module rrtm_vars ! !-------------------- diagnostics fields ------------------------------- - integer :: id_tdt_rad, id_tdt_sw, id_tdt_lw, id_coszen, id_flux_sw, id_flux_lw, id_olr, id_toa_sw, id_toa_sw_down, id_albedo,id_ozone, id_co2, id_fracday, id_half_level_temp, id_full_level_temp + integer :: id_tdt_rad, id_tdt_sw, id_tdt_lw, id_coszen, id_flux_sw, id_flux_lw, id_olr, id_toa_sw, id_albedo,id_ozone, id_co2, id_fracday, id_half_level_temp, id_full_level_temp character(len=14), parameter :: mod_name_rad = 'rrtm_radiation' !s changed parameter name from mod_name to mod_name_rad as compiler objected, presumably because mod_name also defined in idealized_moist_physics.F90 after use rrtm_vars is included. real :: missing_value = -999. @@ -211,7 +210,7 @@ module rrtm_vars &lonstep, do_zm_tracers, do_zm_rad, & &do_precip_albedo, precip_albedo_mode, precip_albedo, precip_lat,& &do_read_co2, co2_file, co2_variable_name, use_dyofyr, solrad, & - &solday, equinox_day,solr_cnst,do_toa_albedo,alb_scaler + &solday, equinox_day,solr_cnst end module rrtm_vars !***************************************************************************************** @@ -306,10 +305,6 @@ subroutine rrtm_radiation_init(axes,Time,ncols,nlay,lonb,latb, Time_step) register_diag_field ( mod_name_rad, 'toa_sw', axes(1:2), Time, & 'Net TOA SW flux', & 'W/m2', missing_value=missing_value ) - id_toa_sw_down = & - register_diag_field ( mod_name_rad, 'toa_sw_down', axes(1:2), Time, & - 'Downward TOA SW flux', & - 'W/m2', missing_value=missing_value ) id_albedo = & register_diag_field ( mod_name_rad, 'rrtm_albedo', axes(1:2), Time, & 'Interactive albedo', & @@ -468,8 +463,6 @@ subroutine rrtm_radiation_init(axes,Time,ncols,nlay,lonb,latb, Time_step) allocate(olr(size(lonb,1)-1,size(latb,2)-1)) if(id_toa_sw > 0) & allocate(toa_sw(size(lonb,1)-1,size(latb,2)-1)) - if(id_toa_sw_down > 0) & - allocate(toa_sw_down(size(lonb,1)-1,size(latb,2)-1)) if(do_precip_albedo)allocate(rrtm_precip(size(lonb,1)-1,size(latb,2)-1)) if(store_intermediate_rad .or. id_tdt_rad > 0)& allocate(tdt_rad(size(lonb,1)-1,size(latb,2)-1,nlay)) @@ -598,11 +591,10 @@ subroutine run_rrtmg(is,js,Time,lat,lon,p_full,p_half,albedo,q,t,t_surf_rad,tdt, real(kind=rb),dimension(size(q,1)/lonstep,size(q,2),size(q,3) ) :: swijk,lwijk real(kind=rb),dimension(size(q,1)/lonstep,size(q,2)) :: swflxijk,lwflxijk real(kind=rb),dimension(ncols_rrt,nlay_rrt+1):: phalf,thalf - real(kind=rb),dimension(ncols_rrt) :: tsrf,cosz_rr,albedo_rr, cloud_coalb_spatial_rr + real(kind=rb),dimension(ncols_rrt) :: tsrf,cosz_rr,albedo_rr real(kind=rb) :: dlon,dlat,dj,di type(time_type) :: Time_loc real(kind=rb),dimension(size(q,1),size(q,2)) :: albedo_loc - real(kind=rb),dimension(size(q,1),size(q,2)) :: cloud_coalb_spatial real(kind=rb),dimension(size(q,1),size(q,2),size(q,3)) :: q_tmp, h2o_vmr real(kind=rb),dimension(size(q,1),size(q,2)) :: fracsun real(kind=rb),dimension(size(q,1),size(q,2)) :: p2 !mp586 addition for annual mean insolation @@ -656,15 +648,8 @@ subroutine run_rrtmg(is,js,Time,lat,lon,p_full,p_half,albedo,q,t,t_surf_rad,tdt, !!!!! mp586 addition for annual mean insolation !!!!! !!!! following https://github.com/sit23/Isca/blob/master/src/atmos_param/socrates/interface/socrates_interface.F90#L888 !!!! - p2 = (1. - 3.*sin(lat(:,:))**2)/4. - if (do_toa_albedo) then - !cloud_coalb_spatial(:,:) = (0.75 + 0.15 * 2. * p2(:,:))*alb_scaler - (alb_scaler - 1.) - cloud_coalb_spatial(:,:) = (0.87 + 0.05 * 2. * p2(:,:))*alb_scaler - (alb_scaler - 1.) - else - cloud_coalb_spatial(:,:) = 1. - endif - - if (frierson_solar_rad) then + if (frierson_solar_rad) then + p2 = (1. - 3.*sin(lat(:,:))**2)/4. coszen = 0.25 * (1.0 + del_sol * p2 + del_sw * sin(lat(:,:))) rrsun = 1 ! needs to be set, set to 1 so that stellar_radiation is unchanged in socrates_interface else @@ -823,7 +808,6 @@ subroutine run_rrtmg(is,js,Time,lat,lon,p_full,p_half,albedo,q,t,t_surf_rad,tdt, cosz_rr = reshape(coszen (1:si:lonstep,:),(/ si*sj/lonstep /)) albedo_rr = reshape(albedo_loc(1:si:lonstep,:),(/ si*sj/lonstep /)) tsrf = reshape(t_surf_rad(1:si:lonstep,:),(/ si*sj/lonstep /)) - cloud_coalb_spatial_rr = reshape(cloud_coalb_spatial (1:si:lonstep,:),(/ si*sj/lonstep /)) !--------------------------------------------------------------------------------------------------------------- ! now actually run RRTM @@ -850,7 +834,7 @@ subroutine run_rrtmg(is,js,Time,lat,lon,p_full,p_half,albedo,q,t,t_surf_rad,tdt, pfull , phalf , tfull , thalf , tsrf , & h2o , o3 , co2 , ch4_val*ones , n2o_val*ones , o2_val*ones , & albedo_rr , albedo_rr, albedo_rr, albedo_rr, & - cosz_rr , cloud_coalb_spatial_rr, solrad , dyofyr , solr_cnst, & + cosz_rr , solrad , dyofyr , solr_cnst, & inflglw , iceflglw , liqflglw , & ! cloud parameters zeros , taucld , sw_zro , sw_zro , sw_zro , & @@ -864,7 +848,7 @@ subroutine run_rrtmg(is,js,Time,lat,lon,p_full,p_half,albedo,q,t,t_surf_rad,tdt, pfull , phalf , tfull , thalf , tsrf , & h2o , o3 , co2 , zeros , zeros, zeros, & albedo_rr , albedo_rr, albedo_rr, albedo_rr, & - cosz_rr , cloud_coalb_spatial_rr, solrad , dyofyr , solr_cnst, & + cosz_rr , solrad , dyofyr , solr_cnst, & inflglw , iceflglw , liqflglw , & ! cloud parameters zeros , taucld , sw_zro , sw_zro , sw_zro , & @@ -1023,21 +1007,6 @@ subroutine run_rrtmg(is,js,Time,lat,lon,p_full,p_half,albedo,q,t,t_surf_rad,tdt, enddo enddo endif - if(id_toa_sw_down > 0)then - swflxijk = reshape(swdflx(:,sk+1),(/ si/lonstep,sj /)) ! net TOA SW flux, +ve down - dlon=1./lonstep - do i=1,size(swijk,1) - i1 = i+1 - ! close toroidally - if(i1 > size(swijk,1)) i1=1 - do ij=1,lonstep - di = (ij-1)*dlon - ij1 = (i-1)*lonstep + ij - toa_sw_down(ij1,:) = di*swflxijk(i1,:) + (1.-di)*swflxijk(i ,:) - enddo - enddo - endif - @@ -1060,8 +1029,7 @@ subroutine write_diag_rrtm(Time,is,js,ozone,cotwo,fracday,albedo_loc, t_full) use rrtm_vars,only: sw_flux,lw_flux,zencos,tdt_rad,tdt_sw_rad,tdt_lw_rad,t_half,& &id_tdt_rad,id_tdt_sw,id_tdt_lw,id_coszen,& &id_flux_sw,id_flux_lw,id_albedo,id_ozone, id_co2, id_fracday,& - &id_olr,id_toa_sw,id_toa_sw_down,olr,toa_sw,toa_sw_down,& - id_half_level_temp, id_full_level_temp + &id_olr,id_toa_sw,olr,toa_sw, id_half_level_temp, id_full_level_temp use diag_manager_mod, only: register_diag_field, send_data use time_manager_mod,only: time_type @@ -1113,10 +1081,6 @@ subroutine write_diag_rrtm(Time,is,js,ozone,cotwo,fracday,albedo_loc, t_full) if ( id_toa_sw > 0 ) then used = send_data ( id_toa_sw, toa_sw, Time) endif -!------- Downward SW toa flux ------------ - if ( id_toa_sw_down > 0 ) then - used = send_data ( id_toa_sw_down, toa_sw_down, Time) - endif !------- Interactive albedo ------------ if ( present(albedo_loc)) then ! used = send_data ( id_albedo, albedo_loc, Time, is, js ) diff --git a/src/atmos_param/rrtm_radiation/rrtmg_lw/gcm_model/modules/parkind.F90 b/src/atmos_param/rrtm_radiation/rrtmg_lw/gcm_model/modules/parkind.f90 similarity index 65% rename from src/atmos_param/rrtm_radiation/rrtmg_lw/gcm_model/modules/parkind.F90 rename to src/atmos_param/rrtm_radiation/rrtmg_lw/gcm_model/modules/parkind.f90 index 8719337f9..5667d460b 100644 --- a/src/atmos_param/rrtm_radiation/rrtmg_lw/gcm_model/modules/parkind.F90 +++ b/src/atmos_param/rrtm_radiation/rrtmg_lw/gcm_model/modules/parkind.f90 @@ -15,11 +15,7 @@ module parkind ! integer kinds ! ------------- ! -#ifdef OVERLOAD_R4 - integer, parameter :: kind_ib = selected_int_kind(6) ! 4 byte integer -#else - integer, parameter :: kind_ib = selected_int_kind(13) ! 8 byte integer -#endif + integer, parameter :: kind_ib = selected_int_kind(13) ! 8 byte integer integer, parameter :: kind_im = selected_int_kind(6) ! 4 byte integer integer, parameter :: kind_in = kind(1) ! native integer @@ -27,11 +23,7 @@ module parkind ! real kinds ! ---------- ! -#ifdef OVERLOAD_R4 - integer, parameter :: kind_rb = selected_real_kind(6) ! 4 byte real -#else - integer, parameter :: kind_rb = selected_real_kind(13) ! 8 byte real -#endif + integer, parameter :: kind_rb = selected_real_kind(12) ! 8 byte real integer, parameter :: kind_rm = selected_real_kind(6) ! 4 byte real integer, parameter :: kind_rn = kind(1.0) ! native real diff --git a/src/atmos_param/rrtm_radiation/rrtmg_lw/gcm_model/src/rrtmg_lw_setcoef.F90 b/src/atmos_param/rrtm_radiation/rrtmg_lw/gcm_model/src/rrtmg_lw_setcoef.f90 similarity index 99% rename from src/atmos_param/rrtm_radiation/rrtmg_lw/gcm_model/src/rrtmg_lw_setcoef.F90 rename to src/atmos_param/rrtm_radiation/rrtmg_lw/gcm_model/src/rrtmg_lw_setcoef.f90 index 3d1afcf8d..cc65c677a 100644 --- a/src/atmos_param/rrtm_radiation/rrtmg_lw/gcm_model/src/rrtmg_lw_setcoef.F90 +++ b/src/atmos_param/rrtm_radiation/rrtmg_lw/gcm_model/src/rrtmg_lw_setcoef.f90 @@ -253,11 +253,8 @@ subroutine setcoef(nlayers, istart, pavel, tavel, tz, tbound, semiss, & ! layer pressure. Store them in JP and JP1. Store in FP the ! fraction of the difference (in ln(pressure)) between these ! two values that the layer pressure lies. -#ifdef OVERLOAD_R4 - plog = alog(pavel(lay)) -#else +! plog = alog(pavel(lay)) plog = dlog(pavel(lay)) -#endif jp(lay) = int(36._rb - 5*(plog+0.04_rb)) if (jp(lay) .lt. 1) then jp(lay) = 1 diff --git a/src/atmos_param/rrtm_radiation/rrtmg_sw/gcm_model/modules/parkind.F90 b/src/atmos_param/rrtm_radiation/rrtmg_sw/gcm_model/modules/parkind.f90 similarity index 65% rename from src/atmos_param/rrtm_radiation/rrtmg_sw/gcm_model/modules/parkind.F90 rename to src/atmos_param/rrtm_radiation/rrtmg_sw/gcm_model/modules/parkind.f90 index 50a2b1f26..5667d460b 100644 --- a/src/atmos_param/rrtm_radiation/rrtmg_sw/gcm_model/modules/parkind.F90 +++ b/src/atmos_param/rrtm_radiation/rrtmg_sw/gcm_model/modules/parkind.f90 @@ -15,12 +15,7 @@ module parkind ! integer kinds ! ------------- ! - -#ifdef OVERLOAD_R4 - integer, parameter :: kind_ib = selected_int_kind(6) ! 4 byte integer -#else - integer, parameter :: kind_ib = selected_int_kind(13) ! 8 byte integer -#endif + integer, parameter :: kind_ib = selected_int_kind(13) ! 8 byte integer integer, parameter :: kind_im = selected_int_kind(6) ! 4 byte integer integer, parameter :: kind_in = kind(1) ! native integer @@ -28,11 +23,7 @@ module parkind ! real kinds ! ---------- ! -#ifdef OVERLOAD_R4 - integer, parameter :: kind_rb = selected_real_kind(6) ! 4 byte real -#else - integer, parameter :: kind_rb = selected_real_kind(13) ! 8 byte real -#endif + integer, parameter :: kind_rb = selected_real_kind(12) ! 8 byte real integer, parameter :: kind_rm = selected_real_kind(6) ! 4 byte real integer, parameter :: kind_rn = kind(1.0) ! native real diff --git a/src/atmos_param/rrtm_radiation/rrtmg_sw/gcm_model/src/rrtmg_sw_rad.f90 b/src/atmos_param/rrtm_radiation/rrtmg_sw/gcm_model/src/rrtmg_sw_rad.f90 index a9071c78d..1ed25daad 100644 --- a/src/atmos_param/rrtm_radiation/rrtmg_sw/gcm_model/src/rrtmg_sw_rad.f90 +++ b/src/atmos_param/rrtm_radiation/rrtmg_sw/gcm_model/src/rrtmg_sw_rad.f90 @@ -82,8 +82,7 @@ subroutine rrtmg_sw & play ,plev ,tlay ,tlev ,tsfc , & h2ovmr , o3vmr ,co2vmr ,ch4vmr ,n2ovmr ,o2vmr , & asdir ,asdif ,aldir ,aldif , & - coszen , cloud_coalb_spatial, & - adjes ,dyofyr ,scon , & + coszen ,adjes ,dyofyr ,scon , & inflgsw ,iceflgsw,liqflgsw,cldfmcl , & taucmcl ,ssacmcl ,asmcmcl ,fsfcmcl , & ciwpmcl ,clwpmcl ,reicmcl ,relqmcl , & @@ -243,7 +242,6 @@ subroutine rrtmg_sw & real(kind=rb), intent(in) :: adjes ! Flux adjustment for Earth/Sun distance real(kind=rb), intent(in) :: coszen(:) ! Cosine of solar zenith angle ! Dimensions: (ncol) - real(kind=rb), intent(in) :: cloud_coalb_spatial(:) ! prescribed coalbedo for clouds real(kind=rb), intent(in) :: scon ! Solar constant (W/m2) integer(kind=im), intent(in) :: inflgsw ! Flag for cloud optical properties @@ -342,7 +340,6 @@ subroutine rrtmg_sw & ! real(kind=rb) :: earth_sun ! function for Earth/Sun distance factor real(kind=rb) :: cossza ! Cosine of solar zenith angle - real(kind=rb) :: cloud_coalb ! co albedo for prescribed clouds real(kind=rb) :: adjflux(jpband) ! adjustment for current Earth/Sun distance real(kind=rb) :: solvar(jpband) ! solar constant scaling factor from rrtmg_sw ! default value of 1368.22 Wm-2 at 1 AU @@ -572,7 +569,6 @@ subroutine rrtmg_sw & ! is below horizon cossza = coszen(iplon) - cloud_coalb = cloud_coalb_spatial(iplon) if (cossza .lt. zepzen) cossza = zepzen @@ -699,7 +695,7 @@ subroutine rrtmg_sw & (nlayers, istart, iend, icpr, idelm, iout, & pavel, tavel, pz, tz, tbound, albdif, albdir, & zcldfmc, ztaucmc, zasycmc, zomgcmc, ztaormc, & - ztaua, zasya, zomga, cossza, cloud_coalb, coldry, wkl, adjflux, & + ztaua, zasya, zomga, cossza, coldry, wkl, adjflux, & laytrop, layswtch, laylow, jp, jt, jt1, & co2mult, colch4, colco2, colh2o, colmol, coln2o, colo2, colo3, & fac00, fac01, fac10, fac11, & diff --git a/src/atmos_param/rrtm_radiation/rrtmg_sw/gcm_model/src/rrtmg_sw_rad.nomcica.f90 b/src/atmos_param/rrtm_radiation/rrtmg_sw/gcm_model/src/rrtmg_sw_rad.nomcica.f90 index 8bf4ddabb..c2d1e68e5 100644 --- a/src/atmos_param/rrtm_radiation/rrtmg_sw/gcm_model/src/rrtmg_sw_rad.nomcica.f90 +++ b/src/atmos_param/rrtm_radiation/rrtmg_sw/gcm_model/src/rrtmg_sw_rad.nomcica.f90 @@ -80,8 +80,7 @@ subroutine rrtmg_sw & play ,plev ,tlay ,tlev ,tsfc , & h2ovmr ,o3vmr ,co2vmr ,ch4vmr ,n2ovmr ,o2vmr, & asdir ,asdif ,aldir ,aldif , & - coszen ,cloud_coalb_spatial, & - adjes ,dyofyr ,scon , & + coszen ,adjes ,dyofyr ,scon , & inflgsw ,iceflgsw,liqflgsw,cldfr , & taucld ,ssacld ,asmcld ,fsfcld , & cicewp ,cliqwp ,reice ,reliq , & @@ -229,7 +228,6 @@ subroutine rrtmg_sw & real(kind=rb), intent(in) :: coszen(:) ! Cosine of solar zenith angle ! Dimensions: (ncol) - real(kind=rb), intent(in) :: cloud_coalb_spatial(:) ! prescribed coalbedo for clouds real(kind=rb), intent(in) :: adjes ! Flux adjustment for Earth/Sun distance integer(kind=im), intent(in) :: dyofyr ! Day of the year (used to get Earth/Sun ! distance if adjflx not provided) @@ -331,7 +329,6 @@ subroutine rrtmg_sw & ! real(kind=rb) :: earth_sun ! function for Earth/Sun distance factor real(kind=rb) :: cossza ! Cosine of solar zenith angle - real(kind=rb) :: cloud_coalb ! co albedo for prescribed clouds real(kind=rb) :: adjflux(jpband) ! adjustment for current Earth/Sun distance real(kind=rb) :: solvar(jpband) ! solar constant scaling factor from rrtmg_sw ! default value of 1368.22 Wm-2 at 1 AU @@ -560,7 +557,6 @@ subroutine rrtmg_sw & ! is below horizon cossza = coszen(iplon) - cloud_coalb = cloud_coalb_spatial(iplon) if (cossza .lt. zepzen) cossza = zepzen @@ -684,7 +680,7 @@ subroutine rrtmg_sw & (nlayers, istart, iend, icpr, idelm, iout, & pavel, tavel, pz, tz, tbound, albdif, albdir, & cldfrac, ztauc, zasyc, zomgc, ztaucorig, & - ztaua, zasya, zomga, cossza, cloud_coalb, coldry, wkl, adjflux, & + ztaua, zasya, zomga, cossza, coldry, wkl, adjflux, & laytrop, layswtch, laylow, jp, jt, jt1, & co2mult, colch4, colco2, colh2o, colmol, coln2o, colo2, colo3, & fac00, fac01, fac10, fac11, & diff --git a/src/atmos_param/rrtm_radiation/rrtmg_sw/gcm_model/src/rrtmg_sw_spcvmc.f90 b/src/atmos_param/rrtm_radiation/rrtmg_sw/gcm_model/src/rrtmg_sw_spcvmc.f90 index d448648d5..749b117f5 100644 --- a/src/atmos_param/rrtm_radiation/rrtmg_sw/gcm_model/src/rrtmg_sw_spcvmc.f90 +++ b/src/atmos_param/rrtm_radiation/rrtmg_sw/gcm_model/src/rrtmg_sw_spcvmc.f90 @@ -35,7 +35,7 @@ subroutine spcvmc_sw & (nlayers, istart, iend, icpr, idelm, iout, & pavel, tavel, pz, tz, tbound, palbd, palbp, & pcldfmc, ptaucmc, pasycmc, pomgcmc, ptaormc, & - ptaua, pasya, pomga, prmu0, cloud_coalb, coldry, wkl, adjflux, & + ptaua, pasya, pomga, prmu0, coldry, wkl, adjflux, & laytrop, layswtch, laylow, jp, jt, jt1, & co2mult, colch4, colco2, colh2o, colmol, coln2o, colo2, colo3, & fac00, fac01, fac10, fac11, & @@ -124,7 +124,6 @@ subroutine spcvmc_sw & real(kind=rb), intent(in) :: palbp(:) ! surface albedo (direct) ! Dimensions: (nbndsw) real(kind=rb), intent(in) :: prmu0 ! cosine of solar zenith angle - real(kind=rb), intent(in) :: cloud_coalb ! prescribed cloud co-albedo real(kind=rb), intent(in) :: pcldfmc(:,:) ! cloud fraction [mcica] ! Dimensions: (nlayers,ngptsw) real(kind=rb), intent(in) :: ptaucmc(:,:) ! cloud optical depth [mcica] @@ -314,7 +313,7 @@ subroutine spcvmc_sw & iw = iw+1 ! Apply adjustment for correct Earth/Sun distance and zenith angle to incoming solar flux - zincflx(iw) = adjflux(jb) * zsflxzen(iw) * prmu0 * cloud_coalb + zincflx(iw) = adjflux(jb) * zsflxzen(iw) * prmu0 ! zincflux = zincflux + adjflux(jb) * zsflxzen(iw) * prmu0 ! inactive ! Compute layer reflectances and transmittances for direct and diffuse sources, diff --git a/src/atmos_param/rrtm_radiation/rrtmg_sw/gcm_model/src/rrtmg_sw_spcvrt.f90 b/src/atmos_param/rrtm_radiation/rrtmg_sw/gcm_model/src/rrtmg_sw_spcvrt.f90 index b92030be0..06c67fb8d 100644 --- a/src/atmos_param/rrtm_radiation/rrtmg_sw/gcm_model/src/rrtmg_sw_spcvrt.f90 +++ b/src/atmos_param/rrtm_radiation/rrtmg_sw/gcm_model/src/rrtmg_sw_spcvrt.f90 @@ -35,7 +35,7 @@ subroutine spcvrt_sw & (nlayers, istart, iend, icpr, idelm, iout, & pavel, tavel, pz, tz, tbound, palbd, palbp, & pclfr, ptauc, pasyc, pomgc, ptaucorig, & - ptaua, pasya, pomga, prmu0, cloud_coalb, coldry, wkl, adjflux, & + ptaua, pasya, pomga, prmu0, coldry, wkl, adjflux, & laytrop, layswtch, laylow, jp, jt, jt1, & co2mult, colch4, colco2, colh2o, colmol, coln2o, colo2, colo3, & fac00, fac01, fac10, fac11, & @@ -121,7 +121,6 @@ subroutine spcvrt_sw & real(kind=rb), intent(in) :: palbp(:) ! surface albedo (direct) ! Dimensions: (nbndsw) real(kind=rb), intent(in) :: prmu0 ! cosine of solar zenith angle - real(kind=rb), intent(in) :: cloud_coalb ! prescribed cloud co-albedo real(kind=rb), intent(in) :: pclfr(:) ! cloud fraction ! Dimensions: (nlayers) real(kind=rb), intent(in) :: ptauc(:,:) ! cloud optical depth @@ -313,7 +312,7 @@ subroutine spcvrt_sw & iw = iw+1 ! Apply adjustments for correct Earth/Sun distance and zenith angle to incoming solar flux - zincflx(iw) = adjflux(jb) * zsflxzen(iw) * prmu0 * cloud_coalb + zincflx(iw) = adjflux(jb) * zsflxzen(iw) * prmu0 ! zincflux = zincflux + adjflux(jb) * zsflxzen(iw) * prmu0 ! inactive ! Compute layer reflectances and transmittances for direct and diffuse sources, diff --git a/src/atmos_param/two_stream_gray_rad/two_stream_gray_rad.F90 b/src/atmos_param/two_stream_gray_rad/two_stream_gray_rad.F90 index def71287e..a9a8d1615 100644 --- a/src/atmos_param/two_stream_gray_rad/two_stream_gray_rad.F90 +++ b/src/atmos_param/two_stream_gray_rad/two_stream_gray_rad.F90 @@ -28,7 +28,7 @@ module two_stream_gray_rad_mod mpp_pe, close_file, error_mesg, & NOTE, FATAL, uppercase - use constants_mod, only: stefan, cp_air, grav, pstd_mks, pstd_mks_earth, seconds_per_sol, orbital_period, radius + use constants_mod, only: stefan, cp_air, grav, pstd_mks, pstd_mks_earth, seconds_per_sol, orbital_period use diag_manager_mod, only: register_diag_field, send_data @@ -85,16 +85,11 @@ module two_stream_gray_rad_mod real :: equinox_day = 0.75 !s Fraction of year [0,1] where NH autumn equinox occurs (only really useful if calendar has defined months). logical :: use_time_average_coszen = .false. !s if .true., then time-averaging is done on coszen so that insolation doesn't depend on timestep real :: dt_rad_avg = -1 -logical :: do_tl = .false. -logical :: do_closein = .false. -logical :: do_normal_integration_method = .true. -real :: R_stellar = 476700000. -real :: d_stellar = 1117512000.0 character(len=32) :: rad_scheme = 'frierson' integer, parameter :: B_GEEN = 1, B_FRIERSON = 2, & - B_BYRNE = 3, B_SCHNEIDER_LIU=4, B_DAVIS=5, B_NEIL=6 + B_BYRNE = 3, B_SCHNEIDER_LIU=4 integer, private :: sw_scheme = B_FRIERSON integer, private :: lw_scheme = B_FRIERSON @@ -124,7 +119,6 @@ module two_stream_gray_rad_mod real :: bog_mu = 1.0 real, allocatable, dimension(:,:) :: insolation, p2, lw_tau_0, sw_tau_0 !s albedo now defined in mixed_layer_init -real, allocatable, dimension(:,:) :: lat_tl, lat_tl_insol real, allocatable, dimension(:,:) :: b_surf, b_surf_gp real, allocatable, dimension(:,:,:) :: b, tdt_rad, tdt_solar real, allocatable, dimension(:,:,:) :: lw_up, lw_down, lw_flux, sw_up, sw_down, sw_flux, rad_flux @@ -156,14 +150,14 @@ module two_stream_gray_rad_mod do_read_co2, co2_file, co2_variable_name, solday, equinox_day, bog_a, bog_b, bog_mu, & use_time_average_coszen, dt_rad_avg,& diabatic_acce, & !Schneider Liu values - do_toa_albedo, do_tl, do_closein, R_stellar, d_stellar, do_normal_integration_method + do_toa_albedo !================================================================================== !-------------------- diagnostics fields ------------------------------- integer :: id_olr, id_swdn_sfc, id_swdn_toa, id_net_lw_surf, id_lwdn_sfc, id_lwup_sfc, & id_tdt_rad, id_tdt_solar, id_flux_rad, id_flux_lw, id_flux_sw, id_coszen, id_fracsun, & - id_lw_dtrans, id_lw_dtrans_win, id_sw_dtrans, id_co2, id_nudge_flux, id_tau, id_sw_tau + id_lw_dtrans, id_lw_dtrans_win, id_sw_dtrans, id_co2, id_nudge_flux character(len=10), parameter :: mod_name = 'two_stream' @@ -229,12 +223,6 @@ subroutine two_stream_gray_rad_init(is, ie, js, je, num_levels, axes, Time, lonb else if(uppercase(trim(rad_scheme)) == 'BYRNE') then lw_scheme = B_BYRNE call error_mesg('two_stream_gray_rad','Using Byrne & OGorman (2013) radiation scheme.', NOTE) -else if(uppercase(trim(rad_scheme)) == 'DAVIS') then - lw_scheme = B_DAVIS - call error_mesg('two_stream_gray_rad','Using Davis and Birner (2019) radiation scheme.', NOTE) -else if(uppercase(trim(rad_scheme)) == 'NEIL') then - lw_scheme = B_NEIL - call error_mesg('two_stream_gray_rad','Using Neils radiation scheme.', NOTE) else if(uppercase(trim(rad_scheme)) == 'SCHNEIDER') then lw_scheme = B_SCHNEIDER_LIU call error_mesg('two_stream_gray_rad','Using Schneider & Liu (2009) radiation scheme for GIANT PLANETS.', NOTE) @@ -280,8 +268,6 @@ subroutine two_stream_gray_rad_init(is, ie, js, je, num_levels, axes, Time, lonb allocate (olr (ie-is+1, je-js+1)) allocate (net_lw_surf (ie-is+1, je-js+1)) allocate (toa_sw_in (ie-is+1, je-js+1)) -allocate (lat_tl (ie-is+1, je-js+1)) -allocate (lat_tl_insol (ie-is+1, je-js+1)) allocate (insolation (ie-is+1, je-js+1)) allocate (p2 (ie-is+1, je-js+1)) @@ -363,23 +349,10 @@ subroutine two_stream_gray_rad_init(is, ie, js, je, num_levels, axes, Time, lonb register_diag_field ( mod_name, 'flux_sw', axes(half), Time, & 'Net shortwave radiative flux (positive up)', & 'W/m^2', missing_value=missing_value ) - - - id_sw_tau = & - register_diag_field ( mod_name, 'sw_tau', axes(half), Time, & - 'shortwave optical depth', & - 'nondim', missing_value=missing_value ) - id_coszen = & register_diag_field ( mod_name, 'coszen', axes(1:2), Time, & 'cosine of zenith angle', & 'none', missing_value=missing_value ) - - id_tau = & - register_diag_field ( mod_name, 'tau', axes(1:2), Time, & - 'surface lw optical depth', & - 'none', missing_value=missing_value ) - id_fracsun = & register_diag_field ( mod_name, 'fracsun', axes(1:2), Time, & 'daylight fraction of time interval', & @@ -415,8 +388,8 @@ end subroutine two_stream_gray_rad_init ! ================================================================================== -subroutine two_stream_gray_rad_down (is, js, Time_diag, lat, lon, p_full, p_half, t, & - net_surf_sw_down, surf_lw_down, albedo, q, const_correct, nudge_out) ! NTL ICE CHANGE +subroutine two_stream_gray_rad_down (is, js, Time_diag, lat, lon, p_half, t, & + net_surf_sw_down, surf_lw_down, albedo, q, const_correct, nudge_out) ! Begin the radiation calculation by computing downward fluxes. ! This part of the calculation does not depend on the surface temperature. @@ -426,7 +399,7 @@ subroutine two_stream_gray_rad_down (is, js, Time_diag, lat, lon, p_full, p_half real, intent(in), dimension(:,:) :: lat, lon, albedo real, intent(out), dimension(:,:) :: net_surf_sw_down real, intent(out), dimension(:,:) :: surf_lw_down -real, intent(in), dimension(:,:,:) :: t, q, p_half, p_full +real, intent(in), dimension(:,:,:) :: t, q, p_half real, intent(in), dimension(:,:) :: const_correct, nudge_out integer :: i, j, k, n, dyofyr @@ -434,19 +407,9 @@ subroutine two_stream_gray_rad_down (is, js, Time_diag, lat, lon, p_full, p_half real :: r_seconds, frac_of_day, frac_of_year, gmt, time_since_ae, rrsun, day_in_s, r_solday, r_total_seconds, r_days, r_dt_rad_avg, dt_rad_radians logical :: used -! all for closein insolation -real :: ci_lat_f, ci_lat_d -real, dimension(size(q,1), size(q,2)) :: ci_rho_dist, ci_J_f, ci_alpha, ci_beta, ci_phi, & - ci_t, ci_l, ci_s, ci_delta1, ci_delta2, ci_J_p - -! for davis ozone -real, dimension(size(q,1), size(q,2)) :: p_center, p_width real ,dimension(size(q,1),size(q,2),size(q,3)) :: co2f -! for neil ozone -real, dimension(size(q,1),size(q,2),size(q,3)) :: dist, Z -real :: H n = size(t,3) @@ -496,58 +459,6 @@ subroutine two_stream_gray_rad_down (is, js, Time_diag, lat, lon, p_full, p_half insolation = insolation * (0.75 + 0.15 * 2. * p2) endif - insolation = insolation - -else if (do_tl) then - - - if (do_closein) then - - lat_tl_insol = asin(cos(lat)*cos(lon)) + pi/2 - - ci_lat_f = acos((radius+R_stellar) / d_stellar) - ci_lat_d = acos((radius-R_stellar) / d_stellar) - - ci_rho_dist = sqrt(d_stellar**2 + radius**2 - 2*d_stellar*radius*cos(lat_tl_insol)) - - ! J_f - ci_J_f = (d_stellar * cos(lat_tl_insol) - radius)*R_stellar**2 / ci_rho_dist**3 - - ! J_p - ci_alpha = acos((d_stellar*cos(lat_tl_insol) - radius) / ci_rho_dist) - ci_beta = pi/2 - ci_alpha - ci_phi = asin(R_stellar / ci_rho_dist) - ci_t = ci_rho_dist * cos(ci_phi) - ci_l = ci_t * sin(ci_phi) - ci_s = ci_t * sin(ci_phi) * cos(ci_alpha) - ci_delta1 = acos(cos(ci_phi)/cos(ci_beta)) - ci_delta2 = acos(tan(ci_beta)/tan(ci_phi)) - ci_J_p = ci_l * ci_s / pi / ci_t**2 * (pi - ci_delta2 + sin(ci_delta2)*cos(ci_delta2)) + 1 / pi * (ci_delta1 - sin(ci_delta1)*cos(ci_delta1)) - - where (lat_tl_insol .le. ci_lat_f) - insolation = solar_constant * (d_stellar / R_stellar)**2 * ci_J_f - elsewhere ((lat_tl_insol .gt. ci_lat_f) .and. (lat_tl_insol .lt. ci_lat_d)) - insolation = solar_constant * (d_stellar / R_stellar)**2 * ci_J_p - elsewhere (lat_tl_insol .ge. ci_lat_d) - insolation = 0. - endwhere - - - - - else - - insolation = solar_constant * cos(lat) * cos(lon) - - where (insolation .le. 0. ) - insolation = 0. - endwhere - endif - - - coszen = insolation / solar_constant - - else if (sw_scheme==B_SCHNEIDER_LIU) then insolation = (solar_constant/pi)*cos(lat) else @@ -557,7 +468,6 @@ subroutine two_stream_gray_rad_down (is, js, Time_diag, lat, lon, p_full, p_half if (do_toa_albedo) then insolation = insolation * (0.75 + 0.15 * 2. * p2) endif - insolation = insolation end if select case(sw_scheme) @@ -581,55 +491,6 @@ subroutine two_stream_gray_rad_down (is, js, Time_diag, lat, lon, p_full, p_half sw_down(:,:,k+1) = sw_down(:,:,k) * sw_dtrans(:,:,k) end do -case(B_DAVIS) - ! Default: Davis and Birner handling of SW radiation - ! SW optical thickness - sw_tau_0 = atm_abs * sqrt(pi) !/ 2. !(1.0 - sw_diff*sin(lat)**2)*atm_abs - p_center = (130. - 90. * cos(0.9*lat)**2.) * 100. - p_width = (20. + 40. * sin(lat)**2.) * 100. - - ! compute optical depths for each model level - ! solar_tau_0(j)*frac_oz*(.25*sqrt(pi)*(erf(((p_half(1,j,k)-140)/60))+1)+& - ! (cos(0.3*lat(1,j))**2)*0.3333*sqrt(pi)*(erf((((p_half(1,j,k)-p_center(j))/p_width(j))))+1)) - - do k = 1, n+1 - sw_tau(:,:,k) = sw_tau_0 * (1./4. * (erf_func((p_half(:,:,k) - 14000.)/6000.) + 1.) + & - 1./3. * cos(0.3*lat(:,:))**2. * (erf_func((p_half(:,:,k) - p_center(:,:))/p_width(:,:)) + 1.)) - end do - - ! compute downward shortwave flux - - do k = 1, n+1 - sw_down(:,:,k) = insolation(:,:) * exp(-sw_tau(:,:,k)) - end do - -case(B_NEIL) - - sw_tau_0 = atm_abs - - H = 287.*300./9.81 - do k = 1, n - Z(:,:,k) = -1 * H *log(p_full(:,:,k) / p_half(:,:,n+1)) - end do - - do k = 1, n - dist(:,:,k) = cos(-(Z(:,:,k) - 3.e4)*pi/3.e4) * (0.4 + 0.6*cos(lat(:,:))**2.) - end do - - where (Z .ge. (3.e4 + 3.e4/2.)) - dist(:,:,:) = 0. - elsewhere (Z .le. (3.e4 - 3.e4/2.)) - dist(:,:,:) = 0. - endwhere - - sw_tau(:,:,1) = 0. - do k = 1, n - sw_tau(:,:,k+1) = sw_tau_0*dist(:,:,k)*(p_half(:,:,k+1) - p_half(:,:,k))/grav + sw_tau(:,:,k) - end do - - do k = 1, n+1 - sw_down(:,:,k) = insolation(:,:) * exp(-sw_tau(:,:,k)) - end do case(B_FRIERSON, B_BYRNE) ! Default: Frierson handling of SW radiation @@ -729,17 +590,9 @@ subroutine two_stream_gray_rad_down (is, js, Time_diag, lat, lon, p_full, p_half lw_down(:,:,k+1) = lw_down(:,:,k)*lw_dtrans(:,:,k) + b(:,:,k)*(1. - lw_dtrans(:,:,k)) end do -case(B_FRIERSON,B_DAVIS,B_NEIL) +case(B_FRIERSON) ! longwave optical thickness function of latitude and pressure - if (do_tl) then - lat_tl = pi/2 - asin(cos(lat)*cos(lon)) - lw_tau_0 = ir_tau_eq + (ir_tau_pole - ir_tau_eq)*sin(lat_tl)**2 - where (lat_tl .ge. (pi/2)) - lw_tau_0 = ir_tau_pole - endwhere - else - lw_tau_0 = ir_tau_eq + (ir_tau_pole - ir_tau_eq)*sin(lat)**2 - endif + lw_tau_0 = ir_tau_eq + (ir_tau_pole - ir_tau_eq)*sin(lat)**2 lw_tau_0 = lw_tau_0 * odp ! scale by optical depth parameter - default 1 ! compute optical depths for each model level @@ -755,17 +608,9 @@ subroutine two_stream_gray_rad_down (is, js, Time_diag, lat, lon, p_full, p_half ! compute downward longwave flux by integrating downward lw_down(:,:,1) = 0. - if (do_normal_integration_method) then - do k = 1, n - lw_down(:,:,k+1) = lw_down(:,:,k)*lw_dtrans(:,:,k) + b(:,:,k)*(1. - lw_dtrans(:,:,k)) - end do - - else - do k = 1, n - lw_down(:,:,k+1) = 2.*b(:,:,k) * (lw_tau(:,:,k+1) - lw_tau(:,:,k))/(2.+ (lw_tau(:,:,k+1) - lw_tau(:,:,k))) + & - lw_down(:,:,k) * (2.- (lw_tau(:,:,k+1) - lw_tau(:,:,k)))/(2.+ (lw_tau(:,:,k+1) - lw_tau(:,:,k))) - end do - endif + do k = 1, n + lw_down(:,:,k+1) = lw_down(:,:,k)*lw_dtrans(:,:,k) + b(:,:,k)*(1. - lw_dtrans(:,:,k)) + end do case(B_SCHNEIDER_LIU) @@ -807,8 +652,6 @@ subroutine two_stream_gray_rad_down (is, js, Time_diag, lat, lon, p_full, p_half used = send_data ( id_lwdn_sfc, surf_lw_down, Time_diag) endif -!surf_lw_down = surf_lw_down + const_correct + nudge_out ! NTL UPDATE NUDGING - !------- incoming sw flux toa ------- if ( id_swdn_toa > 0 ) then used = send_data ( id_swdn_toa, toa_sw_in, Time_diag) @@ -830,14 +673,6 @@ subroutine two_stream_gray_rad_down (is, js, Time_diag, lat, lon, p_full, p_half used = send_data ( id_coszen, coszen, Time_diag) endif -if (id_tau > 0) then - used = send_data ( id_tau, lw_tau_0, Time_diag) -endif - -if (id_sw_tau > 0) then - used = send_data ( id_sw_tau, sw_tau, Time_diag) -endif - !------- carbon dioxide concentration ------------ if ( id_co2 > 0 ) then used = send_data ( id_co2, carbon_conc, Time_diag) @@ -879,20 +714,12 @@ subroutine two_stream_gray_rad_up (is, js, Time_diag, lat, p_half, t_surf, t, td end do lw_up = lw_up + lw_up_win -case(B_FRIERSON, B_BYRNE, B_DAVIS, B_NEIL) +case(B_FRIERSON) ! compute upward longwave flux by integrating upward lw_up(:,:,n+1) = b_surf - if (do_normal_integration_method) then - do k = n, 1, -1 - lw_up(:,:,k) = lw_up(:,:,k+1)*lw_dtrans(:,:,k) + b(:,:,k)*(1.0 - lw_dtrans(:,:,k)) - end do - - else - do k = n, 1, -1 - lw_up(:,:,k) = 2.*b(:,:,k) * -1.*(lw_tau(:,:,k) - lw_tau(:,:,k+1))/ (2.- (lw_tau(:,:,k) - lw_tau(:,:,k+1))) & - + lw_up(:,:,k+1)* (2. + (lw_tau(:,:,k) - lw_tau(:,:,k+1)))/(2. - (lw_tau(:,:,k) - lw_tau(:,:,k+1))) - end do - endif + do k = n, 1, -1 + lw_up(:,:,k) = lw_up(:,:,k+1)*lw_dtrans(:,:,k) + b(:,:,k)*(1.0 - lw_dtrans(:,:,k)) + end do case(B_SCHNEIDER_LIU) ! compute upward longwave flux by integrating upward @@ -1004,19 +831,6 @@ subroutine two_stream_gray_rad_end end subroutine two_stream_gray_rad_end -! =================================================================================== - -function erf_func(x) - - implicit none - - real, intent(in), dimension(:,:) :: x - real, dimension(size(x,1), size(x,2)) :: erf_func - - erf_func=(2/sqrt(pi))*(x/abs(x))*sqrt(1-exp(-(x**2)))*((sqrt(pi)/2)+ (31/200)*exp(-(x**2)) - (341/8000)*(exp(-2*(x**2)))) - -end function erf_func - ! ================================================================================== end module two_stream_gray_rad_mod diff --git a/src/atmos_spectral/driver/solo/idealized_moist_phys.F90 b/src/atmos_spectral/driver/solo/idealized_moist_phys.F90 index acb50da81..93a57c0ae 100644 --- a/src/atmos_spectral/driver/solo/idealized_moist_phys.F90 +++ b/src/atmos_spectral/driver/solo/idealized_moist_phys.F90 @@ -33,8 +33,6 @@ module idealized_moist_phys_mod use dry_convection_mod, only: dry_convection_init, dry_convection -use dry_adj_mod, only: dry_adj_init, dry_adj, dry_adj_end - use diag_manager_mod, only: register_diag_field, send_data use transforms_mod, only: get_grid_domain @@ -104,15 +102,13 @@ module idealized_moist_phys_mod SIMPLE_BETTS_CONV = 1, & FULL_BETTS_MILLER_CONV = 2,& DRY_CONV = 3, & - RAS_CONV = 4, & - DRYADJ_CONV = 5 + RAS_CONV = 4 integer :: r_conv_scheme = UNSET ! the selected convection scheme logical :: lwet_convection = .false. logical :: do_bm = .false. logical :: do_ras = .false. -logical :: do_dryadj = .false. ! Cloud options logical :: do_cloud_simple = .false. ! by default the cloud scheme is off. @@ -151,9 +147,7 @@ module idealized_moist_phys_mod real :: raw_bucket = 0.53 ! default raw coefficient for bucket depth LJJ ! end RG Add bucket -real :: alpha = 1.0 - -namelist / idealized_moist_phys_nml / turb, lwet_convection, do_bm, do_ras, do_dryadj, roughness_heat, & +namelist / idealized_moist_phys_nml / turb, lwet_convection, do_bm, do_ras, roughness_heat, & do_cloud_simple, & two_stream_gray, do_rrtm_radiation, do_damping,& mixed_layer_bc, do_simple, & @@ -377,8 +371,7 @@ subroutine idealized_moist_phys_init(Time, Time_step_in, nhum, rad_lon_2d, rad_l r_conv_scheme = NO_CONV lwet_convection = .false. do_bm = .false. - do_ras = .false. - do_dryadj = .false. + do_ras = .false. call error_mesg('idealized_moist_phys','No convective adjustment scheme used.', NOTE) else if(uppercase(trim(convection_scheme)) == 'SIMPLE_BETTS_MILLER') then @@ -387,7 +380,6 @@ subroutine idealized_moist_phys_init(Time, Time_step_in, nhum, rad_lon_2d, rad_l lwet_convection = .true. do_bm = .false. do_ras = .false. - do_dryadj = .false. else if(uppercase(trim(convection_scheme)) == 'FULL_BETTS_MILLER') then r_conv_scheme = FULL_BETTS_MILLER_CONV @@ -395,7 +387,6 @@ subroutine idealized_moist_phys_init(Time, Time_step_in, nhum, rad_lon_2d, rad_l do_bm = .true. lwet_convection = .false. do_ras = .false. - do_dryadj = .false. else if(uppercase(trim(convection_scheme)) == 'RAS') then r_conv_scheme = RAS_CONV @@ -403,23 +394,13 @@ subroutine idealized_moist_phys_init(Time, Time_step_in, nhum, rad_lon_2d, rad_l do_ras = .true. do_bm = .false. lwet_convection = .false. - do_dryadj = .false. else if(uppercase(trim(convection_scheme)) == 'DRY') then r_conv_scheme = DRY_CONV call error_mesg('idealized_moist_phys','Using dry convection scheme.', NOTE) lwet_convection = .false. do_bm = .false. - do_ras = .false. - do_dryadj = .false. - -else if(uppercase(trim(convection_scheme)) == 'DRYADJ') then - r_conv_scheme = DRYADJ_CONV - call error_mesg('idealized_moist_phys','Using AM3 dry adjustment convection scheme.', NOTE) - lwet_convection = .false. - do_bm = .false. - do_ras = .false. - do_dryadj = .true. + do_ras = .false. else if(uppercase(trim(convection_scheme)) == 'UNSET') then call error_mesg('idealized_moist_phys','determining convection scheme from flags', NOTE) @@ -431,10 +412,6 @@ subroutine idealized_moist_phys_init(Time, Time_step_in, nhum, rad_lon_2d, rad_l r_conv_scheme = FULL_BETTS_MILLER_CONV call error_mesg('idealized_moist_phys','Using Betts-Miller convection scheme.', NOTE) end if - if (do_dryadj) then - r_conv_scheme = DRYADJ_CONV - call error_mesg('idealized_moist_phys','Using AM3 dry adjustment convection scheme.',NOTE) - end if if (do_ras) then r_conv_scheme = RAS_CONV call error_mesg('idealized_moist_phys','Using relaxed Arakawa Schubert convection scheme.', NOTE) @@ -447,20 +424,11 @@ subroutine idealized_moist_phys_init(Time, Time_step_in, nhum, rad_lon_2d, rad_l if(lwet_convection .and. do_bm) & call error_mesg('idealized_moist_phys','lwet_convection and do_bm cannot both be .true.',FATAL) -if(do_dryadj .and. do_bm) & - call error_mesg('idealized_moist_phys','do_dryadj and do_bm cannot both be .true.',FATAL) - if(lwet_convection .and. do_ras) & call error_mesg('idealized_moist_phys','lwet_convection and do_ras cannot both be .true.',FATAL) - -if(lwet_convection .and. do_dryadj) & - call error_mesg('idealized_moist_phys','lwet_convection and do_dryadj cannot both be .true.',FATAL) if(do_bm .and. do_ras) & call error_mesg('idealized_moist_phys','do_bm and do_ras cannot both be .true.',FATAL) - -if(do_dryadj .and. do_ras) & - call error_mesg('idealized_moist_phys', 'do_dryadj and do_ras cannot both be .true.',FATAL) nsphum = nhum Time_step = Time_step_in @@ -765,9 +733,6 @@ subroutine idealized_moist_phys_init(Time, Time_step_in, nhum, rad_lon_2d, rad_l call ras_init (do_strat, axes,Time,tracers_in_ras) -case (DRYADJ_CONV) - call dry_adj_init(alpha) - end select !jp not sure why these diag_fields are fenced when condensation ones above are not... @@ -828,8 +793,6 @@ subroutine idealized_moist_phys_init(Time, Time_step_in, nhum, rad_lon_2d, rad_l id_rh = register_diag_field ( mod_name, 'rh', & axes(1:3), Time, 'relative humidity', 'percent') - - end subroutine idealized_moist_phys_init !================================================================================================================================= subroutine idealized_moist_phys(Time, p_half, p_full, z_half, z_full, ug, vg, tg, grid_tracers, & @@ -926,17 +889,6 @@ subroutine idealized_moist_phys(Time, p_half, p_full, z_half, z_full, ug, vg, tg if(id_cape > 0) used = send_data(id_cape, cape, Time) if(id_cin > 0) used = send_data(id_cin, cin, Time) -case(DRYADJ_CONV) - call dry_adj(tg(:, :, :, previous), & - p_full(:,:,:,previous), p_half(:,:,:,previous), & - conv_dt_tg, delta_t) - - tg_tmp = conv_dt_tg(:,:,:) + tg(:,:,:,previous) - conv_dt_tg = conv_dt_tg / delta_t - conv_dt_qg = 0.0 - qg_tmp = grid_tracers(:,:,:,previous,nsphum) - if(id_conv_dt_tg > 0) used = send_data(id_conv_dt_tg, conv_dt_tg, Time) - case(DRY_CONV) call dry_convection(Time, tg(:, :, :, previous), & p_full(:,:,:,previous), p_half(:,:,:,previous), & @@ -1052,11 +1004,11 @@ subroutine idealized_moist_phys(Time, p_half, p_full, z_half, z_full, ug, vg, tg call two_stream_gray_rad_down(is, js, Time, & rad_lat(:,:), & rad_lon(:,:), & - p_full(:,:,:,current), p_half(:,:,:,current), & + p_half(:,:,:,current), & tg(:,:,:,previous), & net_surf_sw_down(:,:), & surf_lw_down(:,:), albedo, & - grid_tracers(:,:,:,previous,nsphum), const_correct(:,:), nudge_out(:,:)) !NTL_CHANGE + grid_tracers(:,:,:,previous,nsphum), const_correct(:,:), nudge_out(:,:)) end if if(.not.mixed_layer_bc) then @@ -1402,7 +1354,6 @@ subroutine idealized_moist_phys_end if(two_stream_gray) call two_stream_gray_rad_end if(lwet_convection) call qe_moist_convection_end if(do_ras) call ras_end -if(do_dryadj) call dry_adj_end if(turb) then call vert_diff_end diff --git a/src/atmos_spectral/driver/solo/mixed_layer.F90 b/src/atmos_spectral/driver/solo/mixed_layer.F90 index ae7d96ab7..2d9e4df00 100644 --- a/src/atmos_spectral/driver/solo/mixed_layer.F90 +++ b/src/atmos_spectral/driver/solo/mixed_layer.F90 @@ -130,7 +130,7 @@ module mixed_layer_mod character(len=256) :: flux_lhe_anom_field_name = 'flux_lhe_anom' ! THERMO ICE OPTIONS - NTL 01/23 -logical :: do_explicit = .false. ! explicit or implicit timestepping for ice +logical :: do_explicit = .false. ! explicit or implicit timestepping for ocean covered surface logical :: do_thermo_ice = .false. ! do thermodynamic sea ice? real :: thermo_ice_albedo = 0.4 ! from Feldl & Merlis - original XZ is 0.8, M. England 0.6 real :: L_thermo_ice = DENS_ICE * HLF ! latent heat of fusion (J/m^3) @@ -138,19 +138,16 @@ module mixed_layer_mod real, parameter :: thermo_ice_basal_flux_const = 120 ! linear coefficient for heat flux from ocean to ice base (W/m^2/K) real :: t_thermo_ice_base = TFREEZE ! temperature at base of ice, taken to be freshwater freezing point real :: t_surf_freeze = TFREEZE -logical :: do_var_thermo_ice_albedo = .false. real :: thermo_ice_nudge_timescale = 7. * 86400. logical :: do_nudge_ice = .false. logical :: do_nudging_correction = .false. -logical :: correct_nh_only = .false. -logical :: diff_correction = .false. character(len=256) :: nudge_ice_file_name = 'nudge_ice' logical :: time_varying_nudge_ice = .false. logical :: do_prescribe_albedo = .false. character(len=256) :: albedo_file_name = 'albedo' logical :: time_varying_albedo = .false. logical :: prescribe_nh_albedo_only = .false. -logical :: read_const_correct = .true. +logical :: read_const_correct = .true. ! set to false if doing nudging, and want to restart from simulation with no nudging logical :: read_nudge_out = .true. logical :: albedo_from_nudge_file = .false. @@ -172,11 +169,10 @@ module mixed_layer_mod ice_concentration_threshold, & add_latent_heat_flux_anom,flux_lhe_anom_file_name,& flux_lhe_anom_field_name, do_ape_sst, qflux_field_name,& - do_explicit, do_thermo_ice, thermo_ice_albedo, t_surf_freeze, t_thermo_ice_base, & - do_var_thermo_ice_albedo,& ! ICE, NTL 01/23 + do_explicit, do_thermo_ice, thermo_ice_albedo, t_surf_freeze, t_thermo_ice_base, & ! NTL ice do_nudge_ice, nudge_ice_file_name, thermo_ice_nudge_timescale, time_varying_nudge_ice, & - do_nudging_correction, correct_nh_only, do_prescribe_albedo, time_varying_albedo, albedo_file_name, & - prescribe_nh_albedo_only, diff_correction, read_const_correct, read_nudge_out, albedo_from_nudge_file + do_nudging_correction, do_prescribe_albedo, time_varying_albedo, albedo_file_name, & + prescribe_nh_albedo_only, read_const_correct, read_nudge_out, albedo_from_nudge_file !================================================================================================================================= @@ -199,11 +195,6 @@ module mixed_layer_mod id_h_thermo_ice, & ! sea ice thickness id_t_ml, & ! mixed layer temperature id_flux_thermo_ice, & ! conductive heat flux through ice - id_ice_nudge_flux, & - id_ice_nudge_correction, & - id_ice_nudge_flux_avg, & - id_ice_nudge_correction_avg, & - id_ice_nudge_correction_chk !NTL_DELETE real, allocatable, dimension(:,:) :: & ocean_qflux, & ! Q-flux @@ -242,7 +233,6 @@ module mixed_layer_mod flux_thermo_ice, & ! Conductive heat flux through ice corrected_flux_noq real, allocatable, dimension(:,:) :: h_thermo_ice_in, nudge_flux, nudge_correction, albedo_in -real, allocatable, dimension(:,:) :: delta_t_surf_correct, delta_t_surf_nocorrect !NTL_DELETE real :: const_nudge_correction type(interpolate_type),save :: nudge_ice_interp type(interpolate_type),save :: albedo_interp @@ -351,8 +341,6 @@ subroutine mixed_layer_init(is, ie, js, je, num_levels, t_surf, & allocate(albedo_in (is:ie, js:je)) allocate(nudge_flux (is:ie, js:je)) allocate(nudge_correction (is:ie, js:je)) -allocate(delta_t_surf_nocorrect (is:ie, js:je)) !NTL_DELETE -allocate(delta_t_surf_correct (is:ie, js:je)) !NTL_DELETE ! !see if restart file exists for the surface temperature ! @@ -440,9 +428,6 @@ subroutine mixed_layer_init(is, ie, js, je, num_levels, t_surf, & albedo(:,:) = albedo_value where(land) albedo(:,:) = land_albedo_prefactor * albedo_value albedo_initial(:,:) = albedo (:,:) - !where(.not. land) - ! where (t_surf .lt. TFREEZE) albedo(:,:) = thermo_ice_albedo - !endwhere else if (do_prescribe_albedo) then albedo(:,:) = albedo_value where(land) albedo(:,:) = land_albedo_prefactor * albedo_value @@ -479,16 +464,6 @@ subroutine mixed_layer_init(is, ie, js, je, num_levels, t_surf, & axes(1:2), Time, 'mixed layer tempeature','K') id_flux_thermo_ice = register_diag_field(mod_name, 'flux_thermo_ice', & axes(1:2), Time, 'conductive heat flux through sea ice','watts/m2') - id_ice_nudge_flux = register_diag_field(mod_name, 'ice_nudge_flux', & - axes(1:2), Time, 'nudging heat flux','watts/m2') - id_ice_nudge_correction = register_diag_field(mod_name, 'ice_nudge_correction', & - axes(1:2), Time, 'nudging heat flux correction','watts/m2') - id_ice_nudge_correction_chk = register_diag_field(mod_name, 'ice_nudge_correction_chk', & - axes(1:2), Time, 'nudging heat flux correction check','watts/m2') !NTL_DELETE - id_ice_nudge_flux_avg = register_diag_field(mod_name, 'ice_nudge_flux_avg', & - Time, 'nudging heat flux','watts/m2') - id_ice_nudge_correction_avg = register_diag_field(mod_name, 'ice_nudge_correction_avg', & - Time, 'nudging heat flux correction','watts/m2') else id_ice_conc = register_diag_field(mod_name, 'ice_conc', & axes(1:2), Time, 'ice_concentration', 'none') @@ -529,7 +504,8 @@ subroutine mixed_layer_init(is, ie, js, je, num_levels, t_surf, & endif -h_thermo_ice_in = 999. ! large number stops this doing anything if do_nudge_ice is false +h_thermo_ice_in = 999. ! this is used to determine whether nudging should occur + ! setting a large number ensures this will do nothing if not initialised from external file ! load ice if (do_nudge_ice .or. albedo_from_nudge_file) then if (time_varying_nudge_ice .or. albedo_from_nudge_file) then @@ -876,22 +852,7 @@ subroutine mixed_layer ( & t_surf_dependence = beta_t * CP_AIR + beta_lw -!nudge_out = 0.0 -!where ((h_thermo_ice .gt. 0) .and. (h_thermo_ice .le. 0.1)) -!nudge_out = h_thermo_ice / (86400./4.) * L_thermo_ice -!endwhere -!corrected_flux = corrected_flux - nudge_out -!const_nudge_correction = 0.0 -!const_correct = 0.0 -!const_nudge_correction = -area_weighted_global_mean(nudge_out) -!if (.not. (const_nudge_correction .eq. 0.)) then -! where (nudge_out .eq. 0.) -! const_correct(:,:) = const_nudge_correction -! endwhere -! const_correct = const_correct * const_nudge_correction / area_weighted_global_mean(const_correct) -! !const_correct(:,:) = -area_weighted_global_mean(nudge_out) !-area_weighted_global_mean(nudge_flux * L_thermo_ice/dt) -! corrected_flux = corrected_flux - const_correct -!endif + @@ -914,17 +875,10 @@ subroutine mixed_layer ( & if (do_nudging_correction) then const_nudge_correction = -area_weighted_global_mean(nudge_out) if (.not. (const_nudge_correction .eq. 0.)) then - if (correct_nh_only) then - where ((nudge_out .eq. 0.) .and. (rad_lat_2d .gt. 0)) - const_correct(:,:) = const_nudge_correction - endwhere - else - where (nudge_out .eq. 0.) - const_correct(:,:) = const_nudge_correction - endwhere - endif + where (nudge_out .eq. 0.) + const_correct(:,:) = const_nudge_correction + endwhere const_correct = const_correct * const_nudge_correction / area_weighted_global_mean(const_correct) - !const_correct(:,:) = -area_weighted_global_mean(nudge_out) !-area_weighted_global_mean(nudge_flux * L_thermo_ice/dt) corrected_flux = corrected_flux - const_correct endif endif @@ -1025,7 +979,6 @@ subroutine mixed_layer ( & ! = update state = t_ml = t_ml + delta_t_ml h_thermo_ice = h_thermo_ice + delta_h_thermo_ice - !albedo_out(:,:) = albedo_value * land_albedo_prefactor elsewhere ! NTL 01/23 thermodynamic sea ice @@ -1064,10 +1017,6 @@ subroutine mixed_layer ( & elsewhere ! = ice-covered = ( h>0 ) delta_t_ml = - ( thermo_ice_basal_flux_const * ( t_ml - t_thermo_ice_base ) - ocean_qflux ) * dt / land_sea_heat_capacity delta_h_thermo_ice = ( corrected_flux_noq - thermo_ice_basal_flux_const * ( t_ml - t_thermo_ice_base ) ) * dt / L_thermo_ice - !where (h_thermo_ice_in .le. 0 ) ! by default h_thermo_ice_in is zero and this doesn't do anything - ! nudge_flux = - h_thermo_ice_in / thermo_ice_nudge_timescale * dt - ! delta_h_thermo_ice = delta_h_thermo_ice + nudge_flux - !endwhere endwhere where ( t_ml + delta_t_ml .lt. t_surf_freeze ) ! = frazil growth = @@ -1078,15 +1027,6 @@ subroutine mixed_layer ( & where ( ( h_thermo_ice .gt. 0 ) .and. ( h_thermo_ice + delta_h_thermo_ice .le. 0 ) ) ! = complete ablation = delta_t_ml = delta_t_ml - ( h_thermo_ice + delta_h_thermo_ice ) * L_thermo_ice/land_sea_heat_capacity delta_h_thermo_ice = - h_thermo_ice - - ! elsewhere ( ( h_thermo_ice .gt. 0 ) .and. ( h_thermo_ice + delta_h_thermo_ice - h_thermo_ice / thermo_ice_nudge_timescale * dt .le. 0 ) .and. (h_thermo_ice_in .le. 0 )) ! by default h_thermo_ice_in is 999 and this doesn't do anything - ! nudge_flux = (h_thermo_ice + delta_h_thermo_ice) - ! ! no need to update delta_t_ml as dT_ml = dT_ml - (h + dh + n_f), but this = dT_ml - ! delta_h_thermo_ice = - h_thermo_ice - ! elsewhere ( (h_thermo_ice .gt. 0) .and. ( h_thermo_ice + delta_h_thermo_ice - h_thermo_ice / thermo_ice_nudge_timescale * dt .gt. 0 ) .and. (h_thermo_ice_in .le. 0 )) ! by default h_thermo_ice_in is 999 and this doesn't do anything - ! nudge_flux = h_thermo_ice / thermo_ice_nudge_timescale * dt - ! delta_h_thermo_ice = delta_h_thermo_ice - nudge_flux - endwhere ! = update surface temperature = @@ -1117,125 +1057,17 @@ subroutine mixed_layer ( & endwhere - ! delta_t_surf_nocorrect = delta_t_surf !NTL_DELETE - - ! ! do nudging correction to ice free area - ! if (do_nudging_correction) then - ! if(.not. diff_correction) then - ! const_nudge_correction = -area_weighted_global_mean(nudge_flux * L_thermo_ice/land_sea_heat_capacity) - - ! if (.not. (const_nudge_correction .eq. 0.)) then - ! where ( h_thermo_ice .le. 0. ) - ! nudge_correction = const_nudge_correction - ! endwhere - - ! nudge_correction = nudge_correction * const_nudge_correction / area_weighted_global_mean(nudge_correction) - - ! where ( h_thermo_ice .le. 0. ) - ! delta_t_surf = delta_t_surf + nudge_correction - ! t_surf = t_surf + nudge_correction - ! delta_t_ml = delta_t_ml + nudge_correction !NTL_NEWCODE - ! t_ml = t_ml + nudge_correction !NTL_NEWCODE - ! endwhere - ! endif - ! else if (diff_correction) then - ! ! do nudging correction to all locs where there is no nudging - ! !if (do_nudging_correction) then - ! const_nudge_correction = -area_weighted_global_mean(nudge_flux * L_thermo_ice/land_sea_heat_capacity) - - ! if (.not. (const_nudge_correction .eq. 0.)) then - ! where ( nudge_flux .eq. 0. ) - ! nudge_correction = const_nudge_correction - ! endwhere - - ! nudge_correction = nudge_correction * const_nudge_correction / area_weighted_global_mean(nudge_correction) - - ! where ( nudge_flux .eq. 0. ) - ! where ( h_thermo_ice .le. 0. ) - ! delta_t_surf = delta_t_surf + nudge_correction - ! t_surf = t_surf + nudge_correction - ! delta_t_ml = delta_t_ml + nudge_correction !NTL_NEWCODE - ! t_ml = t_ml + nudge_correction !NTL_NEWCODE - ! elsewhere (h_thermo_ice .gt. 0.) - ! h_thermo_ice = h_thermo_ice - delta_h_thermo_ice - ! t_surf = t_surf - delta_t_thermo_ice - ! delta_h_thermo_ice = delta_h_thermo_ice - nudge_correction*depth*RHO_CP/L_thermo_ice - ! flux_thermo_ice = k_thermo_ice / (h_thermo_ice + delta_h_thermo_ice) * ( t_thermo_ice_base - t_surf ) - ! delta_t_thermo_ice = ( - corrected_flux_noq + flux_thermo_ice ) & - ! / ( k_thermo_ice / (h_thermo_ice + delta_h_thermo_ice) + dFdt_surf ) - ! where ( t_surf + delta_t_thermo_ice .gt. t_surf_freeze ) ! surface ablation - ! delta_t_thermo_ice = t_surf_freeze - t_surf - ! endwhere - ! ! surface is ice-covered, so update t_surf as ice surface temperature - ! delta_t_surf = delta_t_thermo_ice - ! t_surf = t_surf + delta_t_thermo_ice - ! h_thermo_ice = h_thermo_ice + delta_h_thermo_ice - ! endwhere - ! endwhere - ! endif - ! endif - ! endif - ! ! end of nudging correction - - ! elsewhere ( ( h_thermo_ice .gt. 0 ) .and. ( h_thermo_ice + delta_h_thermo_ice - h_thermo_ice / thermo_ice_nudge_timescale * dt .le. 0 ) .and. (h_thermo_ice_in .le. 0 )) ! by default h_thermo_ice_in is 999 and this doesn't do anything - ! nudge_flux = (h_thermo_ice + delta_h_thermo_ice) - ! ! no need to update delta_t_ml as dT_ml = dT_ml - (h + dh + n_f), but this = dT_ml - ! delta_h_thermo_ice = - h_thermo_ice - ! elsewhere ( (h_thermo_ice .gt. 0) .and. ( h_thermo_ice + delta_h_thermo_ice - h_thermo_ice / thermo_ice_nudge_timescale * dt .gt. 0 ) .and. (h_thermo_ice_in .le. 0 )) ! by default h_thermo_ice_in is 999 and this doesn't do anything - ! nudge_flux = h_thermo_ice / thermo_ice_nudge_timescale * dt - ! delta_h_thermo_ice = delta_h_thermo_ice - nudge_flux - - - - ! if (do_nudge_ice) then - ! where ((h_thermo_ice .gt. 0) .and. (h_thermo_ice_in .le. 0)) - ! nudge_out = h_thermo_ice / thermo_ice_nudge_timescale * L_thermo_ice - ! elsewhere - ! nudge_out = 0.0 - ! endwhere - ! endif - - - ! if (do_nudging_correction) then - ! const_nudge_correction = -area_weighted_global_mean(nudge_out) - ! where (nudge_out .eq. 0.) - ! const_correct(:,:) = const_nudge_correction - ! elsewhere - ! const_correct(:,:) = 0.0 - ! endwhere - ! const_correct = const_correct * const_nudge_correction / area_weighted_global_mean(const_correct) - ! !const_correct(:,:) = -area_weighted_global_mean(nudge_out) !-area_weighted_global_mean(nudge_flux * L_thermo_ice/dt) - ! else - ! const_correct(:,:) = 0.0 - ! endif - - - !delta_t_surf_correct = delta_t_surf - delta_t_surf_nocorrect !NTL_DELETE - - - ! compute albedo - if (do_var_thermo_ice_albedo) then - where (land_ice_mask) - albedo_out(:,:) = albedo_value * land_albedo_prefactor - elsewhere - where ( h_thermo_ice .gt. 0.0 ) - albedo_out(:,:) = albedo_value + (thermo_ice_albedo - albedo_value) * tanh((h_thermo_ice)/1.) - elsewhere - albedo_out(:,:) = albedo_value - endwhere - endwhere - else - where (land_ice_mask) - albedo_out(:,:) = albedo_value * land_albedo_prefactor + ! compute albedo + where (land_ice_mask) + albedo_out(:,:) = albedo_value * land_albedo_prefactor + elsewhere + where ( h_thermo_ice .gt. 0.0 ) + albedo_out(:,:) = thermo_ice_albedo elsewhere - where ( h_thermo_ice .gt. 0.0 ) - albedo_out(:,:) = thermo_ice_albedo - elsewhere - albedo_out(:,:) = albedo_value - endwhere - endwhere - endif + albedo_out(:,:) = albedo_value + endwhere + endwhere if (do_prescribe_albedo) then @@ -1288,11 +1120,6 @@ subroutine mixed_layer ( & if(id_h_thermo_ice > 0) used = send_data(id_h_thermo_ice, h_thermo_ice, Time_next) if(id_t_ml > 0) used = send_data(id_t_ml, t_ml, Time_next) if(id_flux_thermo_ice > 0) used = send_data(id_flux_thermo_ice, flux_thermo_ice, Time_next) -if(id_ice_nudge_flux > 0) used = send_data(id_ice_nudge_flux, L_thermo_ice * nudge_flux / dt, Time_next) -if(id_ice_nudge_correction > 0) used = send_data(id_ice_nudge_correction, land_sea_heat_capacity * nudge_correction / dt, Time_next) -if(id_ice_nudge_correction_chk > 0) used = send_data(id_ice_nudge_correction_chk, land_sea_heat_capacity * delta_t_surf_correct / dt, Time_next) ! NTL_DELETE -if(id_ice_nudge_flux_avg > 0) used = send_data(id_ice_nudge_flux_avg, area_weighted_global_mean(L_thermo_ice * nudge_flux / dt), Time_next) -if(id_ice_nudge_correction_avg > 0) used = send_data(id_ice_nudge_correction_avg, area_weighted_global_mean(land_sea_heat_capacity * nudge_correction / dt), Time_next) end subroutine mixed_layer diff --git a/src/atmos_spectral/model/spectral_dynamics.F90 b/src/atmos_spectral/model/spectral_dynamics.F90 index e7983c477..ab2e0b9ae 100644 --- a/src/atmos_spectral/model/spectral_dynamics.F90 +++ b/src/atmos_spectral/model/spectral_dynamics.F90 @@ -1568,7 +1568,7 @@ subroutine spectral_diagnostics_init(Time) real,dimension(2) :: vrange character(len=128) :: tname, longname, units -vrange = (/ -10000., 10000. /) +vrange = (/ -400., 400. /) rad_to_deg = 180./pi call get_grid_boundaries(lonb,latb,global=.true.) diff --git a/src/extra/model/grey/path_names b/src/extra/model/grey/path_names index 6d629a670..51c284f52 100644 --- a/src/extra/model/grey/path_names +++ b/src/extra/model/grey/path_names @@ -10,7 +10,6 @@ atmos_param/sea_esf_rad/null/rad_utilities.F90 atmos_param/shallow_conv/shallow_conv.F90 atmos_param/stable_bl_turb/stable_bl_turb.F90 atmos_param/strat_cloud/null/strat_cloud.F90 -atmos_param/dry_adj/dry_adj.F90 atmos_param/cloud_simple/cloud_simple.F90 atmos_param/two_stream_gray_rad/two_stream_gray_rad.F90 atmos_param/qflux/qflux.f90 diff --git a/src/extra/model/isca/path_names b/src/extra/model/isca/path_names index 287cf0c0c..eebdfd2b0 100644 --- a/src/extra/model/isca/path_names +++ b/src/extra/model/isca/path_names @@ -11,7 +11,6 @@ atmos_param/sea_esf_rad/null/rad_utilities.F90 atmos_param/shallow_conv/shallow_conv.F90 atmos_param/stable_bl_turb/stable_bl_turb.F90 atmos_param/strat_cloud/null/strat_cloud.F90 -atmos_param/dry_adj/dry_adj.F90 atmos_param/cloud_simple/cloud_simple.F90 atmos_param/two_stream_gray_rad/two_stream_gray_rad.F90 atmos_param/qflux/qflux.f90 @@ -20,7 +19,7 @@ atmos_param/monin_obukhov/monin_obukhov_kernel.F90 atmos_param/monin_obukhov/monin_obukhov.F90 atmos_param/dry_convection/dry_convection.f90 atmos_param/rayleigh_bottom_drag/rayleigh_bottom_drag.F90 -atmos_param/rrtm_radiation/rrtmg_lw/gcm_model/modules/parkind.F90 +atmos_param/rrtm_radiation/rrtmg_lw/gcm_model/modules/parkind.f90 atmos_param/rrtm_radiation/rrtmg_lw/gcm_model/modules/parrrtm.f90 atmos_param/rrtm_radiation/rrtmg_lw/gcm_model/modules/rrlw_cld.f90 atmos_param/rrtm_radiation/rrtmg_lw/gcm_model/modules/rrlw_con.f90 @@ -50,7 +49,7 @@ atmos_param/rrtm_radiation/rrtmg_lw/gcm_model/src/mcica_random_numbers.f90 atmos_param/rrtm_radiation/rrtmg_lw/gcm_model/src/rrtmg_lw_cldprop.f90 atmos_param/rrtm_radiation/rrtmg_lw/gcm_model/src/rrtmg_lw_rad.nomcica.f90 atmos_param/rrtm_radiation/rrtmg_lw/gcm_model/src/rrtmg_lw_rtrn.f90 -atmos_param/rrtm_radiation/rrtmg_lw/gcm_model/src/rrtmg_lw_setcoef.F90 +atmos_param/rrtm_radiation/rrtmg_lw/gcm_model/src/rrtmg_lw_setcoef.f90 atmos_param/rrtm_radiation/rrtmg_lw/gcm_model/src/mcica_subcol_gen_lw.f90 atmos_param/rrtm_radiation/rrtmg_lw/gcm_model/src/rrtmg_lw_init.f90 atmos_param/rrtm_radiation/rrtmg_lw/gcm_model/src/rrtmg_lw_rtrnmc.f90 @@ -58,7 +57,7 @@ atmos_param/rrtm_radiation/rrtmg_lw/gcm_model/src/rrtmg_lw_taumol.f90 atmos_param/rrtm_radiation/rrtmg_lw/gcm_model/src/rrtmg_lw_cldprmc.f90 atmos_param/rrtm_radiation/rrtmg_lw/gcm_model/src/rrtmg_lw_k_g.f90 atmos_param/rrtm_radiation/rrtmg_lw/gcm_model/src/rrtmg_lw_rtrnmr.f90 -atmos_param/rrtm_radiation/rrtmg_sw/gcm_model/modules/parkind.F90 +atmos_param/rrtm_radiation/rrtmg_sw/gcm_model/modules/parkind.f90 atmos_param/rrtm_radiation/rrtmg_sw/gcm_model/modules/parrrsw.f90 atmos_param/rrtm_radiation/rrtmg_sw/gcm_model/modules/rrsw_aer.f90 atmos_param/rrtm_radiation/rrtmg_sw/gcm_model/modules/rrsw_cld.f90 diff --git a/src/extra/model/socrates/path_names b/src/extra/model/socrates/path_names index 120c82696..c92ad053f 100644 --- a/src/extra/model/socrates/path_names +++ b/src/extra/model/socrates/path_names @@ -11,7 +11,6 @@ atmos_param/sea_esf_rad/null/rad_utilities.F90 atmos_param/shallow_conv/shallow_conv.F90 atmos_param/stable_bl_turb/stable_bl_turb.F90 atmos_param/strat_cloud/null/strat_cloud.F90 -atmos_param/dry_adj/dry_adj.F90 atmos_param/cloud_simple/cloud_simple.F90 atmos_param/two_stream_gray_rad/two_stream_gray_rad.F90 atmos_param/qflux/qflux.f90 diff --git a/src/extra/python/isca/codebase.py b/src/extra/python/isca/codebase.py index ada94b55b..09353e7ce 100644 --- a/src/extra/python/isca/codebase.py +++ b/src/extra/python/isca/codebase.py @@ -34,8 +34,7 @@ def from_repo(cls, repo, commit=None, **kwargs): def from_directory(cls, directory, **kwargs): return cls(directory=directory, **kwargs) - def __init__(self, repo=None, commit=None, directory=None, storedir=P(GFDL_WORK, 'codebase'), safe_mode=False, - compiler='intel', env_name=None): + def __init__(self, repo=None, commit=None, directory=None, storedir=P(GFDL_WORK, 'codebase'), safe_mode=False): """Create a new CodeBase object. A CodeBase can be created with either a git repository or a file directory as it's source. @@ -60,13 +59,9 @@ def __init__(self, repo=None, commit=None, directory=None, storedir=P(GFDL_WORK, self.log.error('Too many sources. Cannot create a CodeBase with both a source directory and a source repository.') raise AttributeError('Either repo= or directory= required to create CodeBase.') - self.compiler = compiler - self.env_name = env_name + self.safe_mode = safe_mode self.storedir = storedir - if compiler == 'gfortran': - self.executable_name = self.executable_name.strip('.x')+'_gfort.x' - #print('HERE!!', self.executable_name) if directory is not None: self.directory = directory @@ -125,13 +120,6 @@ def __init__(self, repo=None, commit=None, directory=None, storedir=P(GFDL_WORK, # read path names from the default file self.path_names = [] self.compile_flags = [] # users can append to this to add additional compiler options - self.precision_compile_flags = ['-DOVERLOAD_C8'] # default double precision - print(compiler, self.executable_name) - if self.compiler == 'intel': - self.other_compile_flags = ['-r8', '-i4'] # default double precision - elif self.compiler == 'gfortran': - self.other_compile_flags = ['-fdefault-real-8', '-fdefault-double-8'] - print('here', compiler, self.executable_name) @property def code_is_available(self): @@ -253,14 +241,10 @@ def _log_line(self, line): @useworkdir @destructive def compile(self, debug=False, optimisation=None): - if self.env_name is not None: - env = get_env_file(self.env_name) - else: - env = get_env_file() + env = get_env_file() mkdir(self.builddir) compile_flags = [] - other_flags = [] # if debug: # compile_flags.append('-g') # compile_flags.append('-traceback') @@ -270,11 +254,7 @@ def compile(self, debug=False, optimisation=None): # compile_flags.append('-O%d' % optimisation) compile_flags.extend(self.compile_flags) - compile_flags.extend(self.precision_compile_flags) compile_flags_str = ' '.join(compile_flags) - - other_flags.extend(self.other_compile_flags) - other_flags_str = ' '.join(other_flags) # get path_names from the directory if not self.path_names: @@ -288,7 +268,6 @@ def compile(self, debug=False, optimisation=None): 'srcdir': self.srcdir, 'workdir': self.workdir, 'compile_flags': compile_flags_str, - 'extra_compile_flags': other_flags_str, 'env_source': env, 'path_names': path_names_str, 'executable_name': self.executable_name, @@ -301,16 +280,6 @@ def compile(self, debug=False, optimisation=None): self._log_line(line) self.log.info('Compilation complete.') - - def use_single_precision(self): - self.log.info('Using single precision') - self.precision_compile_flags = ['-DOVERLOAD_C4', '-DOVERLOAD_R4'] - if self.compiler == 'intel': - self.other_compile_flags = ['-i4'] - elif self.compiler == 'gfortran': - self.other_compile_flags = ['-freal-8-real-4'] - self.executable_name = self.executable_name.strip('.x')+'_single.x' - self.builddir = P(self.workdir, 'build', self.executable_name.split('.')[0]) @@ -424,4 +393,4 @@ class DryCodeBase(GreyCodeBase): # """The Shallow Water Equations. # """ # name = 'shallow' -# executable_name = 'shallow.x' \ No newline at end of file +# executable_name = 'shallow.x' diff --git a/src/extra/python/isca/experiment.py b/src/extra/python/isca/experiment.py index 167078bda..d176033d5 100755 --- a/src/extra/python/isca/experiment.py +++ b/src/extra/python/isca/experiment.py @@ -74,10 +74,7 @@ def __init__(self, name, codebase, safe_mode=False, workbase=GFDL_WORK, database self.restartdir = P(self.datadir, 'restarts') # where restarts will be stored self.template_dir = P(_module_directory, 'templates') - if self.codebase.env_name is not None: - self.env_source = get_env_file(self.codebase.env_name) - else: - self.env_source = get_env_file() + self.env_source = get_env_file() self.templates = Environment(loader=FileSystemLoader(self.template_dir)) diff --git a/src/extra/python/isca/templates/compile.sh b/src/extra/python/isca/templates/compile.sh index 7811a643d..d20377492 100755 --- a/src/extra/python/isca/templates/compile.sh +++ b/src/extra/python/isca/templates/compile.sh @@ -18,9 +18,9 @@ template_debug={{ template_dir }}/mkmf.template.debug execdir={{ execdir }} # where code is compiled and executable is created executable={{ executable_name }} -#netcdf_flags=`nf-config --fflags --flibs` +netcdf_flags=`nf-config --fflags --flibs` -#ulimit -s unlimited # Set stack size to unlimited +ulimit -s unlimited # Set stack size to unlimited export MALLOC_CHECK_=0 # 3. compile the mppncombine tool if it hasn't yet been done. @@ -53,16 +53,14 @@ if [ $debug == True ]; then echo "Compiling in debug mode" # execute mkmf to create makefile -cppDefs="-Duse_libMPI -Duse_netCDF -Duse_LARGEFILE -DINTERNAL_FILE_NML {{compile_flags}}" -othDefs="{{extra_compile_flags}}" -$mkmf -a $sourcedir -t $template_debug -p $executable -c "$cppDefs" -o "$othDefs" $pathnames $sourcedir/shared/include $sourcedir/shared/mpp/include +cppDefs="-Duse_libMPI -Duse_netCDF -Duse_LARGEFILE -DINTERNAL_FILE_NML -DOVERLOAD_C8 {{compile_flags}}" +$mkmf -a $sourcedir -t $template_debug -p $executable -c "$cppDefs" $pathnames $sourcedir/shared/include $sourcedir/shared/mpp/include else # execute mkmf to create makefile -cppDefs="-Duse_libMPI -Duse_netCDF -Duse_LARGEFILE -DINTERNAL_FILE_NML ${CDEFS} {{compile_flags}}" -othDefs="{{extra_compile_flags}}" -$mkmf -a $sourcedir -t $template -p $executable -c "$cppDefs" -o "$othDefs" $pathnames $sourcedir/shared/include $sourcedir/shared/mpp/include +cppDefs="-Duse_libMPI -Duse_netCDF -Duse_LARGEFILE -DINTERNAL_FILE_NML -DOVERLOAD_C8 {{compile_flags}}" +$mkmf -a $sourcedir -t $template -p $executable -c "$cppDefs" $pathnames $sourcedir/shared/include $sourcedir/shared/mpp/include fi @@ -82,4 +80,4 @@ make $executable if [ $? != 0 ]; then echo "ERROR: make failed for $executable" exit 1 -fi \ No newline at end of file +fi diff --git a/src/extra/python/isca/templates/mkmf.template.ia64 b/src/extra/python/isca/templates/mkmf.template.ia64 index e3943e08e..dcb29e579 100644 --- a/src/extra/python/isca/templates/mkmf.template.ia64 +++ b/src/extra/python/isca/templates/mkmf.template.ia64 @@ -18,11 +18,11 @@ NETCDF_LIBS = `nc-config --libs` # -diag-disable 6843: # This suppresses the warning: `warning #6843: A dummy argument with an explicit INTENT(OUT) declaration is not given an explicit value.` of which # there are a lot of instances in the GFDL codebase. -FFLAGS = $(CPPFLAGS) -fpp -stack_temps -safe_cray_ptr -ftz -assume byterecl -shared-intel -g -O2 -diag-disable 6843 -mcmodel large +FFLAGS = $(CPPFLAGS) -fpp -stack_temps -safe_cray_ptr -ftz -assume byterecl -shared-intel -i4 -r8 -g -O2 -diag-disable 6843 -mcmodel large #FFLAGS = $(CPPFLAGS) -fltconsistency -stack_temps -safe_cray_ptr -ftz -shared-intel -assume byterecl -g -O0 -i4 -r8 -check -warn -warn noerrors -debug variable_locations -inline_debug_info -traceback FC = $(F90) LD = $(F90) $(NETCDF_LIBS) #CC = mpicc LDFLAGS = -lnetcdff -lnetcdf -lmpi -shared-intel -CFLAGS = -D__IFC \ No newline at end of file +CFLAGS = -D__IFC From cbdc1c9f6de27cfad52c03fe74e6040e37f135ac Mon Sep 17 00:00:00 2001 From: ntlewis Date: Fri, 31 May 2024 13:14:36 +0100 Subject: [PATCH 14/19] thermodynamic sea ice test case --- .../grey_thermo_ice_test_case.py | 238 ++++++++++++++++++ 1 file changed, 238 insertions(+) create mode 100644 exp/test_cases/thermodynamic_sea_ice/grey_thermo_ice_test_case.py diff --git a/exp/test_cases/thermodynamic_sea_ice/grey_thermo_ice_test_case.py b/exp/test_cases/thermodynamic_sea_ice/grey_thermo_ice_test_case.py new file mode 100644 index 000000000..0e331ad48 --- /dev/null +++ b/exp/test_cases/thermodynamic_sea_ice/grey_thermo_ice_test_case.py @@ -0,0 +1,238 @@ +import os + +import numpy as np + +from isca import GreyCodeBase, DiagTable, Experiment, Namelist, GFDL_BASE + +NCORES =16 + +# Point to code as defined by $GFDL_BASE +cb = GreyCodeBase.from_directory(GFDL_BASE) + +base_dir = os.path.dirname(os.path.realpath(__file__)) + +cb.compile() # compile the source code to working directory $GFDL_WORK/codebase + +# create an Experiment object to handle the configuration of model parameters +# and output diagnostics +exp = Experiment('thermo_ice_test_experiment', codebase=cb) + + +#Tell model how to write diagnostics +diag = DiagTable() +diag.add_file('atmos_monthly', 30, 'days', time_units='days') + +#Tell model which diagnostics to write +diag.add_field('dynamics', 'ps', time_avg=True) +diag.add_field('dynamics', 'bk') +diag.add_field('dynamics', 'pk') +diag.add_field('atmosphere', 'temp_2m', time_avg=True) +diag.add_field('dynamics', 'ucomp', time_avg=True) +diag.add_field('dynamics', 'temp', time_avg=True) +diag.add_field('mixed_layer', 'albedo', time_avg=True) +diag.add_field('mixed_layer', 'h_thermo_ice', time_avg=True) +diag.add_field('mixed_layer', 't_surf', time_avg=True) +diag.add_field('mixed_layer', 'flux_oceanq', time_avg=True) +diag.add_field('two_stream', 'swdn_toa', time_avg=True) + + + +exp.diag_table = diag # register diag table + + +#Empty the run directory ready to run +exp.clear_rundir() + +#Define values for the 'core' namelist +exp.namelist = namelist = Namelist({ + 'main_nml': { + 'days' : 360, # each run lasts one year, and then multiple runs are strung together below (loop on e.g. Line 310) + 'hours' : 0, # a different output file is produced for each run (in this case, each year). Data + 'minutes': 0, # output at the frequency specified in the diag table + 'seconds': 0, + 'dt_atmos':900, + 'current_date' : [1,1,1,0,0,0], + 'calendar' : 'thirty_day' + }, + + 'idealized_moist_phys_nml': { + 'two_stream_gray': True, + 'do_rrtm_radiation': False, + 'convection_scheme': 'SIMPLE_BETTS_MILLER', + 'do_damping': True, + 'turb':True, + 'mixed_layer_bc':True, + 'do_virtual' :True, + 'roughness_mom':5.e-3, + 'roughness_heat':1.e-5, + 'roughness_moist':1.e-5, + 'land_roughness_prefactor':1.0, + 'do_simple':False, + }, + + 'vert_turb_driver_nml': { + 'do_mellor_yamada': False, + 'do_diffusivity': True, + 'do_edt':False, + 'constant_gust': 1.0, + 'use_tau': False, + 'do_entrain':False, + 'do_stable_bl':False, + 'do_shallow_conv':False, + 'do_simple':False + }, + + 'diffusivity_nml': { + 'do_entrain':False, + 'entr_ratio': 0.0, + 'free_atm_diff':False, + 'do_simple': False, + 'parcel_buoy': 0.0, + 'frac_inner': 0.1, + 'fixed_depth': False, + }, + + 'surface_flux_nml': { + 'use_virtual_temp': True, + 'do_simple': False, + 'old_dtaudv': False, + 'gust_const':1.0, + 'land_humidity_prefactor' : 1.0, + 'land_evap_prefactor': 1.0, + }, + + 'atmosphere_nml': { + 'idealized_moist_model': True + }, + + + 'mixed_layer_nml': { + 'depth': 30., + 'albedo_value': 0.1, + 'prescribe_initial_dist':True, + 'tconst' : 305., + 'delta_T': 60., + 'evaporation':True, + 'do_qflux': True, + 'load_qflux':False, + 'time_varying_qflux' : False, + 'do_thermo_ice':True, # turn on thermodynamic ice + 'thermo_ice_albedo':0.55, # ice albedo + 't_thermo_ice_base':273.15-2., # freezing point + 't_surf_freeze':273.15-2., # freezing point + 'do_var_thermo_ice_albedo':False, + 'read_const_correct':False, + 'read_nudge_out':False, + + }, + + 'qflux_nml':{ + 'qflux_amp':30., + }, + + + + 'qe_moist_convection_nml': { + 'rhbm':0.7, + 'tau_bm':7200., + 'Tmin':120., + 'Tmax':360., + 'val_inc': 0.01, + 'precision':1.e-6 + }, + + 'lscale_cond_nml': { + 'do_simple':False, + 'do_evap':False + }, + + 'sat_vapor_pres_nml': { + 'do_simple':False, + + }, + + + + 'damping_driver_nml': { + 'do_rayleigh': True, + 'trayfric': -0.5, # neg. value: time in *days* + 'sponge_pbottom': 2600., + 'do_conserve_energy': True, + }, + + + + + 'two_stream_gray_rad_nml': { + 'rad_scheme': 'frierson', #Select radiation scheme to use + 'do_seasonal': True, + 'use_time_average_coszen':True, + 'dt_rad_avg':86400., + 'atm_abs': 0.22, + 'solar_exponent':2, + 'ir_tau_eq':7.2, + 'ir_tau_pole':3.6,#1.8, + 'del_sol':0.98, + 'solar_constant':1360, + 'linear_tau':0.2, + 'odp':1.4, + 'do_toa_albedo':True, + }, + + + 'spectral_dynamics_nml': { + 'damping_order': 4, # Yields lap^8 damping + 'water_correction_limit': 200.e2, + 'reference_sea_level_press':1.0e5, + 'num_levels':30, + 'valid_range_t':[100.,800.], + 'initial_sphum':[2.e-6], + 'use_virtual_temperature':True, + 'vert_coord_option':'uneven_sigma', + 'robert_coeff':0.03, + # set to T42 resolution + 'lon_max': 128, + 'lat_max': 64, + 'num_fourier': 42, + 'num_spherical':43, + 'surf_res': 0.05, + 'exponent': 3., + 'scale_heights': 5 + + }, + + 'spectral_init_cond_nml':{ + 'initial_temperature':280. + }, + + + + 'diag_manager_nml': { + 'mix_snapshot_average_fields': False # time avg fields are labelled with time in middle of window + }, + + 'fms_nml': { + 'domains_stack_size': 600000 # default: 0 + }, + + 'fms_io_nml': { + 'threading_write': 'single', # default: multi + 'fileset_write': 'single', # default: multi + }, + + + + +}) + +# Lets do a run! +if __name__=="__main__": + + + exp.run(1, use_restart=False, num_cores=NCORES, overwrite_data=False) + + for i in range(1,11): # run for 10 years + exp.run(i, num_cores=NCORES, overwrite_data=False) + + + From 0d3284e75da9162d9db2f9de13ca276c9aac13f3 Mon Sep 17 00:00:00 2001 From: Neil Lewis Date: Fri, 31 May 2024 13:39:55 +0100 Subject: [PATCH 15/19] re-add B_BYRNE (deleted by mistake) --- src/atmos_param/two_stream_gray_rad/two_stream_gray_rad.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/atmos_param/two_stream_gray_rad/two_stream_gray_rad.F90 b/src/atmos_param/two_stream_gray_rad/two_stream_gray_rad.F90 index a9a8d1615..efea06f25 100644 --- a/src/atmos_param/two_stream_gray_rad/two_stream_gray_rad.F90 +++ b/src/atmos_param/two_stream_gray_rad/two_stream_gray_rad.F90 @@ -714,7 +714,7 @@ subroutine two_stream_gray_rad_up (is, js, Time_diag, lat, p_half, t_surf, t, td end do lw_up = lw_up + lw_up_win -case(B_FRIERSON) +case(B_FRIERSON, B_BYRNE) ! compute upward longwave flux by integrating upward lw_up(:,:,n+1) = b_surf do k = n, 1, -1 From e2f9b0dd0bb032498989027ddc8c0e39a4a44961 Mon Sep 17 00:00:00 2001 From: Neil Lewis Date: Fri, 31 May 2024 13:42:32 +0100 Subject: [PATCH 16/19] Delete mkmf.template.ubuntu_conda2 --- .../templates/mkmf.template.ubuntu_conda2 | 28 ------------------- 1 file changed, 28 deletions(-) delete mode 100644 src/extra/python/isca/templates/mkmf.template.ubuntu_conda2 diff --git a/src/extra/python/isca/templates/mkmf.template.ubuntu_conda2 b/src/extra/python/isca/templates/mkmf.template.ubuntu_conda2 deleted file mode 100644 index 65d8eeb68..000000000 --- a/src/extra/python/isca/templates/mkmf.template.ubuntu_conda2 +++ /dev/null @@ -1,28 +0,0 @@ -# template for the gfortran compiler -# typical use with mkmf -# mkmf -t template.ifc -c"-Duse_libMPI -Duse_netCDF" path_names /usr/local/include -CPPFLAGS = `nc-config --cflags` -NC_INC=`nc-config --fflags` -NC_LIB=`nc-config --flibs` - -# FFLAGS: -# -cpp: Use the fortran preprocessor -# -ffree-line-length-none -fno-range-check: Allow arbitrarily long lines -# -fcray-pointer: Cray pointers don't alias other variables. -# -ftz: Denormal numbers are flushed to zero. -# -assume byterecl: Specifies the units for the OPEN statement as bytes. -# -shared-intel: Load intel libraries dynamically -# -i4: 4 byte integers -# -fdefault-real-8: 8 byte reals (compatability for some parts of GFDL code) -# -fdefault-double-8: 8 byte doubles (compat. with RRTM) -# -O2: Level 2 speed optimisations - -FFLAGS = $(CPPFLAGS) $(NC_LIB) -cpp -fcray-pointer \ - -O2 -ffree-line-length-none -fno-range-check \ - -fallow-invalid-boz -fallow-argument-mismatch - -FC = $(F90) -LD = $(F90) - -LDFLAGS = -lnetcdff -lnetcdf -lmpi -CFLAGS = -D__IFC From a898b084267a5e597fcfecf87477b21e0469aa89 Mon Sep 17 00:00:00 2001 From: Neil Lewis Date: Fri, 31 May 2024 13:50:44 +0100 Subject: [PATCH 17/19] fix mixed_layer.F90 syntax error --- src/atmos_spectral/driver/solo/mixed_layer.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/atmos_spectral/driver/solo/mixed_layer.F90 b/src/atmos_spectral/driver/solo/mixed_layer.F90 index 0b1efac1f..55909ee14 100644 --- a/src/atmos_spectral/driver/solo/mixed_layer.F90 +++ b/src/atmos_spectral/driver/solo/mixed_layer.F90 @@ -202,7 +202,7 @@ module mixed_layer_mod ! NTL 01/23 thermodynamic sea ice id_h_thermo_ice, & ! sea ice thickness id_t_ml, & ! mixed layer temperature - id_flux_thermo_ice, & ! conductive heat flux through ice + id_flux_thermo_ice & ! conductive heat flux through ice real, allocatable, dimension(:,:) :: & ocean_qflux, & ! Q-flux From 66268d18f4c905866eb899e1871935242835b0bc Mon Sep 17 00:00:00 2001 From: Neil Lewis Date: Fri, 31 May 2024 14:01:00 +0100 Subject: [PATCH 18/19] try to fix mixed_layer syntax error again --- src/atmos_spectral/driver/solo/mixed_layer.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/atmos_spectral/driver/solo/mixed_layer.F90 b/src/atmos_spectral/driver/solo/mixed_layer.F90 index 55909ee14..ee3df2717 100644 --- a/src/atmos_spectral/driver/solo/mixed_layer.F90 +++ b/src/atmos_spectral/driver/solo/mixed_layer.F90 @@ -200,9 +200,9 @@ module mixed_layer_mod id_ice_conc, & ! st ice concentration id_delta_t_surf, & ! NTL 01/23 thermodynamic sea ice - id_h_thermo_ice, & ! sea ice thickness + id_h_thermo_ice, & ! sea ice thickness id_t_ml, & ! mixed layer temperature - id_flux_thermo_ice & ! conductive heat flux through ice + id_flux_thermo_ice ! conductive heat flux through ice real, allocatable, dimension(:,:) :: & ocean_qflux, & ! Q-flux From a8146a5a0e336f76e078517707922e0759a469b9 Mon Sep 17 00:00:00 2001 From: Neil Lewis Date: Fri, 31 May 2024 14:07:14 +0100 Subject: [PATCH 19/19] deleting code preventing compile that evaded initial tidy --- src/atmos_spectral/driver/solo/mixed_layer.F90 | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/atmos_spectral/driver/solo/mixed_layer.F90 b/src/atmos_spectral/driver/solo/mixed_layer.F90 index ee3df2717..e9ad6c90d 100644 --- a/src/atmos_spectral/driver/solo/mixed_layer.F90 +++ b/src/atmos_spectral/driver/solo/mixed_layer.F90 @@ -908,8 +908,6 @@ subroutine mixed_layer ( & const_nudge_correction = 0. ! for calculation of ice sfc temperature dFdt_surf = dhdt_surf + drdt_surf + HLV * (dedt_surf) ! d(corrected_flux)/d(T_surf) - delta_t_surf_correct = 0. !NTL_DELETE - delta_t_surf_nocorrect = 0. !NTL_DELETE endif !s Surface heat_capacity calculation based on that in MiMA by mj