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