This directory contains a setup for running Aeolus linked to MOM5, ice-sis, land_lad . However, to test and tune Aeolus in the coupled setup, without running the time consuming ocean component, the model is run with prescribed ocean surface conditions, and the ocean and seaice submodels turned off. To start, we copy the configuration from the latest successful coupled run mom4_Aeolus_lad-init-21-mar-clim/run-107 . In input.nml, we set &coupler_nml/do_ocean = .false. ---- outdated, see below --- !---- re-animated, see Thu Mar 10 14:04:31 CET 2016 --- !As regarded by the FMS coupler, the entire ocean is covered by the sea-ice !model. Thus all quantities from the ocean must go through the ice model, and !are then passed on to the coupler. As far as I understand currently, the !coupler expects from the sea-ice these variables: ! !ICE t_surf !ICE rough_mom !ICE rough_heat !ICE rough_moist !ICE albedo # ice albedo for sea-ice-covered areas, 0 for open sea, missing_value for land !ICE albedo_vis_dir # ice albedo for sea-ice-covered areas, ocean albedo for open sea !ICE albedo_nir_dir # ice albedo for sea-ice-covered areas, ocean albedo for open sea !ICE albedo_vis_dif # ice albedo for sea-ice-covered areas, ocean albedo for open sea !ICE albedo_nir_dif # ice albedo for sea-ice-covered areas, ocean albedo for open sea !ICE u_surf !ICE v_surf ! !Also evaporation, which is configured as ``generic tracer'' . ! !Humidity is not delivered directly. !The humidity flux ex_flux_tr(:,sphum) can be overridden from the data_table. !That is calculated in ! surface_flux(... , ! 2=q_atm_in=ex_tr_atm(:,isphum), intent(in) ! 10=q_surf =ex_tr_surf(:,isphum), intent(inout) ! 19=flux_q =ex_flux_tr(:,isphum), intent(out) ! ...) ! !q_surf is taken from land , over ocean it is assumed to be saturation humidity. !surface_flux.F90 ca. line 410: ! ! initilaize surface air humidity according to surface type ! where (land) ! q_surf0 = q_surf ! land calculates it ! elsewhere ! q_surf0 = q_sat ! everything else assumes saturated sfc humidity ! endwhere ! !t_surf, u_surf, v_surf we can most probably take from some ERA dataset. !roughness I dont know. !To obtain usable albedo data appeared somewhat difficult recently, and still !unresolved, in the context of LPJ coupling. !For a first try, I would just record values from a CM2M_coarse_BLING run, and !use that as forcing values. ! ! !First approach: take output from a CM2M_coarse_BLING run and configure to read that in the data_table. !Create a separate subdirectory for those files, so that they are not confused with the ``usual'' INPUT/ files. ! !-------------------- !Mon May 23 12:21:06 CEST 2016 ! !CM2M_coarse_BLING ran for 450 years now. !Created a monthly climatology over years 250-450 of the ocean output. ! !Sigh. it is difficult to get the data for above overrides at the ice-> coupler boundary. !ICE t_surf "TS", "ice_month" !ICE rough_mom "rough_mom", "flux_month" on atmos grid! !ICE rough_heat "rough_heat", "flux_month" on atmos grid! !ICE rough_moist "rough_moist", "flux_month" on atmos grid! !ICE albedo "ALB", "ice_month" !ICE albedo_vis_dir "alb_vis_dir","ice_month" !ICE albedo_nir_dir "alb_nir_dir","ice_month" !ICE albedo_vis_dif "alb_vis_dif","ice_month" !ICE albedo_nir_dif "alb_nir_dif","ice_month" !ICE u_surf "u_surf","ocean_month" ? !ICE v_surf "v_surf","ocean_month" ? ! ---- end outdated --- It might be easier to do the data override at the ocean->ice boundary, which is computed in flux_exchange.F90::flux_ocean_to_ice() In there only these quantities are required from the ocean: call data_override('ICE', 'u', Ocean_Ice_Boundary%u, Time) call data_override('ICE', 'v', Ocean_Ice_Boundary%v, Time) call data_override('ICE', 't', Ocean_Ice_Boundary%t, Time) call data_override('ICE', 's', Ocean_Ice_Boundary%s, Time) call data_override('ICE', 'frazil', Ocean_Ice_Boundary%frazil, Time) call data_override('ICE', 'sea_level', Ocean_Ice_Boundary%sea_level, Time) do n = 1, Ocean_Ice_Boundary%fields%num_bcs !{ do m = 1, Ocean_Ice_Boundary%fields%bc(n)%num_fields !{ call data_override('ICE', Ocean_Ice_Boundary%fields%bc(n)%field(m)%name, & Ocean_Ice_Boundary%fields%bc(n)%field(m)%values, Time) enddo !} m enddo !} n real, pointer, dimension(:,:) :: t_surf =>NULL() ! SST on t-cell (degrees Kelvin) real, pointer, dimension(:,:) :: s_surf =>NULL() ! SSS on t-cell (psu) real, pointer, dimension(:,:) :: u_surf =>NULL() ! i-directed surface ocean velocity on u-cell (m/s) real, pointer, dimension(:,:) :: v_surf =>NULL() ! j-directed surface ocean velocity on u-cell (m/s) real, pointer, dimension(:,:) :: sea_lev =>NULL() ! eta_t + patm/(rho0*grav) - eta_geoid - eta_tide (m) real, pointer, dimension(:,:) :: frazil =>NULL() ! accumulated heating (J/m^2) from ! frazil formation in the ocean ICE u ICE v ICE t ICE s ICE frazil ICE sea_level ICE c14o2_flux_alpha_ocn_ice ICE c14o2_flux_csurf_ocn_ice ICE co2_b_flux_alpha_ocn_ice ICE co2_b_flux_csurf_ocn_ice ICE o2_b_flux_alpha_ocn_ice ICE o2_b_flux_csurf_ocn_ice Sigh. frazil is treated in the ocean model as a generic tracer. What is passed to the ice model is: ocean_public_type%frazil ! accumulated heating (J/m^2) from frazil formation in the ocean It is calculated in ocean_sbc.F90::sum_ocean_sfc() from the diagnostic tracer values: if(index_frazil > 0) then do k=1,nk do j = jsc_bnd,jec_bnd do i = isc_bnd,iec_bnd ii = i + i_shift jj = j + j_shift Ocean_sfc%frazil(i,j) = Ocean_sfc%frazil(i,j) + T_diag(index_frazil)%field(ii,jj,k) enddo enddo enddo endif The summed-up values are then finally time-averaged in ocean_sbc.F90::avg_ocean_sfc(). There some namelist options to choose which qunatities are actually time-averaged, summed up over the ocean-internal sub-timesteps, or taken as instantaneous values. Unfortunately, this quantity is not sent as diagnostic output. The tracer module puts out: ocean_model|frazil|frazil heating|J/m^2|3|T| -1.000000000000000E+020|||xt_ocean,yt_ocean,st_ocean ocean_model|surface_frazil|frazil heating|J/m^2|2|T| -1.000000000000000E+020| -10.0000000000000| 100.000000000000|xt_ocean,yt_ocean ocean_model|frazil_2d|ocn frazil heat flux over time step|W/m^2|2|T| -1.000000000000000E+020| -10000000000.0000| 10000000000.0000|xt_ocean,yt_ocean ocean_model|frazil_3d|ocn frazil heat flux over time step|W/m^2|3|T| -1.000000000000000E+020| -10000000000.0000| 10000000000.0000|xt_ocean,yt_ocean,st_ocean From those, frazil_2d seems to match the ice_model|FRAZIL output most colsely: ice_model|FRAZIL|energy flux of frazil formation|W/m^2|2|T| -9.999999999999999E+033|||xt,yt However, the above-shown ocean_public_type%frazil with units [J/m^2] is matched most closely by surface_frazil. Thus we use for the data override: ICE u "ocean_model","usurf","u_surf","ocean_month" ICE v "ocean_model","vsurf","v_surf","ocean_month" ICE t "ocean_model","surface_temp","sst","ocean_month" ICE s "ocean_model","surface_salt","sss","ocean_month" ICE frazil "ocean_model","surface_frazil","surface_frazil","ocean_month" ICE sea_level "ocean_model","sea_level","sea_level","ocean_month" ICE c14o2_flux_alpha_ocn_ice 0 ICE c14o2_flux_csurf_ocn_ice 0 ICE co2_b_flux_alpha_ocn_ice 0 ICE co2_b_flux_csurf_ocn_ice 0 ICE o2_b_flux_alpha_ocn_ice 0 ICE o2_b_flux_csurf_ocn_ice 0 ------------------ Wed May 25 11:27:20 CEST 2016 To remove doubts, we add diagnostic output to ocean_sbc.F90, sending the actual public data to the diagnostics manager at the end of avg_ocean_sfc() ------------ Fri May 27 15:09:15 CEST 2016 Finally, the CM2M_coarse_BLING model produced years 0452 .. 0502 of monthly data with the 6 new ``public_'' variables as diagnostik output. It turns out that public_frazil == surface_frazil public_sea_lev ~= sea_level (differences on the order of rounding errors?) public_s_surf ~= sss (differences on the order of rounding errors?) public_u_surf ~= u_surf (differences on the order of rounding errors?) public_v_surf ~= v_surf (differences on the order of rounding errors?) public_t_surf might be very similar to t_surf, but must be converted from K to C over ocean, outside the land mask. cdo ymonmean -selname\,public\* 05020101.ocean_month.nc 05020101.ocean_month-publicvars-climatology.nc Sigh. this clobbers the axis information from the MOM diag manager, so that the time_interp_external module cannot read the resulting file lateron. Try to do it with nco. Something like for i in `seq -w 0 11` do ncks -a -O -d time\,$i\,\,12 -v public_\*\,geol\* 05020101.ocean_month.nc o-m-$i.nc rm -f o-m-$i-clim.nc ncra -O o-m-$i.nc o-m-$i-clim.nc done ncrcat -O o-m-$i-clim.nc 05020101.ocean_month-publicvars-climatology.nc Another sigh. The spatial axes look OK, but the time axis does not. For now, just lets use the un-averaged data file. ncks -a -O -v public_\*\,geol\* ../../work/CM2M_coarse_BLING/history-hlrs15-vec-precise-1/05020101.ocean_month.nc 05020101.ocean_month-publicvars.nc cd /p/projects/climber3/petri/POEM/work mkdir Aeolus_lad-init-21-mar-clim cd mom4_Aeolus_lad-init-21-mar-clim/ cp -av INPUT-init-21-mar-clim-RESTART-00520401 ../Aeolus_lad-init-21-mar-clim/ cd run-107 cp -av diag_table* input.nml submit-seq.slurm ../../Aeolus_lad-init-21-mar-clim cd ../../../Aeolus_lad-init-21-mar-clim/ cp -av ../../exp/Aeolus_lad-init-21-mar-clim/data_table . cp -av ../../exp/Aeolus_lad-init-21-mar-clim/field_table . ------- The data override module complains about mismatching dates. Lets try to declare it as climatology: ncatted -a modulo\,time\,c\,c\,' ' 05020101.ocean_month-publicvars.nc Still does not work: FATAL: time_interp_external 2: period of list exceeds modulo period,file=INPUT/05020101.ocean_month-publicvars.nc,field=public_u_surf Next try: let CM2M_coarse_BLING run from 52-01-01 to 53-01-01 with 12 month output. ncks -a -O -v public_\*\,geol\* history/00530101.ocean_month.nc 00530101.ocean_month-publicvars.nc ncatted -a modulo\,time\,c\,c\,' ' 00530101.ocean_month-publicvars.nc And copy the 21-March-init conditions as INPUT directory: cp -av .../work/Aeolus_lad-init-21-mar-clim/INPUT-init-21-mar-clim/INPUT \ ../Aeolus_lad-init-21-mar-clim/INPUT-init-21-mar-clim And change data_table accordingly to point at INPUT/00530101.ocean_month-publicvars.nc Running this for 1 year with -O3 -xHost takes just over 1 hour on 1 CPU. --------- Atmosphere 2115.676064 0.582 atmClock = A-L: update_atmos_model_dow + A-L: update_atmos_model_up Potsdam3_init 0.519297 0.000 Potsdam3_update_up 389.148002 0.107 Potsdam3_update_down 1711.047478 0.471 Potsdam3_diagnostic 305.986022 0.084 Land 113.843636 0.031 Ice 222.424857 0.061 Land-ice-atm coupler 1171.513136 0.322 cplClock = sfcClock + fluxAtmDnClock + regenClock + fluxAtmUpClock SFC boundary layer 928.505936 0.256 sfcClock < cplClock Ice-ocean coupler 6.088290 0.002 cplOcnClock generate_sfc_xgrid 6.670227 0.002 newClock1 flux_ocean_to_ice 5.755799 0.002 newClock2 flux_ice_to_ocean 0.000000 0.000 newClock3 flux_check_stocks 0.000000 0.000 newClock4 ATM 3618.048206 0.996 newClock5 ATM: update_ice_model_slow_up 5.377364 0.001 newClock6 < newClock5 ATM: atmos loop 3465.134553 0.954 newClock7 < newClock5 A-L: atmos_tracer_driver_ga 0.008572 0.000 newClocka < newClock7 A-L: sfc_boundary_layer 928.526885 0.256 newClockb < newClock7 A-L: update_atmos_model_dow 1719.460368 0.473 newClockc < newClock7 | atmClock A-L: flux_down_from_atmos 150.609045 0.041 newClockd < newClock7 A-L: update_land_model_fast 112.910724 0.031 newClocke < newClock7 A-L: update_ice_model_fast 73.158210 0.020 newClockf < newClock7 A-L: flux_up_to_atmos 84.177004 0.023 newClockg < newClock7 A-L: update_atmos_model_up 396.239910 0.109 newClockh < newClock7 | atmClock ATM: update_land_model_slow 0.957953 0.000 newClock8 < newClock5 ATM: flux_land_to_ice 1.580863 0.000 newClock9 < newClock5 ATM: update_ice_model_slow_dn 144.486968 0.040 newClock10 < newClock5 ATM: flux_ice_to_ocean_stocks 0.319596 0.000 newClock11 < newClock5 ------ run-107: 20 years took 201539 seconds, i.e. 2.7 hours per year, with Atmosphere 0.184 Land 0.011 Ice 0.024 Ocean 0.666 Land-ice-atm coupler 0.113 ============ Stefan Rahmstorf recommends to have another go with the ice -> coupler interface. - use t_surf climatology from e.g. ERA - set u_surf and v_surf to 0 because their influence on atmosphere is negligible - set roughness to some small average value - use albedo from observation data In .../Jetstream-Vis/ECMWF we have already downloaded some data: >----- >Thu Mar 10 14:04:31 CET 2016 > >downloaded skin temperature 6hourly (daily) surface level, 1980-1989 > >computed daily and monthly climatology and converted to NetCDF, removed downloaded grib files. > >cdo -W -v -t ecmwf -f nc invertlat -ymonmean -cat skt-6hourly-90S-90N-198\?.grib skt-ymonmean-90S-90N-1980-1989.nc >cdo -W -v -t ecmwf -f nc invertlat -ydaymean -cat skt-6hourly-90S-90N-198\?.grib skt-ydaymean-90S-90N-1980-1989.nc > >also downloaded albedo, sst, t2m >Evaporation seems not to be available on 6hourly (or does it depend on the time step specification???) > >for var in al sst t2m >do > cdo -W -v -t ecmwf -f nc invertlat -ymonmean -cat $var-6hourly-90S-90N-198\?.grib $var-ymonmean-90S-90N-1980-1989.nc > cdo -W -v -t ecmwf -f nc invertlat -ydaymean -cat $var-6hourly-90S-90N-198\?.grib $var-ydaymean-90S-90N-1980-1989.nc >done Maybe we can also get something for roughness From CM2M_coarse_BLING we get: stat if land_mask eq 0 then ROUGH_HEAT ROUGH_HEAT [m] surface roughness for heat Minimum value: 9.319287E-06 Maximum value: 0.00064444 Mean value: 0.00010092 (unweighted average) ROUGH_MOIST [m] surface roughness for moisture Minimum value: 0.000014445 Maximum value: 0.00074363 Mean value: 0.00011408 (unweighted average) Standard deviation: 0.00015717 ROUGH_MOM [m] surface roughness for momentum Minimum value: 0.000034326 Maximum value: 0.0011022 Mean value: 0.00026452 (unweighted average) Standard deviation: 0.00015126 The Climber-2 compatible ocean albedo is where (coszen_radwt > 0.0001) albedo_vis_dir = MIN(0.05/coszen_radwt, 0.20) elsewhere albedo_vis_dir = 0.05 end where albedo_vis_dif = 0.07 albedo_nir_dir = albedo_vis_dir albedo_nir_dif = 0.07 ---- On 06/03/16 21:53, Stefan Petri wrote: > once again Levke asked good questions for debugging the model setup... > Now we have Aeolus+land_lad, MOM and sea-ice turned off. > sst is taken from a monthly climatology file that is delivered with the > GFDL example setups, but not used in our fully-coupled runs. > Originally I tried sst from ERA-Interim, but that is defined on their > land-sea mask, and mismatches the MOM land-sea-mask. > The GFDL-delivered sst climatology is almost identical to ERA-Interim, > but extrapolated across land, so that it is robust against regridding > and differing land-sea-masks. > Albedo is taken from ERA-Interim. > For roughness I have set global constant values, taken from the > sea-surface mean values from CM2M_coarse_BLING output. > surface speeds u and v I have set to just 0. > > This simulated 1 year in 45 minutes. > I have not yet looked into the output. --------- Wed Sep 7 14:49:57 CEST 2016 Levke and Georg looked at output from above run. And also at the input. In the polar regions, the surface temperatures look bad. Also the albedo. It turns out that the SST climatology, from ERA as well as the one delivered by GFDL has temperatures _below_ sea ice. Next try: look at Skin temperature from ERA interim. That was downloaded already on 10 Mar 2016. For non-ice-covered ocean areas, that SKT is identical to the SST, but obviously has top-of-ice temperatures for arctic and antarctic. So we we use that instead of the SST as input for the next run. cp -p /p/projects/petri/Jetstream-Vis/ECMWF/skt-ymonmean-90S-90N-1980-1989.nc INPUT/. (cd INPUT ; ncatted -a calendar\,time\,m\,c\,"365_day" -a modulo\,time\,c\,c\," " skt-ymonmean-90S-90N-1980-1989.nc) ---------------- Wed Sep 28 15:55:54 CEST 2016 The albedo data from ERA interim contains only land surface data. But we need only ocean-surface data for this experiment. Thus we need to construct our own forcing file. Over open ocean, the albedo computed as in ocean_albedo() case 6, which we implemented as inherited from Climber-2. The ocean albedo depends on the radiation-weighted cosine of the zenith angle. Over sea-ice, things are different. ERA-Interim provides data for sea ice thickness / sea ice concentration. We also have the skin temperature data. That should be sufficient to mimic the calculations in ice_optics(). ice_thm.F90 subroutine ice_optics(alb, pen, trn, hs, hi, ts, tfw) real, intent( out) :: alb ! ice surface albedo (0-1) real, intent( out) :: pen ! fraction of down solar penetrating the ice real, intent( out) :: trn ! ratio of down solar at bottom to top of ice real, intent(in ) :: hs ! snow thickness (m-snow) real, intent(in ) :: hi ! ice thickness (m-ice) real, intent(in ) :: ts ! surface temperature real, intent(in ) :: tfw ! seawater freezing temperature called from ice_model.F90::ice_bottom_to_ice_top() The ice model computes first albedo for ice-covered cells, and leaves open-sea-cells at 0. Then it invokes ocean_albedo(), which fills in the open-sea cells with values for albedo_{vis,nir}_di{f,r} (but leaves albedo untouched). The coupler invokes in every time step generate_sfc_xgrid(), which in turn invokes call set_frac_area (Ice%part_size(,,)), which uses the partition sizes as area weights in the exchange grid. This way, the data from the ice partitions is area-weighted averaged in the get_from_xgrid() operation. In the fixed-surface-data case, the restart file ice_model.res.nc contains h_ice, which is used during ice_model_init() to compute part_size. Since no ice_model_update procedures are invoked, the part_size remains fixed throughout the simulation run. However, for each surface cell, the sum over all partitions of part_size is exactly 1. Furthermore, in the data_override procedure, 2-D surface data is just replicated to every partition of a 3-D variable. Thus, we provide albedo data with ocean+sea-ice albedo, for all 5 albedo types, and because sum-over-k of part_size is 1, it does not matter, what particular values are actually contained in part_size. The same holds for the other ice model quantities that are injected via data_override into the coupler (t_surf, rough_mom, rough_heat, rough_moist, u_surf, v_surf, t_surf) In our setup (derived from CM2M_coarse_BLING), the ice model namelist parameters are set to: &ICE_MODEL_NML MOM_ROUGH_ICE = 5.000000000000000E-004, HEAT_ROUGH_ICE = 5.000000000000000E-004, P0 = 27500.0000000000 , C0 = 20.0000000000000 , CDW = 3.240000000000000E-003, WD_TURN = 0.000000000000000E+000, KMELT = 240.000000000000 , ALB_SNO = 0.800000000000000 , ALB_ICE = 0.582600000000000 , PEN_ICE = 0.300000000000000 , OPT_DEP_ICE = 0.670000000000000 , NSTEPS_DYN = 72, NSTEPS_ADV = 1, NUM_PART = 6, ATMOS_WINDS = T, SLAB_ICE = F, SPEC_ICE = F, ICE_BULK_SALIN = 5.000000000000000E-003, LAYOUT = 2*0, DO_ICE_RESTORE = F, DO_ICE_LIMIT = F, MAX_ICE_LIMIT = 4.00000000000000 , ICE_RESTORE_TIMESCALE = 5.00000000000000 , SLP2OCEAN = F, CONSERVATION_CHECK = T, T_RANGE_MELT = 10.0000000000000 , CM2_BUGS = F, KS = 0.310000000000000 , H_LO_LIM = 0.000000000000000E+000, VERBOSE = F, DO_ICEBERGS = F, ADD_DIURNAL_SW = F, IO_LAYOUT = 2*0, CHANNEL_VISCOSITY = 0.000000000000000E+000, SMAG_OCN = 0.150000000000000 , SSH_GRAVITY = 9.81000000000000 , CHAN_CFL_LIMIT = 0.250000000000000 , DO_SUN_ANGLE_FOR_ALB = F, MASK_TABLE = INPUT/ice_mask_table , REPRODUCE_SIENA_201303 = T / The radiation-weighted cosine of the zenith angle we get from running coupled Aeolus with namelist parameter test_astro = T. That writes out 400 days of data. We select just the first year of that. cdo -selvar\,fms_astro_cosz_radwt -seldate\,0001-01-01T00:00:00\,0001-12-31T23:59:59 history/astro_test_daily.nc fms_astro_cosz_radwt-daily.nc It turns out that monthly means of that data yield a not nice fit to the instantaneous data curve, while daily mean data is very close. And in the albedo calculation, several non-linear operations occur (min,max). Also, the data_override module interpolates from the time axis of the given data file to each individual timestep. Thus, we just stay at daily data for now. cancel mode upper use fms_astro_cosz_radwt-daily.nc let ocean_albedo_vis_dir = if fms_astro_cosz_radwt gt 0.0001 then min(0.05/fms_astro_cosz_radwt, 0.2) else 0.05 let ocean_albedo_vis_dif = 0.07 let ocean_albedo_nir_dir = ocean_albedo_vis_dir let ocean_albedo_nir_dif = 0.07 Note that this data has axes time,lat only, and thus needs to be widened to lat dimension as well. That should happen automagically in the subsequent processing steps... from ice_thm.F90: ! albedos are from CSIM4 assumming 0.53 visible and 0.47 near-ir insolation real :: ALB_SNO = 0.85 ! albedo of snow (not melting) real :: ALB_ICE = 0.5826 ! albedo of ice (not melting) real :: T_RANGE_MELT = 1.0 ! melt albedos scaled in below melting T real, parameter :: SI = 1.0 ! salinity of sea ice real, parameter :: MU_TS = 0.054 ! relates freezing temp. to salinity real, parameter :: TFI = -MU_TS*SI ! sea ice freezing temp. = -mu*salinity !! BUT: see namelist above!! subroutine ice_optics(alb, pen, trn, hs, hi, ts, tfw) real, intent( out) :: alb ! ice surface albedo (0-1) !real, intent( out) :: pen ! fraction of down solar penetrating the ice !real, intent( out) :: trn ! ratio of down solar at bottom to top of ice real, intent(in ) :: hs ! snow thickness (m-snow) real, intent(in ) :: hi ! ice thickness (m-ice) real, intent(in ) :: ts ! surface temperature !real, intent(in ) :: tfw ! seawater freezing temperature as = ALB_SNO; ai = ALB_ICE cs = hs/(hs+0.02) ! thin snow partially covers ice fh = min(atan(5.0*hi)/atan(5.0*0.5),1.0) ! use this form from CSIM4 to ! reduce albedo for thin ice if (ts+T_RANGE_MELT > TFI) then ! reduce albedo for melting as in ! CSIM4 assuming 0.53/0.47 vis/ir as = as-0.1235*min((ts+T_RANGE_MELT-TFI)/T_RANGE_MELT,1.0) ai = ai-0.075 *min((ts+T_RANGE_MELT-TFI)/T_RANGE_MELT,1.0) endif ai = fh*ai+(1-fh)*0.06 ! reduce albedo for thin ice alb = cs*as+(1-cs)*ai For this we need hs snow thickness - available from ERA interim, but we need to check if that covers sea ice, or land only. hi ice thickness - ??? ts surface temperature - we use skin temperature skt from ERA interim We also need ocean albedo, as computed above I now found that ERA-Interim offers a variable named ``forecast albedo'' . In contrast to plain ``Albedo'' this has also values for polar sea-ice. Thus our next try: use FRAC_LAND from FMS output, as land/sea mask, and for axes definitions. use sea-ice cover ci from ERA Interim, as sea ice mask use ocean albedo as computed above, for open sea albedo use forecast albedo fal from ERA Interim, for sea-ice albedo first regrid ci, ocean albedo, and fal to FRAC_LAND. let albedo = if ci gt 0 then fal else ocean albedo This will result in ocean albedo also being used over land areas. Which in turn avoids potentially missing values at land-sea boundaries in the coupler. ! non-working ferret script removed sigh. We need to do time interpolation from monthly 1989 to daily climatology 0000 . Also, we need to be careful with calendar attributes. Using 365_day leads to wrong date interpolation. Downloaded 6-hourly ci-data in hope to overcome those problems. ----------- Mon Nov 21 11:05:14 CET 2016 - lets try to use daily climatology of ERA-interim skin temperature instead of monthly climatology, in the hope to avoid the time interpolation overhead in the data_override functions. We have that already: -rw-r--r-- 1 petri climber3 169365364 Mar 10 2016 skt-ydaymean-90S-90N-1980-1989.nc -rw-r--r-- 1 petri climber3 5559656 Mar 10 2016 skt-ymonmean-90S-90N-1980-1989.nc We must adapt the time axis attributes so that FMS data_override recognizes it as climatological ``modulo'' axis, with correct calendar attribute. Also we need to remove the data for 29-feb to adapt it to the noleap 365_day calendar that we use in MOM. cdo del29feb /p/projects/climber3/petri/Jetstream-Vis/ECMWF/skt-ydaymean-90S-90N-1980-1989.nc \ era-interim-skt-6hourly-90S-90N-1980-1989-ydaymean.nc ncatted -a calendar\,time\,m\,c\,365_day -a modulo\,time\,c\,c\," " era-interim-skt-6hourly-90S-90N-1980-1989-ydaymean.nc ERA-Interim daily values for ``forecast albedo'' seem to be available not on ``step 0'' (which appears to be for diagnostics data'', but on ``step 3/6/9/12'' (which is explained in their FAQ as to be for forecast data - which matches the variable's name). And then this forecast data is available for times 00:00:00 and 12:00:00. So, choose time=00/12 and step=6 , which gives us 2 values per day, time-stamped at 06:00 and 18:00 . cdo -W -v -t ecmwf -f nc invertlat -ydaymean -del29feb \ -cat /p/projects/climber3/petri/Jetstream-Vis/ECMWF/fal-6hourly-90S-90N-198\?.grib \ era-interim-fal-12hourly-90S-90N-1980-1989-ydaymean.nc The del29feb operator apparently also updates the calendar attribute to 365_day . ncatted -a modulo\,time\,c\,c\," " era-interim-fal-12hourly-90S-90N-1980-1989-ydaymean.nc cdo -W -v -t ecmwf -f nc invertlat -ydaymean -del29feb \ -cat /p/projects/climber3/petri/Jetstream-Vis/ECMWF/ci-6hourly-90S-90N-198\?.grib \ era-interim-ci-6hourly-90S-90N-1980-1989-ydaymean.nc ncatted -a modulo\,time\,c\,c\," " era-interim-ci-6hourly-90S-90N-1980-1989-ydaymean.nc ferret << EOF can mode up use "history-noocean-noice/00530321.aeolus-monthly.nc" use fms_astro_cosz_radwt-daily.nc use era-interim-ci-6hourly-90S-90N-1980-1989-ydaymean.nc use era-interim-fal-12hourly-90S-90N-1980-1989-ydaymean.nc set mem/size=300 let gridcosz=fms_astro_cosz_radwt[d=2] let gridci=ci[d=3] let ocean_albedo_vis_dir = if fms_astro_cosz_radwt[d=2] gt 0.0001 then min(0.05/fms_astro_cosz_radwt[d=2], 0.2) else 0.05 let ocean_albedo_nir_dir = ocean_albedo_vis_dir ! inherit grid from fms_astro_cosz_radwt let ocean_albedo_vis_dif = 0*fms_astro_cosz_radwt[d=2]+0.07 let ocean_albedo_nir_dif = 0*fms_astro_cosz_radwt[d=2]+0.07 let albedo = if ci[d=3,gt=gridcosz] gt 0 then fal[d=4,gt=gridcosz] else 0 let albedo_vis_dir = if ci[d=3,gt=gridcosz] gt 0 then fal[d=4,gt=gridcosz] else ocean_albedo_vis_dir[gxy=gridci] let albedo_nir_dir = if ci[d=3,gt=gridcosz] gt 0 then fal[d=4,gt=gridcosz] else ocean_albedo_nir_dir[gxy=gridci] let albedo_vis_dif = if ci[d=3,gt=gridcosz] gt 0 then fal[d=4,gt=gridcosz] else ocean_albedo_vis_dif[gxy=gridci] let albedo_nir_dif = if ci[d=3,gt=gridcosz] gt 0 then fal[d=4,gt=gridcosz] else ocean_albedo_nir_dif[gxy=gridci] save/file="ocean_seaice_albedo_erainterim_potsdam2.ydaymean.nc" albedo, albedo_vis_dir, albedo_nir_dir, albedo_vis_dif, albedo_nir_dif EOF ------------- Wed Nov 30 14:19 2016 Run skt_019 [label="eternal spring"]; Explore how to further speed up the coupled model runs to make automatic tuning possible. The idea is to do seasonal runs , instead of the usual annual-cycle-runs. As first test, do an eternal-spring run. That matches our currently used initial condition. Next steps would be to use a 15-Jan and a 15-Jun snapshot from a fully-coupled Aeolus run as initial conditions to test winter and summer season optimisations. That requires to change the simenv skill score function to evaluate only one seasonal time slice instead of monthly values. It is not sufficient to set fix_astr = .true. in the namelist file. The surface forcing climatology data derived from ERA interim still contain the seasonal cycle. Thus we extract only one month from the ymonmean and ydaymean climatologies. #cdo selmonth\,4 ocean_seaice_albedo_erainterim_potsdam2.ydaymean.nc ocean_seaice_albedo_erainterim_potsdam2.ydaymean-april.nc #ncatted -a modulo\,time\,c\,c\," " ocean_seaice_albedo_erainterim_potsdam2.ydaymean-april.nc # #cdo selmonth\,4 era-interim-skt-6hourly-90S-90N-1980-1989-ydaymean.nc era-interim-skt-6hourly-90S-90N-1980-1989-ydaymean-april.nc #ncatted -a modulo\,time\,c\,c\," " era-interim-skt-6hourly-90S-90N-1980-1989-ydaymean-april.nc Sigh. Over the course of April, the daily climatological global mean skin temperature still shows an (almost linear) increase from 287.5 K to 288.5 K. Thus, we should proably convert this to one monthly mean value. cdo ymonmean -selmonth\,3 ocean_seaice_albedo_erainterim_potsdam2.ydaymean.nc ocean_seaice_albedo_erainterim_potsdam2.ymonmean-march.nc ncatted -a modulo\,time\,c\,c\," " ocean_seaice_albedo_erainterim_potsdam2.ymonmean-march.nc cdo ymonmean -selmonth\,3 era-interim-skt-6hourly-90S-90N-1980-1989-ydaymean.nc era-interim-skt-6hourly-90S-90N-1980-1989-ymonmean-march.nc ncatted -a modulo\,time\,c\,c\," " era-interim-skt-6hourly-90S-90N-1980-1989-ymonmean-march.nc cdo ymonmean -selmonth\,4 ocean_seaice_albedo_erainterim_potsdam2.ydaymean.nc ocean_seaice_albedo_erainterim_potsdam2.ymonmean-april.nc ncatted -a modulo\,time\,c\,c\," " ocean_seaice_albedo_erainterim_potsdam2.ymonmean-april.nc cdo ymonmean -selmonth\,4 era-interim-skt-6hourly-90S-90N-1980-1989-ydaymean.nc era-interim-skt-6hourly-90S-90N-1980-1989-ymonmean-april.nc ncatted -a modulo\,time\,c\,c\," " era-interim-skt-6hourly-90S-90N-1980-1989-ymonmean-april.nc Same for January and July: cdo ymonmean -selmonth\,7 ocean_seaice_albedo_erainterim_potsdam2.ydaymean.nc ocean_seaice_albedo_erainterim_potsdam2.ymonmean-july.nc ncatted -a modulo\,time\,c\,c\," " ocean_seaice_albedo_erainterim_potsdam2.ymonmean-july.nc cdo ymonmean -selmonth\,7 era-interim-skt-6hourly-90S-90N-1980-1989-ydaymean.nc era-interim-skt-6hourly-90S-90N-1980-1989-ymonmean-july.nc ncatted -a modulo\,time\,c\,c\," " era-interim-skt-6hourly-90S-90N-1980-1989-ymonmean-july.nc cdo ymonmean -selmonth\,1 ocean_seaice_albedo_erainterim_potsdam2.ydaymean.nc ocean_seaice_albedo_erainterim_potsdam2.ymonmean-january.nc ncatted -a modulo\,time\,c\,c\," " ocean_seaice_albedo_erainterim_potsdam2.ymonmean-january.nc cdo ymonmean -selmonth\,1 era-interim-skt-6hourly-90S-90N-1980-1989-ydaymean.nc era-interim-skt-6hourly-90S-90N-1980-1989-ymonmean-january.nc ncatted -a modulo\,time\,c\,c\," " era-interim-skt-6hourly-90S-90N-1980-1989-ymonmean-january.nc Sigh again. Starting from our usual 21-mar-climatology initial conditions, TS stable state is approached not earlier the 8 years after init. This does not yield the hoped-for acceleration suitable for automated simenv-tuning. Now we run skt-018 for 12 / 15 / 21 months to obtain RESTART conditions for 21-mar, 21-jun, 21-dec . We hope that those are already past the transition phase. --------------- Thu Mar 23 10:18:04 CET 2017 Lapse rates look strange, especially for JJA season over the northern hemisphere, over all tested ranges of parameters aq, Gamma_1, Gamma_2. 2017-03-22 16:15 GMT+03:00 Alexey V. Eliseev : The prescribed skin temperature near the North pole is always below zero. However, surface air temperature in July there is +5C. It's puzzling. However, it is unclear how this temperature difference between the model and the reality (5C) would cause such an outlying behaviour of the lapse rate. Alexey suggests to play with sea ice albedo in the arctic.