!IDEAL:MODEL_LAYER:INITIALIZATION ! This MODULE holds the routines which are used to perform various initializations ! for the individual domains. !----------------------------------------------------------------------- MODULE module_initialize_ideal USE module_domain ! frame/module_domain.F USE module_io_domain ! share USE module_state_description ! frame USE module_model_constants ! share USE module_bc ! share USE module_timing ! frame USE module_configure ! frame USE module_init_utilities ! dyn_em USE module_soil_pre ! share USE module_polarfft USE module_setup_h_grid USE module_setup_v_grid USE module_setup_surface USE module_setup_subsurface USE module_setup_tpz USE module_setup_moist USE module_setup_uv #ifdef DM_PARALLEL USE module_dm #endif CONTAINS !------------------------------------------------------------------- ! this is a wrapper for the solver-specific init_domain routines. ! Also dereferences the grid variables and passes them down as arguments. ! This is crucial, since the lower level routines may do message passing ! and this will get fouled up on machines that insist on passing down ! copies of assumed-shape arrays (by passing down as arguments, the ! data are treated as assumed-size -- ie. f77 -- arrays and the copying ! business is avoided). Fie on the F90 designers. Fie and a pox. SUBROUTINE init_domain ( grid ) IMPLICIT NONE ! Input data. TYPE (domain), POINTER :: grid ! Local data. INTEGER :: idum1, idum2 CALL set_scalar_indices_from_config ( head_grid%id , idum1, idum2 ) CALL init_domain_rk( grid & ! #include ! ) END SUBROUTINE init_domain !------------------------------------------------------------------- SUBROUTINE init_domain_rk ( grid & ! # include ! ) IMPLICIT NONE ! Input data. TYPE (domain), POINTER :: grid # include TYPE (grid_config_rec_type) :: config_flags ! Local data INTEGER :: & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte, & i, j, k INTEGER :: ips , ipe , jps , jpe , kps , kpe INTEGER :: imsx, imex, jmsx, jmex, kmsx, kmex, & ipsx, ipex, jpsx, jpex, kpsx, kpex, & imsy, imey, jmsy, jmey, kmsy, kmey, & ipsy, ipey, jpsy, jpey, kpsy, kpey ! Local data ! Some character strings to report setup conditions and progress CHARACTER(LEN=80) :: description CHARACTER(LEN=80) :: string1, string2, string3, string4, string5 ! See detailed definitions of these variables below INTEGER :: model_setup = 0 INTEGER :: vertical_grid = 1 INTEGER :: surface_data = 0 INTEGER :: subsurface_data = 0 INTEGER :: t_sounding = 1 INTEGER :: iuv = 0 INTEGER :: imoisture = 0 INTEGER :: iscalar = 0 REAL :: T_ref = 0. ! Reference (z=0) temperature (K) REAL :: T_iso = 0. ! Isothermal temperature at the top of the model (K) REAL :: lapse_rate=0. ! Lapse rate from T_ref to T_iso (K/m) REAL :: lat_span = 0. ! Range in degrees of the Y-direction of the grid REAL :: lon_span = 0. ! Range in degrees of the X-direction of the grid REAL :: lat_ref = 0. ! Latitude of the center of the domain REAL :: lon_ref = 0. ! Longitude of the center of the domain REAL :: ht_0 ! Manually set surface height (m) REAL :: angslope_0 ! Manually set surface slope angle (radians) REAL :: azmslope_0 ! Manually set surface slope azimuth (radians) REAL :: albedo_0 ! Manually set surface albedo REAL :: thc_0 ! Manually set surface thermal inertia (MKS units) REAL :: znt_0 ! Manually set surface roughness (m) REAL :: emiss_0 ! Manually set surface infrared emissivity REAL :: co2ice_0 ! Manually set surface CO2 ice cover (kg/m^2) REAL :: h2oice_0 ! Manually set surface H2O ice cover (kg/m^2) LOGICAL :: perturb = .FALSE. ! If true, temperatures are initialized with ! a random perturbation ! WARNING: WRF can be very zonally uniform for ! flat, etc. planets if this is not set true! LOGICAL :: set_base = .FALSE. ! Create a base state to be subtracted off ! in the main WRF code. See notes where it ! is created for how it is used in the code REAL :: lat_np = 0. REAL :: lon_np = 0. REAL :: lon_0 = 0. REAL :: z_scale = 0.5 ! Used with stretched eta-level grid REAL :: p_scale = 1. ! Adjustment factor to make the P0 variable better ! match up with z=0 pressure. This variable will ! multiply all pressures. REAL, DIMENSION(:), ALLOCATABLE :: wsave, fp REAL, DIMENSION(:,:), ALLOCATABLE :: dlat, dlong, htd REAL, DIMENSION(:,:), ALLOCATABLE :: clat_glob REAL, DIMENSION(:,:), ALLOCATABLE :: clong_glob INTEGER :: i_start, i_end, j_start, j_end INTEGER :: nxx, nyy, error, icm, jcm, idealloc INTEGER :: n_filter_arrays REAL :: radius, pi, rpi, massair, column_mass REAL :: p_surf, p_level, pd_surf, qvf1, qvf2, qvf, tperturb REAL :: thtmp, ptmp, temp(3), cof1, cof2 INTEGER :: time_step, time_step_fract_num, time_step_fract_den REAL :: dt REAL :: weight, alb_mean, thc_mean, znt_mean, tsk_mean, tmn_mean REAL, DIMENSION(:), ALLOCATABLE :: tslb_mean CALL get_ijk_from_grid ( grid , & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & ips, ipe, jps, jpe, kps, kpe, & imsx, imex, jmsx, jmex, kmsx, kmex, & ipsx, ipex, jpsx, jpex, kpsx, kpex, & imsy, imey, jmsy, jmey, kmsy, kmey, & ipsy, ipey, jpsy, jpey, kpsy, kpey ) SELECT CASE ( model_data_order ) CASE ( DATA_ORDER_ZXY ) kds = grid%sd31 ; kde = grid%ed31 ; ids = grid%sd32 ; ide = grid%ed32 ; jds = grid%sd33 ; jde = grid%ed33 ; kms = grid%sm31 ; kme = grid%em31 ; ims = grid%sm32 ; ime = grid%em32 ; jms = grid%sm33 ; jme = grid%em33 ; kts = grid%sp31 ; kte = grid%ep31 ; ! note that tile is entire patch its = grid%sp32 ; ite = grid%ep32 ; ! note that tile is entire patch jts = grid%sp33 ; jte = grid%ep33 ; ! note that tile is entire patch CASE ( DATA_ORDER_XYZ ) ids = grid%sd31 ; ide = grid%ed31 ; jds = grid%sd32 ; jde = grid%ed32 ; kds = grid%sd33 ; kde = grid%ed33 ; ims = grid%sm31 ; ime = grid%em31 ; jms = grid%sm32 ; jme = grid%em32 ; kms = grid%sm33 ; kme = grid%em33 ; its = grid%sp31 ; ite = grid%ep31 ; ! note that tile is entire patch jts = grid%sp32 ; jte = grid%ep32 ; ! note that tile is entire patch kts = grid%sp33 ; kte = grid%ep33 ; ! note that tile is entire patch CASE ( DATA_ORDER_XZY ) ids = grid%sd31 ; ide = grid%ed31 ; kds = grid%sd32 ; kde = grid%ed32 ; jds = grid%sd33 ; jde = grid%ed33 ; ims = grid%sm31 ; ime = grid%em31 ; kms = grid%sm32 ; kme = grid%em32 ; jms = grid%sm33 ; jme = grid%em33 ; its = grid%sp31 ; ite = grid%ep31 ; ! note that tile is entire patch kts = grid%sp32 ; kte = grid%ep32 ; ! note that tile is entire patch jts = grid%sp33 ; jte = grid%ep33 ; ! note that tile is entire patch END SELECT CALL model_to_grid_config_rec ( grid%id , model_config_rec , config_flags ) ! here we check to see if the boundary conditions are set properly CALL boundary_condition_check( config_flags, bdyzone, error, grid%id ) grid%itimestep=0 grid%step_number = 0 ! Copied lock, stock and barrel from share/set_timekeeping.F ! ! config_flags%time_step also already has the value of "time_step", so ! we could have accessed that instead of making a call to n[ame]l[ist]_get... CALL nl_get_time_step ( grid%id, time_step ) CALL nl_get_time_step_fract_num( grid%id, time_step_fract_num ) CALL nl_get_time_step_fract_den( grid%id, time_step_fract_den ) dt = REAL(time_step) + REAL(time_step_fract_num) / REAL(time_step_fract_den) ! We are only making "dt" here so it will get written to the ! input file. To conform with what gets written in the main model itself, ! we need to multiply it by P2SI so that it is in real SI units dt = dt * P2SI CALL nl_set_dt( grid%id, dt ) grid%dt = dt #ifdef DM_PARALLEL CALL wrf_dm_bcast_bytes( icm , IWORDSIZE ) CALL wrf_dm_bcast_bytes( jcm , IWORDSIZE ) #endif ! ------------------------------------------------------------------- ! END default WRF stuff ! ------------------------------------------------------------------- ! Other variables needed right away... radius = 1./reradius pi = ACOS(-1.) rpi = 1./pi !************************************************************************* !************************************************************************* !************************************************************************* ! The use of a variable called "model_setup" will provide many of the ! necessary settings for the successful running of this subroutine. ! The different values of model_setup are defined and explained below: ! ! model_setup: ! ! = 1 Default for complete manual modification. ! The user needs to set all variables ! This will probably be especially useful for ! specialized 1D simulations ! ! EARTH ! = 31 Default for a global GCM of Earth ! Vertical levels are linear ! Initial atmosphere is isothermal at 300 K ! ! MARS ! = 41 Default for a global GCM of Mars (3D or 2D) ! Vertical levels are decided by the "vertical_grid" namelist value ! Surface data are found from the surface maps of Mars ! Initial temperature profile is a linear lapse rate from 230 K ! at the surface to 150 at 80 km, and isothermal above that ! = 42 Default for a 1D Martian simulation ! Vertical levels are decided by the "vertical_grid" namelist value ! Surface data are found from the surface maps of Mars ! Initial temperature profile is a linear lapse rate from 230 K ! at the surface to 150 at 80 km, and isothermal above that ! = 43 Default for a global GCM of a flat-planet Mars ! Vertical levels are decided by the "vertical_grid" namelist value ! Surface data are all set to uniform values, and averaged from ! the surface maps of Mars ! Initial temperature profile is a linear lapse rate from 230 K ! at the surface to 150 at 80 km, and isothermal above that ! = 44 Default for a global GCM with a rotated pole (3D) ! Vertical levels are decided by the "vertical_grid" namelist value ! Surface data are found from the surface maps of Mars ! Initial temperature profile is a linear lapse rate from 230 K ! at the surface to 150 at 80 km, and isothermal above that !************************************************************************* ! This is where you can change the default model initialization ! ! Uncomment the following line if you want to override the compilation ! defaults, and be sure to enter the appropriate value for the variable !************************************************************************* !model_setup = ! Here is where choices must be made for use of this subroutine ! Defaults based on what options (-DWRF_PLANET, -DWRF_MARS, -DWRF_TITAN, ! etc.) this file was compiled with are selected. To override those ! options, modify the code below. ! If all goes well, this is all you should ever have to do ... ! ! But if the value above is changed, then ignore the default setttings ! and random help messages... IF (model_setup == 0) THEN #if defined WRF_MARS IF ( (ide == 2) .AND. (jde == 2) ) THEN model_setup = 42 string1 = 'Mars' string2 = '1' string3 = 'N/A' string4 = '1-point interpolation' ELSE IF (config_flags%rotated_pole) THEN model_setup = 44 string1 = 'Mars' string2 = '3' IF ((ide <= 4).OR.(jde <= 4)) string2 = '2' string3 = 'Transverse cylindrical' string4 = 'Variable, location-dependent' ELSE IF (config_flags%flat_planet) THEN model_setup = 43 string1 = 'Mars' string2 = '3' IF ((ide <= 4).OR.(jde <= 4)) string2 = '2' string3 = 'Plate Carree' string4 = 'Uniform/flat, location-independent' ELSE model_setup = 41 string1 = 'Mars' string2 = '3' IF ((ide <= 4).OR.(jde <= 4)) string2 = '2' string3 = 'Plate Carree' string4 = 'Variable, location-dependent' ENDIF #elif defined WRF_EARTH IF (config_flags%rotated_pole) THEN model_setup = 34 string1 = 'Earth' string2 = '3' IF ((ide <= 4).OR.(jde <= 4)) string2 = '2' string3 = 'Transverse cylindrical' string4 = 'Uniform/flat (Held-Suarez), location-independent' ELSE model_setup = 31 string1 = 'Earth' string2 = '3' IF ((ide <= 4).OR.(jde <= 4)) string2 = '2' string3 = 'Plate Carree' string4 = 'Uniform/flat (Held-Suarez), location-independent' ENDIF #else WRITE( wrf_err_message , * ) 'ERROR: Do not use this simple '// & 'cylindrical initialization (module_initialize_simple_cyl.F) '// & 'in non-planetWRF runs' CALL wrf_error_fatal ( TRIM( wrf_err_message ) ) #endif END IF ! However, for finer control of the setup, set "model_setup" above to 0 ! and then just set these sub-variables in the CASE DEFAULT statement below ! ! Sub-variables: ! vertical_grid = 1 Linear grid (default) ! 2 Stretch grid ! this makes for nearly constant spacing in ! physical height near the model top ! 3 Power linear grid (linear grid raised to a power) ! 41 40-level grid ! 43 25-level grid ! surface_data = 0 do nothing: all surface variables will be 0.... ! 11 Earth: use WRFSI stub for data files ! and interpolation ! 41 Mars: use MGS surface maps and ADT code for ! interpolation ! subsurface_data= 0 default: select based on num_soil_layers only ! The following 2 choices require num_soil_layers=12 ! 41 Mars: Use annual average SST file to generate ! initial values of SST's ! t_sounding = 1 Isothermal atmosphere at T_ref ! 2 Constant lapse rate from T_ref to T_iso ! Isothermal at T_iso above that ! T_ref is T at z=0 (geopotential reference, not ! the surface!) ! lat_span = YYY total range in degrees that the model grid will ! span in the latitudinal direction. 180 is the ! default ! 0 Use the "dx" variable in the namelist.input file ! as the model's dx ! lon_span = XXX total range in degrees that the model grid will ! span in the longitudinal direction. 360 is the ! default ! 0 Use the "dy" variable in the namelist.input file ! as the model's dy ! lat_ref = Y Latitude of the center of the grid (-90 < Y < 90) ! lon_ref = X Longitude of the grid center (-180 < X < 180) ! iuv = 0 No initial winds in model ! imoisture = 0 No initial moisture in model ! = 1 All levels 100% relative humidity SELECT CASE (model_setup) !---------------- ! EARTH !---------------- CASE (31) ! Earth, 2D or 3D vertical_grid = 1 surface_data = 0 t_sounding = 1 iuv = 0 imoisture = 0 iscalar = 0 lat_span = 180. lon_span = 360. lat_ref = 0. lon_ref = 0. T_ref = t0 perturb = .TRUE. u_frame = 0. v_frame = 0. CASE (34) ! Earth, 2D or 3D, transverse cylindrical vertical_grid = 1 surface_data = 0 t_sounding = 1 iuv = 0 imoisture = 0 iscalar = 0 lat_span = 180. lon_span = 360. lat_ref = 0. lon_ref = 0. lat_np = 0. lon_np = 90. lon_0 = 0. T_ref = t0 perturb = .TRUE. u_frame = 0. v_frame = 0. !---------------- ! MARS !---------------- CASE (41) ! Mars 2D or 3D SELECT CASE ( kde ) CASE (26) vertical_grid = 43 string5 = '25' CASE (41) vertical_grid = 41 string5 = '40' CASE DEFAULT vertical_grid = 1 WRITE(string5,'(I3)') kde string5 = ADJUSTL(string5) string5 = TRIM(string5)//' (equal mass, linear eta)' END SELECT ! To change choice of vertical grid structure, uncomment the following ! line and enter a value for vertical_grid: !vertical_grid = surface_data = 41 subsurface_data = 41 t_sounding = 2 iuv = 0 imoisture = 0 iscalar = 0 T_ref = 230. T_iso = 150. lapse_rate = (T_ref-T_iso)/80000. lat_span = 180. lon_span = 360. lat_ref = 0. lon_ref = 0. u_frame = 0. v_frame = 0. CASE (42) ! Mars 1D SELECT CASE ( kde ) CASE (26) vertical_grid = 43 string5 = '25' CASE (41) vertical_grid = 41 string5 = '40' CASE DEFAULT vertical_grid = 1 WRITE(string5,'(I3)') kde string5 = ADJUSTL(string5) string5 = TRIM(string5)//' (equal mass, linear eta)' END SELECT ! To change choice of vertical grid structure, uncomment the following ! line and enter a value for vertical_grid: !vertical_grid = surface_data = 41 subsurface_data = 41 t_sounding = 2 iuv = 0 imoisture = 0 iscalar = 0 T_ref = 230. T_iso = 150. lapse_rate = (T_ref-T_iso)/80000. lat_span = 1.*reradius*rpi*180. ! 1 m on the planet lon_span = 1.*reradius*rpi*180. ! 1 m on the planet lat_ref = 0. lon_ref = 0. u_frame = 0. v_frame = 0. CASE (43) ! flat Mars GCM SELECT CASE ( kde ) CASE (26) vertical_grid = 43 string5 = '25' CASE (41) vertical_grid = 41 string5 = '40' CASE DEFAULT vertical_grid = 1 WRITE(string5,'(I3)') kde string5 = ADJUSTL(string5) string5 = TRIM(string5)//' (equal mass, linear eta)' END SELECT ! To change choice of vertical grid structure, uncomment the following ! line and enter a value for vertical_grid: !vertical_grid = surface_data = 41 subsurface_data = 41 t_sounding = 2 iuv = 0 imoisture = 0 iscalar = 0 T_ref = 230. T_iso = 150. lapse_rate = (T_ref-T_iso)/80000. lat_span = 180. lon_span = 360. lat_ref = 0. lon_ref = 0. perturb = .TRUE. u_frame = 0. v_frame = 0. CASE (44) ! Mars 3D with rotated pole SELECT CASE ( kde ) CASE (26) vertical_grid = 43 string5 = '25' CASE (41) vertical_grid = 41 string5 = '40' CASE DEFAULT vertical_grid = 1 WRITE(string5,'(I3)') kde string5 = ADJUSTL(string5) string5 = TRIM(string5)//' (equal mass, linear eta)' END SELECT ! To change choice of vertical grid structure, uncomment the following ! line and enter a value for vertical_grid: !vertical_grid = surface_data = 41 subsurface_data = 41 t_sounding = 2 iuv = 0 imoisture = 0 iscalar = 0 T_ref = 230. T_iso = 150. lapse_rate = (T_ref-T_iso)/80000. lat_span = 180. lon_span = 360. lat_ref = 0. lon_ref = 0. lat_np = 0. lon_np = 90. lon_0 = 0. u_frame = 0. v_frame = 0. CASE DEFAULT vertical_grid = 1 WRITE(string5,'(I3)') kde string5 = ADJUSTL(string5) string5 = TRIM(string5)//' (equal mass, linear eta)' surface_data = 0 subsurface_data = 0 t_sounding = 1 iuv = 0 imoisture = 0 iscalar = 0 lat_span = 0. lon_span = 0. lat_ref = 0. lon_ref = 0. T_ref = t0 ht_0 = 0. angslope_0 = 0. azmslope_0 = 0. albedo_0 = 0. thc_0 = 1000. znt_0 = 0.01 emiss_0 = 1. co2ice_0 = 0. h2oice_0 = 0. z_scale = 0.5 perturb = .FALSE. set_base = .FALSE. p_scale = 0. u_frame = 0. v_frame = 0. END SELECT WRITE( wrf_err_message , * ) '' CALL wrf_message ( TRIM( wrf_err_message ) ) description = 'PLANET:' WRITE(wrf_err_message,*) TRIM(description),' ',TRIM(string1) CALL wrf_message ( TRIM( wrf_err_message ) ) description = 'DIMENSIONS:' WRITE(wrf_err_message,*) TRIM(description),' ' ,TRIM(string2) CALL wrf_message ( TRIM( wrf_err_message ) ) description = 'MAP PROJECTION:' WRITE(wrf_err_message,*) TRIM(description),' ' ,TRIM(string3) CALL wrf_message ( TRIM( wrf_err_message ) ) description = 'SURFACE PROPERTIES:' WRITE(wrf_err_message,*) TRIM(description),' ' ,TRIM(string4) CALL wrf_message ( TRIM( wrf_err_message ) ) description = 'VERTICAL GRID:' WRITE(wrf_err_message,*) TRIM(description),' ' ,TRIM(string5) CALL wrf_message ( TRIM( wrf_err_message ) ) !************************************************************************* !************************************************************************* !************************************************************************* ! The following variables are now namelist.input variables and so ! override the above settings with the namelist values: lat_ref = config_flags%lat_ref lon_ref = config_flags%lon_ref lat_span = config_flags%lat_span lon_span = config_flags%lon_span p_scale = config_flags%p_scale ! ------------------------------------------------------------------- ! Set horizontal grid basic definitions ALLOCATE(dlat(ims:ime,jms:jme)) ALLOCATE(dlong(ims:ime,jms:jme)) ALLOCATE(htd(ims:ime,jms:jme)) CALL SETUP_H_GRID(radius, lat_span, lon_span, lat_ref, lon_ref, & degrad, eomeg, nxx, nyy, grid%dx, grid%dy, & grid%xlat, grid%xlong, & grid%xlat_u, grid%xlong_u, & grid%xlat_v, grid%xlong_v, & dlat, dlong, & grid%e, grid%f, grid%sina, grid%cosa, & grid%msftx, grid%msfty, & grid%msfux, grid%msfuy, & grid%msfvx, grid%msfvy, grid%msfvx_inv, & config_flags%rotated_pole,lat_np,lon_np,lon_0,& grid%clat, grid%clong, & config_flags%stretched_grid, & config_flags%we_fix, config_flags%sn_fix, & config_flags%dx_stretch, & config_flags%dy_stretch, & config_flags%lat_ref_stretch, & config_flags%lon_ref_stretch, & ids, ide, jds, jde, ims, ime, jms, jme, & its, ite, jts, jte ) #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) ) ! WARNING: this might present scaling issues on very large numbers of processors ALLOCATE( clat_glob(ids:ide,jds:jde) ) CALL wrf_patch_to_global_real ( grid%clat, clat_glob, grid%domdesc, 'xy', 'xy', & ids, ide, jds, jde, 1, 1, & ims, ime, jms, jme, 1, 1, & its, ite, jts, jte, 1, 1 ) CALL wrf_dm_bcast_real ( clat_glob , (ide-ids+1)*(jde-jds+1) ) grid%clat_xxx(ipsx:ipex,jpsx:jpex) = clat_glob(ipsx:ipex,jpsx:jpex) DEALLOCATE( clat_glob ) ALLOCATE( clong_glob(ids:ide,jds:jde) ) CALL wrf_patch_to_global_real ( grid%clong, clong_glob, grid%domdesc, 'xy', 'xy', & ids, ide, jds, jde, 1, 1, & ims, ime, jms, jme, 1, 1, & its, ite, jts, jte, 1, 1 ) CALL wrf_dm_bcast_real ( clong_glob , (ide-ids+1)*(jde-jds+1) ) grid%clong_xxx(ipsx:ipex,jpsx:jpex) = clong_glob(ipsx:ipex,jpsx:jpex) DEALLOCATE( clong_glob ) #endif ! ------------------------------------------------------------------- ! Generate L_s for output to file l_s = GET_LS(config_flags%start_day, config_flags%start_hour, & config_flags%start_minute, config_flags%start_second ) description = 'L_s:' WRITE(string1,'(F16.8)') l_s WRITE(wrf_err_message,*) TRIM(description),' ',TRIM(string1) CALL wrf_message ( TRIM( wrf_err_message ) ) ! ------------------------------------------------------------------- ! Namelist definitions grid%rdx = 1./grid%dx grid%rdy = 1./grid%dy description = 'DX:' WRITE(string1,'(F16.8)') grid%dx WRITE(wrf_err_message,*) TRIM(description),' ', & TRIM(string1),' (For use in namelist)' CALL wrf_message ( TRIM( wrf_err_message ) ) description = 'DY:' WRITE(string1,'(F16.8)') grid%dy WRITE(wrf_err_message,*) TRIM(description),' ', & TRIM(string1),' (For use in namelist)' CALL wrf_message ( TRIM( wrf_err_message ) ) WRITE( wrf_err_message , * ) '' CALL wrf_message ( TRIM( wrf_err_message ) ) CALL nl_set_mminlu(1,' ') grid%iswater = 0 grid%cen_lat = lat_ref grid%cen_lon = lon_ref grid%truelat1 = 0. grid%truelat2 = 0. grid%moad_cen_lat = lat_ref grid%stand_lon = lon_ref grid%map_proj = 0 CALL nl_set_cen_lat ( grid%id , grid%cen_lat ) CALL nl_set_cen_lon ( grid%id , grid%cen_lon ) CALL nl_set_truelat1 ( grid%id , grid%truelat1 ) CALL nl_set_truelat2 ( grid%id , grid%truelat2 ) CALL nl_set_moad_cen_lat ( grid%id , grid%moad_cen_lat ) CALL nl_set_stand_lon ( grid%id , grid%stand_lon ) CALL nl_set_map_proj ( grid%id , grid%map_proj ) ! ------------------------------------------------------------------- WRITE( wrf_err_message , * ) 'Doing: Horizontal grid and surface fields' CALL wrf_message ( TRIM( wrf_err_message ) ) ! ------------------------------------------------------------------- ! Do the 2D surface arrays ! DO j=jts,jte DO i=its,ite grid%xland(i,j) = 1 grid%landmask(i,j) = 1 grid%mavail(i,j) = 0. grid%angslope(i,j) = 0. grid%azmslope(i,j) = 0. END DO END DO IF (surface_data == 0) THEN DO j=jts,jte DO i=its,ite grid%ht(i,j) = ht_0 grid%angslope(i,j) = angslope_0 grid%azmslope(i,j) = azmslope_0 grid%albedo(i,j) = albedo_0 grid%thc(i,j) = thc_0 grid%znt(i,j) = znt_0 grid%emiss(i,j) = emiss_0 grid%co2ice(i,j) = co2ice_0 grid%h2oice(i,j) = h2oice_0 END DO END DO ELSE CALL SETUP_SURFACE_ARRAYS(surface_data, nxx, nyy, & grid%map_proj, & grid%cen_lat, grid%cen_lon, & grid%dx, grid%dy, radius, & grid%truelat1, grid%truelat2, & grid%xlat, grid%xlong, dlat, dlong, & grid%msftx, grid%msfty, grid%ht, htd,& grid%albedo, grid%thc, grid%znt, & grid%emiss, grid%co2ice, grid%h2oice,& grid%soil_ice_pc, grid%soil_ice_dp, & grid%mavail, grid%xland, & config_flags%flat_planet, & config_flags%rotated_pole, & lat_np, lon_np, lon_0, & grid%clat, grid%clong, & config_flags%stretched_grid, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) END IF DEALLOCATE(dlat,dlong,stat=idealloc) IF (idealloc > 0) THEN WRITE( wrf_err_message , * ) 'ERROR: There was a problem deallocating', & 'dlat and dlong, STATUS=', idealloc CALL wrf_error_fatal ( TRIM( wrf_err_message ) ) END IF IF ((config_flags%flat_planet).OR.(config_flags%polar_filter)) THEN WRITE( wrf_err_message , * ) 'Doing: Polar filtering or smoothing ', & 'of surface fields' CALL wrf_message ( TRIM( wrf_err_message ) ) END IF ! If we have set the "flat_planet" flag, then we don't need to do ! any polar filtering (it may actually ruin the "flatness"), so put ! both options in one big set of "if"s. IF (config_flags%flat_planet) THEN ! The following averaging code assumes we have access to the ! entire domain on one processor. Check that this is true. IF ( (its /= ids) .OR. & (ite /= ide) .OR. & (jts /= jds) .OR. & (jte /= jde) ) THEN WRITE(wrf_err_message,*) 'ERROR: module_initialize_simple_cyl: ',& '{i,j}t{s,e} /= {i,j}d{s,e}:',its,ids,ite,ide,jts,jds,jte,jde CALL wrf_error_fatal ( TRIM( wrf_err_message ) ) END IF weight = SUM(1./(grid%msftx(ids:ide-1,jds:jde-1)* & grid%msfty(ids:ide-1,jds:jde-1))) alb_mean = SUM(grid%albedo(ids:ide-1,jds:jde-1)/ & grid%msftx(ids:ide-1,jds:jde-1)/ & grid%msfty(ids:ide-1,jds:jde-1)) / weight thc_mean = SUM(grid%thc(ids:ide-1,jds:jde-1)/ & grid%msftx(ids:ide-1,jds:jde-1)/ & grid%msfty(ids:ide-1,jds:jde-1)) / weight znt_mean = SUM(grid%znt(ids:ide-1,jds:jde-1)/ & grid%msftx(ids:ide-1,jds:jde-1)/ & grid%msfty(ids:ide-1,jds:jde-1)) / weight DO j=jts,MIN(jte,jde-1) DO i=its,MIN(ite,ide-1) grid%ht(i,j) = 0. grid%albedo(i,j) = alb_mean grid%thc(i,j) = thc_mean grid%znt(i,j) = znt_mean grid%soil_ice_pc(i,j) = 0. grid%soil_ice_dp(i,j) = -9999. END DO END DO ELSE IF (config_flags%polar_filter) THEN n_filter_arrays = 7 ! Assumes we have at least n_filter_arrays vertical levels. ! Check to make sure. IF ((kte-kts+1) < n_filter_arrays) THEN WRITE (wrf_err_message,*) 'ERROR: module_initialize_simple_cyl: ',& 'vertical levels <',n_filter_arrays,(kte-kts+1),kts,kte CALL wrf_error_fatal ( TRIM( wrf_err_message ) ) END IF ! We stick the arrays in an unused 3d array (unused at this point) DO j = jts, MIN(jte,jde-1) DO k = kts, kte DO i = its, MIN(ite,ide-1) grid%t_init(i,k,j) = 1. END DO END DO DO i = its, MIN(ite,ide-1) grid%t_init(i,1,j) = grid%ht(i,j) grid%t_init(i,2,j) = htd(i,j) grid%t_init(i,3,j) = grid%albedo(i,j) grid%t_init(i,4,j) = grid%thc(i,j) grid%t_init(i,5,j) = grid%znt(i,j) grid%t_init(i,6,j) = grid%co2ice(i,j) grid%t_init(i,7,j) = grid%h2oice(i,j) END DO END DO #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) ) ! Transpose arrays so all latitudes are on a patch # include "XPOSE_POLAR_FILTER_TOPO_z2x.inc" CALL polar_filter_3d(grid%t_xxx, grid%clat_xxx, .FALSE., & config_flags%fft_filter_lat, 0., & ids, ide, jds, jde, kds, kde, & imsx, imex, jmsx, jmex, kmsx, kmex, & ipsx, ipex, jpsx, jpex, kpsx, kpex ) ! Transpose it back to "all z on a patch". # include "XPOSE_POLAR_FILTER_TOPO_x2z.inc" #else ! The following averaging code assumes we have access to the ! entire E/W domain on one processor. Check that this is true. IF ( (its /= ids) .OR. (ite /= ide) ) THEN WRITE (wrf_err_message,*) 'ERROR: module_initialize_simple_cyl: ',& 'it{s,e} /= id{s,e}:',its,ids,ite,ide CALL wrf_error_fatal ( TRIM( wrf_err_message ) ) END IF CALL polar_filter_3d(grid%t_init, grid%clat, .FALSE., & config_flags%fft_filter_lat, 0., & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) #endif ! Get the filtered data back in the original arrays DO j = jts, MIN(jte,jde-1) DO i = its, MIN(ite,ide-1) grid%ht(i,j) = grid%t_init(i,1,j) htd(i,j) = grid%t_init(i,2,j) grid%albedo(i,j) = grid%t_init(i,3,j) grid%thc(i,j) = grid%t_init(i,4,j) grid%znt(i,j) = grid%t_init(i,5,j) grid%co2ice(i,j) = grid%t_init(i,6,j) grid%h2oice(i,j) = grid%t_init(i,7,j) END DO END DO END IF WRITE( wrf_err_message , * ) 'Doing: Topograhics slopes' CALL wrf_message ( TRIM( wrf_err_message ) ) CALL SLOPES(htd, grid%angslope, grid%azmslope, & grid%xlat, grid%xlong, & grid%msftx, grid%msfty, & grid%dx, grid%dy, & grid%cen_lat, grid%cen_lon, & grid%map_proj, & config_flags%rotated_pole, & lat_np, lon_np, lon_0, & config_flags%stretched_grid, & ids , ide , jds , jde , kds , kde , & ims , ime , jms , jme , kms , kme , & its , ite , jts , jte , kts , kte ) IF (config_flags%flat_planet) THEN DO j = jts, MIN(jte,jde-1) DO i = its, MIN(ite,ide-1) grid%angslope(i,j) = 0. grid%azmslope(i,j) = 0. END DO END DO END IF WRITE( wrf_err_message , * ) 'Doing: Vertical grid' CALL wrf_message ( TRIM( wrf_err_message ) ) ! ------------------------------------------------------------------- ! Set up the vertical grid definitions CALL SETUP_V_GRID(vertical_grid, grid%znw, z_scale, kds, kde, kms, kme, kts, kte) DO k=1, kde-1 grid%dnw(k) = grid%znw(k+1) - grid%znw(k) grid%rdnw(k) = 1./grid%dnw(k) grid%znu(k) = 0.5*(grid%znw(k+1)+grid%znw(k)) ENDDO DO k=2, kde-1 grid%dn(k) = 0.5*(grid%dnw(k)+grid%dnw(k-1)) grid%rdn(k) = 1./grid%dn(k) grid%fnp(k) = .5* grid%dnw(k )/grid%dn(k) grid%fnm(k) = .5* grid%dnw(k-1)/grid%dn(k) ENDDO cof1 = (2.*grid%dn(2)+grid%dn(3))/(grid%dn(2)+grid%dn(3))*grid%dnw(1)/grid%dn(2) cof2 = grid%dn(2) /(grid%dn(2)+grid%dn(3))*grid%dnw(1)/grid%dn(3) grid%cf1 = grid%fnp(2) + cof1 grid%cf2 = grid%fnm(2) - cof1 - cof2 grid%cf3 = cof2 grid%cfn = (.5*grid%dnw(kde-1)+grid%dn(kde-1))/grid%dn(kde-1) grid%cfn1 = -.5*grid%dnw(kde-1)/grid%dn(kde-1) ! ------------------------------------------------------------------- ! Do the hard work!!! Set up the initial conditions ! initialize random number generator if necessary IF (perturb) CALL random_seed ! find ptop for the desired ztop (ztop is input from the namelist) grid%p_top = Z2P(t_sounding, config_flags%ztop, T_ref, T_iso, lapse_rate, & g, R_d, p0, p_scale) IF (((ide-ids) == 1) .AND. ((jde-jds) == 1)) THEN ! For the 1D model, all memory (not just domain and/or tile) indices ! have to be filled with the same value. Since memory limits are being ! used, X(:,:) is exactly the same as X(ims:ime,jms:jme) i_start = ims i_end = ime j_start = jms j_end = jme ELSE i_start = its i_end = ite j_start = jts j_end = jte END IF #ifdef WRF_MARS ! Do a quick calculation of total air mass, and adjust it (via p_scale) ! to a "balanced" value, based on a total from the Mars GCM, but ! only if the domain spans the entire planet AND if we have access to ! all those grid points on the same processor IF ((lat_span == 180.) .AND. (lon_span == 360.)) THEN WRITE( wrf_err_message , * ) 'Doing: Mass adjustment (Mars, ', & 'global domain)' CALL wrf_message ( TRIM( wrf_err_message ) ) IF ( (its /= ids) .OR. (ite /= ide) .OR. & (jts /= jds) .OR. (jte /= jde) ) THEN WRITE (wrf_err_message,*) 'ERROR: module_initialize_simple_cyl: ',& 'All grid points are not on the same processor, so summed ', & 'adjustment to global air mass can not be performed.' CALL wrf_message ( TRIM( wrf_err_message ) ) WRITE (wrf_err_message,*) 'ERROR: module_initialize_simple_cyl: ',& '{i,j}t{s,e} /= {i,j}d{s,e}:',its,ids,ite,ide,jts,jds,jte,jde CALL wrf_error_fatal ( TRIM( wrf_err_message ) ) ELSE massair = 0. DO J = j_start, j_end DO I = i_start, i_end p_surf = Z2P(t_sounding, grid%ht(i,j), T_ref, T_iso, lapse_rate, & g, r_d, p0, 1.) !BUG FIX 12 March 2010: incorrect 'p_scale' changed to '1.' column_mass = p_surf-grid%p_top IF (config_flags%stretched_grid) THEN massair = massair + (grid%dx/grid%msftx(i,j)) * & (grid%dy/grid%msfty(i,j)) * & column_mass/g ELSE ! If the grid is not rotated, then clat=xlat and we are OK massair = massair + grid%dy*grid%dx*COS(grid%clat(i,j)*pi/180.)* & column_mass/g END IF END DO END DO ! Total mass should be 2.8582e16 kg to better match VL1/2 pressure ! curve data p_scale = p_scale*2.8582e16/massair CALL wrf_message ( TRIM( wrf_err_message ) ) grid%p_top = Z2P(t_sounding,config_flags%ztop,T_ref,T_iso,lapse_rate, & g, R_d, p0, p_scale) END IF END IF #endif WRITE( wrf_err_message , * ) 'Doing: Meteorological fields' CALL wrf_message ( TRIM( wrf_err_message ) ) ! Values of geopotential (base, perturbation, and at p0) at the surface DO j=j_start,j_end DO i=i_start,i_end ! This code assumes that kts = 1 (see loop from kts+1 to kte below) grid%phb(i,1,j) = grid%ht(i,j)*g grid%php(i,1,j) = 0. ! This is perturbation geopotential ! Since this is an initial condition, there ! should be no perturbation! grid%ph0(i,1,j) = grid%ht(i,j)*g ENDDO ENDDO DO J = j_start, j_end DO I = i_start, i_end p_surf = Z2P(t_sounding, grid%phb(i,1,j)/g, T_ref, T_iso, lapse_rate, & g, r_d, p0, p_scale) grid%mub(i,j) = p_surf-grid%p_top ! given p (coordinate), calculate theta and compute 1/rho from equation ! of state DO K = kts, kte-1 p_level = grid%znu(k)*(p_surf - grid%p_top) + grid%p_top grid%pb(i,k,j) = p_level grid%t_init(i,k,j) = P2T(t_sounding, p_level, T_ref, T_iso, lapse_rate, & g, r_d, p0, p_scale) grid%t_init(i,k,j) = grid%t_init(i,k,j)*(p0/p_level)**rcp grid%t_init(i,k,j) = grid%t_init(i,k,j) - t0 grid%alb(i,k,j)=(r_d/p1000mb)*(grid%t_init(i,k,j)+t0)*(grid%pb(i,k,j)/p1000mb)**cvpm END DO ! calculate hydrostatic balance (alternatively we could interpolate ! the geopotential from the sounding, but this assures that the base ! state is in exact hydrostatic balance with respect to the model eqns. DO k = kts+1, kte grid%phb(i,k,j) = grid%phb(i,k-1,j) - grid%dnw(k-1)*grid%mub(i,j)*grid%alb(i,k-1,j) ENDDO ENDDO ENDDO ! Calculate moisture in the column CALL SETUP_MOIST(imoisture, grid%moist, num_moist, & PARAM_FIRST_SCALAR, P_QV, P_QI, & grid%t_init,t0,grid%pb,p0,rcp, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) ! Now calculate the full (hydrostatically-balanced) state for each column ! We will include moisture DO J = j_start, j_end DO I = i_start, i_end ! At this point p_top is already set. find the DRY mass in the column pd_surf = Z2P(t_sounding, grid%phb(i,1,j)/g, T_ref, T_iso, lapse_rate, & g, r_d, p0, p_scale) ! compute the perturbation mass (mu/mu_1/mu_2) and the full mass grid%mu_1(i,j) = pd_surf-grid%p_top - grid%mub(i,j) grid%mu_2(i,j) = grid%mu_1(i,j) grid%mu0(i,j) = grid%mu_1(i,j) + grid%mub(i,j) ! given the dry pressure and coordinate system, calculate the ! perturbation potential temperature (t/t_1/t_2) DO k = kds, kde-1 p_level = grid%znu(k)*(pd_surf - grid%p_top) + grid%p_top grid%t_1(i,k,j) = P2T(t_sounding,p_level,T_ref,T_iso,lapse_rate,g,r_d,& p0, p_scale) grid%t_1(i,k,j) = grid%t_1(i,k,j)*(p0/p_level)**rcp ! Add a random temperature perturbation if requested ! Maximum amplitude is 0.2% of actual temperature IF (perturb) THEN CALL random_number(tperturb) grid%t_1(i,k,j)=grid%t_1(i,k,j)*(1.0+0.004*(tperturb-0.5)) END IF grid%t_1(i,k,j) = grid%t_1(i,k,j)-t0 grid%t_2(i,k,j) = grid%t_1(i,k,j) END DO ! integrate the hydrostatic equation (from the RHS of the bigstep ! vertical momentum equation) down from the top to get p. ! first from the top of the model to the top pressure k = kte-1 ! top level qvf1 = 0.5*(grid%moist(i,k,j,P_QV)+grid%moist(i,k,j,P_QV)) qvf2 = 1./(1.+qvf1) qvf1 = qvf1*qvf2 grid%p(i,k,j) = - 0.5*(grid%mu_1(i,j)+qvf1*grid%mub(i,j))/grid%rdnw(k)/qvf2 qvf = 1. + rvovrd*grid%moist(i,k,j,P_QV) grid%alt(i,k,j) = (r_d/p1000mb)*(grid%t_1(i,k,j)+t0)*qvf* & (((grid%p(i,k,j)+grid%pb(i,k,j))/p1000mb)**cvpm) grid%al(i,k,j) = grid%alt(i,k,j) - grid%alb(i,k,j) ! down the column do k=kte-2,kts,-1 qvf1 = 0.5*(grid%moist(i,k,j,P_QV)+grid%moist(i,k+1,j,P_QV)) qvf2 = 1./(1.+qvf1) qvf1 = qvf1*qvf2 grid%p(i,k,j) = grid%p(i,k+1,j) - (grid%mu_1(i,j) + qvf1*grid%mub(i,j))/qvf2/grid%rdn(k+1) qvf = 1. + rvovrd*grid%moist(i,k,j,P_QV) grid%alt(i,k,j) = (r_d/p1000mb)*(grid%t_1(i,k,j)+t0)*qvf* & (((grid%p(i,k,j)+grid%pb(i,k,j))/p1000mb)**cvpm) grid%al(i,k,j) = grid%alt(i,k,j) - grid%alb(i,k,j) enddo ! this is the hydrostatic equation used in the model after the ! small timesteps. In the model, al (inverse density) ! is computed from the geopotential. grid%ph_1(i,1,j) = 0. DO k = kts+1,kte grid%ph_1(i,k,j) = grid%ph_1(i,k-1,j) - (1./grid%rdnw(k-1))*( & (grid%mub(i,j)+grid%mu_1(i,j))*grid%al(i,k-1,j)+ & grid%mu_1(i,j)*grid%alb(i,k-1,j) ) grid%ph_2(i,k,j) = grid%ph_1(i,k,j) grid%ph0(i,k,j) = grid%ph_1(i,k,j) + grid%phb(i,k,j) ENDDO END DO END DO ! Now set U & V CALL SETUP_UV(iuv, grid%u_1, grid%u_2, & grid%v_1, grid%v_2, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) ! Initialize w & omega (d-eta/dt) DO j=j_start,j_end DO k=kds,kde DO i=i_start,i_end grid%ww(i,k,j) = 0. grid%w_1(i,k,j) = 0. grid%w_2(i,k,j) = 0. grid%h_diabatic(i,k,j) = 0. END DO END DO END DO ! ------------------------------------------------------------------- ! The "base" variables. DO k=kts,kte grid%t_base(k) = 0. grid%qv_base(k) = 0. grid%u_base(k) = 0. grid%v_base(k) = 0. END DO IF (set_base) THEN DO k=1,kte-1 grid%t_base(k) = grid%t_1(1,k,1) grid%qv_base(k) = grid%moist(1,k,1,P_QV) grid%u_base(k) = grid%u_1(1,k,1) grid%v_base(k) = grid%v_1(1,k,1) END DO END IF WRITE( wrf_err_message , * ) 'Doing: Sub-surface fields' CALL wrf_message ( TRIM( wrf_err_message ) ) ! ------------------------------------------------------------------- ! Do the subsurface arrays SELECT CASE (config_flags%num_soil_layers) CASE (1) ! Only one layer: infinite slab at constant temperature below the surface ! Surface temperature is an infinitely thin "skin" on top of an ! half-infinite slab... ! And the temperature of both the skin and the slab are determined from ! the initial nearest-surface-air-layer temperature... DO J = jts, MIN(jte, jde-1) DO I = its, MIN(ite, ide-1) thtmp = grid%t_2(i,1,j)+t0 ptmp = grid%p(i,1,j)+grid%pb(i,1,j) temp(1) = thtmp * (ptmp/p1000mb)**rcp thtmp = grid%t_2(i,2,j)+t0 ptmp = grid%p(i,2,j)+grid%pb(i,2,j) temp(2) = thtmp * (ptmp/p1000mb)**rcp thtmp = grid%t_2(i,3,j)+t0 ptmp = grid%p(i,3,j)+grid%pb(i,3,j) temp(3) = thtmp * (ptmp/p1000mb)**rcp grid%tsk(I,J)=cf1*temp(1)+cf2*temp(2)+cf3*temp(3) grid%tmn(I,J)=grid%tsk(I,J)-0.5 END DO END DO CASE (12,15) ! 12-(or 15-)layer subsurface model with temperatures initialized ! from a long-term annual average CALL setup_soil_planetary(grid%zs, grid%dzs, grid%tsk, & grid%tmn, grid%tslb, & config_flags%num_soil_layers, & grid%xlat, grid%xlong, & subsurface_data, & config_flags%rotated_pole, & lat_np, lon_np, lon_0, & grid%clat, grid%clong, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte ) IF ((config_flags%flat_planet).OR.(config_flags%polar_filter)) THEN WRITE( wrf_err_message , * ) 'Doing: Polar filtering or smoothing ',& 'of sub-surface fields' CALL wrf_message ( TRIM( wrf_err_message ) ) END IF IF (config_flags%flat_planet) THEN ! The following averaging code assumes we have access to the ! entire domain on one processor. Check that this is true. IF ( (its /= ids) .OR. & (ite /= ide) .OR. & (jts /= jds) .OR. & (jte /= jde) ) THEN WRITE (wrf_err_message,*) 'ERROR: module_initialize_simple_cyl: ',& '{i,j}t{s,e} /= {i,j}d{s,e}:',its,ids,ite,ide,jts,jds,jte,jde CALL wrf_error_fatal ( TRIM( wrf_err_message ) ) END IF weight = SUM(1./(grid%msftx(ids:ide-1,jds:jde-1)* & grid%msfty(ids:ide-1,jds:jde-1))) tsk_mean = SUM(grid%tsk(ids:ide-1,jds:jde-1)/ & grid%msftx(ids:ide-1,jds:jde-1)/ & grid%msfty(ids:ide-1,jds:jde-1)) / weight tmn_mean = SUM(grid%tmn(ids:ide-1,jds:jde-1)/ & grid%msftx(ids:ide-1,jds:jde-1)/ & grid%msfty(ids:ide-1,jds:jde-1)) / weight ALLOCATE(tslb_mean(config_flags%num_soil_layers)) DO k=1,config_flags%num_soil_layers tslb_mean(k) = SUM(grid%tslb(ids:ide-1,k,jds:jde-1)/ & grid%msftx(ids:ide-1,jds:jde-1)/ & grid%msfty(ids:ide-1,jds:jde-1)) / weight END DO DO j=jts,MIN(jte,jde-1) DO i=its,MIN(ite,ide-1) grid%tsk(i,j) = tsk_mean grid%tmn(i,j) = tmn_mean DO k=1,config_flags%num_soil_layers grid%tslb(i,k,j) = tslb_mean(k) END DO END DO END DO ELSE IF (config_flags%polar_filter) THEN n_filter_arrays = config_flags%num_soil_layers+2 ! Assumes we have at least n_filter_arrays vertical levels. ! Check to make sure. IF ((kte-kts+1) < n_filter_arrays) THEN WRITE (wrf_err_message,*) 'ERROR: module_initialize_simple_cyl: ',& 'vertical levels <',n_filter_arrays,(kte-kts+1),kts,kte CALL wrf_error_fatal ( TRIM( wrf_err_message ) ) END IF ! We stick the arrays in an unused 3d array (unused at this point) DO j = jts, MIN(jte,jde-1) DO k = kts, kte DO i = its, MIN(ite,ide-1) grid%ru_m(i,k,j) = 1. END DO END DO DO i = its, MIN(ite,ide-1) grid%ru_m(i,1,j) = grid%tmn(i,j) grid%ru_m(i,2,j) = grid%tsk(i,j) DO k = 1,config_flags%num_soil_layers grid%ru_m(i,k+2,j) = grid%tslb(i,k,j) END DO END DO END DO #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) ) ! Transpose arrays so all latitudes are on a patch # include "XPOSE_POLAR_FILTER_RU_z2x.inc" CALL polar_filter_3d(grid%ru_xxx, grid%clat_xxx, .FALSE., & config_flags%fft_filter_lat, 0., & ids, ide, jds, jde, kds, kde, & imsx, imex, jmsx, jmex, kmsx, kmex, & ipsx, ipex, jpsx, jpex, kpsx, kpex ) ! Transpose it back to "all z on a patch". # include "XPOSE_POLAR_FILTER_RU_x2z.inc" #else ! The following averaging code assumes we have access to the ! entire E/W domain on one processor. Check that this is true. IF ( (its /= ids) .OR. (ite /= ide) ) THEN WRITE (wrf_err_message,*) 'ERROR: module_initialize_simple_cyl: ',& 'it{s,e} /= id{s,e}:',its,ids,ite,ide CALL wrf_error_fatal ( TRIM( wrf_err_message ) ) END IF CALL polar_filter_3d(grid%ru_m, grid%clat, .FALSE., & config_flags%fft_filter_lat, 0., & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) #endif ! Get the filtered data back in the original arrays DO j = jts, MIN(jte,jde-1) DO i = its, MIN(ite,ide-1) grid%tmn(i,j) = grid%ru_m(i,1,j) grid%tsk(i,j) = grid%ru_m(i,2,j) DO k = 1,config_flags%num_soil_layers grid%tslb(i,k,j) = grid%ru_m(i,2+k,j) END DO END DO END DO END IF CASE DEFAULT CALL process_soil_ideal ( grid%xland, grid%xice, & grid%vegfra, grid%snow, & grid%canwat, grid%ivgtyp, & grid%isltyp, grid%tslb, & grid%smois, grid%tsk, grid%tmn, & grid%zs, grid%dzs, & config_flags%num_soil_layers, & config_flags%sf_surface_physics,& ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte ) ! process_soil_ideal resets the values for ivgtyp and xland, ! so we need to reset them here. DO j=jts,jte DO i=its,ite grid%ivgtyp(i,j) = 19 grid%lu_index(i,j) = REAL(grid%ivgtyp(i,j)) grid%xland(i,j) = 1. END DO END DO END SELECT WRITE( wrf_err_message , * ) '' CALL wrf_message ( TRIM( wrf_err_message ) ) RETURN END SUBROUTINE init_domain_rk !--------------------------------------------------------------------- SUBROUTINE init_module_initialize END SUBROUTINE init_module_initialize !--------------------------------------------------------------------- REAL FUNCTION get_ls(iday,ihour,iminute,isecond) RESULT (ls) IMPLICIT NONE INTEGER, INTENT(IN) :: iday, ihour, iminute, isecond ! Local Variables REAL(KIND (0d0)) :: pi REAL(KIND (0d0)) :: deleqn, date_dbl, julian REAL(KIND (0d0)) :: er, qq, e, cd0, ep, em REAL(KIND (0d0)) :: eq, w, als REAL(KIND (0d0)), PARAMETER :: SMALL_VALUE = 1.D-6 PI=ACOS(-1.d0) deleqn = equinox_fraction * REAL(PLANET_YEAR) !-----CALCULATE LONGITUDE OF THE SUN FROM VERNAL EQUINOX: julian=REAL(iday)-0.5 + REAL(ihour)/24. + REAL(iminute)/1440. + & REAL(isecond)/86400. ! DATE = DAYS SINCE LAST PERIHELION PASSAGE date_dbl = julian - zero_date DO WHILE (date_dbl < 0.) date_dbl=date_dbl+REAL(PLANET_YEAR) END DO DO WHILE (date_dbl > REAL(PLANET_YEAR)) date_dbl=date_dbl-REAL(PLANET_YEAR) END DO er = SQRT( (1.0+eccentricity)/(1.0-eccentricity) ) ! qq is the mean anomaly qq = 2.0 * (pi * deleqn / REAL(PLANET_YEAR)) ! determine true anomaly at equinox: eq ! Iteration for eq e = 1.d0 cd0 = 1.d0 DO WHILE (cd0 > SMALL_VALUE) ep = e - (e-eccentricity*SIN(e)-qq) / (1.d0-eccentricity*COS(e)) cd0 = ABS(e-ep) e = ep END DO eq = 2.d0 * ATAN( er * TAN(0.5d0*e) ) ! determine true anomaly at current date: w ! Iteration for w em = 2.d0 * pi * date_dbl / REAL(PLANET_YEAR) e = 1.d0 cd0 = 1.d0 DO WHILE (cd0 > SMALL_VALUE) ep = e-(e-eccentricity*SIN(e)-em) / (1.0-eccentricity*COS(e)) cd0 = ABS(e-ep) e = ep END DO w = 2.d0 * ATAN( er * TAN(0.5d0*e) ) ! Radius vector ( astronomical units: AU ) als= (w - eq)*180.d0/pi !Aerocentric Longitude IF (als < 0.d0) als=als+360.d0 ls = REAL(als) END FUNCTION get_ls END MODULE module_initialize_ideal