PROGRAM tc_data USE module_machine USE module_domain, ONLY : domain, alloc_and_configure_domain, & domain_clock_set, head_grid, program_name, domain_clockprint, & set_current_grid_ptr USE module_io_domain USE module_initialize_real, ONLY : wrfu_initialize USE module_driver_constants USE module_configure, ONLY : grid_config_rec_type, model_config_rec, & initial_config, get_config_as_buffer, set_config_as_buffer USE module_timing USE module_state_description, ONLY: tconly USE module_dm, ONLY: wrf_dm_initialize USE module_symbols_util, ONLY: wrfu_cal_gregorian USE module_utility, ONLY : WRFU_finalize IMPLICIT NONE REAL :: time , bdyfrq INTEGER :: loop , levels_to_process , debug_level TYPE(domain) , POINTER :: null_domain TYPE(domain) , POINTER :: grid , another_grid TYPE(domain) , POINTER :: grid_ptr , grid_ptr2 TYPE (grid_config_rec_type) :: config_flags INTEGER :: number_at_same_level INTEGER :: myproc, nproc INTEGER :: max_dom, domain_id , grid_id , parent_id , parent_id1 , id INTEGER :: e_we , e_sn , i_parent_start , j_parent_start INTEGER :: idum1, idum2 INTEGER :: nbytes INTEGER, PARAMETER :: configbuflen = 4* 65536 INTEGER :: configbuf( configbuflen ) LOGICAL , EXTERNAL :: wrf_dm_on_monitor LOGICAL found_the_id INTEGER :: ids , ide , jds , jde , kds , kde INTEGER :: ims , ime , jms , jme , kms , kme INTEGER :: ips , ipe , jps , jpe , kps , kpe INTEGER :: ijds , ijde , spec_bdy_width INTEGER :: i , j , k , idts, rc INTEGER :: sibling_count , parent_id_hold , dom_loop CHARACTER (LEN=80) :: message CHARACTER(LEN=256) :: some_string INTEGER :: start_year , start_month , start_day , start_hour , start_minute , start_second INTEGER :: end_year , end_month , end_day , end_hour , end_minute , end_second INTEGER :: interval_seconds , real_data_init_type INTEGER :: time_loop_max , time_loop, bogus_id, storm real::t1,t2 real :: latc_loc(max_bogus),lonc_loc(max_bogus),vmax(max_bogus),rmax(max_bogus) real :: rankine_lid INTERFACE SUBROUTINE Setup_Timekeeping( grid ) USE module_domain, ONLY : domain TYPE(domain), POINTER :: grid END SUBROUTINE Setup_Timekeeping END INTERFACE CHARACTER (LEN=*), PARAMETER :: release_version = 'V4.2.1' program_name = "TC_EM " // TRIM(release_version) // " PREPROCESSOR" CALL disable_quilting CALL wrf_debug ( 100 , 'real_em: calling init_modules ' ) CALL init_modules(1) CALL WRFU_Initialize( defaultCalKind=WRFU_CAL_GREGORIAN, rc=rc ) CALL init_modules(2) IF ( wrf_dm_on_monitor () ) THEN CALL wrf_get_myproc (myproc) CALL wrf_get_nproc (nproc) IF ( nproc .GT. 1 ) THEN WRITE (some_string,*) 'My MPI processor number = ',myproc CALL wrf_debug ( 0 , TRIM(some_string) ) WRITE (some_string,*) 'Total MPI requested tasks = ',nproc CALL wrf_debug ( 0 , TRIM(some_string) ) CALL wrf_error_fatal3("",96,& 'MPI TC bogus must run with a single processor only, re-run with num procs set to 1' ) END IF END IF IF ( wrf_dm_on_monitor() ) THEN CALL initial_config END IF CALL get_config_as_buffer( configbuf, configbuflen, nbytes ) CALL wrf_dm_bcast_bytes( configbuf, nbytes ) CALL set_config_as_buffer( configbuf, configbuflen ) CALL wrf_dm_initialize CALL nl_get_debug_level ( 1, debug_level ) CALL set_wrf_debug_level ( debug_level ) CALL wrf_message ( program_name ) CALL nl_set_use_wps_input ( 1 , TCONLY ) NULLIFY( null_domain ) CALL wrf_debug ( 100 , 'real_em: calling alloc_and_configure_domain ' ) CALL alloc_and_configure_domain ( domain_id = 1 , & grid = head_grid , & parent = null_domain , & kid = -1 ) grid => head_grid CALL nl_get_max_dom ( 1 , max_dom ) IF ( model_config_rec%interval_seconds .LE. 0 ) THEN CALL wrf_error_fatal3("",137,& 'namelist value for interval_seconds must be > 0') END IF all_domains : DO domain_id = 1 , max_dom IF ( ( model_config_rec%input_from_file(domain_id) ) .OR. & ( domain_id .EQ. 1 ) ) THEN CALL Setup_Timekeeping ( grid ) CALL set_current_grid_ptr( grid ) CALL domain_clockprint ( 150, grid, & 'DEBUG real: clock after Setup_Timekeeping,' ) CALL domain_clock_set( grid, & time_step_seconds=model_config_rec%interval_seconds ) CALL domain_clockprint ( 150, grid, & 'DEBUG real: clock after timeStep set,' ) CALL wrf_debug ( 100 , 'tc_em: calling set_scalar_indices_from_config ' ) CALL set_scalar_indices_from_config ( grid%id , idum1, idum2 ) CALL wrf_debug ( 100 , 'tc_em: calling model_to_grid_config_rec ' ) lonc_loc(:) = -999. latc_loc(:) = -999. vmax(:) = -999. rmax(:) = -999. CALL model_to_grid_config_rec ( grid%id , model_config_rec , config_flags ) lonc_loc(1) = config_flags%lonc_loc latc_loc(1) = config_flags%latc_loc vmax(1) = config_flags%vmax_meters_per_second rmax(1) = config_flags%rmax rankine_lid = config_flags%rankine_lid do storm = 2,config_flags%num_storm bogus_id = storm CALL model_to_grid_config_rec ( bogus_id , model_config_rec , config_flags ) lonc_loc(storm) = config_flags%lonc_loc latc_loc(storm) = config_flags%latc_loc vmax(storm) = config_flags%vmax_meters_per_second rmax(storm) = config_flags%rmax end do CALL model_to_grid_config_rec ( grid%id , model_config_rec , config_flags ) CALL wrf_debug ( 100 , 'tc_em: calling init_wrfio' ) CALL init_wrfio CALL wrf_debug ( 100 , 'tc_em: re-broadcast the configuration records' ) CALL get_config_as_buffer( configbuf, configbuflen, nbytes ) CALL wrf_dm_bcast_bytes( configbuf, nbytes ) CALL set_config_as_buffer( configbuf, configbuflen ) CALL wrf_debug ( 100 , 'calling tc_med_sidata_input' ) CALL tc_med_sidata_input ( grid , config_flags, latc_loc, lonc_loc, & vmax,rmax,rankine_lid) CALL wrf_debug ( 100 , 'backfrom tc_med_sidata_input' ) ELSE CYCLE all_domains END IF END DO all_domains CALL set_current_grid_ptr( head_grid ) CALL wrf_debug ( 0 , 'tc_em: SUCCESS COMPLETE TC BOGUS' ) CALL wrf_shutdown CALL WRFU_Finalize( rc=rc ) END PROGRAM tc_data SUBROUTINE tc_med_sidata_input ( grid , config_flags, latc_loc, lonc_loc, & vmax, rmax,rankine_lid) USE module_domain USE module_io_domain USE module_configure USE module_bc_time_utilities USE module_optional_input USE module_date_time USE module_utility IMPLICIT NONE INTERFACE SUBROUTINE start_domain ( grid , allowed_to_read ) USE module_domain TYPE (domain) grid LOGICAL, INTENT(IN) :: allowed_to_read END SUBROUTINE start_domain END INTERFACE TYPE(domain) :: grid TYPE (grid_config_rec_type) :: config_flags INTEGER :: time_step_begin_restart INTEGER :: idsi , ierr , myproc, internal_time_loop,iflag INTEGER ::nf_inq CHARACTER (LEN=256) :: si_inpname CHARACTER (LEN=80) :: message CHARACTER(LEN=19) :: start_date_char , end_date_char , current_date_char , next_date_char CHARACTER(LEN=8) :: flag_name INTEGER :: time_loop_max , loop, rc,icnt,itmp INTEGER :: julyr , julday ,metndims, metnvars, metngatts, nunlimdimid,rcode REAL :: gmt real :: t1,t2,t3,t4 real :: latc_loc(max_bogus), lonc_loc(max_bogus) real :: vmax(max_bogus),rmax(max_bogus),rankine_lid grid%input_from_file = .true. grid%input_from_file = .false. CALL tc_compute_si_start ( model_config_rec%start_year (grid%id) , & model_config_rec%start_month (grid%id) , & model_config_rec%start_day (grid%id) , & model_config_rec%start_hour (grid%id) , & model_config_rec%start_minute(grid%id) , & model_config_rec%start_second(grid%id) , & model_config_rec%interval_seconds , & model_config_rec%real_data_init_type , & start_date_char) end_date_char = start_date_char IF ( end_date_char .LT. start_date_char ) THEN CALL wrf_error_fatal3("",288,& 'Ending date in namelist ' // end_date_char // ' prior to beginning date ' // start_date_char ) END IF print *,"the start date char ",start_date_char print *,"the end date char ",end_date_char time_loop_max = 1 CALL domain_clock_set( grid, stop_timestr=end_date_char ) CALL WRFU_ClockStopTimeDisable( grid%domain_clock, rc=rc ) CALL wrf_check_error( WRFU_SUCCESS, rc, & 'WRFU_ClockStopTimeDisable(grid%domain_clock) FAILED', & "tc_em.G" , & 322 ) CALL domain_clockprint ( 150, grid, & 'DEBUG med_sidata_input: clock after stopTime set,' ) current_date_char = start_date_char start_date = start_date_char // '.0000' current_date = start_date CALL nl_set_bdyfrq ( grid%id , REAL(model_config_rec%interval_seconds) ) CALL cpu_time ( t1 ) DO loop = 1 , time_loop_max internal_time_loop = loop IF ( ( grid%id .GT. 1 ) .AND. ( loop .GT. 1 ) .AND. & ( model_config_rec%grid_fdda(grid%id) .EQ. 0 ) .AND. & ( model_config_rec%sst_update .EQ. 0 ) ) EXIT print *,' ' print *,'-----------------------------------------------------------------------------' print *,' ' print '(A,I2,A,A,A,I4,A,I4)' , & ' Domain ',grid%id,': Current date being processed: ',current_date, ', which is loop #',loop,' out of ',time_loop_max CALL geth_julgmt ( config_flags%julyr , config_flags%julday , config_flags%gmt ) print *,'configflags%julyr, %julday, %gmt:',config_flags%julyr, config_flags%julday, config_flags%gmt CALL nl_set_gmt (grid%id, config_flags%gmt) CALL nl_set_julyr (grid%id, config_flags%julyr) CALL nl_set_julday (grid%id, config_flags%julday) CALL cpu_time ( t3 ) WRITE ( wrf_err_message , FMT='(A,A)' )'med_sidata_input: calling open_r_dataset for ', & TRIM(config_flags%auxinput1_inname) CALL wrf_debug ( 100 , wrf_err_message ) IF ( config_flags%auxinput1_inname(1:8) .NE. 'wrf_real' ) THEN CALL construct_filename4a( si_inpname , config_flags%auxinput1_inname , grid%id , 2 , & current_date_char , config_flags%io_form_auxinput1 ) ELSE CALL construct_filename2a( si_inpname , config_flags%auxinput1_inname , grid%id , 2 , & current_date_char ) END IF CALL open_r_dataset ( idsi, TRIM(si_inpname) , grid , config_flags , "DATASET=AUXINPUT1", ierr ) IF ( ierr .NE. 0 ) THEN CALL wrf_error_fatal3("",357,& 'error opening ' // TRIM(si_inpname) // & ' for input; bad date in namelist or file not in directory' ) END IF CALL wrf_debug ( 100 , 'med_sidata_input: calling input_auxinput1' ) CALL input_auxinput1 ( idsi , grid , config_flags , ierr ) WRITE ( wrf_err_message , FMT='(A,I10,A)' ) 'Timing for input ',NINT(t4-t3) ,' s.' CALL wrf_debug( 0, wrf_err_message ) CALL cpu_time ( t3 ) CALL wrf_debug ( 100 , 'med_sidata_input: calling init_module_optional_input' ) CALL init_module_optional_input ( grid , config_flags ) CALL wrf_debug ( 100 , 'med_sidata_input: calling optional_input' ) CALL optional_input ( grid , idsi , config_flags ) flag_name(1:8) = 'SM000010' CALL wrf_get_dom_ti_integer ( idsi, 'FLAG_' // flag_name, itmp, 1, icnt, ierr ) IF ( ierr .EQ. 0 ) THEN grid%flag_sm000010 = 1 end if flag_name(1:8) = 'SM010040' CALL wrf_get_dom_ti_integer ( idsi, 'FLAG_' // flag_name, itmp, 1, icnt, ierr ) IF ( ierr .EQ. 0 ) THEN grid%flag_sm010040 = 1 end if flag_name(1:8) = 'SM040100' CALL wrf_get_dom_ti_integer ( idsi, 'FLAG_' // flag_name, itmp, 1, icnt, ierr ) IF ( ierr .EQ. 0 ) THEN grid%flag_sm040100 = itmp end if flag_name(1:8) = 'SM100200' CALL wrf_get_dom_ti_integer ( idsi, 'FLAG_' // flag_name, itmp, 1, icnt, ierr ) IF ( ierr .EQ. 0 ) THEN grid%flag_sm100200 = itmp end if flag_name(1:8) = 'ST000010' CALL wrf_get_dom_ti_integer ( idsi, 'FLAG_' // flag_name, itmp, 1, icnt, ierr ) IF ( ierr .EQ. 0 ) THEN grid%flag_st000010 = 1 END IF flag_name(1:8) = 'ST010040' CALL wrf_get_dom_ti_integer ( idsi, 'FLAG_' // flag_name, itmp, 1, icnt, ierr ) IF ( ierr .EQ. 0 ) THEN grid%flag_st010040 = 1 END IF flag_name(1:8) = 'ST040100' CALL wrf_get_dom_ti_integer ( idsi, 'FLAG_' // flag_name, itmp, 1, icnt, ierr ) IF ( ierr .EQ. 0 ) THEN grid%flag_st040100 = 1 END IF flag_name(1:8) = 'ST100200' CALL wrf_get_dom_ti_integer ( idsi, 'FLAG_' // flag_name, itmp, 1, icnt, ierr ) IF ( ierr .EQ. 0 ) THEN grid%flag_st100200 = 1 END IF CALL wrf_get_dom_ti_integer ( idsi, 'FLAG_SOIL_LAYERS', itmp, 1, icnt, ierr ) IF ( ierr .EQ. 0 ) THEN grid%flag_soil_layers = 1 END IF CALL close_dataset ( idsi , config_flags , "DATASET=AUXINPUT1" ) CALL cpu_time ( t4 ) CALL cpu_time ( t3 ) already_been_here = .FALSE. CALL model_to_grid_config_rec ( grid%id , model_config_rec , config_flags ) CALL cpu_time ( t3 ) CALL assemble_output ( grid , config_flags , loop , time_loop_max, current_date_char, & latc_loc, lonc_loc, vmax, rmax, rankine_lid,si_inpname) CALL cpu_time ( t4 ) WRITE ( wrf_err_message , FMT='(A,I10,A)' ) 'Timing for output ',NINT(t4-t3) ,' s.' CALL wrf_debug( 0, wrf_err_message ) CALL cpu_time ( t2 ) WRITE ( wrf_err_message , FMT='(A,I4,A,I10,A)' ) 'Timing for loop # ',loop,' = ',NINT(t2-t1) ,' s.' CALL wrf_debug( 0, wrf_err_message ) CALL cpu_time ( t1 ) END DO END SUBROUTINE tc_med_sidata_input SUBROUTINE tc_compute_si_start( & start_year , start_month , start_day , start_hour , start_minute , start_second , & interval_seconds , real_data_init_type , & start_date_char) USE module_date_time IMPLICIT NONE INTEGER :: start_year , start_month , start_day , start_hour , start_minute , start_second INTEGER :: end_year , end_month , end_day , end_hour , end_minute , end_second INTEGER :: interval_seconds , real_data_init_type INTEGER :: time_loop_max , time_loop CHARACTER(LEN=19) :: current_date_char , start_date_char , end_date_char , next_date_char WRITE ( start_date_char , FMT = '(I4.4,"-",I2.2,"-",I2.2,"_",I2.2,":",I2.2,":",I2.2)' ) & start_year,start_month,start_day,start_hour,start_minute,start_second END SUBROUTINE tc_compute_si_start SUBROUTINE assemble_output ( grid , config_flags , loop , time_loop_max,current_date_char, & latc_loc, lonc_loc,vmax,rmax,rankine_lid,si_inpname) USE module_big_step_utilities_em USE module_domain USE module_io_domain USE module_configure USE module_date_time USE module_bc IMPLICIT NONE TYPE(domain) :: grid TYPE (grid_config_rec_type) :: config_flags INTEGER , INTENT(IN) :: loop , time_loop_max real :: vmax(max_bogus),vmax_ratio,rankine_lid real :: rmax(max_bogus),stand_lon,cen_lat,ptop_in_pa real :: latc_loc(max_bogus),lonc_loc(max_bogus) INTEGER :: ijds , ijde , spec_bdy_width INTEGER :: i , j , k , idts,map_proj,remove_only,storms INTEGER :: id1 , interval_seconds , ierr, rc, sst_update, grid_fdda INTEGER , SAVE :: id, id2, id4 CHARACTER (LEN=256) :: tcoutname , bdyname,si_inpname CHARACTER(LEN= 4) :: loop_char CHARACTER(LEN=19) :: current_date_char character *19 :: temp19 character *24 :: temp24 , temp24b real::t1,t2,truelat1,truelat2 spec_bdy_width = model_config_rec%spec_bdy_width interval_seconds = model_config_rec%interval_seconds sst_update = model_config_rec%sst_update grid_fdda = model_config_rec%grid_fdda(grid%id) truelat1 = config_flags%truelat1 truelat2 = config_flags%truelat2 stand_lon = config_flags%stand_lon cen_lat = config_flags%cen_lat map_proj = config_flags%map_proj vmax_ratio = config_flags%vmax_ratio ptop_in_pa = config_flags%p_top_requested remove_only = 0 if(config_flags%remove_storm) then remove_only = 1 end if storms = config_flags%num_storm print *,"number of storms ",config_flags%num_storm call tc_bogus(cen_lat,stand_lon,map_proj,truelat1,truelat2, & grid%dx,grid%e_we,grid%e_sn,grid%num_metgrid_levels,ptop_in_pa, & rankine_lid,latc_loc,lonc_loc,vmax,vmax_ratio,rmax,remove_only, & storms,grid) CALL construct_filename4a( tcoutname , config_flags%auxinput1_outname , grid%id , 2 , & current_date_char , config_flags%io_form_auxinput1 ) print *,"outfile name from construct filename ",tcoutname CALL open_w_dataset ( id1, TRIM(tcoutname) , grid , config_flags ,output_auxinput1,"DATASET=AUXINPUT1",ierr ) IF ( ierr .NE. 0 ) THEN CALL wrf_error_fatal3("",572,& 'tc_em: error opening tc bogus file for writing' ) END IF CALL output_auxinput1( id1, grid , config_flags , ierr ) CALL close_dataset ( id1 , config_flags , "DATASET=AUXINPUT1" ) END SUBROUTINE assemble_output SUBROUTINE tc_bogus(centerlat,stdlon,nproj,truelat1,truelat2,dsm,ew,ns,nz,ptop_in_pa, & rankine_lid,latc_loc,lonc_loc,vmax,vmax_ratio,rmax,remove_only, & storms,grid) USE module_llxy USE module_domain IMPLICIT NONE TYPE(domain) :: grid integer ew,ns,nz integer nproj integer storms,nstrm real :: centerlat,stdlon,conef,truelat1,truelat2,dsm,dx,rankine_lid real :: latc_loc(max_bogus),lonc_loc(max_bogus),vmax(max_bogus),vmax_ratio,rmax(max_bogus) real :: press(ew-1,nz,ns-1),rhmx(nz), vwgt(nz),old_slp(ew-1,ns-1) real, dimension(:,:,:) , allocatable :: u11,v11,t11,rh11,phi11 real, dimension(:,:,:) , allocatable :: u1 , v1 , t1 , rh1 , phi1 real, dimension(ew-1,ns-1) :: lond,terrain,cor,pslx real, dimension(ew,ns-1) :: msfu real, dimension(ew-1,ns) :: msfv real, dimension(ew-1,ns-1) :: msfm CHARACTER*2 jproj LOGICAL :: l_tcbogus real :: r_search,r_vor,beta,devps,humidity_max real :: devpc,const,r_vor2,cnst,alphar,epsilon,vormx , rad , sum_q real :: avg_q ,q_old,ror,q_new,dph,dphx0 real :: rh_max,min_RH_value,ps integer :: vert_variation integer :: i,k,j,kx,remove_only integer :: k00,kfrm ,kto ,k85,n_iter,ew_mvc,ns_mvc,nct,itr integer :: strmci(nz), strmcj(nz) real :: disx,disy,alpha,degran,pie,rovcp,cp REAL :: rho,pprm,phip0,x0,y0,vmx,xico,xjco,xicn,xjcn,p85,xlo,rconst,ew_gcntr,ns_gcntr real :: ptop_in_pa,themax,themin real :: latinc,loninc real :: rtemp,colat0,colat REAL :: q1(ew-1,nz,ns-1), psi1(ew-1,nz,ns-1) TYPE(proj_info) :: proj REAL :: lat1 , lon1 real :: knowni,knownj REAL utcr(ew,nz,ns-1), vtcr(ew-1,nz,ns) REAL utcp(ew,nz,ns-1), vtcp(ew-1,nz,ns) REAL psitc(ew-1,nz,ns-1), psiv(nz) REAL vortc(ew-1,nz,ns-1), vorv(nz) REAL tptc(ew-1,nz,ns-1) REAL phiptc(ew-1,nz,ns-1) REAL uuwork(nz), vvwork(nz), temp2(ew,ns) REAL vort(ew-1,nz,ns-1), div(ew-1,nz,ns-1) REAL vortsv(ew-1,nz,ns-1) REAL theta(ew-1,nz,ns-1), t_reduce(ew-1,nz,ns-1) REAL ug(ew,nz,ns-1), vg(ew-1,nz,ns), vorg(ew-1,nz,ns-1) REAL delpx(ew-1,ns-1) REAL outold(ew-1,ns-1) REAL rd(ew-1,ns-1), ff(ew-1,ns-1) REAL tmp1(ew-1,ns-1), tmp2(ew-1,ns-1) REAL , DIMENSION (ew-1,nz,ns-1) :: t0, t00, rh0, q0, phi0, psi0, chi REAL , DIMENSION (ew-1,nz,ns-1) :: psipos, tpos, psi ,phipos, phip REAL u2(ew,nz,ns-1), v2(ew-1,nz,ns) REAL t2(ew-1,nz,ns-1),z2(ew-1,nz,ns-1) REAL phi2(ew-1,nz,ns-1),rh2(ew-1,nz,ns-1) print *,"the dimensions: north-south = ",ns," east-west =",ew IF (nproj .EQ. 1) THEN jproj = 'LC' print *,"Lambert Conformal projection" ELSE IF (nproj .EQ. 2) THEN jproj = 'ST' ELSE IF (nproj .EQ. 3) THEN jproj = 'ME' print *,"A mercator projection" END IF knowni = 1. knownj = 1. pie = 3.141592653589793 degran = pie/180. rconst = 287.04 min_RH_value = 5.0 cp = 1004.0 rovcp = rconst/cp r_search = 400000.0 r_vor = 300000.0 r_vor2 = r_vor * 4 beta = 0.5 devpc= 40.0 vert_variation = 1 humidity_max = 95.0 alphar = 1.8 latinc = -999. loninc = -999. if(remove_only .eq. 1) then do nstrm=1,storms vmax(nstrm) = 0.1 end do end if dx = dsm lat1 = grid%xlat_gc(1,1) lon1 = grid%xlong_gc(1,1) IF( jproj .EQ. 'ME' )THEN IF ( lon1 .LT. -180. ) lon1 = lon1 + 360. IF ( lon1 .GT. 180. ) lon1 = lon1 - 360. IF ( stdlon .LT. -180. ) stdlon = stdlon + 360. IF ( stdlon .GT. 180. ) stdlon = stdlon - 360. CALL map_set ( proj_merc, proj, lat1, lon1, lat1, lon1, knowni, knownj, dx, & latinc,loninc,stdlon , truelat1 , truelat2) conef = 0. ELSE IF ( jproj .EQ. 'LC' ) THEN if((truelat1 .eq. 0.0) .and. (truelat2 .eq. 0.0)) then print *,"Truelat1 and Truelat2 are both 0" stop end if CALL map_set (proj_lc,proj, lat1, lon1, lat1, lon1, knowni, knownj, dx, & latinc,loninc,stdlon , truelat1 , truelat2) conef = proj%cone ELSE IF ( jproj .EQ. 'ST' ) THEN conef = 1. CALL map_set ( proj_ps,proj,lat1, lon1, lat1, lon1, knowni, knownj, dx, & latinc,loninc,stdlon , truelat1 , truelat2) END IF kx = nz do j = 1,ns-1 do k = 1,nz do i = 1,ew-1 press(i,k,j) = grid%p_gc(i,k,j)*0.01 end do end do end do IF ( ( ptop_in_pa .EQ. 40000. ) .OR. ( ptop_in_pa .EQ. 60000. ) ) THEN PRINT '(A)','Hold on pardner, your value for PTOP is gonna cause problems for the TC bogus option.' PRINT '(A)','Make it higher up than 400 mb.' STOP 'ptop_woes_for_tc_bogus' END IF IF ( vert_variation .EQ. 1 ) THEN DO k=1,kx IF ( press(1,k,1) .GT. 400. ) THEN rhmx(k) = humidity_max ELSE rhmx(k) = humidity_max * MAX( 0.1 , (press(1,k,1) - ptop_in_pa/100.)/(400.-ptop_in_pa/100.) ) END IF IF ( press(1,k,1) .GT. 600. ) THEN vwgt(k) = 1.0 ELSE IF ( press(1,k,1) .LE. 100. ) THEN vwgt(k) = 0.0001 ELSE vwgt(k) = MAX ( 0.0001 , (press(1,k,1)-ptop_in_pa/100.)/(600.-ptop_in_pa/100.) ) END IF END DO ELSE IF ( vert_variation .EQ. 2 ) THEN IF ( kx .eq. 24 ) THEN rhmx = (/ 95., 95., 95., 95., 95., 95., 95., 95., & 95., 95., 95., 95., 95., 90., 85., 80., 75., & 70., 66., 60., 39., 10., 10., 10./) vwgt = (/ 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 0.9850, & 0.9680, 0.9500, 0.9290, 0.9060, 0.8810, 0.8500, 0.7580, 0.6500, 0.5100, & 0.3500, 0.2120, 0.0500, 0.0270, 0.0001, 0.0001, 0.0001/) ELSE PRINT '(A)','Number of vertical levels assumed to be 24 for AFWA TC bogus option' STOP 'AFWA_TC_BOGUS_LEVEL_ERROR' END IF END IF allocate(u11 (1:ew, 1:nz, 1:ns-1)) allocate(u1 (1:ew, 1:nz, 1:ns-1)) allocate(v11 (1:ew-1, 1:nz, 1:ns)) allocate(v1 (1:ew-1, 1:nz, 1:ns)) do j = 1,ns-1 do k = 1,nz do i = 1,ew u11(i,k,j) = grid%u_gc(i,k,j) u1(i,k,j) = grid%u_gc(i,k,j) msfu(i,j) = grid%msfu(i,j) end do end do end do do j = 1,ns do k = 1,nz do i = 1,ew-1 v11(i,k,j) = grid%v_gc(i,k,j) v1(i,k,j) = grid%v_gc(i,k,j) msfv(i,j) = grid%msfv(i,j) end do end do end do allocate(t11 (1:ew-1, 1:nz, 1:ns-1)) allocate(t1 (1:ew-1, 1:nz, 1:ns-1)) allocate(rh11 (1:ew-1, 1:nz, 1:ns-1)) allocate(rh1 (1:ew-1, 1:nz, 1:ns-1)) allocate(phi11(1:ew-1, 1:nz, 1:ns-1)) allocate(phi1 (1:ew-1, 1:nz, 1:ns-1)) do j = 1,ns-1 do k = 1,nz do i = 1,ew-1 t11(i,k,j) = grid%t_gc(i,k,j) t1(i,k,j) = grid%t_gc(i,k,j) rh11(i,k,j) = grid%rh_gc(i,k,j) rh1(i,k,j) = grid%rh_gc(i,k,j) msfm(i,j) = grid%msft(i,j) if(k .eq. 1)then phi11(i,k,j) = grid%ht_gc(i,j) phi1(i,k,j) = grid%ht_gc(i,j) * 9.81 else phi11(i,k,j) = grid%ght_gc(i,k,j) phi1(i,k,j) = grid%ght_gc(i,k,j) * 9.81 end if end do end do end do do j = 1,ns-1 do i = 1,ew-1 pslx(i,j) = grid%pslv_gc(i,j) * 0.01 cor(i,j) = grid%f(i,j) lond(i,j) = grid%xlong_gc(i,j) terrain(i,j) = grid%ht_gc(i,j) old_slp(i,j) = grid%pslv_gc(i,j) end do end do l_tcbogus = .FALSE. all_storms : DO nstrm=1,storms if(rmax(nstrm) .eq. -999.) then print *,"Please enter a value for rmax in the namelist" stop end if k00 = 2 kfrm = k00 p85 = 850. kto = kfrm DO k=kfrm+1,kx IF ( press(1,k,1) .GE. p85 ) THEN kto = kto + 1 END IF END DO k85 = kto rho = 1.2 pprm = devpc*100. phip0= pprm/rho CALL latlon_to_ij ( proj , latc_loc(nstrm) , lonc_loc(nstrm) , x0 , y0 ) IF ( ( x0 .LT. 1. ) .OR. ( x0 .GT. REAL(ew-1) ) .OR. & ( y0 .LT. 1. ) .OR. ( y0 .GT. REAL(ns-1) ) ) THEN PRINT '(A,I3,A,A,A)',' Storm position is outside the computational domain.' PRINT '(A,2F6.2,A)' ,' Storm postion: (x,y) = ',x0,y0,'.' stop END IF l_tcbogus = .TRUE. vmx = vmax(nstrm) * vmax_ratio IF ( latc_loc(nstrm) .LT. 0. ) THEN vmx = -vmx END IF IF ( vmax(nstrm) .LE. 0. ) THEN vmx = SQRT( 2.*(1-beta)*ABS(phip0) ) END IF ew_gcntr = x0 ns_gcntr = y0 ew_gcntr = x0 + 0.5 ns_gcntr = y0 + 0.5 n_iter = 1 PRINT '(/,A,I3,A,A,A)' ,'---> TC: Processing storm number= ',nstrm PRINT '(A,F6.2,A,F7.2,A)' ,' Storm center lat= ',latc_loc(nstrm),' lon= ',lonc_loc(nstrm),'.' PRINT '(A,2F6.2,A)' ,' Storm center grid position (x,y)= ',ew_gcntr,ns_gcntr,'.' PRINT '(A,F5.2,F9.2,A)' ,' Storm max wind (m/s) and max radius (m)= ',vmx,rmax(nstrm),'.' PRINT '(A,F5.2,A)' ,' Estimated central press dev (mb)= ',devpc,'.' DO k=1,kx strmci(k) = 1 strmcj(k) = 1 END DO utcp(:,:,:) = 0.0 vtcp(:,:,:) = 0.0 print *,"nstrm ",rmax(nstrm),ew_gcntr,ns_gcntr DO j=1,ns-1 DO i=1,ew-1 disx = REAL(i) - ew_gcntr disy = REAL(j) - ns_gcntr CALL rankine(disx,disy,dx,kx,vwgt,rmax(nstrm),vmx,uuwork,vvwork,psiv,vorv) DO k=1,kx utcp(i,k,j) = uuwork(k) vtcp(i,k,j) = vvwork(k) psitc(i,k,j) = psiv(k) vortc(i,k,j) = vorv(k) END DO END DO END DO call stagger_rankine_winds(utcp,vtcp,ew,ns,nz) utcr(:,:,:) = 0.0 vtcr(:,:,:) = 0.0 DO j=1,ns-1 DO i=2,ew-1 xlo = stdlon-grid%xlong_u(i,j) IF ( xlo .GT. 180.)xlo = xlo-360. IF ( xlo .LT.-180.)xlo = xlo+360. alpha = xlo*conef*degran*SIGN(1.,centerlat) DO k=1,kx utcr(i,k,j) = (vtcp(i-1,k,j)+vtcp(i,k,j)+vtcp(i,k,j+1)+vtcp(i-1,k,j+1))/4 *SIN(alpha)+utcp(i,k,j)*COS(alpha) if(utcr(i,k,j) .gt. 300.) then print *,i,k,j,"a very bad value of utcr" stop end if END DO END DO END DO DO j=2,ns-1 DO i=1,ew-1 xlo = stdlon-grid%xlong_v(i,j) IF ( xlo .GT. 180.)xlo = xlo-360. IF ( xlo .LT.-180.)xlo = xlo+360. alpha = xlo*conef*degran*SIGN(1.,centerlat) DO k=1,kx vtcr(i,k,j) = vtcp(i,k,j)*COS(alpha)-(utcp(i,k,j-1)+utcp(i+1,k,j-1)+utcp(i+1,k,j)+utcp(i,k,j))/4*SIN(alpha) if(vtcr(i,k,j) .gt. 300.) then print *,i,k,j,"a very bad value of vtcr" stop end if END DO END DO END DO do j = 1,ns-1 utcr(1,:,j) = utcr(2,:,j) utcr(ew,:,j) = utcr(ew-1,:,j) end do do i = 1,ew-1 vtcr(i,:,1) = vtcr(i,:,2) vtcr(i,:,ns) = vtcr(i,:,ns-1) end do CALL vor(u1,v1,msfu,msfv,msfm,ew,ns,kx,dx,vort) CALL diverg(u1,v1,msfu,msfv,msfm,ew,ns,kx,dx,div) CALL mxratprs(rh1,t1,press*100.,ew,ns,kx,q1,min_RH_value) q1(:,1,:) = q1(:,2,:) vortsv = vort q0 = q1 DO k=1,kx DO j=1,ns-1 DO i=1,ew-1 ff(i,j) = vort(i,k,j) tmp1(i,j)= 0.0 END DO END DO epsilon = 1.E-2 CALL relax(tmp1,ff,rd,ew,ns,dx,epsilon,alphar) DO j=1,ns-1 DO i=1,ew-1 psi1(i,k,j) = tmp1(i,j) END DO END DO END DO DO k=1,kx IF ( latc_loc(nstrm) .GE. 0. ) THEN vormx = -1.e10 ELSE vormx = 1.e10 END IF ew_mvc = 1 ns_mvc = 1 DO j=1,ns-1 DO i=1,ew-1 rad = SQRT((REAL(i)-ew_gcntr)**2.+(REAL(j)-ns_gcntr)**2.)*dx IF ( rad .LE. r_search ) THEN IF ( latc_loc(nstrm) .GE. 0. ) THEN IF ( vortsv(i,k,j) .GT. vormx ) THEN vormx = vortsv(i,k,j) ew_mvc = i ns_mvc = j END IF ELSE IF (latc_loc(nstrm) .LT. 0. ) THEN IF ( vortsv(i,k,j) .LT. vormx ) THEN vormx = vortsv(i,k,j) ew_mvc = i ns_mvc = j END IF END IF END IF END DO END DO strmci(k) = ew_mvc strmcj(k) = ns_mvc DO j=1,ns-1 DO i=1,ew-1 rad = SQRT(REAL((i-ew_mvc)**2.+(j-ns_mvc)**2.))*dx IF ( rad .GT. r_vor ) THEN vort(i,k,j) = 0. div(i,k,j) = 0. END IF END DO END DO DO itr=1,n_iter sum_q = 0. nct = 0 DO j=1,ns-1 DO i=1,ew-1 rad = SQRT(REAL(i-ew_mvc)**2.+REAL(j-ns_mvc)**2.)*dx IF ( (rad .LT. r_vor2).AND.(rad .GE. 0.8*r_vor2) ) THEN sum_q = sum_q + q0(i,k,j) nct = nct + 1 END IF END DO END DO avg_q = sum_q/MAX(REAL(nct),1.) DO j=1,ns-1 DO i=1,ew-1 q_old = q0(i,k,j) rad = SQRT(REAL(i-ew_mvc)**2.+REAL(j-ns_mvc)**2.)*dx IF ( rad .LT. r_vor2 ) THEN ror = rad/r_vor2 q_new = ((1.-ror)*avg_q) + (ror*q_old) q0(i,k,j) = q_new END IF END DO END DO END DO END DO DO k=1,kx DO j=1,ns-1 DO i=1,ew-1 ff(i,j) = div(i,k,j) tmp1(i,j)= 0.0 END DO END DO epsilon = 1.e-2 CALL relax(tmp1,ff,rd,ew,ns,dx,epsilon,alphar) DO j=1,ns-1 DO i=1,ew-1 chi(i,k,j) = tmp1(i,j) END DO END DO END DO DO k=1,kx DO j=1,ns-1 DO i=1,ew-1 ff(i,j)=vort(i,k,j) tmp1(i,j)=0.0 END DO END DO epsilon = 1.e-2 CALL relax(tmp1,ff,rd,ew,ns,dx,epsilon,alphar) DO j=1,ns-1 DO i=1,ew-1 psi(i,k,j)=tmp1(i,j) END DO END DO END DO call final_ew_velocity(u2,u1,chi,psi,utcr,dx,ew,ns,nz) call final_ns_velocity(v2,v1,chi,psi,vtcr,dx,ew,ns,nz) DO k=1,kx DO j=1,ns-1 DO i=1,ew-1 psi0(i,k,j) = psi1(i,k,j)-psi(i,k,j) END DO END DO END DO DO k=k00,kx DO j=1,ns-1 DO i=1,ew-1 psipos(i,k,j)=psi(i,k,j) END DO END DO END DO CALL geowind(phi1,ew,ns,kx,dx,ug,vg) CALL vor(ug,vg,msfu,msfv,msfm,ew,ns,kx,dx,vorg) DO k=1,kx ew_mvc = strmci(k) ns_mvc = strmcj(k) DO j=1,ns-1 DO i=1,ew-1 rad = SQRT(REAL(i-ew_mvc)**2.+REAL(j-ns_mvc)**2.)*dx IF ( rad .GT. r_vor ) THEN vorg(i,k,j) = 0. END IF END DO END DO END DO DO k=k00,kx DO j=1,ns-1 DO i=1,ew-1 ff(i,j) = vorg(i,k,j) tmp1(i,j)= 0.0 END DO END DO epsilon = 1.e-3 CALL relax(tmp1,ff,rd,ew,ns,dx,epsilon,alphar) DO j=1,ns-1 DO i=1,ew-1 phip(i,k,j) = tmp1(i,j) END DO END DO END DO DO k=k00,kx DO j=1,ns-1 DO i=1,ew-1 phi0(i,k,j) = phi1(i,k,j) - phip(i,k,j) END DO END DO END DO DO k=k00,kx DO j=1,ns-1 DO i=1,ew-1 IF( k .EQ. 2 ) THEN tpos(i,k,j) = (-1./rconst)*(phip(i,k+1,j)-phip(i,k,j ))/LOG(press(i,k+1,j)/press(i,k,j)) ELSE IF ( k .EQ. kx ) THEN tpos(i,k,j) = (-1./rconst)*(phip(i,k ,j)-phip(i,k-1,j))/LOG(press(i,k,j )/press(i,k-1,j)) ELSE tpos(i,k,j) = (-1./rconst)*(phip(i,k+1,j)-phip(i,k-1,j))/LOG(press(i,k+1,j)/press(i,k-1,j)) END IF t0(i,k,j) = t1(i,k,j)-tpos(i,k,j) t00(i,k,j) = t0(i,k,j) if(t0(i,k,j) .gt. 400) then print *,"interesting temperature ",t0(i,k,j)," at ",i,j,k stop end if END DO END DO END DO CALL qvtorh (q0,t0,press*100.,k00,ew,ns,kx,rh0,min_RH_value) call final_RH(rh2,rh0,rhmx,strmci,strmcj,rmax(nstrm),ew,ns,nz,k00,dx,ew_gcntr,ns_gcntr,r_vor2) DO k=k00,kx DO j=1,ns-1 DO i=1,ew-1 theta(i,k,j) = t1(i,k,j)*(1000./press(i,k,j))**rovcp END DO END DO END DO ew_mvc = strmci(k00) ns_mvc = strmcj(k00) DO k=kfrm,kto DO j=1,ns-1 DO i=1,ew-1 rad = SQRT(REAL(i-ew_mvc)**2.+REAL(j-ns_mvc)**2.)*dx IF ( rad .LT. r_vor2 ) THEN t_reduce(i,k,j) = theta(i,k85,j)-0.03*(press(i,k,j)-press(i,k85,j)) t0(i,k,j) = t00(i,k,j)*(rad/r_vor2) + (((press(i,k,j)/1000.)**rovcp)*t_reduce(i,k,j))*(1.-(rad/r_vor2)) END IF END DO END DO END DO DO k=1,kx DO j=1,ns-1 DO i=1,ew-1 tmp1(i,j)=psitc(i,k,j) END DO END DO CALL balance(cor,tmp1,ew,ns,dx,outold) DO j=1,ns-1 DO i=1,ew-1 ff(i,j)=outold(i,j) tmp1(i,j)=0.0 END DO END DO epsilon = 1.e-3 CALL relax (tmp1,ff,rd,ew,ns,dx,epsilon,alphar) DO j=1,ns-1 DO i=1,ew-1 phiptc(i,k,j) = tmp1(i,j) END DO END DO END DO DO j=1,ns-1 DO k=1,kx DO i=1,ew-1 phi2(i,k,j) = phi0(i,k,j) + phiptc(i,k,j) END DO END DO END DO DO j=1,ns-1 DO k=k00,kx DO i=1,ew-1 IF( k .EQ. 2 ) THEN tptc(i,k,j)=(-1./rconst)*(phiptc(i,k+1,j)-phiptc(i,k,j ))/LOG(press(i,k+1,j)/press(i,k,j)) ELSE IF ( k .EQ. kx ) THEN tptc(i,k,j)=(-1./rconst)*(phiptc(i,k,j )-phiptc(i,k-1,j))/LOG(press(i,k,j)/press(i,k-1,j)) ELSE tptc(i,k,j)=(-1./rconst)*(phiptc(i,k+1,j)-phiptc(i,k-1,j))/LOG(press(i,k+1,j)/press(i,k-1,j)) END IF t2(i,k,j) = t0(i,k,j) + tptc(i,k,j) if(t2(i,k,j) .gt. 400) then print *,"interesting temperature " print *,t2(i,k,j),i,k,j,tptc(i,k,j) stop end if END DO END DO END DO DO j=1,ns-1 DO i=1,ew-1 dph = phi2(i,k00,j)-phi1(i,k00,j) delpx(i,j) = rho*dph*0.01 END DO END DO DO j=1,ns-1 DO i=1,ew-1 pslx(i,j) = pslx(i,j)+delpx(i,j) grid%pslv_gc(i,j) = pslx(i,j) * 100. END DO END DO DO j=1,ns-1 DO i=1,ew-1 z2(i,1,j) = terrain(i,j) END DO END DO DO j=1,ns-1 DO k=k00,kx DO i=1,ew-1 z2(i,k,j) = phi2(i,k,j)/9.81 END DO END DO END DO DO j=1,ns-1 DO i=1,ew-1 ps = pslx(i,j) t2(i,1,j) = t2(i,k00,j)*((ps/1000.)**rovcp) if(t2(i,1,j) .gt. 400) then print *,"Interesting surface temperature" print *,t2(i,1,j),t2(i,k00,j),ps,i,j stop end if END DO END DO DO j=1,ns-1 DO i=1,ew-1 rh2(i,1,j) = rh2(i,k00,j) END DO END DO PRINT '(A,I3,A)' ,' Bogus storm number ',nstrm,' completed.' do j = 1,ns-1 do k = 1,nz do i = 1,ew u1(i,k,j) = u2(i,k,j) grid%u_gc(i,k,j) = u2(i,k,j) end do end do end do do j = 1,ns do k = 1,nz do i = 1,ew-1 v1(i,k,j) = v2(i,k,j) grid%v_gc(i,k,j) = v2(i,k,j) end do end do end do do j = 1,ns-1 do k = 1,nz do i = 1,ew-1 t1(i,k,j) = t2(i,k,j) grid%t_gc(i,k,j) = t2(i,k,j) rh1(i,k,j) = rh2(i,k,j) grid%rh_gc(i,k,j) = rh2(i,k,j) phi1(i,k,j) = phi2(i,k,j) grid%ght_gc(i,k,j) = z2(i,k,j) END DO END DO END DO END DO all_storms deallocate(u11) deallocate(v11) deallocate(t11) deallocate(rh11) deallocate(phi11) deallocate(u1) deallocate(v1) deallocate(t1) deallocate(rh1) deallocate(phi1) do j = 1,ns-1 do i = 1,ew-1 if(grid%ht_gc(i,j) .gt. 1) then grid%p_gc(i,1,j) = grid%p_gc(i,1,j) + (pslx(i,j) * 100. - old_slp(i,j)) grid%psfc(i,j) = grid%psfc(i,j) + (pslx(i,j) * 100. - old_slp(i,j)) else grid%p_gc(i,1,j) = pslx(i,j) * 100. grid%psfc(i,j) = pslx(i,j) * 100. end if end do end do END SUBROUTINE tc_bogus SUBROUTINE rankine(dx,dy,ds,nlvl,vwgt,rmax,vmax,uu,vv,psi,vor) IMPLICIT NONE INTEGER nlvl REAL , DIMENSION(nlvl) :: uu, vv, psi, vor REAL , DIMENSION(nlvl) :: vwgt REAL :: dx,dy,ds,rmax,vmax REAL , PARAMETER :: alpha1= 1. REAL , PARAMETER :: alpha2= -0.75 real :: pi INTEGER :: k REAL :: vr , ang , rr , term1 , bb , term2 , alpha pi = 3.141592653589793 DO k=1,nlvl rr = SQRT(dx**2+dy**2)*ds IF ( rr .LT. rmax ) THEN alpha = 1. ELSE IF ( rr .GE. rmax ) THEN alpha = alpha2 END IF vr = vmax * (rr/rmax)**(alpha) IF ( dx.GE.0. ) THEN ang = (pi/2.) - ATAN2(dy,MAX(dx,1.e-6)) uu(k) = vwgt(k)*(-vr*COS(ang)) vv(k) = vwgt(k)*( vr*SIN(ang)) ELSE IF ( dx.LT.0. ) THEN ang = ((3.*pi)/2.) + ATAN2(dy,dx) uu(k) = vwgt(k)*(-vr*COS(ang)) vv(k) = vwgt(k)*(-vr*SIN(ang)) END IF END DO DO k=1,nlvl rr = SQRT(dx**2+dy**2)*ds IF ( rr .LT. rmax ) THEN psi(k) = vwgt(k) * (vmax*rr*rr)/(2.*rmax) ELSE IF ( rr .GE. rmax ) THEN IF (alpha1.EQ.1.0 .AND. alpha2.eq.-1.0) THEN psi(k) = vwgt(k) * vmax*rmax*(0.5+LOG(rr/rmax)) ELSE IF (alpha1.EQ.1.0 .AND. alpha2.NE.-1.0) THEN term1 = vmax/(rmax**alpha1)*(rmax**(alpha1+1)/(alpha1+1)) bb = (rr**(alpha2+1)/(alpha2+1))-(rmax**(alpha2+1))/(alpha2+1) term2 = vmax/(rmax**alpha2)*bb psi(k) = vwgt(k) * (term1 + term2) END IF END IF END DO DO k=1,nlvl rr = SQRT(dx**2+dy**2)*ds IF ( rr .LT. rmax ) THEN vor(k) = vwgt(k) * (2.*vmax)/rmax ELSE IF ( rr .GE. rmax ) THEN vor(k) = vwgt(k) * ( (vmax/rmax**alpha2)*(rr**(alpha2-1.))*(1.+alpha2) ) END IF END DO END SUBROUTINE rankine SUBROUTINE vor(uin,vin,msfu,msfv,msfm,ew,ns,nz,ds,vort) IMPLICIT NONE INTEGER :: jp1,jm1,ip1,im1,i,j,k INTEGER :: ns, ew, nz, k1 REAL , DIMENSION(ew,nz,ns-1) :: uin REAL , DIMENSION(ew-1,nz,ns) :: vin REAL , DIMENSION(ew-1,nz,ns-1) :: vort REAL , DIMENSION(ew,ns-1) :: msfu REAL , DIMENSION(ew-1,ns) :: msfv REAL , DIMENSION(ew-1,ns-1) :: msfm real :: u(ew,ns-1),v(ew-1,ns) REAL :: ds REAL :: dsx,dsy , u1 , u2 , u3 , u4 , v1 , v2 , v3 , v4 real :: dudy,dvdx,mm vort(:,:,:) = -999. do k = 1,nz do j = 1,ns-1 do i = 1,ew u(i,j) = uin(i,k,j) end do end do do j = 1,ns do i = 1,ew-1 v(i,j) = vin(i,k,j) end do end do do j = 2,ns-2 do i = 2,ew-2 mm = msfm(i,j) * msfm(i,j) u1 = u(i ,j-1)/msfu(i ,j-1) u2 = u(i+1,j-1)/msfu(i+1,j-1) u3 = u(i+1,j+1)/msfu(i+1,j+1) u4 = u(i ,j+1)/msfu(i ,j+1) dudy = mm * (u4 + u3 -(u1 + u2)) /(4*ds) v1 = v(i-1,j )/msfv(i-1,j) v2 = v(i+1,j )/msfv(i+1,j) v3 = v(i-1 ,j+1)/msfv(i-1,j+1) v4 = v(i+1,j+1)/msfv(i+1,j+1) dvdx = mm * (v4 + v2 - (v1 + v3))/(4*ds) vort(i,k,j) = dvdx - dudy end do end do do i = 2,ew-2 vort(i,k,1) = vort(i,k,2) vort(i,k,ns-1) = vort(i,k,ns-2) end do do j = 1,ns-1 vort(ew-1,k,j) = vort(ew-2,k,j) vort(1,k,j) = vort(2,k,j) end do end do END SUBROUTINE SUBROUTINE diverg(uin,vin,msfu,msfv,msfm,ew,ns,nz,ds,div) IMPLICIT NONE INTEGER :: jp1,jm1,ip1,im1,i,j,k INTEGER :: ns, ew, nz, k1 REAL , DIMENSION(ew,nz,ns-1) :: uin REAL , DIMENSION(ew-1,nz,ns) :: vin REAL , DIMENSION(ew-1,nz,ns-1) :: div REAL , DIMENSION(ew,ns-1) :: msfu REAL , DIMENSION(ew-1,ns) :: msfv REAL , DIMENSION(ew-1,ns-1) :: msfm real :: u(ew,ns-1),v(ew-1,ns) REAL :: ds REAL :: dsr,u1,u2,v1,v2 real :: dudx,dvdy,mm,arg1,arg2 dsr = 1/ds do k = 1,nz do j = 1,ns-1 do i = 1,ew u(i,j) = uin(i,k,j) end do end do do j = 1,ns do i = 1,ew-1 v(i,j) = vin(i,k,j) end do end do do j = 2,ns-2 do i = 2,ew-2 mm = msfm(i,j) * msfm(i,j) u1 = u(i+1,j)/msfu(i+1,j) u2 = u(i ,j)/msfu(i ,j) v1 = v(i,j+1)/msfv(i,j+1) v2 = v(i,j) /msfv(i,j) div(i,k,j) = mm * (u1 - u2 + v1 - v2) * dsr end do end do do i = 2,ew-2 div(i,k,1) = div(i,k,2) div(i,k,ns-1) = div(i,k,ns-2) end do do j = 1,ns-1 div(ew-1,k,j) = div(ew-2,k,j) div(1,k,j) = div(2,k,j) end do end do END SUBROUTINE diverg SUBROUTINE mxratprs (rh, t, ppa, ew, ns, nz, q, min_RH_value) IMPLICIT NONE INTEGER :: i , ew , j , ns , k , nz REAL :: min_RH_value REAL :: ppa(ew-1,nz,ns-1) REAL :: p( ew-1,nz,ns-1 ) REAL :: q (ew-1,nz,ns-1),rh(ew-1,nz,ns-1),t(ew-1,nz,ns-1) REAL :: es REAL :: qs REAL :: cp = 1004.0 REAL :: svp1,svp2,svp3 REAL :: celkel REAL :: eps p = ppa * 0.01 DO j = 1, ns - 1 DO k = 1, nz DO i = 1, ew - 1 rh(i,k,j) = MIN ( MAX ( rh(i,k,j) ,min_RH_value ) , 100. ) END DO END DO END DO svp3 = 29.65 svp1 = 0.6112 svp2 = 17.67 celkel = 273.15 eps = 0.622 DO j = 1, ns-1 DO k = 1, nz DO i = 1,ew-1 es = svp1 * 10. * EXP(svp2 * (t(i,k,j) - celkel ) / (t(i,k,j) - svp3 )) qs = eps * es / (p(i,k,j) - es) q(i,k,j) = MAX(0.01 * rh(i,k,j) * qs,0.0) END DO END DO END DO END SUBROUTINE mxratprs SUBROUTINE mass2_Ustag(field,dim1,dim2,dim3) IMPLICIT NONE INTEGER :: dim1 , dim2 , dim3 REAL , DIMENSION(dim1,dim2,dim3) :: field,dummy dummy = 0.0 dummy(:,2:dim2-1,:) = ( field(:,1:dim2-2,:) + & field(:,2:dim2-1,:) ) * 0.5 dummy(:,1,:) = field(:,1,:) dummy(:,dim2,:) = field(:,dim2-1,:) field = dummy END SUBROUTINE mass2_Ustag SUBROUTINE mass2_Vstag(field,dim1,dim2,dim3) IMPLICIT NONE INTEGER :: dim1 , dim2 , dim3 REAL , DIMENSION(dim1,dim2,dim3) :: field,dummy dummy = 0.0 dummy(2:dim1-1,:,:) = ( field(1:dim1-2,:,:) + & field(2:dim1-1,:,:) ) * 0.5 dummy(1,:,:) = field(1,:,:) dummy(dim1,:,:) = field(dim1-1,:,:) field = dummy END SUBROUTINE mass2_Vstag SUBROUTINE relax (chi, ff, rd, ew, ns, ds, smallres, alpha) IMPLICIT NONE INTEGER, PARAMETER :: mm = 20000 INTEGER :: i INTEGER :: ie INTEGER :: ew INTEGER :: iter INTEGER :: j INTEGER :: je INTEGER :: jm INTEGER :: ns INTEGER :: mi REAL :: alpha REAL :: alphaov4 REAL :: chi(ew-1,ns-1) REAL :: chimx(ns-1) REAL :: ds REAL :: epx REAL :: fac REAL :: ff(ew-1,ns-1) REAL :: rd(ew-1,ns-1) REAL :: rdmax(ns-1) REAL :: smallres LOGICAL :: converged = .FALSE. fac = ds * ds alphaov4 = alpha * 0.25 ie=ew-2 je=ns-2 DO j = 1, ns-1 DO i = 1, ew-1 ff(i,j) = fac * ff(i,j) rd(i,j) = 0.0 END DO END DO iter_loop : DO iter = 1, mm mi = iter chimx = 0.0 DO j = 2, ns-1 DO i = 2, ew-1 chimx(j) = MAX(ABS(chi(i,j)),chimx(j)) END DO END DO epx = MAXVAL(chimx) * SMALLRES * 4.0 / alpha DO j = 2, ns-2 DO i = 2, ew-2 rd(i,j) = chi(i,j+1) + chi(i,j-1) + chi(i+1,j) + chi(i-1,j) - 4.0 * chi(i,j) - ff(i,j) chi(i,j) = chi(i,j) + rd(i,j) * alphaov4 END DO END DO rdmax = 0.0 DO j = 2, ns-2 DO i = 2, ew-2 rdmax(j) = MAX(ABS(rd(i,j)),rdmax(j)) END DO END DO IF (MAXVAL(rdmax) .lt. epx) THEN converged = .TRUE. EXIT iter_loop END IF END DO iter_loop IF (converged ) THEN ELSE PRINT '(A,I5,A)','Relaxation did not converge in',mm,' iterations.' STOP 'no_converge' END IF do i = 2,ew-2 chi(i,ns-1) = chi(i,ns-2) chi(i,1) = chi(i,2) end do do j = 2,ns-2 chi(ew-1,j) = chi(ew-2,j) chi(1,j) = chi(2,j) end do chi(1,1) = chi(2,1) chi(ew-1,1) = chi(ew-2,1) chi(1,ns-1) = chi(2,ns-1) chi(ew-1,ns-1) = chi(ew-2,ns-1) END SUBROUTINE relax SUBROUTINE geowind(height,ew,ns,nz,ds,ug,vg) IMPLICIT NONE INTEGER :: ew , ns , nz REAL :: ds REAL , DIMENSION(ew-1,nz,ns-1) :: height REAL , DIMENSION(ew,nz,ns-1) :: ug REAL , DIMENSION(ew-1,nz,ns) :: vg REAL :: ds2r , h1 , h2 , h3 , h4, ds4r INTEGER :: i , j , k ds4r=1./(4.*ds) ug(:,:,:) = -999. do j=2,ns-2 do k=1,nz do i=2,ew-1 h1 = height(i,k,j+1) h2 = height(i-1,k,j+1) h3 = height(i ,k,j-1) h4 = height(i-1,k,j-1) ug(i,k,j) = -( (h1 + h2) - ( h3 + h4) ) * ds4r end do end do end do do i = 2,ew-1 ug(i,:,1) = ug(i,:,2) ug(i,:,ns-1) = ug(i,:,ns-2) end do do j = 2,ns-2 ug(1,:,j) = ug(2,:,j) ug(ew,:,j) = ug(ew-1,:,j) end do ug(1,:,1) = ug(2,:,1) ug(1,:,ns-1) = ug(2,:,ns-1) ug(ew,:,1) = ug(ew-1,:,1) ug(ew,:,ns-1) = ug(ew-1,:,ns-1) vg(:,:,:) = -999. DO j=2,ns-1 DO k=1,nz DO i=2,ew-2 h1 = height(i+1,k,j) h2 = height(i-1,k,j) h3 = height(i+1,k,j-1) h4 = height(i-1,k,j-1) vg(i,k,j) = ( (h1 + h3) - ( h2 + h4) ) * ds4r end do end do end do do i = 2,ew-2 vg(i,:,1) = vg(i,:,2) vg(i,:,ns) = vg(i,:,ns-1) end do do j = 2,ns-1 vg(1,:,j) = vg(2,:,j) vg(ew-1,:,j) = vg(ew-2,:,j) end do vg(1,:,1) = vg(2,:,1) vg(1,:,ns) = vg(2,:,ns) vg(ew-1,:,1) = vg(ew-2,:,1) vg(ew-1,:,ns) = vg(ew-2,:,ns) END SUBROUTINE geowind SUBROUTINE balance (f,psi,ew,ns,ds,out) IMPLICIT NONE INTEGER :: ew , ns,nslast,ewlast,ifill REAL , DIMENSION(ew-1,ns-1) :: f,psi,out REAL :: ds REAL :: psixx , psiyy , psiy , psix, psixy REAL :: dssq , ds2 , dssq4,arg1,arg2,arg3,arg4 INTEGER :: i , j dssq = ds * ds ds2 = ds * 2. dssq4 = ds * ds * 4. out(:,:) = -999.0 DO j=2,ns-2 DO i=2,ew-2 psiyy = ( psi(i,j+1) + psi(i,j-1) - 2.*psi(i,j) ) / dssq psixx = ( psi(i+1,j) + psi(i-1,j) - 2.*psi(i,j) ) / dssq psiy = ( psi(i,j+1) - psi(i,j-1) ) / ds2 psix = ( psi(i+1,j) - psi(i-1,j) ) / ds2 psixy = ( psi(i+1,j+1)+psi(i-1,j-1)-psi(i-1,j+1)-psi(i+1,j-1)) / dssq4 arg1 = f(i,j)* (psixx+psiyy) arg2 = ( ( f(i,j+1) - f(i,j-1)) / ds2 ) * psiy arg3 = ( ( f(i+1,j) - f(i-1,j)) / ds2 ) * psix arg4 = 2 *(psixy*psixy-psixx*psiyy) out(i,j)= arg1 + arg2 + arg3 - arg4 END DO END DO do i = 2,ew-2 out(i,ns-1) = out(i,ns-2) out(i,1) = out(i,2) end do do j = 2,ns-2 out(ew-1,j) = out(ew-2,j) out(1,j) = out(2,j) end do out(1,1) = out(2,1) out(ew-1,1) = out(ew-2,1) out(1,ns-1) = out(2,ns-1) out(ew-1,ns-1) = out(ew-2,ns-1) END SUBROUTINE balance SUBROUTINE qvtorh ( q , t , p , k00, ew , ns , nz , rh, min_RH_value ) IMPLICIT NONE INTEGER , INTENT(IN) :: ew , ns , nz , k00 REAL , INTENT(IN) , DIMENSION(ew-1,nz,ns-1) :: q ,t, p REAL , INTENT(OUT) , DIMENSION(ew-1,nz,ns-1) :: rh real min_RH_value INTEGER :: i , j , k,fill REAL :: es REAL :: qs REAL :: cp = 1004.0 REAL :: svp1,svp2,svp3 REAL :: celkel REAL :: eps svp3 = 29.65 svp1 = 0.6112 svp2 = 17.67 celkel = 273.15 eps = 0.622 DO j = 1 , ns - 1 DO k = k00 , nz DO i = 1 , ew -1 es = svp1 * 10. * EXP(svp2 * (t(i,k,j) - celkel ) / (t(i,k,j) - svp3 )) qs = eps*es/(0.01*p(i,k,j) - es) rh(i,k,j) = MIN ( 100. , MAX ( 100.*q(i,k,j)/qs , min_RH_value ) ) END DO END DO END DO END SUBROUTINE qvtorh SUBROUTINE stagger_rankine_winds(utcp,vtcp,ew,ns,nz) INTEGER :: ew, ns, nz, i,k,j REAL utcp(ew,nz,ns-1), vtcp(ew-1,nz,ns) DO j=1,ns-1 DO i=2,ew-1 DO k=1,nz utcp(i,k,j) = ( utcp(i-1,k,j) + utcp(i,k,j) ) /2 end do end do end do do j = 1,ns-1 utcp(1,:,j) = utcp(2,:,j) utcp(ew,:,j) = utcp(ew-1,:,j) end do DO j=2,ns-1 DO i=1,ew-1 DO k=1,nz vtcp(i,k,j) = ( vtcp(i,k,j+1) + vtcp(i,k,j-1) ) /2 end do end do end do do i = 1,ew-1 vtcp(i,:,1) = vtcp(i,:,2) vtcp(i,:,ns) = vtcp(i,:,ns-1) end do END SUBROUTINE stagger_rankine_winds subroutine final_ew_velocity(u2,u1,chi,psi,utcr,dx,ew,ns,nz) integer :: ew,ns,nz,i,j,k real :: u1(ew,nz,ns-1),utcr(ew,nz,ns-1) real :: psi(ew-1,nz,ns-1),chi(ew-1,nz,ns-1) real :: u2(ew,nz,ns-1) real upos(ew,nz,ns-1),u0(ew,nz,ns-1),uchi(ew,nz,ns-1) real :: dx,arg1,arg2 uchi(:,:,:) = -999. DO k=1,nz DO j=1,ns-1 DO i=2,ew-1 uchi(i,k,j) = ( chi(i,k,j) - chi(i-1,k,j) )/dx END DO END DO do j = 1,ns-1 uchi(1,k,j) = uchi(2,k,j) uchi(ew,k,j) = uchi(ew-1,k,j) end do end do upos(:,:,:) = -999. DO k=1,nz DO j=2,ns-2 DO i=2,ew-1 arg1 = psi(i,k,j+1) + psi(i-1,k,j+1) arg2 = psi(i,k,j-1) + psi(i-1,k,j-1) upos(i,k,j) = -( arg1 - arg2 )/(4.*dx) END DO END DO do i = 2,ew-1 upos(i,k,1) = upos(i,k,2) upos(i,k,ns-1) = upos(i,k,ns-2) end do do j = 1,ns-2 upos(1,k,j) = upos(2,k,j) upos(ew,k,j) = upos(ew-1,k,j) end do upos(1,k,1) = upos(2,k,1) upos(1,k,ns-1) = upos(2,k,ns-1) upos(ew,k,1) = upos(ew-1,k,1) upos(ew,k,ns-1) = upos(ew-1,k,ns-1) end do do j=1,ns-1 do k=1,nz do i=1,ew u0(i,k,j) = u1(i,k,j)-(upos(i,k,j)+uchi(i,k,j)) end do end do end do do j=1,ns-1 do k=1,nz do i=1,ew u2(i,k,j) = u0(i,k,j)+utcr(i,k,j) end do end do end do end subroutine final_ew_velocity subroutine final_ns_velocity(v2,v1,chi,psi,vtcr,dx,ew,ns,nz) integer :: ew,ns,nz,i,j,k real :: v1(ew-1,nz,ns),vtcr(ew-1,nz,ns) real :: psi(ew-1,nz,ns-1),chi(ew-1,nz,ns-1) real :: v2(ew-1,nz,ns) real vpos(ew-1,nz,ns),v0(ew-1,nz,ns),vchi(ew-1,nz,ns) real :: dx,arg1,arg2 vchi(:,:,:) = -999.0 do k = 1,nz DO j=2,ns-1 DO i=1,ew-1 vchi(i,k,j) = ( chi(i,k,j) - chi(i,k,j-1))/dx END DO END DO do i = 1,ew-1 vchi(i,k,1) = vchi(i,k,2) vchi(i,k,ns) = vchi(i,k,ns-1) end do end do vpos(:,:,:) = -999. DO k=1,nz DO j=2,ns-1 DO i=2,ew-2 arg1 = psi(i+1,k,j) + psi(i+1,k,j-1) arg2 = psi(i-1,k,j) + psi(i-1,k,j-1) vpos(i,k,j) = ( arg1 - arg2 )/(4.*dx) END DO END DO do i = 2,ew-2 vpos(i,k,1) = vpos(i,k,2) vpos(i,k,ns) = vpos(i,k,ns-1) end do do j = 1,ns vpos(1,k,j) = vpos(2,k,j) vpos(ew-1,k,j) = vpos(ew-2,k,j) end do vpos(1,k,1) = vpos(2,k,1) vpos(1,k,ns) = vpos(2,k,ns) vpos(ew-1,k,1) = vpos(ew-2,k,1) vpos(ew-1,k,ns) = vpos(ew-2,k,ns) END DO do j=1,ns do k=1,nz do i=1,ew-1 v0(i,k,j) = v1(i,k,j)-(vpos(i,k,j)+vchi(i,k,j)) if( v0(i,k,j) .gt. 100.) then print *,vchi(i,k,j),i,k,j stop end if end do end do end do do j=1,ns do k=1,nz do i=1,ew-1 v2(i,k,j) = v0(i,k,j)+vtcr(i,k,j) end do end do end do end subroutine final_ns_velocity subroutine final_RH(rh2,rh0,rhmx,strmci,strmcj,rmax_nstrm,ew,ns,nz,k00, & dx,ew_gcntr,ns_gcntr,r_vor2) integer :: ew,ns,nz real :: rh2(ew-1,nz,ns-1) real :: rh0(ew-1,nz,ns-1) real :: rhmx(nz) real :: ew_gcntr real :: ns_gcntr real :: dx real :: rmax_nstrm real :: sum_rh,avg_rh,rh_min,rhbkg,rhbog,r_ratio real :: rad real :: rhtc(ew-1,nz,ns-1) integer :: nct,k00,i,j,k,ew_mvc,ns_mvc integer :: strmci(nz), strmcj(nz) DO k=k00,nz rh_max= rhmx(k) ew_mvc = strmci(k) ns_mvc = strmcj(k) sum_rh = 0. nct = 0 DO j=1,ns-1 DO i=1,ew-1 rad = SQRT(REAL(i-ew_mvc)**2.+REAL(j-ns_mvc)**2.)*dx IF ( (rad .LT. r_vor2).AND.(rad .GE. 0.8*r_vor2) ) THEN sum_rh = sum_rh + rh0(i,k,j) nct = nct + 1 END IF END DO END DO avg_rh = sum_rh/MAX(REAL(nct),1.) DO j=1,ns-1 DO i=1,ew-1 rh_min = avg_rh rad = SQRT((REAL(i)-ew_gcntr)**2.+(REAL(j)-ns_gcntr)**2.)*dx IF ( rad .LE. rmax_nstrm ) THEN rhtc(i,k,j) = rh_max ELSE rhtc(i,k,j) = (rmax_nstrm/rad)*rh_max+(1.-(rmax_nstrm/rad))*rh_min END IF END DO END DO END DO DO j=1,ns-1 DO k=k00,nz DO i=1,ew-1 rhbkg = rh0(i,k,j) rhbog = rhtc(i,k,j) rad = SQRT((REAL(i)-ew_mvc)**2.+(REAL(j)-ns_mvc)**2.)*dx IF ( (rad.GT.rmax_nstrm) .AND. (rad.LE.r_vor2) ) THEN r_ratio = (rad-rmax_nstrm)/(r_vor2-rmax_nstrm) rh2(i,k,j) = ((1.-r_ratio)*rhbog) + (r_ratio*rhbkg) ELSE IF (rad .LE. rmax_nstrm ) THEN rh2(i,k,j) = rhbog ELSE rh2(i,k,j) = rhbkg END IF END DO END DO END DO end subroutine final_RH