wrfv2_fire/main/tc_em.F

   1 !  Create an initial data set for the WRF model based on real data.  This
   2 !  program is specifically set up for the Eulerian, mass-based coordinate.
   3 PROGRAM tc_data
   4    USE module_machine
   5    USE module_domain, ONLY : domain, alloc_and_configure_domain, &
   6         domain_clock_set, head_grid, program_name, domain_clockprint, &
   7         set_current_grid_ptr
   8    USE module_io_domain
   9    USE module_initialize_real, ONLY : wrfu_initialize
  10    USE module_driver_constants
  11    USE module_configure, ONLY : grid_config_rec_type, model_config_rec, &
  12         initial_config, get_config_as_buffer, set_config_as_buffer
  13    USE module_timing
  14    USE module_state_description, ONLY : realonly
  15 #ifdef NO_LEAP_CALENDAR
  16    USE module_symbols_util, ONLY: wrfu_cal_noleap
  17 #else
  18    USE module_symbols_util, ONLY: wrfu_cal_gregorian
  19 #endif
  20    USE module_utility, ONLY : WRFU_finalize
  21
  22    IMPLICIT NONE
  23
  24
  25    REAL    :: time , bdyfrq
  26
  27    INTEGER :: loop , levels_to_process , debug_level
  28
  29
  30    TYPE(domain) , POINTER :: null_domain
  31    TYPE(domain) , POINTER :: grid , another_grid
  32    TYPE(domain) , POINTER :: grid_ptr , grid_ptr2
  33    TYPE (grid_config_rec_type)              :: config_flags
  34    INTEGER                :: number_at_same_level
  35
  36    INTEGER :: max_dom, domain_id , grid_id , parent_id , parent_id1 , id
  37    INTEGER :: e_we , e_sn , i_parent_start , j_parent_start
  38    INTEGER :: idum1, idum2
  39 #ifdef DM_PARALLEL
  40    INTEGER                 :: nbytes
  41    INTEGER, PARAMETER      :: configbuflen = 4* CONFIG_BUF_LEN
  42    INTEGER                 :: configbuf( configbuflen )
  43    LOGICAL , EXTERNAL      :: wrf_dm_on_monitor
  44 #endif
  45    LOGICAL found_the_id
  46
  47    INTEGER :: ids , ide , jds , jde , kds , kde
  48    INTEGER :: ims , ime , jms , jme , kms , kme
  49    INTEGER :: ips , ipe , jps , jpe , kps , kpe
  50    INTEGER :: ijds , ijde , spec_bdy_width
  51    INTEGER :: i , j , k , idts, rc
  52    INTEGER :: sibling_count , parent_id_hold , dom_loop
  53
  54    CHARACTER (LEN=80)     :: message
  55
  56    INTEGER :: start_year , start_month , start_day , start_hour , start_minute , start_second
  57    INTEGER ::   end_year ,   end_month ,   end_day ,   end_hour ,   end_minute ,   end_second
  58    INTEGER :: interval_seconds , real_data_init_type
  59    INTEGER :: time_loop_max , time_loop, bogus_id, storm
  60    real::t1,t2
  61    real    :: latc_loc(max_bogus),lonc_loc(max_bogus),vmax(max_bogus),rmax(max_bogus)
  62    real    :: rankine_lid
  63    INTERFACE
  64      SUBROUTINE Setup_Timekeeping( grid )
  65       USE module_domain, ONLY : domain
  66       TYPE(domain), POINTER :: grid
  67      END SUBROUTINE Setup_Timekeeping
  68    END INTERFACE
  69
  70 #include "version_decl"
  71
  72    !  Define the name of this program (program_name defined in module_domain)
  73
  74    program_name = "TC_EM " // TRIM(release_version) // " PREPROCESSOR"
  75
  76 !  The TC bogus algorithm assumes that the user defines a central point, and then
  77 !  allows the program to remove a typhoon based on a distance in km.  This is
  78 !  implemented on a single processor only.
  79
  80 #ifdef DM_PARALLEL
  81    IF ( .NOT. wrf_dm_on_monitor() ) THEN
  82       CALL wrf_error_fatal( 'TC bogus must run with a single processor only, re-run with num procs set to 1' )
  83    END IF
  84 #endif
  85
  86 #ifdef DM_PARALLEL
  87    CALL disable_quilting
  88 #endif
  89
  90    !  Initialize the modules used by the WRF system.  Many of the CALLs made from the
  91    !  init_modules routine are NO-OPs.  Typical initializations are: the size of a
  92    !  REAL, setting the file handles to a pre-use value, defining moisture and
  93    !  chemistry indices, etc.
  94
  95    CALL       wrf_debug ( 100 , 'real_em: calling init_modules ' )
  96    CALL init_modules(1)   ! Phase 1 returns after MPI_INIT() (if it is called)
  97 #ifdef NO_LEAP_CALENDAR
  98    CALL WRFU_Initialize( defaultCalKind=WRFU_CAL_NOLEAP, rc=rc )
  99 #else
 100    CALL WRFU_Initialize( defaultCalKind=WRFU_CAL_GREGORIAN, rc=rc )
 101 #endif
 102    CALL init_modules(2)   ! Phase 2 resumes after MPI_INIT() (if it is called)
 103
 104    !  The configuration switches mostly come from the NAMELIST input.
 105
 106 #ifdef DM_PARALLEL
 107    IF ( wrf_dm_on_monitor() ) THEN
 108       CALL initial_config
 109    END IF
 110    CALL get_config_as_buffer( configbuf, configbuflen, nbytes )
 111    CALL wrf_dm_bcast_bytes( configbuf, nbytes )
 112    CALL set_config_as_buffer( configbuf, configbuflen )
 113 !   CALL wrf_dm_initialize
 114 #else
 115    CALL initial_config
 116 #endif
 117
 118
 119    CALL nl_get_debug_level ( 1, debug_level )
 120    CALL set_wrf_debug_level ( debug_level )
 121
 122    CALL  wrf_message ( program_name )
 123
 124    !  There are variables in the Registry that are only required for the real
 125    !  program, fields that come from the WPS package.  We define the run-time
 126    !  flag that says to allocate space for these input-from-WPS-only arrays.
 127
 128    CALL nl_set_use_wps_input ( 1 , REALONLY )
 129
 130    !  Allocate the space for the mother of all domains.
 131
 132    NULLIFY( null_domain )
 133    CALL       wrf_debug ( 100 , 'real_em: calling alloc_and_configure_domain ' )
 134    CALL alloc_and_configure_domain ( domain_id  = 1           , &
 135                                      grid       = head_grid   , &
 136                                      parent     = null_domain , &
 137                                      kid        = -1            )
 138
 139    grid => head_grid
 140    CALL nl_get_max_dom ( 1 , max_dom )
 141
 142    IF ( model_config_rec%interval_seconds .LE. 0 ) THEN
 143      CALL wrf_error_fatal( 'namelist value for interval_seconds must be > 0')
 144    END IF
 145
 146    all_domains : DO domain_id = 1 , max_dom
 147
 148       IF ( ( model_config_rec%input_from_file(domain_id) ) .OR. &
 149            ( domain_id .EQ. 1 ) ) THEN
 150
 151          CALL Setup_Timekeeping ( grid )
 152          CALL set_current_grid_ptr( grid )
 153          CALL domain_clockprint ( 150, grid, &
 154                 'DEBUG real:  clock after Setup_Timekeeping,' )
 155          CALL domain_clock_set( grid, &
 156                                 time_step_seconds=model_config_rec%interval_seconds )
 157          CALL domain_clockprint ( 150, grid, &
 158                 'DEBUG real:  clock after timeStep set,' )
 159
 160
 161          CALL       wrf_debug ( 100 , 'tc_em: calling set_scalar_indices_from_config ' )
 162          CALL set_scalar_indices_from_config ( grid%id , idum1, idum2 )
 163
 164 !This is goofy but we need to loop through the number of storms to get
 165 !the namelist variables for the tc_bogus.  But then we need to
 166 !call model_to_grid_config_rec with the grid%id = to 1 in order to
 167 !reset to the correct information.
 168          CALL       wrf_debug ( 100 , 'tc_em: calling model_to_grid_config_rec ' )
 169          lonc_loc(:) = -999.
 170          latc_loc(:) = -999.
 171          vmax(:)     = -999.
 172          rmax(:)     = -999.
 173          CALL model_to_grid_config_rec ( grid%id , model_config_rec , config_flags )
 174          lonc_loc(1) = config_flags%lonc_loc
 175          latc_loc(1) = config_flags%latc_loc
 176          vmax(1)     = config_flags%vmax_meters_per_second
 177          rmax(1)     = config_flags%rmax
 178          rankine_lid = config_flags%rankine_lid
 179          do storm = 2,config_flags%num_storm
 180              bogus_id = storm
 181              CALL model_to_grid_config_rec ( bogus_id , model_config_rec , config_flags )
 182              lonc_loc(storm) = config_flags%lonc_loc
 183              latc_loc(storm) = config_flags%latc_loc
 184              vmax(storm)     = config_flags%vmax_meters_per_second
 185              rmax(storm)     = config_flags%rmax
 186 !             print *,"in loop ",storm,lonc_loc(storm),latc_loc(storm),vmax(storm),rmax(storm)
 187          end do
 188          CALL model_to_grid_config_rec ( grid%id , model_config_rec , config_flags )
 189
 190          !  Initialize the WRF IO: open files, init file handles, etc.
 191
 192          CALL       wrf_debug ( 100 , 'tc_em: calling init_wrfio' )
 193          CALL init_wrfio
 194
 195          !  Some of the configuration values may have been modified from the initial READ
 196          !  of the NAMELIST, so we re-broadcast the configuration records.
 197
 198 #ifdef DM_PARALLEL
 199          CALL       wrf_debug ( 100 , 'tc_em: re-broadcast the configuration records' )
 200          CALL get_config_as_buffer( configbuf, configbuflen, nbytes )
 201          CALL wrf_dm_bcast_bytes( configbuf, nbytes )
 202          CALL set_config_as_buffer( configbuf, configbuflen )
 203 #endif
 204
 205          !   No looping in this layer.
 206
 207          CALL       wrf_debug ( 100 , 'calling tc_med_sidata_input' )
 208          CALL tc_med_sidata_input ( grid , config_flags, latc_loc, lonc_loc, &
 209                                     vmax,rmax,rankine_lid)
 210          CALL       wrf_debug ( 100 , 'backfrom tc_med_sidata_input' )
 211
 212       ELSE
 213          CYCLE all_domains
 214       END IF
 215
 216    END DO all_domains
 217
 218    CALL set_current_grid_ptr( head_grid )
 219
 220    !  We are done.
 221
 222    CALL       wrf_debug (   0 , 'tc_em: SUCCESS COMPLETE TC BOGUS' )
 223
 224    CALL wrf_shutdown
 225
 226    CALL WRFU_Finalize( rc=rc )
 227
 228
 229 END PROGRAM tc_data
 230
 231
 232 !-----------------------------------------------------------------
 233 SUBROUTINE tc_med_sidata_input ( grid , config_flags, latc_loc, lonc_loc, &
 234                                  vmax, rmax,rankine_lid)
 235   ! Driver layer
 236    USE module_domain
 237    USE module_io_domain
 238   ! Model layer
 239    USE module_configure
 240    USE module_bc_time_utilities
 241    USE module_optional_input
 242
 243    USE module_date_time
 244    USE module_utility
 245
 246    IMPLICIT NONE
 247
 248
 249   ! Interface
 250    INTERFACE
 251      SUBROUTINE start_domain ( grid , allowed_to_read )  ! comes from module_start in appropriate dyn_ directory
 252        USE module_domain
 253        TYPE (domain) grid
 254        LOGICAL, INTENT(IN) :: allowed_to_read
 255      END SUBROUTINE start_domain
 256    END INTERFACE
 257
 258   ! Arguments
 259    TYPE(domain)                :: grid
 260    TYPE (grid_config_rec_type) :: config_flags
 261   ! Local
 262    INTEGER                :: time_step_begin_restart
 263    INTEGER                :: idsi , ierr , myproc, internal_time_loop,iflag
 264 ! Declarations for the netcdf routines.
 265    INTEGER                ::nf_inq
 266 !
 267    CHARACTER (LEN=80)     :: si_inpname
 268    CHARACTER (LEN=80)     :: message
 269
 270    CHARACTER(LEN=19) :: start_date_char , end_date_char , current_date_char , next_date_char
 271    CHARACTER(LEN=8)  :: flag_name
 272
 273    INTEGER :: time_loop_max , loop, rc,icnt,itmp
 274    INTEGER :: julyr , julday ,metndims, metnvars, metngatts, nunlimdimid,rcode
 275    REAL    :: gmt
 276    real    :: t1,t2,t3,t4
 277    real    :: latc_loc(max_bogus), lonc_loc(max_bogus)
 278    real    :: vmax(max_bogus),rmax(max_bogus),rankine_lid
 279
 280    grid%input_from_file = .true.
 281    grid%input_from_file = .false.
 282
 283    CALL tc_compute_si_start ( model_config_rec%start_year  (grid%id) , &
 284                                    model_config_rec%start_month (grid%id) , &
 285                                    model_config_rec%start_day   (grid%id) , &
 286                                    model_config_rec%start_hour  (grid%id) , &
 287                                    model_config_rec%start_minute(grid%id) , &
 288                                    model_config_rec%start_second(grid%id) , &
 289                                    model_config_rec%interval_seconds      , &
 290                                    model_config_rec%real_data_init_type   , &
 291                                    start_date_char)
 292
 293    end_date_char = start_date_char
 294    IF ( end_date_char .LT. start_date_char ) THEN
 295       CALL wrf_error_fatal( 'Ending date in namelist ' // end_date_char // ' prior to beginning date ' // start_date_char )
 296    END IF
 297    print *,"the start date char ",start_date_char
 298    print *,"the end date char ",end_date_char
 299
 300    time_loop_max = 1
 301    !  Override stop time with value computed above.
 302    CALL domain_clock_set( grid, stop_timestr=end_date_char )
 303
 304    ! TBH:  for now, turn off stop time and let it run data-driven
 305    CALL WRFU_ClockStopTimeDisable( grid%domain_clock, rc=rc )
 306    CALL wrf_check_error( WRFU_SUCCESS, rc, &
 307                          'WRFU_ClockStopTimeDisable(grid%domain_clock) FAILED', &
 308                          __FILE__ , &
 309                          __LINE__  )
 310    CALL domain_clockprint ( 150, grid, &
 311           'DEBUG med_sidata_input:  clock after stopTime set,' )
 312
 313    !  Here we define the initial time to process, for later use by the code.
 314
 315    current_date_char = start_date_char
 316    start_date = start_date_char // '.0000'
 317    current_date = start_date
 318
 319    CALL nl_set_bdyfrq ( grid%id , REAL(model_config_rec%interval_seconds) )
 320
 321
 322    CALL cpu_time ( t1 )
 323    DO loop = 1 , time_loop_max
 324
 325       internal_time_loop = loop
 326       IF ( ( grid%id .GT. 1 ) .AND. ( loop .GT. 1 ) .AND. &
 327            ( model_config_rec%grid_fdda(grid%id) .EQ. 0 ) .AND. &
 328            ( model_config_rec%sst_update .EQ. 0 ) ) EXIT
 329
 330       print *,' '
 331       print *,'-----------------------------------------------------------------------------'
 332       print *,' '
 333       print '(A,I2,A,A,A,I4,A,I4)' , &
 334       ' Domain ',grid%id,': Current date being processed: ',current_date, ', which is loop #',loop,' out of ',time_loop_max
 335
 336       !  After current_date has been set, fill in the julgmt stuff.
 337
 338       CALL geth_julgmt ( config_flags%julyr , config_flags%julday , config_flags%gmt )
 339
 340         print *,'configflags%julyr, %julday, %gmt:',config_flags%julyr, config_flags%julday, config_flags%gmt
 341       !  Now that the specific Julian info is available, save these in the model config record.
 342
 343       CALL nl_set_gmt (grid%id, config_flags%gmt)
 344       CALL nl_set_julyr (grid%id, config_flags%julyr)
 345       CALL nl_set_julday (grid%id, config_flags%julday)
 346
 347       !  Open the input file for tc stuff.  Either the "new" one or the "old" one.  The "new" one could have
 348       !  a suffix for the type of the data format.  Check to see if either is around.
 349
 350       CALL cpu_time ( t3 )
 351       WRITE ( wrf_err_message , FMT='(A,A)' )'med_sidata_input: calling open_r_dataset for ', &
 352                                              TRIM(config_flags%auxinput1_inname)
 353       CALL wrf_debug ( 100 , wrf_err_message )
 354       IF ( config_flags%auxinput1_inname(1:8) .NE. 'wrf_real' ) THEN
 355          CALL construct_filename4a( si_inpname , config_flags%auxinput1_inname , grid%id , 2 , &
 356                                     current_date_char , config_flags%io_form_auxinput1 )
 357       ELSE
 358          CALL construct_filename2a( si_inpname , config_flags%auxinput1_inname , grid%id , 2 , &
 359                                     current_date_char )
 360       END IF
 361       CALL open_r_dataset ( idsi, TRIM(si_inpname) , grid , config_flags , "DATASET=AUXINPUT1", ierr )
 362       IF ( ierr .NE. 0 ) THEN
 363          CALL wrf_error_fatal( 'error opening ' // TRIM(si_inpname) // &
 364                                ' for input; bad date in namelist or file not in directory' )
 365       END IF
 366
 367       !  Input data.
 368
 369       CALL wrf_debug ( 100 , 'med_sidata_input: calling input_auxinput1' )
 370       CALL input_auxinput1 ( idsi ,   grid , config_flags , ierr )
 371       WRITE ( wrf_err_message , FMT='(A,I10,A)' ) 'Timing for input ',NINT(t4-t3) ,' s.'
 372       CALL wrf_debug( 0, wrf_err_message )
 373
 374       !  Possible optional SI input.  This sets flags used by init_domain.
 375
 376       CALL cpu_time ( t3 )
 377       CALL       wrf_debug ( 100 , 'med_sidata_input: calling init_module_optional_input' )
 378       CALL init_module_optional_input ( grid , config_flags )
 379       CALL       wrf_debug ( 100 , 'med_sidata_input: calling optional_input' )
 380       CALL optional_input ( grid , idsi , config_flags )
 381
 382 !Here we check the flags yet again.  The flags are checked in optional_input but
 383 !the grid% flags are not set.
 384       flag_name(1:8) = 'SM000010'
 385       CALL wrf_get_dom_ti_integer ( idsi, 'FLAG_' // flag_name, itmp, 1, icnt, ierr )
 386       IF ( ierr .EQ. 0 ) THEN
 387           grid%flag_sm000010 = 1
 388       end if
 389
 390        flag_name(1:8) = 'SM010040'
 391        CALL wrf_get_dom_ti_integer ( idsi, 'FLAG_' // flag_name, itmp, 1, icnt, ierr )
 392        IF ( ierr .EQ. 0 ) THEN
 393           grid%flag_sm010040 = 1
 394        end if
 395
 396        flag_name(1:8) = 'SM040100'
 397        CALL wrf_get_dom_ti_integer ( idsi, 'FLAG_' // flag_name, itmp, 1, icnt, ierr )
 398        IF ( ierr .EQ. 0 ) THEN
 399             grid%flag_sm040100 = itmp
 400        end if
 401
 402
 403        flag_name(1:8) = 'SM100200'
 404        CALL wrf_get_dom_ti_integer ( idsi, 'FLAG_' // flag_name, itmp, 1, icnt, ierr )
 405        IF ( ierr .EQ. 0 ) THEN
 406             grid%flag_sm100200 = itmp
 407        end if
 408
 409 !       flag_name(1:8) = 'SM010200'
 410 !       CALL wrf_get_dom_ti_integer ( idsi, 'FLAG_' // flag_name, itmp, 1, icnt, ierr )
 411 !       IF ( ierr .EQ. 0 ) THEN
 412 !            config_flags%flag_sm010200 = itmp
 413 !            print *,"found the flag_sm010200 "
 414 !       end if
 415
 416 !Now the soil temperature flags
 417         flag_name(1:8) = 'ST000010'
 418         CALL wrf_get_dom_ti_integer ( idsi, 'FLAG_' // flag_name, itmp, 1, icnt, ierr )
 419         IF ( ierr .EQ. 0 ) THEN
 420             grid%flag_st000010 = 1
 421         END IF
 422
 423
 424          flag_name(1:8) = 'ST010040'
 425          CALL wrf_get_dom_ti_integer ( idsi, 'FLAG_' // flag_name, itmp, 1, icnt, ierr )
 426          IF ( ierr .EQ. 0 ) THEN
 427             grid%flag_st010040 = 1
 428          END IF
 429
 430          flag_name(1:8) = 'ST040100'
 431          CALL wrf_get_dom_ti_integer ( idsi, 'FLAG_' // flag_name, itmp, 1, icnt, ierr )
 432          IF ( ierr .EQ. 0 ) THEN
 433             grid%flag_st040100 = 1
 434          END IF
 435
 436
 437          flag_name(1:8) = 'ST100200'
 438          CALL wrf_get_dom_ti_integer ( idsi, 'FLAG_' // flag_name, itmp, 1, icnt, ierr )
 439          IF ( ierr .EQ. 0 ) THEN
 440             grid%flag_st100200 = 1
 441          END IF
 442
 443
 444
 445       CALL close_dataset ( idsi , config_flags , "DATASET=AUXINPUT1" )
 446       CALL cpu_time ( t4 )
 447
 448       !  Possible optional SI input.  This sets flags used by init_domain.
 449       !  We need to call the optional input routines to get the flags that
 450       !  are in the metgrid output file so they can be put in the tc bogus
 451       !  output file for real to read.
 452       CALL cpu_time ( t3 )
 453       already_been_here = .FALSE.
 454       CALL model_to_grid_config_rec ( grid%id , model_config_rec , config_flags )
 455
 456
 457       CALL cpu_time ( t3 )
 458
 459       CALL assemble_output ( grid , config_flags , loop , time_loop_max, current_date_char, &
 460                              latc_loc, lonc_loc, vmax, rmax, rankine_lid,si_inpname)
 461       CALL cpu_time ( t4 )
 462       WRITE ( wrf_err_message , FMT='(A,I10,A)' ) 'Timing for output ',NINT(t4-t3) ,' s.'
 463       CALL wrf_debug( 0, wrf_err_message )
 464       CALL cpu_time ( t2 )
 465       WRITE ( wrf_err_message , FMT='(A,I4,A,I10,A)' ) 'Timing for loop # ',loop,' = ',NINT(t2-t1) ,' s.'
 466       CALL wrf_debug( 0, wrf_err_message )
 467
 468       CALL cpu_time ( t1 )
 469    END DO
 470
 471 END SUBROUTINE tc_med_sidata_input
 472
 473
 474 !-------------------------------------------------------------------------------------
 475 SUBROUTINE tc_compute_si_start(  &
 476    start_year , start_month , start_day , start_hour , start_minute , start_second , &
 477    interval_seconds , real_data_init_type , &
 478    start_date_char)
 479
 480    USE module_date_time
 481
 482    IMPLICIT NONE
 483
 484    INTEGER :: start_year , start_month , start_day , start_hour , start_minute , start_second
 485    INTEGER ::   end_year ,   end_month ,   end_day ,   end_hour ,   end_minute ,   end_second
 486    INTEGER :: interval_seconds , real_data_init_type
 487    INTEGER :: time_loop_max , time_loop
 488
 489    CHARACTER(LEN=19) :: current_date_char , start_date_char , end_date_char , next_date_char
 490
 491 #ifdef PLANET
 492    WRITE ( start_date_char , FMT = '(I4.4,"-",I5.5,"_",I2.2,":",I2.2,":",I2.2)' ) &
 493            start_year,start_day,start_hour,start_minute,start_second
 494 #else
 495    WRITE ( start_date_char , FMT = '(I4.4,"-",I2.2,"-",I2.2,"_",I2.2,":",I2.2,":",I2.2)' ) &
 496            start_year,start_month,start_day,start_hour,start_minute,start_second
 497 #endif
 498
 499
 500 END SUBROUTINE tc_compute_si_start
 501
 502 !-----------------------------------------------------------------------
 503 SUBROUTINE assemble_output ( grid , config_flags , loop , time_loop_max,current_date_char, &
 504                              latc_loc, lonc_loc,vmax,rmax,rankine_lid,si_inpname)
 505
 506    USE module_big_step_utilities_em
 507    USE module_domain
 508    USE module_io_domain
 509    USE module_configure
 510    USE module_date_time
 511    USE module_bc
 512    IMPLICIT NONE
 513
 514    TYPE(domain)                 :: grid
 515    TYPE (grid_config_rec_type)  :: config_flags
 516
 517    INTEGER , INTENT(IN)         :: loop , time_loop_max
 518
 519 !These values are in the name list and are avaiable from
 520 !from the config_flags.
 521    real    :: vmax(max_bogus),vmax_ratio,rankine_lid
 522    real    :: rmax(max_bogus),stand_lon,cen_lat,ptop_in_pa
 523    real    :: latc_loc(max_bogus),lonc_loc(max_bogus)
 524
 525    INTEGER :: ijds , ijde , spec_bdy_width
 526    INTEGER :: i , j , k , idts,map_proj,remove_only,storms
 527
 528    INTEGER :: id1 , interval_seconds , ierr, rc, sst_update, grid_fdda
 529    INTEGER , SAVE :: id, id2,  id4
 530    CHARACTER (LEN=80) :: tcoutname , bdyname,si_inpname
 531    CHARACTER(LEN= 4) :: loop_char
 532    CHARACTER(LEN=19) ::  current_date_char
 533
 534 character *19 :: temp19
 535 character *24 :: temp24 , temp24b
 536
 537 real::t1,t2,truelat1,truelat2
 538
 539
 540    !  Boundary width, scalar value.
 541
 542    spec_bdy_width = model_config_rec%spec_bdy_width
 543    interval_seconds = model_config_rec%interval_seconds
 544    sst_update = model_config_rec%sst_update
 545    grid_fdda = model_config_rec%grid_fdda(grid%id)
 546    truelat1  = config_flags%truelat1
 547    truelat2  = config_flags%truelat2
 548
 549    stand_lon = config_flags%stand_lon
 550    cen_lat   = config_flags%cen_lat
 551    map_proj  = config_flags%map_proj
 552
 553    vmax_ratio = config_flags%vmax_ratio
 554    ptop_in_pa = config_flags%p_top_requested
 555    remove_only = 0
 556    if(config_flags%remove_storm) then
 557       remove_only = 1
 558    end if
 559
 560    storms = config_flags%num_storm
 561    print *,"number of storms ",config_flags%num_storm
 562    call tc_bogus(cen_lat,stand_lon,map_proj,truelat1,truelat2, &
 563                  grid%dx,grid%e_we,grid%e_sn,grid%num_metgrid_levels,ptop_in_pa, &
 564                  rankine_lid,latc_loc,lonc_loc,vmax,vmax_ratio,rmax,remove_only, &
 565                  storms,grid)
 566
 567
 568
 569    !  Open the tc bogused output file. cd
 570    CALL construct_filename4a( tcoutname , config_flags%auxinput1_outname , grid%id , 2 , &
 571                                     current_date_char , config_flags%io_form_auxinput1 )
 572
 573    print *,"outfile name from construct filename ",tcoutname
 574    CALL open_w_dataset ( id1, TRIM(tcoutname) , grid , config_flags ,output_auxinput1,"DATASET=AUXINPUT1",ierr )
 575    IF ( ierr .NE. 0 ) THEN
 576         CALL wrf_error_fatal( 'tc_em: error opening tc bogus file for writing' )
 577    END IF
 578    CALL output_auxinput1( id1, grid , config_flags , ierr )
 579    CALL close_dataset ( id1 , config_flags , "DATASET=AUXINPUT1" )
 580
 581
 582 END SUBROUTINE assemble_output
 583
 584 !----------------------------------------------------------------------------------------------
 585
 586 SUBROUTINE tc_bogus(centerlat,stdlon,nproj,truelat1,truelat2,dsm,ew,ns,nz,ptop_in_pa, &
 587                     rankine_lid,latc_loc,lonc_loc,vmax,vmax_ratio,rmax,remove_only, &
 588                     storms,grid)
 589
 590 !!Original Author Dave Gill.  Modified by Sherrie Fredrick
 591 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 592 !These are read in from the netcdf file.
 593 !centerlat  The center latitude from the global attributes in the netcdf file.
 594 !stdlon     The center longitude from the global attributes in the netcdf file.
 595 !nproj      The map projection from the global attributes in the netcdf file.
 596 !dsm        The spacing in meters from the global attributes in the netcdf file.
 597 !ew         The west_east_stag from the dimensions in the netcdf file..
 598 !ns         The south_north_stag from the dimensions in the netcdf file. .
 599 !nz         The number of metgrid levels from the dimensions in the netcdf file.
 600
 601 !ptop_in_pa This is part of the namelist.input file under the &domains section.
 602
 603 !These values are part of the namelist.input file under the &tc section specifically
 604 !for the tc bogus code.
 605 !NOTES: There can be up to five bogus storms.  The variable max_bogus is set in
 606 !the WRF subroutine called module_driver_constants.F in the ./WRFV3/frame directory.
 607
 608 !latc_loc    The center latitude of the bogus strorm. This is an array dimensioned max_bogus.
 609
 610 !lonc_loc    The center longitude of the bogus strorm. This is an array dimensioned max_bogus.
 611
 612 !vmax        The max vortex in meters/second it comes from the namelist entry.
 613 !             This is an array dimensioned max_bogus.
 614
 615 !vmax_ratio  This comes from the namelist entry.
 616
 617 !rmax        The maximum radius this comes from the namelist entry.
 618 !             This is an array dimensioned max_bogus
 619
 620 !remove_only If this is set to true in the namelist.input file a value of 0.1
 621 !             is automatically assigned to vmax.
 622
 623 !rankine_lid This comes from the namelist entry.  It can be used to determine
 624 !            what model levels the bogus storm affects.
 625
 626 !storms      The number of bogus storms.
 627
 628 !grid        This is a Fortran structure which holds all of the field data values
 629 !            for the netcdf that was read in.
 630 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 631
 632
 633 !module_llxy resides in the share directory.
 634   USE module_llxy
 635 !This is for the large structure (grid)
 636   USE module_domain
 637
 638
 639
 640   IMPLICIT NONE
 641   TYPE(domain)                 :: grid
 642   integer ew,ns,nz
 643   integer nproj
 644   integer storms,nstrm
 645   real :: centerlat,stdlon,conef,truelat1,truelat2,dsm,dx,rankine_lid
 646   real :: latc_loc(max_bogus),lonc_loc(max_bogus),vmax(max_bogus),vmax_ratio,rmax(max_bogus)
 647
 648   real :: press(ew-1,nz,ns-1),rhmx(nz), vwgt(nz),old_slp(ew-1,ns-1)
 649   real, dimension(:,:,:) , allocatable :: u11,v11,t11,rh11,phi11
 650   real, dimension(:,:,:) , allocatable :: u1 , v1 , t1 , rh1 , phi1
 651   real, dimension(ew-1,ns-1) :: lond,terrain,cor,pslx
 652
 653
 654 !The map scale factors.
 655   real, dimension(ew,ns-1)    :: msfu   !The mapscale factor for the ew wind staggered grid
 656   real, dimension(ew-1,ns)    :: msfv   !The mapscale factor for the ns wind staggered grid
 657   real, dimension(ew-1,ns-1)  :: msfm   !The mapscale factor for the unstaggered grid.
 658
 659   CHARACTER*2  jproj
 660   LOGICAL :: l_tcbogus
 661
 662
 663   real :: r_search,r_vor,beta,devps,humidity_max
 664   real :: devpc,const,r_vor2,cnst,alphar,epsilon,vormx , rad , sum_q
 665   real :: avg_q ,q_old,ror,q_new,dph,dphx0
 666   real :: rh_max,min_RH_value,ps
 667   integer :: vert_variation
 668   integer :: i,k,j,kx,remove_only
 669   integer :: k00,kfrm ,kto ,k85,n_iter,ew_mvc,ns_mvc,nct,itr
 670   integer :: strmci(nz), strmcj(nz)
 671   real :: disx,disy,alpha,degran,pie,rovcp,cp
 672   REAL :: rho,pprm,phip0,x0,y0,vmx,xico,xjco,xicn,xjcn,p85,xlo,rconst,ew_gcntr,ns_gcntr
 673   real :: ptop_in_pa,themax,themin
 674   real :: latinc,loninc
 675   real :: rtemp,colat0,colat
 676   REAL :: q1(ew-1,nz,ns-1), psi1(ew-1,nz,ns-1)
 677
 678 ! This is the entire map projection enchilada.
 679   TYPE(proj_info) :: proj
 680
 681
 682
 683   REAL :: lat1 , lon1
 684 ! These values are read in from the data set.
 685    real :: knowni,knownj
 686
 687 !  TC bogus
 688    REAL utcr(ew,nz,ns-1),  vtcr(ew-1,nz,ns)
 689    REAL utcp(ew,nz,ns-1),  vtcp(ew-1,nz,ns)
 690    REAL psitc(ew-1,nz,ns-1), psiv(nz)
 691    REAL vortc(ew-1,nz,ns-1), vorv(nz)
 692    REAL tptc(ew-1,nz,ns-1)
 693    REAL phiptc(ew-1,nz,ns-1)
 694
 695 !  Work arrays
 696    REAL uuwork(nz), vvwork(nz), temp2(ew,ns)
 697    REAL vort(ew-1,nz,ns-1), div(ew-1,nz,ns-1)
 698    REAL vortsv(ew-1,nz,ns-1)
 699    REAL theta(ew-1,nz,ns-1), t_reduce(ew-1,nz,ns-1)
 700    REAL ug(ew,nz,ns-1),   vg(ew-1,nz,ns),  vorg(ew-1,nz,ns-1)
 701    REAL delpx(ew-1,ns-1)
 702
 703 !subroutines for relaxation
 704    REAL outold(ew-1,ns-1)
 705    REAL rd(ew-1,ns-1),     ff(ew-1,ns-1)
 706    REAL tmp1(ew-1,ns-1),   tmp2(ew-1,ns-1)
 707
 708 !  Background fields.
 709    REAL , DIMENSION (ew-1,nz,ns-1) :: t0, t00, rh0, q0, phi0, psi0, chi
 710
 711 !  Perturbations
 712    REAL , DIMENSION (ew-1,nz,ns-1) :: psipos, tpos, psi ,phipos, phip
 713
 714 !  Final fields.
 715    REAL  u2(ew,nz,ns-1),  v2(ew-1,nz,ns)
 716    REAL  t2(ew-1,nz,ns-1),z2(ew-1,nz,ns-1)
 717    REAL  phi2(ew-1,nz,ns-1),rh2(ew-1,nz,ns-1)
 718
 719    print *,"the dimensions: north-south = ",ns," east-west =",ew
 720    IF (nproj .EQ. 1) THEN
 721         jproj = 'LC'
 722         print *,"Lambert Conformal projection"
 723    ELSE IF (nproj .EQ. 2) THEN
 724         jproj = 'ST'
 725    ELSE IF (nproj .EQ. 3) THEN
 726         jproj = 'ME'
 727         print *,"A mercator projection"
 728    END IF
 729
 730
 731   knowni = 1.
 732   knownj = 1.
 733   pie     = 3.141592653589793
 734   degran = pie/180.
 735   rconst = 287.04
 736   min_RH_value = 5.0
 737   cp = 1004.0
 738   rovcp = rconst/cp
 739
 740    r_search = 400000.0
 741    r_vor = 300000.0
 742    r_vor2 = r_vor * 4
 743    beta = 0.5
 744    devpc= 40.0
 745    vert_variation = 1
 746    humidity_max   = 95.0
 747    alphar         = 1.8
 748    latinc        = -999.
 749    loninc        = -999.
 750
 751    if(remove_only .eq. 1) then
 752      do nstrm=1,storms
 753          vmax(nstrm) = 0.1
 754      end do
 755    end if
 756
 757   !  Set up initializations for map projection so that the lat/lon
 758   !  of the tropical storm can be put into model (i,j) space.  This needs to be done once per
 759   !  map projection definition.  Since this is the domain that we are "GOING TO", it is a once
 760   !  per regridder requirement.  If the user somehow ends up calling this routine for several
 761   !  time periods, there is no problemos, just a bit of overhead with redundant calls.
 762
 763    dx = dsm
 764    lat1 = grid%xlat_gc(1,1)
 765    lon1 = grid%xlong_gc(1,1)
 766    IF( jproj .EQ. 'ME' )THEN
 767        IF ( lon1  .LT. -180. ) lon1  = lon1  + 360.
 768        IF ( lon1  .GT.  180. ) lon1  = lon1  - 360.
 769        IF ( stdlon .LT. -180. ) stdlon = stdlon + 360.
 770        IF ( stdlon .GT.  180. ) stdlon = stdlon - 360.
 771        CALL map_set ( proj_merc, proj, lat1, lon1, lat1, lon1, knowni, knownj, dx, &
 772                       latinc,loninc,stdlon , truelat1 , truelat2)
 773        conef = 0.
 774    ELSE IF ( jproj .EQ. 'LC' ) THEN
 775         if((truelat1 .eq. 0.0)  .and. (truelat2 .eq. 0.0)) then
 776             print *,"Truelat1 and Truelat2 are both 0"
 777             stop
 778          end if
 779         CALL map_set (proj_lc,proj, lat1, lon1, lat1, lon1, knowni, knownj, dx, &
 780                        latinc,loninc,stdlon , truelat1 , truelat2)
 781        conef = proj%cone
 782    ELSE IF ( jproj .EQ. 'ST' ) THEN
 783         conef = 1.
 784         CALL map_set ( proj_ps,proj,lat1, lon1, lat1, lon1, knowni, knownj, dx, &
 785                       latinc,loninc,stdlon , truelat1 , truelat2)
 786    END IF
 787
 788 ! Load the pressure array.
 789  kx = nz
 790  do j = 1,ns-1
 791     do k = 1,nz
 792        do i = 1,ew-1
 793            press(i,k,j) = grid%p_gc(i,k,j)*0.01
 794        end do
 795     end do
 796  end do
 797
 798
 799 !  Initialize the vertical profiles for humidity and weighting.
 800 !The ptop variable will be read in from the namelist
 801    IF ( ( ptop_in_pa .EQ. 40000. ) .OR. ( ptop_in_pa .EQ. 60000. ) ) THEN
 802          PRINT '(A)','Hold on pardner, your value for PTOP is gonna cause problems for the TC bogus option.'
 803          PRINT '(A)','Make it higher up than 400 mb.'
 804          STOP 'ptop_woes_for_tc_bogus'
 805    END IF
 806
 807  IF ( vert_variation .EQ. 1 ) THEN
 808     DO k=1,kx
 809        IF ( press(1,k,1) .GT. 400. ) THEN
 810                rhmx(k) = humidity_max
 811        ELSE
 812                rhmx(k) = humidity_max * MAX( 0.1 , (press(1,k,1) - ptop_in_pa/100.)/(400.-ptop_in_pa/100.) )
 813        END IF
 814
 815         IF ( press(1,k,1) .GT. 600. ) THEN
 816              vwgt(k) = 1.0
 817         ELSE IF ( press(1,k,1) .LE. 100. ) THEN
 818              vwgt(k) = 0.0001
 819         ELSE
 820              vwgt(k) = MAX ( 0.0001 , (press(1,k,1)-ptop_in_pa/100.)/(600.-ptop_in_pa/100.) )
 821         END IF
 822       END DO
 823
 824  ELSE IF ( vert_variation .EQ. 2 ) THEN
 825          IF ( kx .eq. 24 ) THEN
 826             rhmx = (/ 95.,       95., 95., 95., 95., 95., 95., 95.,      &
 827                       95., 95.,  95., 95., 95., 90., 85., 80., 75.,      &
 828                       70., 66.,  60., 39., 10., 10., 10./)
 829             vwgt = (/ 1.0000,         1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 0.9850,      &
 830                       0.9680, 0.9500, 0.9290, 0.9060, 0.8810, 0.8500, 0.7580, 0.6500, 0.5100,      &
 831                       0.3500, 0.2120, 0.0500, 0.0270, 0.0001, 0.0001, 0.0001/)
 832          ELSE
 833             PRINT '(A)','Number of vertical levels assumed to be 24 for AFWA TC bogus option'
 834             STOP 'AFWA_TC_BOGUS_LEVEL_ERROR'
 835          END IF
 836  END IF
 837
 838 !Remember that ns = the north south staggered. This is one more than the ns mass point grid.
 839 !              ew = the east west staggered. This is one more than the ew mass point grid.
 840
 841
 842 !Put the U and V into the new arrays.
 843 !Remember that the WRF ordering is ew,vert level,ns
 844 !Vorticity and Divergence calculatins are done on
 845 !the staggered grids so the winds are not destaggered
 846  allocate(u11 (1:ew, 1:nz, 1:ns-1))
 847  allocate(u1  (1:ew, 1:nz, 1:ns-1))
 848  allocate(v11 (1:ew-1, 1:nz, 1:ns))
 849  allocate(v1  (1:ew-1, 1:nz, 1:ns))
 850  do j = 1,ns-1
 851     do k = 1,nz
 852        do i = 1,ew
 853             u11(i,k,j) = grid%u_gc(i,k,j)
 854              u1(i,k,j) = grid%u_gc(i,k,j)
 855              msfu(i,j) = grid%msfu(i,j) !map scale factor on the U staggered grid
 856        end do
 857     end do
 858  end do
 859
 860
 861  do j = 1,ns
 862     do k = 1,nz
 863        do i = 1,ew-1
 864             v11(i,k,j) = grid%v_gc(i,k,j)
 865              v1(i,k,j) = grid%v_gc(i,k,j)
 866            msfv(i,j)   = grid%msfv(i,j)  !map scale factor on the V staggered grid
 867        end do
 868     end do
 869  end do
 870
 871
 872 !Put the temperature, relative humidity and height fields
 873 !into arrays.  Save the initial fields also.
 874 !These arrays are on the WRF mass points
 875  allocate(t11  (1:ew-1, 1:nz, 1:ns-1))
 876  allocate(t1   (1:ew-1, 1:nz, 1:ns-1))
 877  allocate(rh11 (1:ew-1, 1:nz, 1:ns-1))
 878  allocate(rh1  (1:ew-1, 1:nz, 1:ns-1))
 879  allocate(phi11(1:ew-1, 1:nz, 1:ns-1))
 880  allocate(phi1 (1:ew-1, 1:nz, 1:ns-1))
 881  do j = 1,ns-1
 882     do k = 1,nz
 883        do i = 1,ew-1
 884              t11(i,k,j)  =  grid%t_gc(i,k,j)
 885               t1(i,k,j)  =  grid%t_gc(i,k,j)
 886             rh11(i,k,j)  =  grid%rh_gc(i,k,j)
 887              rh1(i,k,j)  =  grid%rh_gc(i,k,j)
 888               msfm(i,j)  = grid%msft(i,j)
 889             if(k .eq. 1)then
 890                phi11(i,k,j) =  grid%ht_gc(i,j)
 891                phi1(i,k,j)  =  grid%ht_gc(i,j) * 9.81
 892             else
 893                phi11(i,k,j) =  grid%ght_gc(i,k,j)
 894                phi1(i,k,j)  =  grid%ght_gc(i,k,j) * 9.81
 895             end if
 896        end do
 897     end do
 898  end do
 899
 900 !The two D fields
 901 !The terrain soil height is from ght at level 1
 902  do j = 1,ns-1
 903     do i = 1,ew-1
 904        pslx(i,j)    = grid%pslv_gc(i,j) * 0.01
 905        cor(i,j)     = grid%f(i,j)               !coreolous
 906        lond(i,j)    = grid%xlong_gc(i,j)
 907        terrain(i,j) = grid%ht_gc(i,j)
 908        old_slp(i,j) = grid%pslv_gc(i,j)
 909     end do
 910  end do
 911
 912
 913
 914 !  Loop over the number of storms to process.
 915
 916  l_tcbogus = .FALSE.
 917  all_storms : DO nstrm=1,storms
 918
 919
 920 !Make sure the user has defined the rmax variable
 921  if(rmax(nstrm) .eq. -999.) then
 922     print *,"Please enter a value for rmax in the namelist"
 923     stop
 924  end if
 925
 926
 927  k00  = 2
 928  kfrm = k00
 929  p85  = 850.
 930
 931  kto  = kfrm
 932  DO k=kfrm+1,kx
 933      IF ( press(1,k,1) .GE. p85 ) THEN
 934            kto = kto + 1
 935      END IF
 936  END DO
 937  k85 = kto
 938
 939
 940 !  Parameters for max wind
 941  rho  = 1.2
 942  pprm = devpc*100.
 943  phip0= pprm/rho
 944
 945
 946 !latc_loc and lonc_loc come in from the namelist.
 947 !These x0 and y0 points are relative to the mass points.
 948  CALL latlon_to_ij ( proj , latc_loc(nstrm) , lonc_loc(nstrm) , x0 , y0 )
 949  IF ( ( x0 .LT. 1. ) .OR. ( x0 .GT. REAL(ew-1) ) .OR. &
 950               ( y0 .LT. 1. ) .OR. ( y0 .GT. REAL(ns-1) ) ) THEN
 951          PRINT '(A,I3,A,A,A)','         Storm position is outside the computational domain.'
 952          PRINT '(A,2F6.2,A)' ,'         Storm postion: (x,y) = ',x0,y0,'.'
 953          stop
 954  END IF
 955
 956  l_tcbogus = .TRUE.
 957 !  Bogus vortex specifications, vmax (m/s); rmax (m);
 958  vmx = vmax(nstrm)  * vmax_ratio
 959
 960  IF (  latc_loc(nstrm) .LT. 0.  ) THEN
 961        vmx = -vmx
 962  END IF
 963
 964  IF (  vmax(nstrm)  .LE. 0.  ) THEN
 965        vmx = SQRT( 2.*(1-beta)*ABS(phip0) )
 966  END IF
 967
 968  ew_gcntr    = x0  !ew center grid location
 969  ns_gcntr    = y0  !ns center grid location
 970 !For right now we are adding 0.5 to the grid location this
 971 !makes the output of the wrf tc_bogus scheme analogous to the
 972 !ouput of the MM5 tc_bogus scheme.
 973  ew_gcntr    = x0 + 0.5
 974  ns_gcntr    = y0 + 0.5
 975
 976  n_iter  = 1
 977
 978 !  Start computing.
 979
 980  PRINT '(/,A,I3,A,A,A)'     ,'---> TC: Processing storm number= ',nstrm
 981  PRINT '(A,F6.2,A,F7.2,A)'  ,'         Storm center lat= ',latc_loc(nstrm),' lon= ',lonc_loc(nstrm),'.'
 982  PRINT '(A,2F6.2,A)'        ,'         Storm center grid position (x,y)= ',ew_gcntr,ns_gcntr,'.'
 983  PRINT '(A,F5.2,F9.2,A)'    ,'         Storm max wind (m/s) and max radius (m)= ',vmx,rmax(nstrm),'.'
 984  PRINT '(A,F5.2,A)'         ,'         Estimated central press dev (mb)= ',devpc,'.'
 985
 986
 987 !  Initialize storm center to (1,1)
 988
 989   DO k=1,kx
 990      strmci(k) = 1
 991      strmcj(k) = 1
 992   END DO
 993
 994 !  Define complete field of bogus storm
 995 !Note dx is spacing in meters.
 996 !The output arrays from the rankine subroutine vvwork,uuwork,psiv and vorv
 997 !are defined on the WRF mass points.
 998   utcp(:,:,:) = 0.0
 999   vtcp(:,:,:) = 0.0
1000   print *,"nstrm  ",rmax(nstrm),ew_gcntr,ns_gcntr
1001   DO j=1,ns-1
1002      DO i=1,ew-1
1003         disx = REAL(i) - ew_gcntr
1004         disy = REAL(j) - ns_gcntr
1005         CALL rankine(disx,disy,dx,kx,vwgt,rmax(nstrm),vmx,uuwork,vvwork,psiv,vorv)
1006         DO k=1,kx
1007             utcp(i,k,j) = uuwork(k)
1008             vtcp(i,k,j) = vvwork(k)
1009            psitc(i,k,j) = psiv(k)
1010            vortc(i,k,j) = vorv(k)
1011         END DO
1012      END DO
1013   END DO
1014   call stagger_rankine_winds(utcp,vtcp,ew,ns,nz)
1015
1016
1017   utcr(:,:,:) = 0.0
1018   vtcr(:,:,:) = 0.0
1019 ! dave Rotate wind to map proj, on the correct staggering
1020   DO j=1,ns-1
1021      DO i=2,ew-1
1022         xlo = stdlon-grid%xlong_u(i,j)
1023         IF ( xlo .GT. 180.)xlo = xlo-360.
1024         IF ( xlo .LT.-180.)xlo = xlo+360.
1025
1026         alpha = xlo*conef*degran*SIGN(1.,centerlat)
1027         DO k=1,kx
1028            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)
1029            if(utcr(i,k,j) .gt. 300.) then
1030               print *,i,k,j,"a very bad value of utcr"
1031               stop
1032            end if
1033         END DO
1034      END DO
1035   END DO
1036
1037
1038   DO j=2,ns-1
1039      DO i=1,ew-1
1040         xlo = stdlon-grid%xlong_v(i,j)
1041         IF ( xlo .GT. 180.)xlo = xlo-360.
1042         IF ( xlo .LT.-180.)xlo = xlo+360.
1043
1044         alpha = xlo*conef*degran*SIGN(1.,centerlat)
1045         DO k=1,kx
1046            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)
1047            if(vtcr(i,k,j) .gt. 300.) then
1048               print *,i,k,j,"a very bad value of vtcr"
1049               stop
1050            end if
1051         END DO
1052      END DO
1053   END DO
1054
1055
1056 !Fill in UTCR's along the left and right side.
1057   do j = 1,ns-1
1058      utcr(1,:,j)  = utcr(2,:,j)
1059      utcr(ew,:,j) = utcr(ew-1,:,j)
1060  end do
1061
1062 !Fill in V's along the bottom and top.
1063   do i = 1,ew-1
1064      vtcr(i,:,1)  = vtcr(i,:,2)
1065      vtcr(i,:,ns) = vtcr(i,:,ns-1)
1066   end do
1067
1068
1069 !  Compute vorticity of FG.  This is the vorticity of the original winds
1070 !  on the staggered grid.  The vorticity and divergence are defined at
1071 !  the mass points when done.
1072    CALL vor(u1,v1,msfu,msfv,msfm,ew,ns,kx,dx,vort)
1073
1074
1075 !  Compute divergence of FG
1076    CALL diverg(u1,v1,msfu,msfv,msfm,ew,ns,kx,dx,div)
1077
1078
1079 !  Compute mixing ratio of FG
1080    CALL mxratprs(rh1,t1,press*100.,ew,ns,kx,q1,min_RH_value)
1081    q1(:,1,:) = q1(:,2,:)
1082
1083
1084 !  Compute initial streamfunction - PSI1
1085    vortsv = vort
1086    q0 = q1
1087
1088
1089 !  Solve for streamfunction.
1090    DO k=1,kx
1091       DO j=1,ns-1
1092          DO i=1,ew-1
1093             ff(i,j) = vort(i,k,j)
1094             tmp1(i,j)= 0.0
1095          END DO
1096       END DO
1097       epsilon = 1.E-2
1098       CALL relax(tmp1,ff,rd,ew,ns,dx,epsilon,alphar)
1099       DO j=1,ns-1
1100          DO i=1,ew-1
1101             psi1(i,k,j) = tmp1(i,j)
1102          END DO
1103       END DO
1104    END DO
1105
1106
1107    DO k=1,kx  !start of the k loop
1108       IF ( latc_loc(nstrm) .GE. 0. ) THEN
1109            vormx = -1.e10
1110       ELSE
1111            vormx =  1.e10
1112       END IF
1113
1114       ew_mvc = 1
1115       ns_mvc = 1
1116
1117       DO j=1,ns-1
1118          DO i=1,ew-1
1119             rad = SQRT((REAL(i)-ew_gcntr)**2.+(REAL(j)-ns_gcntr)**2.)*dx
1120             IF ( rad .LE. r_search ) THEN
1121                IF ( latc_loc(nstrm) .GE. 0. ) THEN
1122                    IF ( vortsv(i,k,j) .GT. vormx ) THEN
1123                         vormx = vortsv(i,k,j)
1124                         ew_mvc = i
1125                         ns_mvc = j
1126                     END IF
1127                ELSE IF (latc_loc(nstrm) .LT. 0. ) THEN
1128                     IF ( vortsv(i,k,j) .LT. vormx ) THEN
1129                          vormx = vortsv(i,k,j)
1130                          ew_mvc = i
1131                          ns_mvc = j
1132                     END IF
1133                END IF
1134             END IF
1135          END DO
1136       END DO
1137
1138       strmci(k) = ew_mvc
1139       strmcj(k) = ns_mvc
1140
1141       DO j=1,ns-1
1142          DO i=1,ew-1
1143             rad = SQRT(REAL((i-ew_mvc)**2.+(j-ns_mvc)**2.))*dx
1144             IF ( rad .GT. r_vor ) THEN
1145                  vort(i,k,j) = 0.
1146                  div(i,k,j)  = 0.
1147             END IF
1148          END DO
1149       END DO
1150
1151       DO itr=1,n_iter
1152          sum_q = 0.
1153          nct = 0
1154          DO j=1,ns-1
1155             DO i=1,ew-1
1156                rad = SQRT(REAL(i-ew_mvc)**2.+REAL(j-ns_mvc)**2.)*dx
1157                IF ( (rad .LT. r_vor2).AND.(rad .GE. 0.8*r_vor2) ) THEN
1158                      sum_q = sum_q + q0(i,k,j)
1159                      nct = nct + 1
1160                END IF
1161              END DO
1162           END DO
1163           avg_q = sum_q/MAX(REAL(nct),1.)
1164
1165           DO j=1,ns-1
1166              DO i=1,ew-1
1167                  q_old = q0(i,k,j)
1168                  rad = SQRT(REAL(i-ew_mvc)**2.+REAL(j-ns_mvc)**2.)*dx
1169                  IF ( rad .LT. r_vor2 ) THEN
1170                       ror = rad/r_vor2
1171                       q_new = ((1.-ror)*avg_q) + (ror*q_old)
1172                       q0(i,k,j) = q_new
1173                  END IF
1174               END DO
1175            END DO
1176      END DO !end of itr loop
1177  END DO !of the k loop
1178
1179
1180 !  Compute divergent wind (chi) at the mass points
1181    DO k=1,kx
1182       DO j=1,ns-1
1183          DO i=1,ew-1
1184             ff(i,j) = div(i,k,j)
1185             tmp1(i,j)= 0.0
1186          END DO
1187       END DO
1188
1189       epsilon = 1.e-2
1190       CALL relax(tmp1,ff,rd,ew,ns,dx,epsilon,alphar)
1191       DO j=1,ns-1
1192          DO i=1,ew-1
1193             chi(i,k,j) = tmp1(i,j)
1194          END DO
1195       END DO
1196     END DO !of the k loop for divergent winds
1197
1198
1199
1200 !  Compute background streamfunction (PSI0) and perturbation field (PSI)
1201 !     print *,"perturbation field (PSI) relax three"
1202      DO k=1,kx
1203          DO j=1,ns-1
1204             DO i=1,ew-1
1205                ff(i,j)=vort(i,k,j)
1206                tmp1(i,j)=0.0
1207             END DO
1208          END DO
1209          epsilon = 1.e-2
1210          CALL relax(tmp1,ff,rd,ew,ns,dx,epsilon,alphar)
1211          DO j=1,ns-1
1212             DO i=1,ew-1
1213                psi(i,k,j)=tmp1(i,j)
1214             END DO
1215          END DO
1216      END DO
1217
1218
1219  !We can now calculate the final wind fields.
1220    call final_ew_velocity(u2,u1,chi,psi,utcr,dx,ew,ns,nz)
1221    call final_ns_velocity(v2,v1,chi,psi,vtcr,dx,ew,ns,nz)
1222
1223      DO k=1,kx
1224         DO j=1,ns-1
1225            DO i=1,ew-1
1226               psi0(i,k,j) = psi1(i,k,j)-psi(i,k,j)
1227            END DO
1228         END DO
1229      END DO
1230
1231      DO k=k00,kx
1232         DO j=1,ns-1
1233            DO i=1,ew-1
1234               psipos(i,k,j)=psi(i,k,j)
1235            END DO
1236         END DO
1237      END DO
1238
1239
1240 !  Geostrophic vorticity.
1241 !We calculate the ug and vg on the wrf U and V staggered grids
1242 !since this is where the vorticity subroutine expects them.
1243
1244      CALL geowind(phi1,ew,ns,kx,dx,ug,vg)
1245      CALL vor(ug,vg,msfu,msfv,msfm,ew,ns,kx,dx,vorg)
1246
1247      DO k=1,kx
1248         ew_mvc = strmci(k)
1249         ns_mvc = strmcj(k)
1250
1251          DO j=1,ns-1
1252            DO i=1,ew-1
1253                rad = SQRT(REAL(i-ew_mvc)**2.+REAL(j-ns_mvc)**2.)*dx
1254                IF ( rad .GT. r_vor ) THEN
1255                     vorg(i,k,j) = 0.
1256                END IF
1257            END DO
1258          END DO
1259      END DO
1260
1261       DO k=k00,kx
1262          DO j=1,ns-1
1263             DO i=1,ew-1
1264                ff(i,j) = vorg(i,k,j)
1265                tmp1(i,j)= 0.0
1266             END DO
1267          END DO
1268          epsilon = 1.e-3
1269          CALL relax(tmp1,ff,rd,ew,ns,dx,epsilon,alphar)
1270          DO j=1,ns-1
1271             DO i=1,ew-1
1272                phip(i,k,j) = tmp1(i,j)
1273             END DO
1274          END DO
1275      END DO
1276
1277
1278      !  Background geopotential.
1279      DO k=k00,kx
1280          DO j=1,ns-1
1281             DO i=1,ew-1
1282                phi0(i,k,j) = phi1(i,k,j) - phip(i,k,j)
1283             END DO
1284          END DO
1285      END DO
1286
1287
1288      !  Background temperature
1289      DO k=k00,kx
1290         DO j=1,ns-1
1291            DO i=1,ew-1
1292               IF( k .EQ.  2 ) THEN
1293                   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))
1294               ELSE IF ( k .EQ. kx ) THEN
1295                   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))
1296               ELSE
1297                   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))
1298               END IF
1299               t0(i,k,j) = t1(i,k,j)-tpos(i,k,j)
1300               t00(i,k,j) = t0(i,k,j)
1301               if(t0(i,k,j) .gt. 400) then
1302                  print *,"interesting temperature ",t0(i,k,j)," at ",i,j,k
1303                  stop
1304               end if
1305            END DO
1306         END DO
1307      END DO
1308
1309      !  New RH.
1310      CALL qvtorh (q0,t0,press*100.,k00,ew,ns,kx,rh0,min_RH_value)
1311      call final_RH(rh2,rh0,rhmx,strmci,strmcj,rmax(nstrm),ew,ns,nz,k00,dx,ew_gcntr,ns_gcntr,r_vor2)
1312
1313
1314
1315      ! adjust T0
1316      DO k=k00,kx
1317         DO j=1,ns-1
1318            DO i=1,ew-1
1319               theta(i,k,j) = t1(i,k,j)*(1000./press(i,k,j))**rovcp
1320            END DO
1321         END DO
1322      END DO
1323
1324
1325      ew_mvc = strmci(k00)
1326      ns_mvc = strmcj(k00)
1327      DO k=kfrm,kto
1328         DO j=1,ns-1
1329            DO i=1,ew-1
1330               rad = SQRT(REAL(i-ew_mvc)**2.+REAL(j-ns_mvc)**2.)*dx
1331               IF ( rad .LT. r_vor2 ) THEN
1332                   t_reduce(i,k,j) = theta(i,k85,j)-0.03*(press(i,k,j)-press(i,k85,j))
1333                   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))
1334               END IF
1335            END DO
1336         END DO
1337      END DO
1338
1339     !  Geopotential perturbation
1340     DO k=1,kx
1341        DO j=1,ns-1
1342           DO i=1,ew-1
1343               tmp1(i,j)=psitc(i,k,j)
1344           END DO
1345        END DO
1346        CALL balance(cor,tmp1,ew,ns,dx,outold)
1347        DO j=1,ns-1
1348           DO i=1,ew-1
1349              ff(i,j)=outold(i,j)
1350              tmp1(i,j)=0.0
1351           END DO
1352        END DO
1353        epsilon = 1.e-3
1354        CALL relax (tmp1,ff,rd,ew,ns,dx,epsilon,alphar)
1355        DO j=1,ns-1
1356           DO i=1,ew-1
1357              phiptc(i,k,j) = tmp1(i,j)
1358           END DO
1359        END DO
1360     END DO
1361
1362
1363 !  New geopotential field.
1364    DO j=1,ns-1
1365       DO k=1,kx
1366          DO i=1,ew-1
1367             phi2(i,k,j)  = phi0(i,k,j) + phiptc(i,k,j)
1368          END DO
1369       END DO
1370    END DO
1371
1372
1373    !  New temperature field.
1374     DO j=1,ns-1
1375        DO k=k00,kx
1376           DO i=1,ew-1
1377              IF( k .EQ.  2 ) THEN
1378                  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))
1379              ELSE IF ( k .EQ. kx ) THEN
1380                  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))
1381              ELSE
1382                  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))
1383              END IF
1384              t2(i,k,j) = t0(i,k,j) + tptc(i,k,j)
1385              if(t2(i,k,j) .gt. 400) then
1386                 print *,"interesting temperature "
1387                 print *,t2(i,k,j),i,k,j,tptc(i,k,j)
1388                 stop
1389              end if
1390            END DO
1391         END DO
1392     END DO
1393
1394
1395    !  Sea level pressure change.
1396       DO j=1,ns-1
1397          DO i=1,ew-1
1398             dph = phi2(i,k00,j)-phi1(i,k00,j)
1399             delpx(i,j) = rho*dph*0.01
1400          END DO
1401       END DO
1402
1403
1404     !  New SLP.
1405 !      print *,"new slp",nstrm
1406       DO j=1,ns-1
1407          DO i=1,ew-1
1408             pslx(i,j) = pslx(i,j)+delpx(i,j)
1409             grid%pslv_gc(i,j) = pslx(i,j) * 100.
1410 !            print *,pslx(i,j)
1411          END DO
1412       END DO
1413
1414   !  Set new geopotential at surface to terrain elevation.
1415      DO j=1,ns-1
1416         DO i=1,ew-1
1417            z2(i,1,j) = terrain(i,j)
1418         END DO
1419      END DO
1420
1421   !  Geopotential back to height.
1422
1423      DO j=1,ns-1
1424         DO k=k00,kx
1425            DO i=1,ew-1
1426                z2(i,k,j) = phi2(i,k,j)/9.81
1427             END DO
1428          END DO
1429      END DO
1430
1431
1432      !  New surface temperature, assuming same theta as from 1000 mb.
1433 !     print *,"new surface temperature"
1434      DO j=1,ns-1
1435         DO i=1,ew-1
1436            ps = pslx(i,j)
1437            t2(i,1,j) = t2(i,k00,j)*((ps/1000.)**rovcp)
1438            if(t2(i,1,j) .gt. 400) then
1439               print *,"Interesting surface temperature"
1440               print *,t2(i,1,j),t2(i,k00,j),ps,i,j
1441               stop
1442            end if
1443         END DO
1444      END DO
1445
1446
1447      !  Set surface RH to the value from 1000 mb.
1448      DO j=1,ns-1
1449         DO i=1,ew-1
1450            rh2(i,1,j) = rh2(i,k00,j)
1451         END DO
1452      END DO
1453
1454     !  Modification of tropical storm complete.
1455     PRINT '(A,I3,A)'       ,'         Bogus storm number ',nstrm,' completed.'
1456
1457    do j = 1,ns-1
1458       do k = 1,nz
1459          do i = 1,ew
1460             u1(i,k,j) =  u2(i,k,j)
1461             grid%u_gc(i,k,j) = u2(i,k,j)
1462          end do
1463       end do
1464    end do
1465
1466    do j = 1,ns
1467       do k = 1,nz
1468          do i = 1,ew-1
1469             v1(i,k,j)   = v2(i,k,j)
1470             grid%v_gc(i,k,j) = v2(i,k,j)
1471          end do
1472       end do
1473    end do
1474
1475     do j = 1,ns-1
1476       do k = 1,nz
1477          do i = 1,ew-1
1478             t1(i,k,j)   = t2(i,k,j)
1479             grid%t_gc(i,k,j) = t2(i,k,j)
1480             rh1(i,k,j)  = rh2(i,k,j)
1481             grid%rh_gc(i,k,j)  = rh2(i,k,j)
1482             phi1(i,k,j) = phi2(i,k,j)
1483             grid%ght_gc(i,k,j) = z2(i,k,j)
1484          END DO
1485       END DO
1486    END DO
1487
1488
1489 END DO all_storms
1490  deallocate(u11)
1491  deallocate(v11)
1492  deallocate(t11)
1493  deallocate(rh11)
1494  deallocate(phi11)
1495  deallocate(u1)
1496  deallocate(v1)
1497  deallocate(t1)
1498  deallocate(rh1)
1499  deallocate(phi1)
1500
1501  do j = 1,ns-1
1502     do i = 1,ew-1
1503        if(grid%ht_gc(i,j) .gt. 1) then
1504          grid%p_gc(i,1,j)  = grid%p_gc(i,1,j)  + (pslx(i,j) * 100. - old_slp(i,j))
1505          grid%psfc(i,j) = grid%psfc(i,j) + (pslx(i,j) * 100. - old_slp(i,j))
1506        else
1507          grid%p_gc(i,1,j)  = pslx(i,j) * 100.
1508          grid%psfc(i,j) = pslx(i,j) * 100.
1509        end if
1510     end do
1511  end do
1512
1513 END SUBROUTINE tc_bogus
1514
1515 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1516 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1517
1518    SUBROUTINE rankine(dx,dy,ds,nlvl,vwgt,rmax,vmax,uu,vv,psi,vor)
1519
1520    !  Define analytical bogus vortex
1521
1522       IMPLICIT NONE
1523
1524       INTEGER nlvl
1525       REAL , DIMENSION(nlvl) :: uu, vv, psi, vor
1526       REAL , DIMENSION(nlvl) :: vwgt
1527       REAL :: dx,dy,ds,rmax,vmax
1528
1529       REAL , PARAMETER :: alpha1= 1.
1530       REAL , PARAMETER :: alpha2= -0.75
1531       real :: pi
1532
1533
1534       INTEGER :: k
1535       REAL :: vr , ang , rr , term1 , bb , term2 , alpha
1536
1537
1538       pi = 3.141592653589793
1539       !  Wind component
1540
1541       DO k=1,nlvl
1542          rr = SQRT(dx**2+dy**2)*ds
1543          IF ( rr .LT. rmax ) THEN
1544             alpha = 1.
1545          ELSE IF ( rr .GE. rmax ) THEN
1546             alpha = alpha2
1547          END IF
1548          vr = vmax * (rr/rmax)**(alpha)
1549          IF ( dx.GE.0. ) THEN
1550             ang = (pi/2.) - ATAN2(dy,MAX(dx,1.e-6))
1551             uu(k) = vwgt(k)*(-vr*COS(ang))
1552             vv(k) = vwgt(k)*( vr*SIN(ang))
1553          ELSE IF ( dx.LT.0. ) THEN
1554             ang = ((3.*pi)/2.) + ATAN2(dy,dx)
1555             uu(k) = vwgt(k)*(-vr*COS(ang))
1556             vv(k) = vwgt(k)*(-vr*SIN(ang))
1557          END IF
1558       END DO
1559
1560       !  psi
1561
1562       DO k=1,nlvl
1563          rr = SQRT(dx**2+dy**2)*ds
1564          IF ( rr .LT. rmax ) THEN
1565             psi(k) = vwgt(k) * (vmax*rr*rr)/(2.*rmax)
1566          ELSE IF ( rr .GE. rmax ) THEN
1567             IF (alpha1.EQ.1.0 .AND. alpha2.eq.-1.0) THEN
1568                psi(k) = vwgt(k) * vmax*rmax*(0.5+LOG(rr/rmax))
1569             ELSE IF (alpha1.EQ.1.0 .AND. alpha2.NE.-1.0) THEN
1570                term1 = vmax/(rmax**alpha1)*(rmax**(alpha1+1)/(alpha1+1))
1571                bb    = (rr**(alpha2+1)/(alpha2+1))-(rmax**(alpha2+1))/(alpha2+1)
1572                term2 = vmax/(rmax**alpha2)*bb
1573                psi(k) = vwgt(k) * (term1 + term2)
1574             END IF
1575          END IF
1576       END DO
1577
1578       ! vort
1579
1580       DO k=1,nlvl
1581          rr = SQRT(dx**2+dy**2)*ds
1582          IF ( rr .LT. rmax ) THEN
1583             vor(k) = vwgt(k) * (2.*vmax)/rmax
1584          ELSE IF ( rr .GE. rmax ) THEN
1585             vor(k) = vwgt(k) * ( (vmax/rmax**alpha2)*(rr**(alpha2-1.))*(1.+alpha2) )
1586          END IF
1587       END DO
1588
1589    END SUBROUTINE rankine
1590
1591 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1592 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1593
1594    SUBROUTINE vor(uin,vin,msfu,msfv,msfm,ew,ns,nz,ds,vort)
1595
1596 !Here we assume that the U and V's are still on the WRF staggered grid.
1597 !The vorticity is then calculated at the mass points on the WRF grid.
1598
1599
1600       IMPLICIT NONE
1601
1602       INTEGER :: jp1,jm1,ip1,im1,i,j,k
1603       INTEGER :: ns, ew, nz, k1
1604
1605       REAL , DIMENSION(ew,nz,ns-1)   :: uin   !u values on unstaggered U grid
1606       REAL , DIMENSION(ew-1,nz,ns)   :: vin   !v values on unstaggered V grid
1607       REAL , DIMENSION(ew-1,nz,ns-1) :: vort  !vort is defined on the mass points
1608
1609       REAL , DIMENSION(ew,ns-1)    :: msfu  !map scale factors on U staggered grid
1610       REAL , DIMENSION(ew-1,ns)    :: msfv  !map scale factors on V staggered grid
1611       REAL , DIMENSION(ew-1,ns-1)  :: msfm  !map scale factors on unstaggered grid
1612
1613       real :: u(ew,ns-1),v(ew-1,ns)
1614
1615
1616       REAL :: ds
1617
1618       REAL :: dsx,dsy , u1 , u2 , u3 , u4 , v1 , v2 , v3 , v4
1619       real :: dudy,dvdx,mm
1620
1621
1622       vort(:,:,:) = -999.
1623       do k = 1,nz
1624
1625          do j = 1,ns-1
1626             do i = 1,ew
1627                u(i,j) = uin(i,k,j)
1628             end do
1629          end do
1630
1631
1632          do j = 1,ns
1633             do i = 1,ew-1
1634                v(i,j) = vin(i,k,j)
1635             end do
1636          end do
1637
1638 !Our indicies are from 2 to ns-2 and ew-2.  This is because out
1639 !map scale factors are not defined for the entire grid.
1640          do j = 2,ns-2
1641             do i = 2,ew-2
1642                mm = msfm(i,j) * msfm(i,j)
1643                u1 = u(i  ,j-1)/msfu(i  ,j-1)
1644                u2 = u(i+1,j-1)/msfu(i+1,j-1)
1645                u3 = u(i+1,j+1)/msfu(i+1,j+1)
1646                u4 = u(i  ,j+1)/msfu(i  ,j+1)
1647                dudy = mm * (u4 + u3 -(u1 + u2)) /(4*ds)
1648
1649                v1 = v(i-1,j  )/msfv(i-1,j)
1650                v2 = v(i+1,j  )/msfv(i+1,j)
1651                v3 = v(i-1 ,j+1)/msfv(i-1,j+1)
1652                v4 = v(i+1,j+1)/msfv(i+1,j+1)
1653                dvdx = mm * (v4 + v2 - (v1 + v3))/(4*ds)
1654
1655                vort(i,k,j) = dvdx - dudy
1656             end do
1657          end do
1658 !Our vorticity array goes out to ew-1 and ns-1 which is the
1659 !mass point grid dimensions.
1660          do i = 2,ew-2
1661             vort(i,k,1)    = vort(i,k,2)    !bottom not corners
1662             vort(i,k,ns-1) = vort(i,k,ns-2) !top not corners
1663          end do
1664
1665          do j = 1,ns-1
1666             vort(ew-1,k,j) = vort(ew-2,k,j) !right side including corners
1667             vort(1,k,j)    = vort(2,k,j)    !left side including corners
1668          end do
1669
1670      end do ! this is the k loop end
1671
1672    END SUBROUTINE
1673
1674 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1675 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1676
1677    SUBROUTINE diverg(uin,vin,msfu,msfv,msfm,ew,ns,nz,ds,div)
1678
1679    !  Computes divergence on unstaggered grid.  The divergence is calculated
1680    !  at the mass points on the WRF grid.
1681    !  div = m*m (du/dx + dv/dy)
1682
1683       IMPLICIT NONE
1684
1685       INTEGER :: jp1,jm1,ip1,im1,i,j,k
1686       INTEGER :: ns, ew, nz, k1
1687
1688       REAL , DIMENSION(ew,nz,ns-1)   :: uin   !u values on unstaggered U grid
1689       REAL , DIMENSION(ew-1,nz,ns)   :: vin   !v values on unstaggered V grid
1690       REAL , DIMENSION(ew-1,nz,ns-1) :: div   !divergence is calculate on the mass points
1691       REAL , DIMENSION(ew,ns-1)    :: msfu  !map scale factors on U staggered grid
1692       REAL , DIMENSION(ew-1,ns)    :: msfv  !map scale factors on V staggered grid
1693       REAL , DIMENSION(ew-1,ns-1)  :: msfm  !map scale factors on unstaggered grid
1694
1695       real :: u(ew,ns-1),v(ew-1,ns)
1696
1697
1698       REAL :: ds
1699
1700       REAL :: dsr,u1,u2,v1,v2
1701       real :: dudx,dvdy,mm,arg1,arg2
1702
1703       dsr = 1/ds
1704       do k = 1,nz
1705
1706          do j = 1,ns-1
1707             do i = 1,ew
1708                u(i,j) = uin(i,k,j)
1709             end do
1710          end do
1711
1712
1713          do j = 1,ns
1714             do i = 1,ew-1
1715                v(i,j) = vin(i,k,j)
1716             end do
1717          end do
1718 !Our indicies are from 2 to ns-2 and ew-2.  This is because out
1719 !map scale factors are not defined for the entire grid.
1720          do j = 2,ns-2
1721             do i = 2,ew-2
1722                mm = msfm(i,j) * msfm(i,j)
1723                u1 = u(i+1,j)/msfu(i+1,j)
1724                u2 = u(i  ,j)/msfu(i  ,j)
1725
1726                v1 = v(i,j+1)/msfv(i,j+1)
1727                v2 = v(i,j)  /msfv(i,j)
1728
1729                div(i,k,j) = mm * (u1 - u2 + v1 - v2) * dsr
1730             end do
1731           end do
1732
1733 !Our divergence array is defined on the mass points.
1734          do i = 2,ew-2
1735             div(i,k,1)    = div(i,k,2)    !bottom not corners
1736             div(i,k,ns-1) = div(i,k,ns-2) !top not corners
1737          end do
1738
1739          do j = 1,ns-1
1740             div(ew-1,k,j) = div(ew-2,k,j) !right side including corners
1741             div(1,k,j)    = div(2,k,j)    !left side including corners
1742          end do
1743
1744      end do !end for the k loop
1745
1746    END SUBROUTINE diverg
1747
1748 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1749 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1750
1751    SUBROUTINE mxratprs (rh, t, ppa, ew, ns, nz, q, min_RH_value)
1752
1753
1754       IMPLICIT NONE
1755
1756       INTEGER   :: i , ew , j , ns , k , nz
1757
1758
1759       REAL      :: min_RH_value
1760       REAL      :: ppa(ew-1,nz,ns-1)
1761       REAL      :: p( ew-1,nz,ns-1 )
1762       REAL      :: q (ew-1,nz,ns-1),rh(ew-1,nz,ns-1),t(ew-1,nz,ns-1)
1763
1764       REAL      :: es
1765       REAL      :: qs
1766       REAL      :: cp              = 1004.0
1767       REAL      :: svp1,svp2,svp3
1768       REAL      :: celkel
1769       REAL      :: eps
1770
1771
1772       !  This function is designed to compute (q) from basic variables
1773       !  p (mb), t(K) and rh(0-100%) to give (q) in (kg/kg).
1774
1775
1776       p = ppa * 0.01
1777
1778       DO j = 1, ns - 1
1779          DO k = 1, nz
1780             DO i = 1, ew - 1
1781                   rh(i,k,j) = MIN ( MAX ( rh(i,k,j) ,min_RH_value ) , 100. )
1782             END DO
1783         END DO
1784      END DO
1785
1786       svp3   =  29.65
1787       svp1   =  0.6112
1788       svp2   =  17.67
1789       celkel =  273.15
1790          eps =  0.622
1791
1792       DO j = 1, ns-1
1793          DO k = 1, nz
1794             DO i = 1,ew-1
1795                es = svp1 * 10. * EXP(svp2 * (t(i,k,j) - celkel ) / (t(i,k,j) - svp3 ))
1796                qs = eps * es / (p(i,k,j) - es)
1797                q(i,k,j) = MAX(0.01 * rh(i,k,j) * qs,0.0)
1798             END DO
1799          END DO
1800       END DO
1801
1802    END SUBROUTINE mxratprs
1803
1804 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1805 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1806 SUBROUTINE mass2_Ustag(field,dim1,dim2,dim3)
1807
1808    IMPLICIT NONE
1809
1810    INTEGER :: dim1 , dim2 , dim3
1811    REAL , DIMENSION(dim1,dim2,dim3) :: field,dummy
1812
1813    dummy = 0.0
1814    dummy(:,2:dim2-1,:)         = ( field(:,1:dim2-2,:) + &
1815                                    field(:,2:dim2-1,:) ) * 0.5
1816    dummy(:,1,:)                = field(:,1,:)
1817    dummy(:,dim2,:)             = field(:,dim2-1,:)
1818
1819    field                       =   dummy
1820
1821 END SUBROUTINE mass2_Ustag
1822
1823 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1824 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1825 SUBROUTINE mass2_Vstag(field,dim1,dim2,dim3)
1826
1827    IMPLICIT NONE
1828
1829    INTEGER :: dim1 , dim2 , dim3
1830    REAL , DIMENSION(dim1,dim2,dim3) :: field,dummy
1831
1832    dummy = 0.0
1833    dummy(2:dim1-1,:,:)         = ( field(1:dim1-2,:,:) + &
1834                                    field(2:dim1-1,:,:) ) * 0.5
1835    dummy(1,:,:)                = field(1,:,:)
1836    dummy(dim1,:,:)             = field(dim1-1,:,:)
1837
1838    field                       =   dummy
1839
1840 END SUBROUTINE mass2_Vstag
1841
1842
1843 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1844 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1845
1846    SUBROUTINE relax (chi, ff, rd, ew, ns, ds, smallres, alpha)
1847
1848       IMPLICIT NONE
1849
1850       INTEGER, PARAMETER    :: mm = 20000
1851
1852       INTEGER               :: i
1853       INTEGER               :: ie
1854       INTEGER               :: ew  !ew direction
1855       INTEGER               :: iter
1856       INTEGER               :: j
1857       INTEGER               :: je
1858       INTEGER               :: jm
1859       INTEGER               :: ns  !ns direction
1860       INTEGER               :: mi
1861
1862       REAL                  :: alpha
1863       REAL                  :: alphaov4
1864       REAL                  :: chi(ew-1,ns-1)
1865       REAL                  :: chimx(ns-1)
1866       REAL                  :: ds
1867       REAL                  :: epx
1868       REAL                  :: fac
1869       REAL                  :: ff(ew-1,ns-1)
1870       REAL                  :: rd(ew-1,ns-1)
1871       REAL                  :: rdmax(ns-1)
1872       REAL                  :: smallres
1873
1874       LOGICAL               :: converged = .FALSE.
1875
1876       fac = ds * ds
1877       alphaov4 = alpha * 0.25
1878
1879       ie=ew-2
1880       je=ns-2
1881
1882       DO j = 1, ns-1
1883          DO i = 1, ew-1
1884             ff(i,j) = fac * ff(i,j)
1885             rd(i,j) = 0.0
1886          END DO
1887       END DO
1888
1889       iter_loop : DO iter = 1, mm
1890          mi = iter
1891          chimx = 0.0
1892
1893
1894          DO j = 2, ns-1
1895             DO i = 2, ew-1
1896                chimx(j) = MAX(ABS(chi(i,j)),chimx(j))
1897             END DO
1898          END DO
1899
1900          epx = MAXVAL(chimx) * SMALLRES * 4.0 / alpha
1901
1902          DO j = 2, ns-2
1903             DO i = 2, ew-2
1904                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)
1905                chi(i,j) = chi(i,j) + rd(i,j) * alphaov4
1906             END DO
1907          END DO
1908
1909          rdmax = 0.0
1910
1911          DO j = 2, ns-2
1912             DO i = 2, ew-2
1913                rdmax(j) = MAX(ABS(rd(i,j)),rdmax(j))
1914             END DO
1915          END DO
1916
1917
1918          IF (MAXVAL(rdmax) .lt. epx) THEN
1919             converged = .TRUE.
1920             EXIT iter_loop
1921          END IF
1922
1923       END DO iter_loop
1924
1925       IF (converged ) THEN
1926 !        PRINT '(A,I5,A)','Relaxation converged in ',mi,' iterations.'
1927       ELSE
1928          PRINT '(A,I5,A)','Relaxation did not converge in',mm,' iterations.'
1929          STOP 'no_converge'
1930       END IF
1931
1932
1933       do i = 2,ew-2
1934             chi(i,ns-1) = chi(i,ns-2) !top not including corners
1935             chi(i,1)    = chi(i,2)    !bottom not including corners
1936       end do
1937
1938       do j = 2,ns-2
1939             chi(ew-1,j) = chi(ew-2,j) !right side not including corners
1940             chi(1,j)    = chi(2,j)    !left side not including corners
1941       end do
1942
1943  !Fill in the corners
1944       chi(1,1)       = chi(2,1)
1945       chi(ew-1,1)    = chi(ew-2,1)
1946       chi(1,ns-1)    = chi(2,ns-1)
1947       chi(ew-1,ns-1) = chi(ew-2,ns-1)
1948
1949
1950
1951    END SUBROUTINE relax
1952 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1953 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1954    SUBROUTINE geowind(height,ew,ns,nz,ds,ug,vg)
1955
1956       IMPLICIT NONE
1957
1958       !     input       height   geopotential on wrf mass grid points
1959       !                 ns       wrf staggered V dimension n-s
1960       !                 ew       wrf staggered U dimension e-w
1961       !                 nz       number of vertical levels
1962       !
1963       !     output      ug       u component of geo wind on wrf staggered V points
1964       !                 vg       v component of geo wind on wrf staggered U points
1965
1966       INTEGER :: ew , ns , nz
1967       REAL :: ds
1968       REAL , DIMENSION(ew-1,nz,ns-1) :: height
1969       REAL , DIMENSION(ew,nz,ns-1) :: ug
1970       REAL , DIMENSION(ew-1,nz,ns) :: vg
1971
1972       REAL :: ds2r , h1 , h2 , h3 , h4, ds4r
1973       INTEGER :: i , j , k
1974
1975       ds4r=1./(4.*ds)
1976
1977 ! The height field comes in on the WRF mass points.
1978
1979
1980
1981 ! ug is the derivative of height in the ns direction  ug = -dheight/dy
1982       ug(:,:,:) = -999.
1983       do j=2,ns-2
1984          do k=1,nz
1985             do i=2,ew-1
1986               h1 = height(i,k,j+1)
1987               h2 = height(i-1,k,j+1)
1988               h3 = height(i  ,k,j-1)
1989               h4 = height(i-1,k,j-1)
1990               ug(i,k,j) = -( (h1 + h2) - ( h3 + h4) ) * ds4r
1991            end do
1992         end do
1993       end do
1994
1995        do i = 2,ew-1
1996           ug(i,:,1)    = ug(i,:,2)    !bottom not including corner points
1997           ug(i,:,ns-1) = ug(i,:,ns-2) !top not including corner points
1998        end do
1999
2000        do j = 2,ns-2
2001           ug(1,:,j)  = ug(2,:,j)    !left side
2002           ug(ew,:,j) = ug(ew-1,:,j) !right side
2003        end do
2004
2005        ug(1,:,1)     = ug(2,:,1)         !Lower left hand corner
2006        ug(1,:,ns-1)  = ug(2,:,ns-1)      !upper left hand corner
2007        ug(ew,:,1)    = ug(ew-1,:,1)      !Lower right hand corner
2008        ug(ew,:,ns-1) = ug(ew-1,:,ns-1)   !upper right hand corner
2009
2010
2011 ! ug is the derivative of height in the ns direction  vg = dheight/dx
2012     vg(:,:,:) = -999.
2013     DO j=2,ns-1
2014        DO k=1,nz
2015           DO i=2,ew-2
2016               h1 = height(i+1,k,j)
2017               h2 = height(i-1,k,j)
2018               h3 = height(i+1,k,j-1)
2019               h4 = height(i-1,k,j-1)
2020               vg(i,k,j) = ( (h1 + h3) - ( h2 + h4) ) * ds4r
2021           end do
2022        end do
2023     end do
2024
2025     do i = 2,ew-2
2026        vg(i,:,1)  = vg(i,:,2)    !bottom not including corner points
2027        vg(i,:,ns) = vg(i,:,ns-1) !top not including corner points
2028     end do
2029
2030     do j = 2,ns-1
2031        vg(1,:,j)    = vg(2,:,j)    !left side not including corner points
2032        vg(ew-1,:,j) = vg(ew-2,:,j) !right side not including corner points
2033    end do
2034
2035    vg(1,:,1)     = vg(2,:,1)        !Lower left hand corner
2036    vg(1,:,ns)    = vg(2,:,ns)       !upper left hand corner
2037    vg(ew-1,:,1)  = vg(ew-2,:,1)     !Lower right hand corner
2038    vg(ew-1,:,ns) = vg(ew-2,:,ns)    !upper right hand corner
2039
2040
2041    END SUBROUTINE geowind
2042 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2043 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2044
2045    SUBROUTINE balance (f,psi,ew,ns,ds,out)
2046
2047    !  Calculates the forcing terms in balance equation
2048
2049    IMPLICIT NONE
2050
2051       !  f       coriolis force
2052       !  psi     stream function
2053       !  ew, ns  grid points in east west, north south direction, respectively
2054       !  ds      grid distance
2055       !  out     output array
2056
2057       INTEGER :: ew , ns,nslast,ewlast,ifill
2058       REAL , DIMENSION(ew-1,ns-1) :: f,psi,out
2059       REAL :: ds
2060
2061       REAL :: psixx , psiyy , psiy , psix, psixy
2062       REAL :: dssq , ds2 , dssq4,arg1,arg2,arg3,arg4
2063
2064       INTEGER :: i , j
2065
2066       dssq  = ds * ds
2067       ds2   = ds * 2.
2068       dssq4 = ds * ds * 4.
2069
2070 !The forcing terms are calculated on the WRF mass points.
2071       out(:,:) = -999.0
2072       DO j=2,ns-2
2073          DO i=2,ew-2
2074             psiyy = ( psi(i,j+1) + psi(i,j-1) - 2.*psi(i,j) ) / dssq
2075             psixx = ( psi(i+1,j) + psi(i-1,j) - 2.*psi(i,j) ) / dssq
2076             psiy  = ( psi(i,j+1) - psi(i,j-1) ) / ds2
2077             psix  = ( psi(i+1,j) - psi(i-1,j) ) / ds2
2078             psixy = ( psi(i+1,j+1)+psi(i-1,j-1)-psi(i-1,j+1)-psi(i+1,j-1)) / dssq4
2079
2080             arg1  = f(i,j)* (psixx+psiyy)
2081             arg2  = ( ( f(i,j+1) - f(i,j-1)) / ds2 ) * psiy
2082             arg3  = ( ( f(i+1,j) - f(i-1,j)) / ds2 ) * psix
2083             arg4  = 2 *(psixy*psixy-psixx*psiyy)
2084             out(i,j)= arg1 + arg2  + arg3 - arg4
2085          END DO
2086       END DO
2087
2088       do i = 2,ew-2
2089             out(i,ns-1) = out(i,ns-2) !top not including corners
2090             out(i,1)    = out(i,2)    !bottom not including corners
2091       end do
2092
2093       do j = 2,ns-2
2094             out(ew-1,j) = out(ew-2,j) !right side not including corners
2095             out(1,j)    = out(2,j)    !left side not including corners
2096       end do
2097
2098  !Fill in the corners
2099       out(1,1)       = out(2,1)
2100       out(ew-1,1)    = out(ew-2,1)
2101       out(1,ns-1)    = out(2,ns-1)
2102       out(ew-1,ns-1) = out(ew-2,ns-1)
2103
2104    END SUBROUTINE balance
2105
2106 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2107 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2108
2109    SUBROUTINE qvtorh ( q , t , p , k00, ew , ns , nz , rh, min_RH_value   )
2110
2111       IMPLICIT NONE
2112
2113       INTEGER , INTENT(IN) :: ew , ns , nz , k00
2114       REAL , INTENT(IN) ,  DIMENSION(ew-1,nz,ns-1) :: q ,t, p
2115       REAL , INTENT(OUT) , DIMENSION(ew-1,nz,ns-1) :: rh
2116
2117       real    min_RH_value
2118
2119       !  Local variables.
2120
2121       INTEGER :: i , j , k,fill
2122       REAL      :: es
2123       REAL      :: qs
2124       REAL      :: cp              = 1004.0
2125       REAL      :: svp1,svp2,svp3
2126       REAL      :: celkel
2127       REAL      :: eps
2128       svp3   =  29.65
2129       svp1   =  0.6112
2130       svp2   =  17.67
2131       celkel =  273.15
2132          eps =  0.622
2133
2134       DO j = 1 , ns - 1
2135          DO k = k00 , nz
2136             DO i = 1 , ew -1
2137                es = svp1 * 10. * EXP(svp2 * (t(i,k,j) - celkel ) / (t(i,k,j) - svp3 ))
2138                qs = eps*es/(0.01*p(i,k,j) - es)
2139                rh(i,k,j) = MIN ( 100. , MAX ( 100.*q(i,k,j)/qs , min_RH_value ) )
2140             END DO
2141          END DO
2142       END DO
2143
2144    END SUBROUTINE qvtorh
2145
2146 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2147 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2148
2149    SUBROUTINE stagger_rankine_winds(utcp,vtcp,ew,ns,nz)
2150
2151
2152 !utcp and vtcp are the output winds of the rankine subroutine
2153 !The winds are calculated on the mass points of the WRF grid
2154 !so they need to be staggered out to the WRF staggering.
2155 !The utcp is placed on the staggered ew wind grid.
2156 !The vtcp is placed on the staggered ns wind grid.
2157
2158 !ew is the full grid dimension in the ew direction.
2159 !ns is the full grid dimension in the ns direction.
2160
2161 !nz is the number of vertical levels.
2162
2163  INTEGER :: ew, ns, nz, i,k,j
2164  REAL utcp(ew,nz,ns-1),  vtcp(ew-1,nz,ns)
2165
2166 !----------------------------------------------------
2167 !Stagger UTCP
2168   DO j=1,ns-1
2169      DO i=2,ew-1
2170         DO k=1,nz
2171            utcp(i,k,j)  = ( utcp(i-1,k,j) + utcp(i,k,j) ) /2
2172         end do
2173     end do
2174   end do
2175
2176 !Fill in U's along the left and right side.
2177  do j = 1,ns
2178     utcp(1,:,j)  = utcp(2,:,j)
2179     utcp(ew,:,j) = utcp(ew-1,:,j)
2180  end do
2181
2182
2183 !Stagger VTCP
2184   DO j=2,ns-1
2185      DO i=1,ew-1
2186         DO k=1,nz
2187            vtcp(i,k,j)  = ( vtcp(i,k,j+1) + vtcp(i,k,j-1) ) /2
2188         end do
2189     end do
2190   end do
2191
2192 !Fill in V's along the bottom and bottom.
2193   do i = 1,ew
2194      vtcp(i,:,1)  = vtcp(i,:,2)
2195      vtcp(i,:,ns) = vtcp(i,:,ns-1)
2196   end do
2197
2198
2199    END SUBROUTINE stagger_rankine_winds
2200
2201 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2202 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2203
2204   subroutine final_ew_velocity(u2,u1,chi,psi,utcr,dx,ew,ns,nz)
2205
2206
2207   integer :: ew,ns,nz,i,j,k
2208   real :: u1(ew,nz,ns-1),utcr(ew,nz,ns-1)
2209   real :: psi(ew-1,nz,ns-1),chi(ew-1,nz,ns-1)
2210 ! input arrays:
2211 !       u1 is the original wind coming in from the metgrid file.
2212 !       utcr is the rankine winds rotated to the map projection put on u WRF staggered grid points.
2213
2214 ! psi comes in on the WRF mass points.  psi is the perturbation field
2215 ! calculated from the relaxation of the vorticity.
2216
2217 ! chi is the relaxation of the divergent winds on WRF mass points.
2218
2219
2220 ! ew is the grid dimension of the WRF ew staggered wind
2221 ! ns is the grid dimension of the WRF ns staggered wind
2222 ! nz is the number of vertical levels
2223 ! dx is the grid spacing
2224 !-------------------------------------------------------------------------------------------
2225
2226   real :: u2(ew,nz,ns-1)
2227 ! output array u2 is the new wind in the ew direction. Is is on WRF staggering.
2228 !-------------------------------------------------------------------------------------------
2229
2230   real upos(ew,nz,ns-1),u0(ew,nz,ns-1),uchi(ew,nz,ns-1)
2231 ! upos is the derivative of psi in the ns direction  u = -dpsi/dy
2232 ! u0 is the background ew velocity
2233 ! uchi is the array chi on the u staggered grid.
2234
2235   real    :: dx,arg1,arg2
2236
2237 !-------------------------------------------------------------
2238 !Take the derivative of chi in the ew direction.
2239    uchi(:,:,:) = -999.
2240    DO k=1,nz !start of k loop
2241       DO j=1,ns-1
2242          DO i=2,ew-1
2243             uchi(i,k,j) = ( chi(i,k,j) - chi(i-1,k,j) )/dx
2244          END DO
2245       END DO
2246
2247       do j = 1,ns-1
2248        uchi(1,k,j)    = uchi(2,k,j)    !fill in the left side
2249        uchi(ew,k,j)   = uchi(ew-1,k,j) !fill in the right side
2250       end do
2251    end do !k loop
2252
2253 !-----------------------------------------------------------------------------------------
2254 ! Take the derivative of psi in the ns direction.
2255     upos(:,:,:) = -999.
2256     DO k=1,nz
2257
2258        DO j=2,ns-2
2259           DO i=2,ew-1
2260               arg1 = psi(i,k,j+1) + psi(i-1,k,j+1)
2261               arg2 = psi(i,k,j-1) + psi(i-1,k,j-1)
2262               upos(i,k,j) = -( arg1 - arg2 )/(4.*dx)
2263           END DO
2264        END DO
2265
2266        do i = 2,ew-1
2267           upos(i,k,1)    = upos(i,k,2)    !bottom not including corner points
2268           upos(i,k,ns-1) = upos(i,k,ns-2) !top not including corner points
2269        end do
2270
2271        do j = 1,ns-2
2272           upos(1,k,j)  = upos(2,k,j)    !left side not including corners
2273           upos(ew,k,j) = upos(ew-1,k,j) !right side not including corners
2274        end do
2275
2276
2277        upos(1,k,1)     = upos(2,k,1)         !Lower left hand corner
2278        upos(1,k,ns-1)  = upos(2,k,ns-1)      !upper left hand corner
2279        upos(ew,k,1)    = upos(ew-1,k,1)      !Lower right hand corner
2280        upos(ew,k,ns-1) = upos(ew-1,k,ns-1)   !upper right hand corner
2281
2282     end do  !k loop for dspi/dy
2283
2284
2285
2286 !-----------------------------------------------------------------------------------------
2287
2288 !  Background u field.
2289 !  Subtract the first quess wind field from the original winds.
2290    do j=1,ns-1
2291       do k=1,nz
2292          do i=1,ew
2293             u0(i,k,j) = u1(i,k,j)-(upos(i,k,j)+uchi(i,k,j))
2294          end do
2295       end do
2296    end do
2297
2298
2299 !   Calculate the final velocity
2300     do j=1,ns-1
2301        do k=1,nz
2302           do i=1,ew
2303              u2(i,k,j) = u0(i,k,j)+utcr(i,k,j)
2304           end do
2305        end do
2306     end do
2307
2308  end subroutine final_ew_velocity
2309
2310 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2311 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2312
2313   subroutine final_ns_velocity(v2,v1,chi,psi,vtcr,dx,ew,ns,nz)
2314
2315
2316   integer :: ew,ns,nz,i,j,k
2317   real :: v1(ew-1,nz,ns),vtcr(ew-1,nz,ns)
2318   real :: psi(ew-1,nz,ns-1),chi(ew-1,nz,ns-1)
2319 ! input arrays:
2320 !       v1 is the original wind coming in from the metgrid file.
2321 !       vtcr is the is the rankine winds rotated to the map projection put on v WRF staggered grid points.
2322
2323 ! psi comes on the WRF mass points.  psi is the perturbation field
2324 ! calculated from the relaxation of the vorticity.
2325
2326 ! chi is the relaxation of the divergent winds on WRF mass points.
2327
2328 ! ew is the grid dimension of the WRF ew staggered wind
2329 ! ns is the grid dimension of the WRF ns staggered wind
2330 ! nz is the number of vertical levels
2331
2332
2333   real :: v2(ew-1,nz,ns)
2334 ! output array v2 is the new wind in the ns direction. Is is on WRF staggering.
2335
2336
2337   real vpos(ew-1,nz,ns),v0(ew-1,nz,ns),vchi(ew-1,nz,ns)
2338 ! vpos is the derivative of psi in the ew direction  v = dpsi/dx
2339 ! v0 is the background ns velocity
2340 ! vchi is the relaxation of the divergent wind put on v WRF staggered grid points.
2341
2342   real    :: dx,arg1,arg2
2343
2344
2345 !-----------------------------------------------------------------------------------------
2346  vchi(:,:,:) = -999.0
2347 !The derivative dchi in the ns direction.
2348     do k = 1,nz
2349        DO j=2,ns-1
2350           DO i=1,ew-1
2351               vchi(i,k,j) = ( chi(i,k,j) - chi(i,k,j-1))/dx
2352           END DO
2353        END DO
2354
2355     do i = 1,ew-1
2356        vchi(i,k,1)  = vchi(i,k,2)
2357        vchi(i,k,ns) = vchi(i,k,ns-1)
2358     end do
2359
2360     end do !end of k loop
2361
2362 !-----------------------------------------------------------------------------------------
2363 !Take the derivative of psi in the ew direction.
2364     vpos(:,:,:) = -999.
2365
2366     DO k=1,nz
2367        DO j=2,ns-1
2368           DO i=2,ew-2
2369               arg1 = psi(i+1,k,j) + psi(i+1,k,j-1)
2370               arg2 = psi(i-1,k,j) + psi(i-1,k,j-1)
2371               vpos(i,k,j) =  ( arg1 - arg2 )/(4.*dx)
2372           END DO
2373        END DO
2374
2375        do i = 2,ew-2
2376           vpos(i,k,1)  = vpos(i,k,2)    !bottom not including corner points
2377           vpos(i,k,ns) = vpos(i,k,ns-1) !top not including corner points
2378       end do
2379
2380        do j = 1,ns
2381           vpos(1,k,j)    = vpos(2,k,j)    !left side not including corner points
2382           vpos(ew-1,k,j) = vpos(ew-2,k,j) !right side not including corner points
2383       end do
2384
2385
2386       vpos(1,k,1)     = vpos(2,k,1)        !Lower left hand corner
2387       vpos(1,k,ns)    = vpos(2,k,ns)       !upper left hand corner
2388       vpos(ew-1,k,1)  = vpos(ew-2,k,1)     !Lower right hand corner
2389       vpos(ew-1,k,ns) = vpos(ew-2,k,ns)    !upper right hand corner
2390
2391     END DO!k loop for dspi/dx
2392
2393
2394     do j=1,ns
2395        do k=1,nz
2396           do i=1,ew-1
2397               v0(i,k,j) = v1(i,k,j)-(vpos(i,k,j)+vchi(i,k,j))
2398               if( v0(i,k,j) .gt. 100.) then
2399                 print *,vchi(i,k,j),i,k,j
2400                 stop
2401               end if
2402           end do
2403        end do
2404     end do
2405
2406
2407 !   Calculate the final velocity
2408     do j=1,ns
2409        do k=1,nz
2410           do i=1,ew-1
2411              v2(i,k,j) = v0(i,k,j)+vtcr(i,k,j)
2412           end do
2413        end do
2414     end do
2415
2416     end subroutine final_ns_velocity
2417 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2418 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2419 subroutine final_RH(rh2,rh0,rhmx,strmci,strmcj,rmax_nstrm,ew,ns,nz,k00, &
2420                     dx,ew_gcntr,ns_gcntr,r_vor2)
2421
2422
2423
2424      integer :: ew,ns,nz
2425      real :: rh2(ew-1,nz,ns-1)  !The final output relative humidity.
2426      real :: rh0(ew-1,nz,ns-1)  !First quess rh read from the metgrid input file.
2427      real :: rhmx(nz)
2428      real :: ew_gcntr !ew grid center as returned from the map projection routines.
2429      real :: ns_gcntr !ns grid center as returned from the map projection routines.
2430      real :: dx       !grid spacing
2431      real :: rmax_nstrm
2432
2433
2434 !Local real variables
2435      real :: sum_rh,avg_rh,rh_min,rhbkg,rhbog,r_ratio
2436      real :: rad
2437      real :: rhtc(ew-1,nz,ns-1)
2438
2439      integer :: nct,k00,i,j,k,ew_mvc,ns_mvc
2440      integer :: strmci(nz), strmcj(nz)
2441
2442
2443 !-----------------------------------------------------------------------
2444      DO k=k00,nz
2445         rh_max= rhmx(k)
2446         ew_mvc = strmci(k)
2447         ns_mvc = strmcj(k)
2448
2449
2450         sum_rh = 0.
2451         nct = 0
2452         DO j=1,ns-1
2453            DO i=1,ew-1
2454               rad = SQRT(REAL(i-ew_mvc)**2.+REAL(j-ns_mvc)**2.)*dx
2455               IF ( (rad .LT. r_vor2).AND.(rad .GE. 0.8*r_vor2) ) THEN
2456                   sum_rh = sum_rh + rh0(i,k,j)
2457                   nct = nct + 1
2458               END IF
2459            END DO
2460         END DO
2461         avg_rh = sum_rh/MAX(REAL(nct),1.)
2462
2463         DO j=1,ns-1
2464             DO i=1,ew-1
2465                rh_min = avg_rh
2466                rad = SQRT((REAL(i)-ew_gcntr)**2.+(REAL(j)-ns_gcntr)**2.)*dx
2467                IF ( rad .LE. rmax_nstrm ) THEN
2468                   rhtc(i,k,j) = rh_max
2469                ELSE
2470                   rhtc(i,k,j) = (rmax_nstrm/rad)*rh_max+(1.-(rmax_nstrm/rad))*rh_min
2471                END IF
2472             END DO
2473          END DO
2474      END DO
2475
2476
2477      !  New RH.
2478      DO j=1,ns-1
2479         DO k=k00,nz
2480            DO i=1,ew-1
2481               rhbkg = rh0(i,k,j)
2482               rhbog = rhtc(i,k,j)
2483               rad = SQRT((REAL(i)-ew_mvc)**2.+(REAL(j)-ns_mvc)**2.)*dx
2484                IF ( (rad.GT.rmax_nstrm) .AND. (rad.LE.r_vor2) ) THEN
2485                     r_ratio = (rad-rmax_nstrm)/(r_vor2-rmax_nstrm)
2486                     rh2(i,k,j) = ((1.-r_ratio)*rhbog) + (r_ratio*rhbkg)
2487               ELSE IF (rad .LE. rmax_nstrm ) THEN
2488                     rh2(i,k,j) = rhbog
2489               ELSE
2490                     rh2(i,k,j) = rhbkg
2491               END IF
2492
2493           END DO
2494         END DO
2495     END DO
2496
2497
2498
2499     end subroutine final_RH
2500
2501 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2502 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!