main/module_wrf_top.F

   1 !WRF:DRIVER_LAYER:TOP
   2 !
   3
   4 !TBH:  $$$  move this to ../frame?
   5
   6 MODULE module_wrf_top
   7 !<DESCRIPTION>
   8 ! This module defines top-level wrf_init(), wrf_adtl_check(), wrf_run(), and wrf_finalize()
   9 ! routines.
  10 !</DESCRIPTION>
  11
  12    USE module_machine
  13    USE module_domain
  14    USE module_integrate
  15    USE module_driver_constants
  16    USE module_configure
  17    USE module_check_a_mundo
  18
  19    USE module_timing
  20    USE module_wrf_error
  21 #if ( WRFPLUS == 1 )
  22    USE module_adtl_grid_utilities, ONLY : allocate_grid, deallocate_grid, copy_grid_to_s, copy_grid_to_b, &
  23                                           restore_grid, inner_dot_g, add_grid, inner_dot, init_domain_size, &
  24                                           copy_s_to_g_adjtest, copy_g_to_b_adjtest, inner_dot_g_adjtest, &
  25                                           copy_g_to_a_adjtest, inner_dot_a_b_adjtest, zero_out_ad, zero_out_tl, &
  26                                           get_forc_locations, allocate_locations, gen_scenario_matrix, &
  27                                           spot_force_ad, spot_force_tl, spot_force_nl, &
  28                                           extract_ad_der, extract_tl_der, extract_nl_den, extract_nl_num
  29    USE mediation_pertmod_io, ONLY : xtraj_pointer, xtraj_head, xtraj_tail, xtraj_io_initialize, adtl_initialize
  30    USE module_io_domain, ONLY: close_dataset, wrf_inquire_opened
  31 #if (EM_CORE==1)
  32    USE module_state_description, ONLY : SURFDRAGSCHEME, PARAM_FIRST_SCALAR, num_moist, num_tracer
  33 #endif
  34 #endif
  35
  36    USE module_nesting
  37
  38 #ifdef DM_PARALLEL
  39    USE module_dm, ONLY : wrf_dm_initialize,wrf_get_hostid,domain_active_this_task,mpi_comm_allcompute
  40 #if ( WRFPLUS == 1 )
  41    USE module_dm, ONLY : wrf_dm_sum_real
  42 #endif
  43 #else
  44    USE module_dm, ONLY : domain_active_this_task
  45 #endif
  46
  47 #if ( WRFPLUS == 1 )
  48    USE module_date_time, ONLY : current_date, start_date
  49    USE module_io_domain, ONLY : open_w_dataset, output_boundary
  50 #endif
  51
  52    USE module_cpl, ONLY : coupler_on, cpl_finalize, cpl_defdomain
  53
  54    IMPLICIT NONE
  55
  56 #if ( WRFPLUS == 1 )
  57 #include "wrf_io_flags.h"
  58 #endif
  59
  60    REAL    :: time
  61
  62    INTEGER :: loop , &
  63               levels_to_process
  64
  65    TYPE (domain) , POINTER :: keep_grid, grid_ptr, null_domain
  66 #if ( WRFPLUS == 1 )
  67    TYPE (domain) , POINTER :: grid
  68 #endif
  69    TYPE (domain) , pointer :: parent_grid, new_nest
  70    LOGICAL                                :: a_nest_was_opened
  71    TYPE (grid_config_rec_type), SAVE :: config_flags
  72    INTEGER        :: kid, nestid
  73    INTEGER                 :: number_at_same_level
  74    INTEGER                 :: time_step_begin_restart
  75
  76    INTEGER :: max_dom , domain_id , fid , oid , idum1 , idum2 , ierr
  77    INTEGER :: debug_level
  78    LOGICAL :: input_from_file
  79
  80 #ifdef DM_PARALLEL
  81    INTEGER                 :: nbytes
  82    INTEGER, PARAMETER      :: configbuflen = 4* CONFIG_BUF_LEN
  83    INTEGER                 :: configbuf( configbuflen )
  84    LOGICAL , EXTERNAL      :: wrf_dm_on_monitor
  85 #endif
  86
  87    CHARACTER (LEN=256)     :: rstname
  88    CHARACTER (LEN=80)      :: message
  89 #if ( WRFPLUS == 1 )
  90    INTEGER                 :: alarmid
  91    LOGICAL                 :: gradient_out = .FALSE.
  92 #endif
  93    CHARACTER (LEN=256) , PRIVATE :: a_message
  94
  95    INTERFACE
  96      SUBROUTINE Setup_Timekeeping( grid )
  97       USE module_domain
  98       TYPE(domain), POINTER :: grid
  99      END SUBROUTINE Setup_Timekeeping
 100
 101 ! #if (EM_CORE == 1)
 102      SUBROUTINE wrf_dfi_write_initialized_state( )
 103      END SUBROUTINE wrf_dfi_write_initialized_state
 104
 105      SUBROUTINE wrf_dfi_startfwd_init( )
 106      END SUBROUTINE wrf_dfi_startfwd_init
 107
 108      SUBROUTINE wrf_dfi_startbck_init( )
 109      END SUBROUTINE wrf_dfi_startbck_init
 110
 111      SUBROUTINE wrf_dfi_bck_init( )
 112      END SUBROUTINE wrf_dfi_bck_init
 113
 114      SUBROUTINE wrf_dfi_fwd_init( )
 115      END SUBROUTINE wrf_dfi_fwd_init
 116
 117      SUBROUTINE wrf_dfi_fst_init( )
 118      END SUBROUTINE wrf_dfi_fst_init
 119
 120      SUBROUTINE wrf_dfi_array_reset ( )
 121      END SUBROUTINE wrf_dfi_array_reset
 122 ! #endif
 123
 124      SUBROUTINE med_nest_initial ( parent , grid , config_flags )
 125        USE module_domain
 126        USE module_configure
 127        TYPE (domain), POINTER ::  grid , parent
 128        TYPE (grid_config_rec_type) config_flags
 129      END SUBROUTINE med_nest_initial
 130
 131    END INTERFACE
 132
 133
 134 CONTAINS
 135
 136
 137    SUBROUTINE wrf_init( no_init1 )
 138 !<DESCRIPTION>
 139 !     WRF initialization routine.
 140 !</DESCRIPTION>
 141 #ifdef _OPENMP
 142      use omp_lib
 143 #endif
 144 #ifdef _ACCEL
 145      use accel_lib
 146 #endif
 147      LOGICAL, OPTIONAL, INTENT(IN) :: no_init1
 148      INTEGER i, myproc, nproc, hostid, loccomm, ierr, buddcounter, mydevice, save_comm
 149      INTEGER, ALLOCATABLE :: hostids(:), budds(:)
 150      CHARACTER*512 hostname
 151      CHARACTER*512 mminlu_loc
 152 #ifdef _ACCEL
 153      INTEGER :: it, nt, in, devnum
 154 #endif
 155 #if defined(DM_PARALLEL) && !defined(STUBMPI) && ( defined(RUN_ON_GPU) || defined(_ACCEL))
 156      include "mpif.h"
 157 #endif
 158 #include "version_decl"
 159 #include "commit_decl"
 160
 161
 162 !<DESCRIPTION>
 163 ! Program_name, a global variable defined in frame/module_domain.F, is
 164 ! set, then a routine <a href=init_modules.html>init_modules</a> is
 165 ! called. This calls all the init programs that are provided by the
 166 ! modules that are linked into WRF.  These include initialization of
 167 ! external I/O packages.   Also, some key initializations for
 168 ! distributed-memory parallelism occur here if DM_PARALLEL is specified
 169 ! in the compile: setting up I/O quilt processes to act as I/O servers
 170 ! and dividing up MPI communicators among those as well as initializing
 171 ! external communication packages such as RSL or RSL_LITE.
 172 !
 173 !</DESCRIPTION>
 174
 175    program_name = "WRF " // TRIM(release_version) // " MODEL"
 176
 177    ! Initialize WRF modules:
 178    ! Phase 1 returns after MPI_INIT() (if it is called)
 179    CALL init_modules(1)
 180    IF ( .NOT. PRESENT( no_init1 ) ) THEN
 181      ! Initialize utilities (time manager, etc.)
 182 #ifdef NO_LEAP_CALENDAR
 183      CALL WRFU_Initialize( defaultCalKind=WRFU_CAL_NOLEAP )
 184 #else
 185      CALL WRFU_Initialize( defaultCalKind=WRFU_CAL_GREGORIAN )
 186 #endif
 187    ENDIF
 188    ! Phase 2 resumes after MPI_INIT() (if it is called)
 189    CALL init_modules(2)
 190
 191 !<DESCRIPTION>
 192 ! The wrf namelist.input file is read and stored in the USE associated
 193 ! structure model_config_rec, defined in frame/module_configure.F, by the
 194 ! call to <a href=initial_config.html>initial_config</a>.  On distributed
 195 ! memory parallel runs this is done only on one processor, and then
 196 ! broadcast as a buffer.  For distributed-memory, the broadcast of the
 197 ! configuration information is accomplished by first putting the
 198 ! configuration information into a buffer (<a
 199 ! href=get_config_as_buffer.html>get_config_as_buffer</a>), broadcasting
 200 ! the buffer, then setting the configuration information (<a
 201 ! href=set_config_as_buffer.html>set_config_as_buffer</a>).
 202 !
 203 !</DESCRIPTION>
 204
 205 #ifdef DM_PARALLEL
 206    CALL wrf_get_dm_communicator( save_comm )
 207    CALL wrf_set_dm_communicator( mpi_comm_allcompute )
 208    IF ( wrf_dm_on_monitor() ) THEN
 209      CALL initial_config
 210    ENDIF
 211    CALL get_config_as_buffer( configbuf, configbuflen, nbytes )
 212    CALL wrf_dm_bcast_bytes( configbuf, nbytes )
 213    CALL set_config_as_buffer( configbuf, configbuflen )
 214    CALL wrf_dm_initialize
 215    CALL wrf_set_dm_communicator( save_comm )
 216 #else
 217    CALL initial_config
 218 #endif
 219
 220    CALL setup_physics_suite
 221    CALL set_derived_rconfigs
 222    CALL check_nml_consistency
 223    CALL set_physics_rconfigs
 224
 225 #ifdef _ACCEL
 226    buddcounter = 1
 227    mydevice = 0
 228 # if defined(DM_PARALLEL) && !defined(STUBMPI)
 229    CALL wrf_get_myproc( myproc )
 230    CALL wrf_get_nproc( nproc )
 231    CALL wrf_get_hostid ( hostid )
 232    CALL wrf_get_dm_communicator ( loccomm )
 233
 234    ALLOCATE( hostids(nproc) )
 235    ALLOCATE( budds(nproc) )
 236    CALL mpi_allgather( hostid, 1, MPI_INTEGER, hostids, 1, MPI_INTEGER, loccomm, ierr )
 237    if ( ierr .NE. 0 ) print * ,'error in mpi_allgather ',ierr
 238    budds = -1
 239    buddcounter = 0
 240    ! mark the ones i am on the same node with
 241    DO i = 1, nproc
 242       IF ( hostid .EQ. hostids(i) ) THEN
 243          budds(i) = buddcounter
 244          buddcounter = buddcounter + 1
 245       ENDIF
 246    ENDDO
 247    mydevice = budds(myproc+1)
 248    DEALLOCATE( hostids )
 249    DEALLOCATE( budds )
 250 # endif
 251    in = acc_get_num_devices(acc_device_nvidia)
 252    if (in .le. 0) print *, 'error:  No GPUS present: ',in
 253 # ifdef _OPENMP
 254    !$OMP PARALLEL SHARED(mydevice,in) PRIVATE(it,nt,devnum)
 255    it = omp_get_thread_num()
 256    nt = omp_get_num_threads()
 257    devnum = mod(mod(mydevice*nt,in) + it, in)
 258 # ifdef _ACCEL_PROF
 259    print *, "Process, Thread, Device: ",mydevice, it, devnum
 260 # endif
 261    call acc_set_device_num(devnum, acc_device_nvidia)
 262
 263    !$OMP END PARALLEL
 264 # else
 265    it = 0
 266    nt = 1
 267    devnum = mod(mod(mydevice*nt,in) + it, in)
 268 #  ifdef _ACCEL_PROF
 269    print *, "Process, Thread, Device: ",mydevice, it, devnum
 270 #  endif
 271    call acc_set_device_num(devnum, acc_device_nvidia)
 272 # endif
 273 #endif
 274
 275 #ifdef RUN_ON_GPU
 276    CALL wrf_get_myproc( myproc )
 277    CALL wrf_get_nproc( nproc )
 278 # ifdef DM_PARALLEL
 279    CALL wrf_get_hostid ( hostid )
 280    CALL wrf_get_dm_communicator ( loccomm )
 281    ALLOCATE( hostids(nproc) )
 282    ALLOCATE( budds(nproc) )
 283    CALL mpi_allgather( hostid, 1, MPI_INTEGER, hostids, 1, MPI_INTEGER, loccomm, ierr )
 284    IF ( ierr .NE. 0 ) THEN
 285       write(a_message,*)__FILE__,__LINE__,'error in mpi_allgather ',ierr
 286       CALL wrf_message ( a_message )
 287    END IF
 288    budds = -1
 289    buddcounter = 0
 290    ! mark the ones i am on the same node with
 291    DO i = 1, nproc
 292       IF ( hostid .EQ. hostids(i) ) THEN
 293          budds(i) = buddcounter
 294          buddcounter = buddcounter + 1
 295       ENDIF
 296    ENDDO
 297    mydevice = budds(myproc+1)
 298    DEALLOCATE( hostids )
 299    DEALLOCATE( budds )
 300 # else
 301    mydevice = 0
 302 # endif
 303    CALL wsm5_gpu_init( myproc, nproc, mydevice )
 304 #endif
 305
 306 !<DESCRIPTION>
 307 ! Among the configuration variables read from the namelist is
 308 ! debug_level. This is retrieved using nl_get_debug_level (Registry
 309 ! generated and defined in frame/module_configure.F).  The value is then
 310 ! used to set the debug-print information level for use by <a
 311 ! href=wrf_debug.html>wrf_debug</a> throughout the code. Debug_level
 312 ! of zero (the default) causes no information to be printed when the
 313 ! model runs. The higher the number (up to 1000) the more information is
 314 ! printed.
 315 !
 316 !</DESCRIPTION>
 317
 318    CALL nl_get_debug_level ( 1, debug_level )
 319    CALL set_wrf_debug_level ( debug_level )
 320
 321    ! allocated and configure the mother domain
 322
 323    NULLIFY( null_domain )
 324
 325 !<DESCRIPTION>
 326 ! RSL is required for WRF nesting options.
 327 ! The non-MPI build that allows nesting is only supported on machines
 328 ! with the -DSTUBMPI option.  Check to see if the WRF model is being asked
 329 ! for a for a multi-domain run (max_dom > 1, from the namelist).  If so,
 330 ! then we check to make sure that we are under the parallel
 331 ! run option or we are on an acceptable machine.
 332 !</DESCRIPTION>
 333
 334    CALL nl_get_max_dom( 1, max_dom )
 335    IF ( max_dom > 1 ) THEN
 336 #if ( ! defined(DM_PARALLEL)  &&   ! defined(STUBMPI) )
 337    CALL wrf_error_fatal( &
 338      'nesting requires either an MPI build or use of the -DSTUBMPI option' )
 339 #endif
 340    END IF
 341
 342 !<DESCRIPTION>
 343 ! The top-most domain in the simulation is then allocated and configured
 344 ! by calling <a href=alloc_and_configure_domain.html>alloc_and_configure_domain</a>.
 345 ! Here, in the case of this root domain, the routine is passed the
 346 ! globally accessible pointer to TYPE(domain), head_grid, defined in
 347 ! frame/module_domain.F.  The parent is null and the child index is given
 348 ! as negative, signifying none.  Afterwards, because the call to
 349 ! alloc_and_configure_domain may modify the model's configuration data
 350 ! stored in model_config_rec, the configuration information is again
 351 ! repacked into a buffer, broadcast, and unpacked on each task (for
 352 ! DM_PARALLEL compiles). The call to <a
 353 ! href=setup_timekeeping.html>setup_timekeeping</a> for head_grid relies
 354 ! on this configuration information, and it must occur after the second
 355 ! broadcast of the configuration information.
 356 !
 357 !</DESCRIPTION>
 358    CALL       wrf_message ( program_name )
 359    CALL       wrf_message ( commit_version )
 360    CALL       wrf_debug ( 100 , 'wrf: calling alloc_and_configure_domain ' )
 361    CALL alloc_and_configure_domain ( domain_id  = 1 ,                  &
 362                                active_this_task = domain_active_this_task(1), &
 363                                      grid       = head_grid ,          &
 364                                      parent     = null_domain ,        &
 365                                      kid        = -1                   )
 366
 367    CALL       wrf_debug ( 100 , 'wrf: calling model_to_grid_config_rec ' )
 368    CALL model_to_grid_config_rec ( head_grid%id , model_config_rec , config_flags )
 369    CALL       wrf_debug ( 100 , 'wrf: calling set_scalar_indices_from_config ' )
 370    CALL set_scalar_indices_from_config ( head_grid%id , idum1, idum2 )
 371    CALL       wrf_debug ( 100 , 'wrf: calling init_wrfio' )
 372    CALL init_wrfio
 373
 374 #ifdef DM_PARALLEL
 375    CALL wrf_get_dm_communicator( save_comm )
 376    CALL wrf_set_dm_communicator( mpi_comm_allcompute )
 377    CALL get_config_as_buffer( configbuf, configbuflen, nbytes )
 378    CALL wrf_dm_bcast_bytes( configbuf, nbytes )
 379    CALL set_config_as_buffer( configbuf, configbuflen )
 380    CALL wrf_set_dm_communicator( save_comm )
 381 #endif
 382
 383 ! #if (EM_CORE == 1)
 384    ! In case we are doing digital filter initialization, set dfi_stage = DFI_SETUP
 385    !   to indicate in Setup_Timekeeping that we want forecast start and
 386    !   end times at this point
 387    IF ( head_grid%dfi_opt .NE. DFI_NODFI ) head_grid%dfi_stage = DFI_SETUP
 388 ! #endif
 389
 390    CALL Setup_Timekeeping (head_grid)
 391
 392 !<DESCRIPTION>
 393 ! The head grid is initialized with read-in data through the call to <a
 394 ! href=med_initialdata_input.html>med_initialdata_input</a>, which is
 395 ! passed the pointer head_grid and a locally declared configuration data
 396 ! structure, config_flags, that is set by a call to <a
 397 ! href=model_to_grid_config_rec.html>model_to_grid_config_rec</a>.  It is
 398 ! also necessary that the indices into the 4d tracer arrays such as
 399 ! moisture be set with a call to <a
 400 ! href=set_scalar_indices_from_config.html>set_scalar_indices_from_config</a>
 401 ! prior to the call to initialize the domain.  Both of these calls are
 402 ! told which domain they are setting up for by passing in the integer id
 403 ! of the head domain as <tt>head_grid%id</tt>, which is 1 for the
 404 ! top-most domain.
 405 !
 406 ! In the case that write_restart_at_0h is set to true in the namelist,
 407 ! the model simply generates a restart file using the just read-in data
 408 ! and then shuts down. This is used for ensemble breeding, and is not
 409 ! typically enabled.
 410 !
 411 !</DESCRIPTION>
 412
 413  IF ( domain_active_this_task(1) ) THEN
 414    CALL med_initialdata_input( head_grid , config_flags )
 415
 416    IF ( config_flags%write_restart_at_0h ) THEN
 417       CALL med_restart_out ( head_grid, config_flags )
 418 #ifndef AUTODOC_BUILD
 419 ! prevent this from showing up before the call to integrate in the autogenerated call tree
 420       CALL wrf_debug ( 0 , ' 0 h restart only wrf: SUCCESS COMPLETE WRF' )
 421 ! TBH:  $$$ Unscramble this later...
 422 ! TBH:  $$$ Need to add state to avoid calling wrf_finalize() twice when ESMF
 423 ! TBH:  $$$ library is used.  Maybe just set clock stop_time=start_time and
 424 ! TBH:  $$$ do not call wrf_finalize here...
 425       CALL wrf_finalize( )
 426 #endif
 427    END IF
 428   ENDIF  ! domain_active_this_task
 429
 430    ! set default values for subtimes
 431    head_grid%start_subtime = domain_get_start_time ( head_grid )
 432    head_grid%stop_subtime = domain_get_stop_time ( head_grid )
 433
 434  IF ( domain_active_this_task(1) ) THEN
 435    !  For EM (but not DA), if this is a DFI run, we can allocate some space.  We are
 436    !  not allowing anyting tricky for nested DFI.  If there are any nested domains,
 437    !  they all need to start at the same time.  Otherwise, why even do the DFI?  If
 438    !  the domains do not all start at the same time, then there will be inconsistencies,
 439    !  which is what DFI is supposed to address.
 440
 441 #if (EM_CORE == 1)
 442    IF ( head_grid%dfi_opt .NE. DFI_NODFI ) THEN
 443       CALL alloc_doms_for_dfi ( head_grid )
 444    END IF
 445 #endif
 446
 447    IF (coupler_on) CALL cpl_defdomain( head_grid )
 448   ENDIF  ! domain_active_this_task
 449
 450    END SUBROUTINE wrf_init
 451
 452 #if ( WRFPLUS == 1 )
 453 !--------------AD/TL code begin---------------------------------------
 454    SUBROUTINE wrf_adtl_check( )
 455
 456    IF ( config_flags%scenario_type .EQ. 0 ) THEN
 457       CALL wrf_adtl_check_sum
 458    ELSE
 459       CALL wrf_adtl_check_spot
 460    END IF
 461
 462    END SUBROUTINE wrf_adtl_check
 463
 464
 465    SUBROUTINE wrf_adtl_check_sum( )
 466 !<DESCRIPTION>
 467 !     WRF adjoint and tangent linear code check routine.
 468 !</DESCRIPTION>
 469
 470 !<DESCRIPTION>
 471 ! Once the top-level domain has been allocated, configured, and
 472 ! initialized, the model time integration is ready to proceed.
 473 !
 474 !</DESCRIPTION>
 475
 476    REAL :: alpha_m, factor, val_l, val_a, save_l, val_n, coef
 477    INTEGER :: nt, ij, time_step, rc
 478    CHARACTER*256 :: timestr
 479
 480    ! Return immediately if not dyn_em_check
 481    IF ( config_flags%dyn_opt .NE. dyn_em_check ) RETURN
 482
 483    ! Force to turn off history output in this case
 484    CALL WRFU_AlarmRingerOff( head_grid%alarms( HISTORY_ALARM ), rc=rc )
 485
 486    ! Force to read the lateral boundary condition here.
 487    CALL med_before_solve_io ( head_grid, config_flags )
 488
 489    ! Close boundary file, as it will be read again from start later
 490    CALL close_dataset ( head_grid%lbc_fid , config_flags , "DATASET=BOUNDARY" )
 491
 492    CALL init_domain_size ( head_grid, config_flags )
 493
 494    ! Release linked list for trajectory
 495    call xtraj_io_initialize
 496
 497 IF ( config_flags%check_TL .or. config_flags%check_AD ) THEN
 498
 499    ! Save the initial condition and boundary condition, x
 500    CALL allocate_grid ( )
 501
 502    CALL copy_grid_to_s ( head_grid , &
 503                          head_grid%i_start(1), head_grid%i_end(1),                &
 504                          head_grid%j_start(1), head_grid%j_end(1)                 )
 505
 506    CALL       wrf_message ( "wrf: calling nonlinear integrate" )
 507    model_config_rec%dyn_opt = dyn_em
 508
 509    ! Set up basic states output
 510    CALL nl_get_time_step ( head_grid%id, time_step )
 511    CALL nl_set_auxhist6_interval_s ( head_grid%id, time_step )
 512    CALL nl_set_io_form_auxhist6 ( head_grid%id, 2 )
 513    CALL nl_set_frames_per_auxhist6 ( head_grid%id, 1 )
 514    IF ( ASSOCIATED( head_grid%domain_clock ) ) THEN
 515      IF ( head_grid%domain_clock_created ) THEN
 516        CALL WRFU_ClockDestroy( head_grid%domain_clock )
 517        head_grid%domain_clock_created = .FALSE.
 518      ENDIF
 519    ENDIF
 520    IF ( ASSOCIATED( head_grid%alarms ) .AND. &
 521         ASSOCIATED( head_grid%alarms_created ) ) THEN
 522      DO alarmid = 1, MAX_WRF_ALARMS
 523        IF ( head_grid%alarms_created( alarmid ) ) THEN
 524          CALL WRFU_AlarmDestroy( head_grid%alarms( alarmid ) )
 525          head_grid%alarms_created( alarmid ) = .FALSE.
 526        ENDIF
 527      ENDDO
 528    ENDIF
 529    CALL Setup_Timekeeping ( head_grid )
 530
 531    ! Force to turn off history output in this case
 532    CALL WRFU_AlarmRingerOff( head_grid%alarms( HISTORY_ALARM ), rc=rc )
 533    IF ( config_flags%check_TL ) THEN
 534       IF ( config_flags%mp_physics .NE. 0 .and. config_flags%mp_physics .NE. 99 .and. &
 535            config_flags%mp_physics .NE. 98 ) CALL nl_set_mp_physics (head_grid%id, config_flags%mp_physics_ad)
 536       IF ( config_flags%bl_pbl_physics .GT. 0 ) CALL nl_set_bl_pbl_physics (head_grid%id, 98)
 537       IF ( config_flags%cu_physics .GT. 0 ) THEN
 538          CALL nl_set_cu_physics (head_grid%id, 98)
 539          head_grid%cudt = 0
 540       ENDIF
 541       CALL nl_set_ra_lw_physics (head_grid%id, 0)
 542       CALL nl_set_ra_sw_physics (head_grid%id, 0)
 543       CALL nl_set_sf_sfclay_physics (head_grid%id, 0)
 544
 545       ! Force to turn off boundary input as we can only perturb the initial boundary condition.
 546       CALL nl_set_io_form_boundary( head_grid%id, 0 )
 547    ENDIF
 548
 549    CALL wrf_run
 550
 551    ! Turn off basic states output
 552    CALL nl_set_io_form_auxhist6 ( head_grid%id, 0 )
 553    CALL nl_set_auxhist6_interval_s ( head_grid%id, 0 )
 554    if ( .not. head_grid%trajectory_io ) CALL nl_set_inputout_interval_s ( head_grid%id, 0 )
 555    IF ( ASSOCIATED( head_grid%domain_clock ) ) THEN
 556      IF ( head_grid%domain_clock_created ) THEN
 557        CALL WRFU_ClockDestroy( head_grid%domain_clock )
 558        head_grid%domain_clock_created = .FALSE.
 559      ENDIF
 560    ENDIF
 561    IF ( ASSOCIATED( head_grid%alarms ) .AND. &
 562         ASSOCIATED( head_grid%alarms_created ) ) THEN
 563      DO alarmid = 1, MAX_WRF_ALARMS
 564        IF ( head_grid%alarms_created( alarmid ) ) THEN
 565          CALL WRFU_AlarmDestroy( head_grid%alarms( alarmid ) )
 566          head_grid%alarms_created( alarmid ) = .FALSE.
 567        ENDIF
 568      ENDDO
 569    ENDIF
 570    CALL Setup_Timekeeping ( head_grid )
 571
 572    CALL       wrf_message ( "wrf: back from nonlinear integrate" )
 573
 574    IF ( config_flags%check_TL ) THEN
 575
 576    wrf_err_message = '=============== Tangent Linear Check ==================='
 577    CALL wrf_message(TRIM(wrf_err_message))
 578
 579    wrf_err_message = 'check==== U === V === W == PH === T == MU == MOIST == TRACER ==='
 580    CALL wrf_message(TRIM(wrf_err_message))
 581
 582    WRITE(wrf_err_message,'(A,8(1x,L5))') 'check',head_grid%check_u, head_grid%check_v,  &
 583                                head_grid%check_w, head_grid%check_ph, &
 584                                head_grid%check_t, head_grid%check_mu, &
 585                                head_grid%check_moist, head_grid%check_tracer
 586    CALL wrf_message(TRIM(wrf_err_message))
 587
 588    ! Save the f(x)
 589    CALL copy_grid_to_b ( head_grid )
 590
 591    ! Restore the x and assign the \delta x
 592    CALL restore_grid ( head_grid )
 593    CALL copy_s_to_g_adjtest ( head_grid, 1.0 )
 594
 595    ! Set up basic states reading
 596    model_config_rec%auxinput6_inname = "auxhist6_d<domain>_<date>"
 597    CALL nl_get_time_step ( head_grid%id, time_step )
 598    CALL nl_set_auxinput6_interval_s (head_grid%id, time_step )
 599    CALL nl_set_io_form_auxinput6 ( head_grid%id, 2 )
 600    CALL nl_set_frames_per_auxinput6 ( head_grid%id, 1 )
 601
 602    CALL wrf_run_tl
 603
 604    ! Turn off auxinput6 reading
 605    CALL nl_set_auxinput6_interval_s (head_grid%id, 0 )
 606    IF ( ASSOCIATED( head_grid%domain_clock ) ) THEN
 607      IF ( head_grid%domain_clock_created ) THEN
 608        CALL WRFU_ClockDestroy( head_grid%domain_clock )
 609        head_grid%domain_clock_created = .FALSE.
 610      ENDIF
 611    ENDIF
 612    IF ( ASSOCIATED( head_grid%alarms ) .AND. &
 613         ASSOCIATED( head_grid%alarms_created ) ) THEN
 614      DO alarmid = 1, MAX_WRF_ALARMS
 615        IF ( head_grid%alarms_created( alarmid ) ) THEN
 616          CALL WRFU_AlarmDestroy( head_grid%alarms( alarmid ) )
 617          head_grid%alarms_created( alarmid ) = .FALSE.
 618        ENDIF
 619      ENDDO
 620    ENDIF
 621    CALL Setup_Timekeeping ( head_grid )
 622
 623    save_l = 0.0
 624    ! Calculate the inner dot.
 625    !$OMP PARALLEL DO    &
 626    !$OMP DEFAULT (SHARED) PRIVATE ( ij ) &
 627    !$OMP REDUCTION (+:save_l)
 628    DO ij = 1 , head_grid%num_tiles
 629      CALL inner_dot_g ( head_grid , save_l,                                        &
 630                         head_grid%i_start(ij), head_grid%i_end(ij),                &
 631                         head_grid%j_start(ij), head_grid%j_end(ij)                 )
 632    END DO
 633    !$OMP END PARALLEL DO
 634 #ifdef DM_PARALLEL
 635    save_l = wrf_dm_sum_real ( save_l )
 636 #endif
 637
 638    alpha_m = 1.
 639    tangentLinearCheck : DO nt = 1 , 11
 640
 641       alpha_m = 0.1 * alpha_m
 642       factor = 1.0 + alpha_m
 643
 644       CALL add_grid ( head_grid,  factor )
 645
 646       CALL       wrf_message ( "wrf: calling nonlinear integrate" )
 647       model_config_rec%dyn_opt = dyn_em
 648       CALL domain_clock_get( head_grid, start_timestr=timestr )
 649       CALL domain_clock_set( head_grid, current_timestr=timestr )
 650
 651       ! Force to turn off history output in this case
 652       CALL WRFU_AlarmRingerOff( head_grid%alarms( HISTORY_ALARM ), rc=rc )
 653       ! Force to turn off boundary input as it was perturbated in add_grid.
 654       CALL nl_set_io_form_boundary( head_grid%id, 0 )
 655
 656       head_grid%dtbc = 0
 657       head_grid%itimestep = 0
 658       CALL wrf_run
 659       CALL       wrf_message ( "wrf: back from nonlinear integrate" )
 660
 661       val_n = 0.0
 662       !$OMP PARALLEL DO    &
 663       !$OMP PRIVATE ( ij ) &
 664       !$OMP REDUCTION (+:val_n)
 665       DO ij = 1 , head_grid%num_tiles
 666         CALL inner_dot ( head_grid, val_n,                                          &
 667                          head_grid%i_start(ij), head_grid%i_end(ij),                &
 668                          head_grid%j_start(ij), head_grid%j_end(ij)                 )
 669       END DO
 670       !$OMP END PARALLEL DO
 671 #ifdef DM_PARALLEL
 672       val_n = wrf_dm_sum_real ( val_n )
 673 #endif
 674       val_l=save_l*alpha_m**2
 675       coef=val_n/val_l
 676       WRITE(wrf_err_message, FMT='(A,E9.4,A,E23.14,A,E14.7,A,E14.7)') &
 677            'tl_check: alpha_m=',alpha_m,'  coef=',coef, &
 678            '  val_n=',val_n,'  val_l=',val_l
 679       CALL wrf_message(TRIM(wrf_err_message))
 680
 681    ENDDO  tangentLinearCheck
 682
 683    END IF !end of Tangent linear check
 684
 685    IF ( config_flags%check_AD ) THEN
 686
 687    wrf_err_message = '==================== Adjoint check ====================='
 688    CALL wrf_message(TRIM(wrf_err_message))
 689
 690    wrf_err_message = 'check==== U === V === W == PH === T == MU == MOIST == TRACER ==='
 691    CALL wrf_message(TRIM(wrf_err_message))
 692
 693    WRITE(wrf_err_message,'(A,8(1x,L5))') 'check',head_grid%check_u, head_grid%check_v,  &
 694                                head_grid%check_w, head_grid%check_ph, &
 695                                head_grid%check_t, head_grid%check_mu, &
 696                                head_grid%check_moist, head_grid%check_tracer
 697
 698    CALL wrf_message(TRIM(wrf_err_message))
 699
 700    CALL restore_grid ( head_grid )
 701    CALL copy_s_to_g_adjtest ( head_grid, 0.01 )
 702
 703    CALL copy_g_to_b_adjtest ( head_grid )
 704
 705    ! Set up basic states reading
 706    model_config_rec%auxinput6_inname = "auxhist6_d<domain>_<date>"
 707    CALL nl_get_time_step ( head_grid%id, time_step )
 708    CALL nl_set_auxinput6_interval_s (head_grid%id, time_step )
 709    CALL nl_set_io_form_auxinput6 ( head_grid%id, 2 )
 710    CALL nl_set_frames_per_auxinput6 ( head_grid%id, 1 )
 711
 712    IF ( config_flags%mp_physics .NE. 0 .and. config_flags%mp_physics .NE. 99 .and. &
 713         config_flags%mp_physics .NE. 98 ) CALL nl_set_mp_physics (head_grid%id, config_flags%mp_physics_ad)
 714    IF ( config_flags%bl_pbl_physics .GT. 0 ) CALL nl_set_bl_pbl_physics (head_grid%id, 98)
 715    IF ( config_flags%cu_physics .GT. 0 ) THEN
 716       CALL nl_set_cu_physics (head_grid%id, 98)
 717       head_grid%cudt = 0
 718    ENDIF
 719    CALL nl_set_ra_lw_physics (head_grid%id, 0)
 720    CALL nl_set_ra_sw_physics (head_grid%id, 0)
 721    CALL nl_set_sf_sfclay_physics (head_grid%id, 0)
 722
 723    CALL wrf_run_tl
 724
 725    ! Turn off auxinput6 reading
 726    CALL nl_set_auxinput6_interval_s (head_grid%id, 0 )
 727    IF ( ASSOCIATED( head_grid%domain_clock ) ) THEN
 728      IF ( head_grid%domain_clock_created ) THEN
 729        CALL WRFU_ClockDestroy( head_grid%domain_clock )
 730        head_grid%domain_clock_created = .FALSE.
 731      ENDIF
 732    ENDIF
 733    IF ( ASSOCIATED( head_grid%alarms ) .AND. &
 734         ASSOCIATED( head_grid%alarms_created ) ) THEN
 735      DO alarmid = 1, MAX_WRF_ALARMS
 736        IF ( head_grid%alarms_created( alarmid ) ) THEN
 737          CALL WRFU_AlarmDestroy( head_grid%alarms( alarmid ) )
 738          head_grid%alarms_created( alarmid ) = .FALSE.
 739        ENDIF
 740      ENDDO
 741    ENDIF
 742    CALL Setup_Timekeeping ( head_grid )
 743
 744    val_l = 0.0
 745    !$OMP PARALLEL DO    &
 746    !$OMP PRIVATE ( ij ) &
 747    !$OMP REDUCTION (+:val_l)
 748    DO ij = 1 , head_grid%num_tiles
 749      CALL inner_dot_g_adjtest ( head_grid , val_l,                                         &
 750                         head_grid%i_start(ij), head_grid%i_end(ij),                &
 751                         head_grid%j_start(ij), head_grid%j_end(ij)                 )
 752    END DO
 753    !$OMP END PARALLEL DO
 754 #ifdef DM_PARALLEL
 755    val_l =  wrf_dm_sum_real ( val_l )
 756 #endif
 757
 758 !  CALL restore_grid ( head_grid )
 759    !$OMP PARALLEL DO    &
 760    !$OMP DEFAULT (SHARED) PRIVATE ( ij )
 761    DO ij = 1 , head_grid%num_tiles
 762       CALL copy_g_to_a_adjtest ( head_grid,                                           &
 763                         head_grid%i_start(ij), head_grid%i_end(ij),                &
 764                         head_grid%j_start(ij), head_grid%j_end(ij)                 )
 765    END DO
 766    !$OMP END PARALLEL DO
 767
 768    ! Set the file names and interval for reading basic states.
 769    model_config_rec%auxinput6_inname = "auxhist6_d<domain>_<date>"
 770    CALL nl_get_time_step ( head_grid%id, time_step )
 771    CALL nl_set_auxinput6_interval_s (head_grid%id, time_step )
 772    CALL nl_set_io_form_auxinput6 ( head_grid%id, 2 )
 773    CALL nl_set_frames_per_auxinput6 ( head_grid%id, 1 )
 774
 775    CALL wrf_run_ad
 776
 777    ! Turn off auxinput6 reading
 778    CALL nl_set_auxinput6_interval_s (head_grid%id, 0 )
 779    IF ( ASSOCIATED( head_grid%domain_clock ) ) THEN
 780      IF ( head_grid%domain_clock_created ) THEN
 781        CALL WRFU_ClockDestroy( head_grid%domain_clock )
 782        head_grid%domain_clock_created = .FALSE.
 783      ENDIF
 784    ENDIF
 785    IF ( ASSOCIATED( head_grid%alarms ) .AND. &
 786         ASSOCIATED( head_grid%alarms_created ) ) THEN
 787      DO alarmid = 1, MAX_WRF_ALARMS
 788        IF ( head_grid%alarms_created( alarmid ) ) THEN
 789          CALL WRFU_AlarmDestroy( head_grid%alarms( alarmid ) )
 790          head_grid%alarms_created( alarmid ) = .FALSE.
 791        ENDIF
 792      ENDDO
 793    ENDIF
 794    CALL Setup_Timekeeping ( head_grid )
 795
 796    val_a = 0.0
 797    !$OMP PARALLEL DO    &
 798    !$OMP PRIVATE ( ij ) &
 799    !$OMP REDUCTION (+:val_a)
 800    DO ij = 1 , head_grid%num_tiles
 801      CALL inner_dot_a_b_adjtest ( head_grid, val_a,                                  &
 802                           head_grid%i_start(ij), head_grid%i_end(ij),                &
 803                           head_grid%j_start(ij), head_grid%j_end(ij)                 )
 804    END DO
 805    !$OMP END PARALLEL DO
 806 #ifdef DM_PARALLEL
 807    val_a =  wrf_dm_sum_real ( val_a )
 808 #endif
 809
 810    WRITE(wrf_err_message, FMT='(A,E23.14)') 'ad_check: VAL_TL: ', val_l
 811    CALL wrf_message(TRIM(wrf_err_message))
 812    WRITE(wrf_err_message, FMT='(A,E23.14)') 'ad_check: VAL_AD: ', val_a
 813    CALL wrf_message(TRIM(wrf_err_message))
 814
 815    END IF !ADJ test
 816
 817    ! Deallocate all working space
 818    CALL deallocate_grid()
 819
 820 ENDIF
 821
 822    ! Release linked list for trajectory
 823    call xtraj_io_initialize
 824
 825    ! WRF model clean-up.  This calls MPI_FINALIZE() for DM parallel runs.
 826    CALL wrf_finalize
 827    STOP
 828    END SUBROUTINE wrf_adtl_check_sum
 829
 830    SUBROUTINE wrf_adtl_check_spot( )
 831 !<DESCRIPTION>
 832 !     WRF adjoint and tangent linear code check routine.
 833 !</DESCRIPTION>
 834
 835 !<DESCRIPTION>
 836 ! Once the top-level domain has been allocated, configured, and
 837 ! initialized, the model time integration is ready to proceed.
 838 !
 839 !</DESCRIPTION>
 840
 841    REAL :: alpha_m, val_l, val_a, save_l, val_n, coef, nl_num, nl_den
 842    REAL, DIMENSION(:,:), ALLOCATABLE :: factors
 843    REAL, DIMENSION(:,:,:), ALLOCATABLE :: ad_derivative, tl_derivative
 844    REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: nl_derivative
 845    INTEGER, DIMENSION(:,:), ALLOCATABLE :: locations_f, locations_i, scenario_matrix
 846    INTEGER :: ni, nf, nsc, nvar, ij, time_step, rc, &
 847               ninverse, nforward, firatio, iter, check_type, psign, &
 848               iloc_f, jloc_f, kloc_f, iloc_i, jloc_i, kloc_i, &
 849               vars_count, nnumer, ndenom
 850    CHARACTER*10, DIMENSION(:), ALLOCATABLE :: check_name !large for expectation of many future chem variables
 851    CHARACTER*256 :: timestr
 852
 853    ! Return immediately if not dyn_em_check
 854    IF ( config_flags%dyn_opt .NE. dyn_em_check ) RETURN
 855
 856    ! Force to turn off history output in this case
 857    CALL WRFU_AlarmRingerOff( head_grid%alarms( HISTORY_ALARM ), rc=rc )
 858
 859    ! Force to read the lateral boundary condition here.
 860    CALL med_before_solve_io ( head_grid, config_flags )
 861
 862    ! Close boundary file, as it will be read again from start later
 863    CALL close_dataset ( head_grid%lbc_fid , config_flags , &
 864                                             "DATASET=BOUNDARY" )
 865
 866    CALL init_domain_size ( head_grid, config_flags )
 867
 868    ! Release linked list for trajectory
 869    call xtraj_io_initialize
 870
 871 IF ( config_flags%check_TL .or. config_flags%check_AD ) THEN
 872
 873    ! Save the initial condition and boundary condition, x
 874    CALL allocate_grid ( )
 875
 876    !$OMP PARALLEL DO    &
 877    !$OMP DEFAULT (SHARED) PRIVATE ( ij )
 878    DO ij = 1 , head_grid%num_tiles
 879       CALL copy_grid_to_s ( head_grid , &
 880                             head_grid%i_start(ij), head_grid%i_end(ij),                &
 881                             head_grid%j_start(ij), head_grid%j_end(ij)                 )
 882    ENDDO
 883    !$OMP END PARALLEL DO
 884
 885    CALL       wrf_message ( "wrf: calling nonlinear integrate" )
 886    model_config_rec%dyn_opt = dyn_em
 887
 888    ! Set up basic states output
 889    CALL nl_get_time_step ( head_grid%id, time_step )
 890    CALL nl_set_auxhist6_interval_s ( head_grid%id, time_step )
 891    CALL nl_set_io_form_auxhist6 ( head_grid%id, 2 )
 892    CALL nl_set_frames_per_auxhist6 ( head_grid%id, 1 )
 893    IF ( ASSOCIATED( head_grid%domain_clock ) ) THEN
 894      IF ( head_grid%domain_clock_created ) THEN
 895        CALL WRFU_ClockDestroy( head_grid%domain_clock )
 896        head_grid%domain_clock_created = .FALSE.
 897      ENDIF
 898    ENDIF
 899    IF ( ASSOCIATED( head_grid%alarms ) .AND. &
 900         ASSOCIATED( head_grid%alarms_created ) ) THEN
 901      DO alarmid = 1, MAX_WRF_ALARMS
 902        IF ( head_grid%alarms_created( alarmid ) ) THEN
 903          CALL WRFU_AlarmDestroy( head_grid%alarms( alarmid ) )
 904          head_grid%alarms_created( alarmid ) = .FALSE.
 905        ENDIF
 906      ENDDO
 907    ENDIF
 908    CALL Setup_Timekeeping ( head_grid )
 909
 910    ! Force to turn off history output in this case
 911    CALL WRFU_AlarmRingerOff( head_grid%alarms( HISTORY_ALARM ), rc=rc )
 912    IF ( config_flags%mp_physics .NE. 0 .and. config_flags%mp_physics .NE. 99 .and. &
 913         config_flags%mp_physics .NE. 98 ) CALL nl_set_mp_physics (head_grid%id, config_flags%mp_physics_ad)
 914    IF ( config_flags%bl_pbl_physics .GT. 0 ) CALL nl_set_bl_pbl_physics (head_grid%id, 98)
 915    IF ( config_flags%cu_physics .GT. 0 ) THEN
 916       CALL nl_set_cu_physics (head_grid%id, 98)
 917       head_grid%cudt = 0
 918    ENDIF
 919    CALL nl_set_ra_lw_physics (head_grid%id, 0)
 920    CALL nl_set_ra_sw_physics (head_grid%id, 0)
 921    CALL nl_set_sf_sfclay_physics (head_grid%id, 0)
 922
 923    ! Force to turn off boundary input as we can only perturb the initial boundary condition.
 924    CALL nl_set_io_form_boundary( head_grid%id, 0 )
 925
 926    CALL wrf_run
 927
 928    ! Turn off basic states output
 929    CALL nl_set_io_form_auxhist6 ( head_grid%id, 0 )
 930    CALL nl_set_auxhist6_interval_s ( head_grid%id, 0 )
 931    if ( .not. head_grid%trajectory_io ) CALL nl_set_inputout_interval_s ( head_grid%id, 0 )
 932    IF ( ASSOCIATED( head_grid%domain_clock ) ) THEN
 933      IF ( head_grid%domain_clock_created ) THEN
 934        CALL WRFU_ClockDestroy( head_grid%domain_clock )
 935        head_grid%domain_clock_created = .FALSE.
 936      ENDIF
 937    ENDIF
 938    IF ( ASSOCIATED( head_grid%alarms ) .AND. &
 939         ASSOCIATED( head_grid%alarms_created ) ) THEN
 940      DO alarmid = 1, MAX_WRF_ALARMS
 941        IF ( head_grid%alarms_created( alarmid ) ) THEN
 942          CALL WRFU_AlarmDestroy( head_grid%alarms( alarmid ) )
 943          head_grid%alarms_created( alarmid ) = .FALSE.
 944        ENDIF
 945      ENDDO
 946    ENDIF
 947    CALL Setup_Timekeeping ( head_grid )
 948
 949    CALL       wrf_message ( "wrf: back from nonlinear integrate" )
 950
 951    CALL copy_grid_to_b ( head_grid )
 952
 953    wrf_err_message = '============== Adjoint - Tangent Linear Check ==============='
 954    CALL wrf_message(TRIM(wrf_err_message))
 955
 956    ! Get the locations of adjoint and tangent linear forcings
 957    CALL allocate_locations( 'locations_f', locations_f, ninverse, ierr)
 958    IF( ierr .GT. 0) THEN
 959       CALL wrf_error_fatal( 'adtl_check: error opening locations_f for reading' )
 960    ENDIF
 961    CALL allocate_locations( 'locations_i', locations_i, nforward, ierr)
 962    IF( ierr .GT. 0) THEN
 963       CALL wrf_error_fatal(&
 964                    'adtl_check: error opening locations_i for reading' )
 965    ENDIF
 966
 967    firatio = nforward/ninverse
 968    ierr = 0
 969    CALL get_forc_locations( 'locations_f', locations_f, ninverse, ierr)
 970    CALL get_forc_locations( 'locations_i', locations_i, nforward, ierr)
 971
 972
 973
 974    vars_count = 6
 975    DO nsc = 1, num_moist
 976       vars_count = vars_count + 1
 977    ENDDO
 978    DO nsc = 1, num_tracer
 979       vars_count = vars_count + 1
 980    ENDDO
 981    ALLOCATE(factors(vars_count,2))
 982    ALLOCATE(scenario_matrix(vars_count,vars_count))
 983    scenario_matrix = 0
 984    CALL gen_scenario_matrix(config_flags, scenario_matrix, vars_count, &
 985                             model_config_rec%numer_vars, model_config_rec%denom_vars)
 986
 987    ALLOCATE(check_name(vars_count))
 988    WRITE(check_name(1),FMT='(A)') 'U'
 989    WRITE(check_name(2),FMT='(A)') 'V'
 990    WRITE(check_name(3),FMT='(A)') 'W'
 991    WRITE(check_name(4),FMT='(A)') 'T'
 992    WRITE(check_name(5),FMT='(A)') 'PH'
 993    WRITE(check_name(6),FMT='(A)') 'MU'
 994    DO nsc = 1, num_moist
 995       IF( nsc .LT. PARAM_FIRST_SCALAR ) THEN
 996          WRITE(check_name(6+nsc),FMT='(A)')'DUMMY'
 997       ELSE
 998          WRITE(check_name(6+nsc),FMT='(A,I2.1)')'MOIST',nsc-PARAM_FIRST_SCALAR+1
 999       ENDIF
1000    ENDDO
1001    DO nsc = 1, num_tracer
1002       IF( nsc .LT. PARAM_FIRST_SCALAR ) THEN
1003          WRITE(check_name(6+nsc+num_moist),FMT='(A)')'DUMMY'
1004       ELSE
1005          WRITE(check_name(6+nsc+num_moist),FMT='(A,I2.1)')'TRACER',nsc-PARAM_FIRST_SCALAR+1
1006       ENDIF
1007    ENDDO
1008
1009    CALL wrf_message("AVAILABLE CONTROL VARIABLES FOR VALIDATION")
1010    DO nvar = 1, vars_count
1011       WRITE(wrf_err_message, FMT='(I2,A,A)') nvar,'  ',check_name(nvar)
1012       CALL wrf_message(TRIM(wrf_err_message))
1013    ENDDO
1014
1015    CALL wrf_message("")
1016    CALL wrf_message("SELECTED CONTROL VARIABLES")
1017    DO nnumer = 0, vars_count
1018       WRITE(wrf_err_message, FMT='(I2,A)') nnumer, ' ,  '
1019       DO ndenom = 1, vars_count
1020          IF ( nnumer .EQ. 0 ) THEN
1021             WRITE(wrf_err_message, FMT='(A,I2,A)') TRIM(wrf_err_message),&
1022                                 ndenom, ' ,  '
1023          ELSE
1024             WRITE(wrf_err_message, FMT='(A,I2,A)') TRIM(wrf_err_message),&
1025                                 scenario_matrix(nnumer,ndenom), ' ,  '
1026          ENDIF
1027       ENDDO
1028       CALL wrf_message(TRIM(wrf_err_message))
1029    ENDDO
1030
1031    ALLOCATE(ad_derivative(nforward,vars_count,vars_count))
1032    ALLOCATE(tl_derivative(nforward,vars_count,vars_count))
1033    ALLOCATE(nl_derivative(nforward,3,vars_count,vars_count))
1034
1035    ad_derivative = 0.0
1036    tl_derivative = 0.0
1037    nl_derivative = 0.0
1038
1039    !Find all the adjoint sensitivities
1040    numer_vars_Loop1: DO nnumer = 1, vars_count
1041
1042       IF( ALL(scenario_matrix(nnumer,:) .EQ. 0) ) CYCLE numer_vars_Loop1
1043       factors(:,1) = 0.0
1044       factors(nnumer,1) = 1.0
1045
1046
1047       iter = 0
1048
1049       DO nf = 1, ninverse
1050          iloc_f = locations_f(nf,1)
1051          jloc_f = locations_f(nf,2)
1052          kloc_f = locations_f(nf,3)
1053
1054          CALL zero_out_ad( head_grid )
1055          !$OMP PARALLEL DO    &
1056          !$OMP PRIVATE ( ij )
1057          adforce: DO ij = 1 , head_grid%num_tiles
1058             IF( head_grid%i_start(ij) .le. iloc_f .AND. &
1059                 head_grid%i_end(ij)  .ge. iloc_f .AND. &
1060                 head_grid%j_start(ij) .le. jloc_f .AND. &
1061                 head_grid%j_end(ij)  .ge. jloc_f ) THEN
1062                 CALL spot_force_ad( head_grid, iloc_f, jloc_f , kloc_f, factors(:,1), vars_count )
1063                 exit adforce
1064             ENDIF
1065          ENDDO adforce
1066          !$OMP END PARALLEL DO
1067
1068          ! Set the file names and interval for reading basic states.
1069          model_config_rec%auxinput6_inname = &
1070                                           "auxhist6_d<domain>_<date>"
1071          CALL nl_get_time_step ( head_grid%id, time_step )
1072          CALL nl_set_auxinput6_interval_s (head_grid%id, time_step )
1073          CALL nl_set_io_form_auxinput6 ( head_grid%id, 2 )
1074          CALL nl_set_frames_per_auxinput6 ( head_grid%id, 1 )
1075
1076          CALL wrf_run_ad
1077
1078
1079          ! Turn off auxinput6 reading
1080          CALL nl_set_auxinput6_interval_s (head_grid%id, 0 )
1081          IF ( ASSOCIATED( head_grid%domain_clock ) ) THEN
1082            IF ( head_grid%domain_clock_created ) THEN
1083              CALL WRFU_ClockDestroy( head_grid%domain_clock )
1084              head_grid%domain_clock_created = .FALSE.
1085            ENDIF
1086          ENDIF
1087          IF ( ASSOCIATED( head_grid%alarms ) .AND. &
1088               ASSOCIATED( head_grid%alarms_created ) ) THEN
1089            DO alarmid = 1, MAX_WRF_ALARMS
1090              IF ( head_grid%alarms_created( alarmid ) ) THEN
1091                CALL WRFU_AlarmDestroy( head_grid%alarms( alarmid ) )
1092                head_grid%alarms_created( alarmid ) = .FALSE.
1093              ENDIF
1094            ENDDO
1095          ENDIF
1096          CALL Setup_Timekeeping ( head_grid )
1097
1098          DO ni = 1,firatio
1099             iter = iter + 1
1100             iloc_i = locations_i(iter,1)
1101             jloc_i = locations_i(iter,2)
1102             kloc_i = locations_i(iter,3)
1103
1104             denom_vars_Loop1: DO ndenom = 1,vars_count
1105                IF ( scenario_matrix(nnumer,ndenom) .EQ. 0) CYCLE denom_vars_Loop1
1106                factors(:,2)=0.0
1107                factors(ndenom,2) = 1.0
1108                WRITE(wrf_err_message, FMT='(2(A,I3),5(A))') 'nf = ',nf, &
1109                             ',  ni = ',ni, ',  check_type = d[', &
1110                             TRIM(check_name(nnumer)),']/d[',TRIM(check_name(ndenom)),']'
1111                CALL wrf_message(TRIM(wrf_err_message))
1112
1113                CALL wrf_message("Extracting the AD sensitivity")
1114                !Store the AD derivative
1115                !$OMP PARALLEL DO    &
1116                !$OMP PRIVATE ( ij ) &
1117                !$OMP REDUCTION (+: ad_derivative(iter,nnumer,ndenom))
1118                adextract: DO ij = 1 , head_grid%num_tiles
1119                   IF( head_grid%i_start(ij) .le. iloc_i .AND. &
1120                      head_grid%i_end(ij)  .ge. iloc_i .AND. &
1121                      head_grid%j_start(ij) .le. jloc_i .AND. &
1122                      head_grid%j_end(ij)  .ge. jloc_i ) THEN
1123                      CALL extract_ad_der( head_grid, ad_derivative(iter,nnumer,ndenom), iloc_i, jloc_i, kloc_i, factors (:,2), vars_count)
1124                      exit adextract
1125                   ENDIF
1126                END DO adextract
1127                !$OMP END PARALLEL DO
1128 #ifdef DM_PARALLEL
1129                ad_derivative(iter,nnumer,ndenom) = wrf_dm_sum_real ( ad_derivative(iter,nnumer,ndenom) )
1130 #endif
1131
1132             ENDDO denom_vars_Loop1
1133          ENDDO
1134       ENDDO
1135    ENDDO numer_vars_Loop1
1136
1137
1138    !Find all the nonlinear and TL sensitivities
1139    iter = 0
1140    DO nf = 1, ninverse
1141       iloc_f = locations_f(nf,1)
1142       jloc_f = locations_f(nf,2)
1143       kloc_f = locations_f(nf,3)
1144       DO ni = 1,firatio
1145          iter = iter + 1
1146          iloc_i = locations_i(iter,1)
1147          jloc_i = locations_i(iter,2)
1148          kloc_i = locations_i(iter,3)
1149
1150          denom_vars_Loop2: DO ndenom = 1,vars_count
1151             IF ( ALL(scenario_matrix(:,ndenom) .EQ. 0) ) CYCLE denom_vars_Loop2
1152             factors(:,2)=0.0
1153             factors(ndenom,2) = 1.0
1154
1155             ! Do Finite Difference test of nonlinear model
1156             IF ( head_grid%check_AD .AND. head_grid%check_TL ) THEN
1157                nonlinearLoop: DO psign=-1,1,2
1158                   CALL restore_grid ( head_grid )
1159
1160                   !$OMP PARALLEL DO    &
1161                   !$OMP PRIVATE ( ij )
1162                   nlforce: DO ij = 1 , head_grid%num_tiles
1163                      IF( head_grid%i_start(ij) .le. iloc_i .AND. &
1164                          head_grid%i_end(ij)  .ge. iloc_i .AND. &
1165                          head_grid%j_start(ij) .le. jloc_i .AND. &
1166                          head_grid%j_end(ij)  .ge. jloc_i ) THEN
1167                          CALL spot_force_nl( head_grid, iloc_i, jloc_i, kloc_i, &
1168                                              factors(:,2), vars_count, REAL(psign) * config_flags%nl_pert)
1169                          exit nlforce
1170                      ENDIF
1171                   ENDDO nlforce
1172                   !$OMP END PARALLEL DO
1173
1174                   CALL wrf_message ( "wrf: calling nonlinear integrate" )
1175                   model_config_rec%dyn_opt = dyn_em
1176                   CALL domain_clock_get( head_grid, start_timestr=timestr )
1177                   CALL domain_clock_set( head_grid, current_timestr=timestr )
1178
1179                   ! Force to turn off history output in this case
1180                   CALL WRFU_AlarmRingerOff( head_grid%alarms( HISTORY_ALARM ), rc=rc )
1181                   ! Force to turn off boundary input as it was perturbated in add_grid.
1182                   CALL nl_set_io_form_boundary( head_grid%id, 0 )
1183
1184                   head_grid%dtbc = 0
1185                   head_grid%itimestep = 0
1186                   CALL wrf_run
1187                   CALL wrf_message( "wrf: back from nonlinear integrate" )
1188                   numer_vars_Loop2: DO nnumer = 1, vars_count
1189                      nl_num = 0.0
1190                      nl_den = 0.0
1191
1192                      IF( scenario_matrix(nnumer,ndenom) .EQ. 0 ) CYCLE numer_vars_Loop2
1193                      factors(:,1) = 0.0
1194                      factors(nnumer,1) = 1.0
1195                      WRITE(wrf_err_message, FMT='(2(A,I3),5(A))')'nf = ',nf, &
1196                                   ',  ni = ',ni, ',  check_type = d[', &
1197                                   TRIM(check_name(nnumer)),']/d[',TRIM(check_name(ndenom)),']'
1198                      CALL wrf_message(TRIM(wrf_err_message))
1199
1200                      CALL wrf_message("Extracting the FD sensitivity")
1201                      !Store the NL derivative
1202                      !$OMP PARALLEL DO    &
1203                      !$OMP PRIVATE ( ij ) &
1204                      !$OMP REDUCTION (+: nl_num)
1205                      numextract: DO ij = 1 , head_grid%num_tiles
1206                         IF( head_grid%i_start(ij) .le. iloc_f .AND. &
1207                            head_grid%i_end(ij)  .ge. iloc_f .AND. &
1208                            head_grid%j_start(ij) .le. jloc_f .AND. &
1209                            head_grid%j_end(ij)  .ge. jloc_f ) THEN
1210                            CALL extract_nl_num( head_grid, nl_num, &
1211                                                 iloc_f, jloc_f, kloc_f, &
1212                                                 factors(:,1), vars_count)
1213                            exit numextract
1214                         ENDIF
1215                      END DO numextract
1216                      !$OMP END PARALLEL DO
1217 #ifdef DM_PARALLEL
1218                      nl_num = wrf_dm_sum_real ( nl_num )
1219 #endif
1220
1221                      !$OMP PARALLEL DO    &
1222                      !$OMP PRIVATE ( ij ) &
1223                      !$OMP REDUCTION (+: nl_den)
1224                      denextract: DO ij = 1 , head_grid%num_tiles
1225                         IF( head_grid%i_start(ij) .le. iloc_i .AND. &
1226                            head_grid%i_end(ij)  .ge. iloc_i .AND. &
1227                            head_grid%j_start(ij) .le. jloc_i .AND. &
1228                            head_grid%j_end(ij)  .ge. jloc_i ) THEN
1229                            CALL extract_nl_den( nl_den, &
1230                                                 iloc_i, jloc_i, kloc_i, &
1231                                                 factors(:,2), vars_count, REAL(psign)*config_flags%nl_pert)
1232                            exit denextract
1233                         ENDIF
1234                      END DO denextract
1235                      !$OMP END PARALLEL DO
1236 #ifdef DM_PARALLEL
1237                      nl_den = wrf_dm_sum_real ( nl_den )
1238 #endif
1239
1240                      nl_derivative(iter,2+psign,nnumer,ndenom) = nl_num / nl_den
1241                   ENDDO numer_vars_Loop2
1242
1243                ENDDO nonlinearLoop
1244                nl_derivative(iter,2,nnumer,ndenom) = (nl_derivative(iter,1,nnumer,ndenom) + nl_derivative(iter,3,nnumer,ndenom)) * 0.5
1245             ENDIF
1246
1247             IF ( ASSOCIATED( head_grid%domain_clock ) ) THEN
1248               IF ( head_grid%domain_clock_created ) THEN
1249                 CALL WRFU_ClockDestroy( head_grid%domain_clock )
1250                 head_grid%domain_clock_created = .FALSE.
1251               ENDIF
1252             ENDIF
1253             IF ( ASSOCIATED( head_grid%alarms ) .AND. &
1254                  ASSOCIATED( head_grid%alarms_created ) ) THEN
1255               DO alarmid = 1, MAX_WRF_ALARMS
1256                 IF ( head_grid%alarms_created( alarmid ) ) THEN
1257                   CALL WRFU_AlarmDestroy( head_grid%alarms( alarmid ) )
1258                   head_grid%alarms_created( alarmid ) = .FALSE.
1259                 ENDIF
1260               ENDDO
1261             ENDIF
1262             CALL Setup_Timekeeping ( head_grid )
1263
1264             !Set the single forcing value for each tangent linear run
1265             CALL zero_out_tl( head_grid )
1266
1267             !$OMP PARALLEL DO    &
1268             !$OMP PRIVATE ( ij )
1269             tlforce: DO ij = 1 , head_grid%num_tiles
1270                IF( head_grid%i_start(ij) .le. iloc_i .AND. &
1271                   head_grid%i_end(ij)  .ge. iloc_i .AND. &
1272                   head_grid%j_start(ij) .le. jloc_i .AND. &
1273                   head_grid%j_end(ij)  .ge. jloc_i ) THEN
1274                   CALL spot_force_tl( head_grid, iloc_i, jloc_i, kloc_i, factors(:,2), vars_count)
1275                   exit tlforce
1276                ENDIF
1277             ENDDO tlforce
1278             !$OMP END PARALLEL DO
1279
1280             ! Set up basic states reading
1281             model_config_rec%auxinput6_inname = &
1282                                        "auxhist6_d<domain>_<date>"
1283             CALL nl_get_time_step ( head_grid%id, time_step )
1284             CALL nl_set_auxinput6_interval_s (head_grid%id, time_step )
1285             CALL nl_set_io_form_auxinput6 ( head_grid%id, 2 )
1286             CALL nl_set_frames_per_auxinput6 ( head_grid%id, 1 )
1287
1288             CALL wrf_run_tl
1289
1290             ! Turn off auxinput6 reading
1291             CALL nl_set_auxinput6_interval_s (head_grid%id, 0 )
1292             IF ( ASSOCIATED( head_grid%domain_clock ) ) THEN
1293               IF ( head_grid%domain_clock_created ) THEN
1294                 CALL WRFU_ClockDestroy( head_grid%domain_clock )
1295                 head_grid%domain_clock_created = .FALSE.
1296               ENDIF
1297             ENDIF
1298             IF ( ASSOCIATED( head_grid%alarms ) .AND. &
1299                  ASSOCIATED( head_grid%alarms_created ) ) THEN
1300               DO alarmid = 1, MAX_WRF_ALARMS
1301                 IF ( head_grid%alarms_created( alarmid ) ) THEN
1302                   CALL WRFU_AlarmDestroy( head_grid%alarms( alarmid ) )
1303                   head_grid%alarms_created( alarmid ) = .FALSE.
1304                 ENDIF
1305               ENDDO
1306             ENDIF
1307             CALL Setup_Timekeeping ( head_grid )
1308             numer_vars_Loop3: DO nnumer = 1, vars_count
1309                IF( scenario_matrix(nnumer,ndenom) .EQ. 0 ) CYCLE numer_vars_Loop3
1310                factors(:,1) = 0.0
1311                factors(nnumer,1) = 1.0
1312                WRITE(wrf_err_message, FMT='(2(A,I3),5(A))') 'nf = ',nf, &
1313                             ',  ni = ',ni, ',  check_type = d[', &
1314                             TRIM(check_name(nnumer)),']/d[',TRIM(check_name(ndenom)),']'
1315                CALL wrf_message(TRIM(wrf_err_message))
1316
1317                CALL wrf_message("Extracting the TL sensitivity")
1318                !Store the TL derivative
1319                !$OMP PARALLEL DO    &
1320                !$OMP PRIVATE ( ij ) &
1321                !$OMP REDUCTION (+: tl_derivative(iter,nnumer,ndenom))
1322                tlextract: DO ij = 1 , head_grid%num_tiles
1323                   IF( head_grid%i_start(ij) .le. iloc_f .AND. &
1324                      head_grid%i_end(ij)  .ge. iloc_f .AND. &
1325                      head_grid%j_start(ij) .le. jloc_f .AND. &
1326                      head_grid%j_end(ij)  .ge. jloc_f ) THEN
1327                      CALL extract_tl_der( head_grid, tl_derivative(iter,nnumer,ndenom), iloc_f, jloc_f, kloc_f, factors(:,1), vars_count)
1328                      exit tlextract
1329                   ENDIF
1330                END DO tlextract
1331                !$OMP END PARALLEL DO
1332 #ifdef DM_PARALLEL
1333                tl_derivative(iter,nnumer,ndenom) = wrf_dm_sum_real ( tl_derivative(iter,nnumer,ndenom) )
1334 #endif
1335             ENDDO numer_vars_Loop3
1336          ENDDO denom_vars_Loop2
1337       ENDDO
1338    ENDDO
1339
1340    ! Print out the sensitivities
1341    CALL wrf_message(&
1342         "====================== Results =======================")
1343    numer_vars_Loop4: DO nnumer = 1, vars_count
1344       IF( ALL(scenario_matrix(nnumer,:) .EQ. 0) ) CYCLE numer_vars_Loop4
1345       iter = 0
1346       DO nf = 1, ninverse
1347          DO ni = 1, firatio
1348             iter = iter+1
1349             iloc_i = locations_i(iter,1)
1350             jloc_i = locations_i(iter,2)
1351             kloc_i = locations_i(iter,3)
1352             iloc_f = locations_f(nf,1)
1353             jloc_f = locations_f(nf,2)
1354             kloc_f = locations_f(nf,3)
1355
1356             WRITE(wrf_err_message, FMT='(A,6(1x,I8))')&
1357                          'REPORT:',&
1358                          iloc_f,&
1359                          jloc_f,&
1360                          kloc_f,&
1361                          iloc_i,&
1362                          jloc_i,&
1363                          kloc_i
1364             CALL wrf_message(TRIM(wrf_err_message))
1365          ENDDO
1366       ENDDO
1367
1368       DO ndenom = 1, vars_count
1369          IF( scenario_matrix(nnumer,ndenom) .EQ. 0) CYCLE
1370          WRITE(wrf_err_message, FMT='(5(A))') 'REPORT: ------check_type = d[',&
1371                            TRIM(check_name(nnumer)),']/d[',TRIM(check_name(ndenom)),']------'
1372          CALL wrf_message(TRIM(wrf_err_message))
1373
1374          IF ( head_grid%check_TL .AND. head_grid%check_AD ) THEN
1375             WRITE(wrf_err_message, FMT='(A,4(1x,A17))') 'REPORT:','TL','AD','NL_BD','NL_FD'
1376             CALL wrf_message(TRIM(wrf_err_message))
1377
1378             iter = 0
1379             DO nf = 1, ninverse
1380                DO ni = 1, firatio
1381                   iter = iter+1
1382
1383                   WRITE(wrf_err_message, FMT='(A,4(E18.9))')&
1384                                'REPORT:',&
1385                                tl_derivative(iter,nnumer,ndenom), &
1386                                ad_derivative(iter,nnumer,ndenom), &
1387                                nl_derivative(iter,1,nnumer,ndenom), &
1388                                nl_derivative(iter,3,nnumer,ndenom)
1389                   CALL wrf_message(TRIM(wrf_err_message))
1390                ENDDO
1391             ENDDO
1392          ELSE
1393             WRITE(wrf_err_message, FMT='(A,2(1x,A17))') 'REPORT:  ','TL','AD'
1394             CALL wrf_message(TRIM(wrf_err_message))
1395
1396             iter = 0
1397             DO nf = 1, ninverse
1398                DO ni = 1, firatio
1399                   iter = iter+1
1400                   WRITE(wrf_err_message, FMT='(A,2(E18.9))')&
1401                                'REPORT:',&
1402                                tl_derivative(iter,nnumer,ndenom), &
1403                                ad_derivative(iter,nnumer,ndenom)
1404                   CALL wrf_message(TRIM(wrf_err_message))
1405                ENDDO
1406             ENDDO
1407          ENDIF
1408       ENDDO
1409    ENDDO numer_vars_Loop4
1410
1411    DEALLOCATE(ad_derivative)
1412    DEALLOCATE(tl_derivative)
1413    DEALLOCATE(nl_derivative)
1414
1415    DEALLOCATE(check_name)
1416    DEALLOCATE(scenario_matrix)
1417    DEALLOCATE(factors)
1418    DEALLOCATE(locations_f)
1419    DEALLOCATE(locations_i)
1420
1421    ! Deallocate all working space
1422    CALL deallocate_grid()
1423
1424 ENDIF
1425
1426    ! Release linked list for trajectory
1427    call xtraj_io_initialize
1428
1429    ! WRF model clean-up.  This calls MPI_FINALIZE() for DM parallel runs.
1430    CALL wrf_finalize
1431    STOP
1432    END SUBROUTINE wrf_adtl_check_spot
1433
1434
1435    SUBROUTINE wrf_run_tl( )
1436 !<DESCRIPTION>
1437 !     WRF tangent linear run routine.
1438 !</DESCRIPTION>
1439
1440 !<DESCRIPTION>
1441 ! Once the top-level domain has been allocated, configured, and
1442 ! initialized, the model time integration is ready to proceed.  The start
1443 ! and stop times for the domain are set to the start and stop time of the
1444 ! model run, and then <a href=integrate.html>integrate</a> is called to
1445 ! advance the domain forward through that specified time interval.  On
1446 ! return, the simulation is completed.
1447 !
1448 !</DESCRIPTION>
1449
1450    CHARACTER*256 :: timestr
1451    INTEGER       :: rc
1452    INTEGER       :: io_auxh7
1453    CHARACTER (LEN=80)      :: bdyname
1454    INTEGER                 :: open_status
1455
1456    !  The forecast integration for the most coarse grid is now started.  The
1457    !  integration is from the first step (1) to the last step of the simulation.
1458
1459    CALL       wrf_message ( "wrf: calling tangent linear integrate" )
1460    model_config_rec%dyn_opt = dyn_em_tl
1461    IF ( ASSOCIATED( head_grid%domain_clock ) ) THEN
1462      IF ( head_grid%domain_clock_created ) THEN
1463        CALL WRFU_ClockDestroy( head_grid%domain_clock )
1464        head_grid%domain_clock_created = .FALSE.
1465      ENDIF
1466    ENDIF
1467    IF ( ASSOCIATED( head_grid%alarms ) .AND. &
1468         ASSOCIATED( head_grid%alarms_created ) ) THEN
1469      DO alarmid = 1, MAX_WRF_ALARMS
1470        IF ( head_grid%alarms_created( alarmid ) ) THEN
1471          CALL WRFU_AlarmDestroy( head_grid%alarms( alarmid ) )
1472          head_grid%alarms_created( alarmid ) = .FALSE.
1473        ENDIF
1474      ENDDO
1475    ENDIF
1476    CALL Setup_Timekeeping ( head_grid )
1477    CALL domain_clock_get( head_grid, start_timestr=timestr )
1478    CALL domain_clock_set( head_grid, current_timestr=timestr )
1479
1480    ! Force to turn off history output in this case
1481    CALL WRFU_AlarmRingerOff( head_grid%alarms( HISTORY_ALARM ), rc=rc )
1482
1483    IF ( model_config_rec%bl_pbl_physics(head_grid%id) .EQ. SURFDRAGSCHEME ) THEN
1484       head_grid%rublten = 0.0
1485       head_grid%rvblten = 0.0
1486       head_grid%rthblten = 0.0
1487       head_grid%rqvblten = 0.0
1488    ENDIF
1489    IF ( model_config_rec%mp_physics(head_grid%id) .EQ. LSCONDSCHEME ) THEN
1490       head_grid%h_diabatic = 0.0
1491       head_grid%qv_diabatic = 0.0
1492       head_grid%qc_diabatic = 0.0
1493       head_grid%rainnc = 0.0
1494       head_grid%rainncv = 0.0
1495    ENDIF
1496    IF ( model_config_rec%mp_physics(head_grid%id) .EQ. MKESSLERSCHEME ) THEN
1497       head_grid%h_diabatic = 0.0
1498       head_grid%qv_diabatic = 0.0
1499       head_grid%qc_diabatic = 0.0
1500       head_grid%rainnc = 0.0
1501       head_grid%rainncv = 0.0
1502    ENDIF
1503    IF ( model_config_rec%cu_physics(head_grid%id) .EQ. DUCUSCHEME ) THEN
1504       head_grid%rthcuten = 0.0
1505       head_grid%rqvcuten = 0.0
1506       head_grid%rainc = 0.0
1507       head_grid%raincv = 0.0
1508       head_grid%pratec = 0.0
1509    ENDIF
1510
1511    IF ( .NOT. config_flags%trajectory_io ) THEN
1512       CALL nl_get_io_form_auxhist7( head_grid%id, io_auxh7 )
1513       CALL nl_set_io_form_auxhist7( head_grid%id, 0 )
1514    ENDIF
1515
1516    head_grid%itimestep = 0
1517    CALL start_domain ( head_grid , .TRUE. )
1518
1519    IF ( ASSOCIATED(xtraj_tail) ) xtraj_pointer => xtraj_tail
1520    CALL integrate ( head_grid )
1521
1522    IF ( .NOT. config_flags%trajectory_io ) THEN
1523       CALL nl_set_io_form_auxhist7( head_grid%id, io_auxh7 )
1524    ENDIF
1525
1526    ! Close boundary file, as it will be read again from start
1527    CALL construct_filename2a ( bdyname , config_flags%bdy_inname , head_grid%id , 2 , 'dummydate' )
1528    CALL wrf_inquire_opened(head_grid%lbc_fid , TRIM(bdyname) , open_status , ierr )
1529    IF ( open_status .EQ. WRF_FILE_OPENED_FOR_READ ) THEN
1530       CALL close_dataset ( head_grid%lbc_fid , config_flags , "DATASET=BOUNDARY" )
1531    ENDIF
1532
1533    CALL       wrf_message ( "wrf: back from tangent linear integrate" )
1534
1535    END SUBROUTINE wrf_run_tl
1536
1537    SUBROUTINE wrf_run_tl_standalone( )
1538 !<DESCRIPTION>
1539 !     WRF tangent linear code standalone run
1540 !</DESCRIPTION>
1541
1542 !</DESCRIPTION>
1543
1544    INTEGER :: rc, time_step, id, ierr
1545    CHARACTER*256 :: timestr
1546
1547    ! Return immediately if not dyn_em_tl
1548    IF ( config_flags%dyn_opt .NE. dyn_em_tl ) RETURN
1549
1550    ! Force to turn off history output in this case
1551    CALL WRFU_AlarmRingerOff( head_grid%alarms( HISTORY_ALARM ), rc=rc )
1552
1553    IF ( head_grid%trajectory_io ) THEN
1554
1555    ! Force to read the lateral boundary condition here.
1556    !CALL med_before_solve_io ( head_grid, config_flags )
1557
1558    ! Close boundary file, as it will be read again from start later
1559    !CALL close_dataset ( head_grid%lbc_fid , config_flags , "DATASET=BOUNDARY" )
1560
1561    CALL init_domain_size ( head_grid, config_flags )
1562
1563    ! Release linked list for trajectory
1564    call xtraj_io_initialize
1565
1566    ! Initialize linked list for adjoint forcing and tl. pertbation
1567    call adtl_initialize
1568
1569    CALL       wrf_message ( "wrf: calling nonlinear integrate" )
1570
1571    ! Set up basic states output
1572    CALL nl_get_time_step ( head_grid%id, time_step )
1573    CALL nl_set_auxhist6_interval_s ( head_grid%id, time_step )
1574    CALL nl_set_io_form_auxhist6 ( head_grid%id, 2 )
1575    CALL nl_set_frames_per_auxhist6 ( head_grid%id, 1 )
1576    IF ( ASSOCIATED( head_grid%domain_clock ) ) THEN
1577      IF ( head_grid%domain_clock_created ) THEN
1578        CALL WRFU_ClockDestroy( head_grid%domain_clock )
1579        head_grid%domain_clock_created = .FALSE.
1580      ENDIF
1581    ENDIF
1582    IF ( ASSOCIATED( head_grid%alarms ) .AND. &
1583         ASSOCIATED( head_grid%alarms_created ) ) THEN
1584      DO alarmid = 1, MAX_WRF_ALARMS
1585        IF ( head_grid%alarms_created( alarmid ) ) THEN
1586          CALL WRFU_AlarmDestroy( head_grid%alarms( alarmid ) )
1587          head_grid%alarms_created( alarmid ) = .FALSE.
1588        ENDIF
1589      ENDDO
1590    ENDIF
1591    CALL Setup_Timekeeping ( head_grid )
1592
1593    CALL wrf_run
1594
1595    ! Turn off basic states output
1596    CALL nl_set_io_form_auxhist6 ( head_grid%id, 0 )
1597    CALL nl_set_auxhist6_interval_s ( head_grid%id, 0 )
1598    IF ( ASSOCIATED( head_grid%domain_clock ) ) THEN
1599      IF ( head_grid%domain_clock_created ) THEN
1600        CALL WRFU_ClockDestroy( head_grid%domain_clock )
1601        head_grid%domain_clock_created = .FALSE.
1602      ENDIF
1603    ENDIF
1604    IF ( ASSOCIATED( head_grid%alarms ) .AND. &
1605         ASSOCIATED( head_grid%alarms_created ) ) THEN
1606      DO alarmid = 1, MAX_WRF_ALARMS
1607        IF ( head_grid%alarms_created( alarmid ) ) THEN
1608          CALL WRFU_AlarmDestroy( head_grid%alarms( alarmid ) )
1609          head_grid%alarms_created( alarmid ) = .FALSE.
1610        ENDIF
1611      ENDDO
1612    ENDIF
1613    CALL Setup_Timekeeping ( head_grid )
1614
1615    CALL       wrf_message ( "wrf: back from nonlinear integrate" )
1616
1617    ENDIF
1618
1619    ! Set the file names and interval for reading basic states.
1620    model_config_rec%auxinput6_inname = "auxhist6_d<domain>_<date>"
1621    CALL nl_get_time_step ( head_grid%id, time_step )
1622    CALL nl_set_auxinput6_interval_s (head_grid%id, time_step )
1623    CALL nl_set_io_form_auxinput6 ( head_grid%id, 2 )
1624    CALL nl_set_frames_per_auxinput6 ( head_grid%id, 1 )
1625
1626    IF ( config_flags%mp_physics .NE. 0 .and. config_flags%mp_physics .NE. 99 .and. &
1627         config_flags%mp_physics .NE. 98 ) CALL nl_set_mp_physics (head_grid%id, config_flags%mp_physics_ad)
1628    IF ( config_flags%bl_pbl_physics .GT. 0 ) CALL nl_set_bl_pbl_physics (head_grid%id, 98)
1629    IF ( config_flags%cu_physics .GT. 0 ) THEN
1630       CALL nl_set_cu_physics (head_grid%id, 98)
1631       head_grid%cudt = 0
1632    ENDIF
1633    CALL nl_set_ra_lw_physics (head_grid%id, 0)
1634    CALL nl_set_ra_sw_physics (head_grid%id, 0)
1635    CALL nl_set_sf_sfclay_physics (head_grid%id, 0)
1636
1637    head_grid%g_u_1 = 0.0
1638    head_grid%g_v_1 = 0.0
1639    head_grid%g_w_1 = 0.0
1640    head_grid%g_t_1 = 0.0
1641    head_grid%g_ph_1 = 0.0
1642    head_grid%g_mu_1 = 0.0
1643
1644    head_grid%g_u_2 = 0.0
1645    head_grid%g_v_2 = 0.0
1646    head_grid%g_w_2 = 0.0
1647    head_grid%g_t_2 = 0.0
1648    head_grid%g_ph_2 = 0.0
1649    head_grid%g_mu_2 = 0.0
1650
1651    head_grid%g_p = 0.0
1652
1653    head_grid%g_moist = 0.0
1654    head_grid%g_tracer = 0.0
1655    head_grid%g_scalar = 0.0
1656    head_grid%g_rainnc  = 0.0
1657    head_grid%g_rainncv = 0.0
1658    head_grid%g_rainc  = 0.0
1659    head_grid%g_raincv = 0.0
1660
1661    CALL domain_clock_get( head_grid, start_timestr=timestr )
1662    CALL domain_clock_set( head_grid, current_timestr=timestr )
1663    config_flags%auxinput9_inname = "init_pert_d01"
1664    config_flags%io_form_auxinput9 = 2
1665    CALL nl_set_io_form_auxinput9 ( head_grid%id, 2 )
1666    config_flags%frames_per_auxinput9 = 1
1667    CALL med_auxinput_in ( head_grid, auxinput9_alarm, config_flags )
1668    CALL  wrf_debug ( 0 , 'Read in initial perturbation' )
1669    config_flags%io_form_auxinput9 = 0
1670
1671    CALL wrf_run_tl
1672
1673    ! Release linked list for trajectory
1674    call xtraj_io_initialize
1675
1676    END SUBROUTINE wrf_run_tl_standalone
1677
1678    SUBROUTINE wrf_run_ad( )
1679 !<DESCRIPTION>
1680 !     WRF adjoint run routine.
1681 !</DESCRIPTION>
1682
1683 !<DESCRIPTION>
1684 ! Once the top-level domain has been allocated, configured, and
1685 ! initialized, the model time integration is ready to proceed.  The start
1686 ! and stop times for the domain are set to the start and stop time of the
1687 ! model run, and then <a href=integrate.html>integrate</a> is called to
1688 ! advance the domain forward through that specified time interval.  On
1689 ! return, the simulation is completed.
1690 !
1691 !</DESCRIPTION>
1692
1693    CHARACTER*256 :: timestr, timestr1
1694    INTEGER :: start_year,start_month,start_day,start_hour,start_minute,start_second
1695    INTEGER :: end_year,end_month,end_day,end_hour,end_minute,end_second
1696    INTEGER :: rc, runad
1697    REAL    :: timestep
1698    TYPE(WRFU_TimeInterval) :: run_interval
1699    INTEGER       :: io_auxh8
1700    CHARACTER (LEN=80)      :: bdyname
1701    INTEGER                 :: open_status
1702
1703    !  The forecast integration for the most coarse grid is now started.  The
1704    !  integration is from the first step (1) to the last step of the simulation.
1705
1706    CALL       wrf_message ( "wrf: calling adjoint integrate" )
1707
1708    ! Seeting for AD model
1709    model_config_rec%dyn_opt = dyn_em_ad
1710
1711    model_config_rec%run_days = -1 * model_config_rec%run_days
1712    model_config_rec%run_hours = -1 * model_config_rec%run_hours
1713    model_config_rec%run_minutes = -1 * model_config_rec%run_minutes
1714    model_config_rec%run_seconds = -1 * model_config_rec%run_seconds
1715
1716    start_year = model_config_rec%start_year(head_grid%id)
1717    start_month = model_config_rec%start_month(head_grid%id)
1718    start_day = model_config_rec%start_day(head_grid%id)
1719    start_hour = model_config_rec%start_hour(head_grid%id)
1720    start_minute = model_config_rec%start_minute(head_grid%id)
1721    start_second = model_config_rec%start_second(head_grid%id)
1722
1723    end_year = model_config_rec%end_year(head_grid%id)
1724    end_month = model_config_rec%end_month(head_grid%id)
1725    end_day = model_config_rec%end_day(head_grid%id)
1726    end_hour = model_config_rec%end_hour(head_grid%id)
1727    end_minute = model_config_rec%end_minute(head_grid%id)
1728    end_second = model_config_rec%end_second(head_grid%id)
1729
1730    model_config_rec%start_year(head_grid%id)  = end_year
1731    model_config_rec%start_month(head_grid%id)  = end_month
1732    model_config_rec%start_day(head_grid%id)  = end_day
1733    model_config_rec%start_hour(head_grid%id)  = end_hour
1734    model_config_rec%start_minute(head_grid%id)  = end_minute
1735    model_config_rec%start_second(head_grid%id)  = end_second
1736
1737    model_config_rec%end_year(head_grid%id)    = start_year
1738    model_config_rec%end_month(head_grid%id)    = start_month
1739    model_config_rec%end_day(head_grid%id)    = start_day
1740    model_config_rec%end_hour(head_grid%id)    = start_hour
1741    model_config_rec%end_minute(head_grid%id)    = start_minute
1742    model_config_rec%end_second(head_grid%id)    = start_second
1743
1744    IF ( ASSOCIATED( head_grid%domain_clock ) ) THEN
1745      IF ( head_grid%domain_clock_created ) THEN
1746        CALL WRFU_ClockDestroy( head_grid%domain_clock )
1747        head_grid%domain_clock_created = .FALSE.
1748      ENDIF
1749    ENDIF
1750    IF ( ASSOCIATED( head_grid%alarms ) .AND. &
1751         ASSOCIATED( head_grid%alarms_created ) ) THEN
1752      DO alarmid = 1, MAX_WRF_ALARMS
1753        IF ( head_grid%alarms_created( alarmid ) ) THEN
1754          CALL WRFU_AlarmDestroy( head_grid%alarms( alarmid ) )
1755          head_grid%alarms_created( alarmid ) = .FALSE.
1756        ENDIF
1757      ENDDO
1758    ENDIF
1759    CALL Setup_Timekeeping ( head_grid )
1760    head_grid%start_subtime = domain_get_start_time ( head_grid )
1761    head_grid%stop_subtime = domain_get_stop_time ( head_grid )
1762
1763    ! Force to turn off history output in this case
1764    CALL WRFU_AlarmRingerOff( head_grid%alarms( HISTORY_ALARM ), rc=rc )
1765
1766    CALL domain_clock_get( head_grid, start_timestr=timestr )
1767    CALL domain_clock_set( head_grid, current_timestr=timestr )
1768    CALL domain_clock_set( head_grid, time_step_seconds=-1*model_config_rec%time_step )
1769
1770    IF ( ASSOCIATED(xtraj_head) ) xtraj_pointer => xtraj_head%next
1771
1772    ! head_grid%itimestep should be the last
1773    IF ( head_grid%itimestep .EQ. 0 ) THEN
1774       timestep=abs(head_grid%dt)
1775       run_interval = head_grid%stop_subtime - head_grid%start_subtime
1776
1777       CALL WRFU_TimeIntervalGet( run_interval, S=runad, rc=rc )
1778       runad = abs(runad)
1779
1780       head_grid%itimestep = ceiling(real(runad)/timestep)
1781    ENDIF
1782
1783    IF ( .NOT. config_flags%trajectory_io ) THEN
1784       CALL nl_get_io_form_auxhist8( head_grid%id, io_auxh8 )
1785       CALL nl_set_io_form_auxhist8( head_grid%id, 0 )
1786    ENDIF
1787
1788    CALL integrate ( head_grid )
1789
1790    IF ( .NOT. config_flags%trajectory_io ) THEN
1791       CALL nl_set_io_form_auxhist8( head_grid%id, io_auxh8 )
1792    ENDIF
1793
1794    CALL start_domain ( head_grid , .TRUE. )
1795
1796    IF ( .NOT. config_flags%trajectory_io .OR. gradient_out ) THEN
1797       config_flags%auxhist7_outname = "gradient_wrfplus_d<domain>_<date>"
1798       config_flags%io_form_auxhist7 = 2
1799       CALL nl_set_io_form_auxhist7 ( head_grid%id, 2 )
1800       CALL  med_hist_out ( head_grid , AUXHIST7_ALARM , config_flags )
1801       CALL  wrf_debug ( 0 , 'Output gradient in WRFPLUS (not including LBC gradient)' )
1802    ENDIF
1803
1804 !  Recover to NL model
1805    model_config_rec%dyn_opt = dyn_em
1806
1807    model_config_rec%run_days = -1 * model_config_rec%run_days
1808    model_config_rec%run_hours = -1 * model_config_rec%run_hours
1809    model_config_rec%run_minutes = -1 * model_config_rec%run_minutes
1810    model_config_rec%run_seconds = -1 * model_config_rec%run_seconds
1811
1812    model_config_rec%start_year(head_grid%id)  = start_year
1813    model_config_rec%start_month(head_grid%id)  = start_month
1814    model_config_rec%start_day(head_grid%id)  = start_day
1815    model_config_rec%start_hour(head_grid%id)  = start_hour
1816    model_config_rec%start_minute(head_grid%id)  = start_minute
1817    model_config_rec%start_second(head_grid%id)  = start_second
1818
1819    model_config_rec%end_year(head_grid%id)    = end_year
1820    model_config_rec%end_month(head_grid%id)    = end_month
1821    model_config_rec%end_day(head_grid%id)    = end_day
1822    model_config_rec%end_hour(head_grid%id)    = end_hour
1823    model_config_rec%end_minute(head_grid%id)    = end_minute
1824    model_config_rec%end_second(head_grid%id)    = end_second
1825
1826    IF ( ASSOCIATED( head_grid%domain_clock ) ) THEN
1827      IF ( head_grid%domain_clock_created ) THEN
1828        CALL WRFU_ClockDestroy( head_grid%domain_clock )
1829        head_grid%domain_clock_created = .FALSE.
1830      ENDIF
1831    ENDIF
1832    IF ( ASSOCIATED( head_grid%alarms ) .AND. &
1833         ASSOCIATED( head_grid%alarms_created ) ) THEN
1834      DO alarmid = 1, MAX_WRF_ALARMS
1835        IF ( head_grid%alarms_created( alarmid ) ) THEN
1836          CALL WRFU_AlarmDestroy( head_grid%alarms( alarmid ) )
1837          head_grid%alarms_created( alarmid ) = .FALSE.
1838        ENDIF
1839      ENDDO
1840    ENDIF
1841    CALL Setup_Timekeeping ( head_grid )
1842    head_grid%start_subtime = domain_get_start_time ( head_grid )
1843    head_grid%stop_subtime = domain_get_stop_time ( head_grid )
1844    CALL domain_clock_set( head_grid, time_step_seconds=model_config_rec%time_step )
1845
1846    ! Close boundary file, as it will be read again from start
1847    CALL construct_filename2a ( bdyname , config_flags%bdy_inname , head_grid%id , 2 , 'dummydate' )
1848    CALL wrf_inquire_opened(head_grid%lbc_fid , TRIM(bdyname) , open_status , ierr )
1849    IF ( open_status .EQ. WRF_FILE_OPENED_FOR_READ ) THEN
1850       CALL close_dataset ( head_grid%lbc_fid , config_flags , "DATASET=BOUNDARY" )
1851    ENDIF
1852
1853    CALL       wrf_message ( "wrf: back from adjoint integrate" )
1854
1855    END SUBROUTINE wrf_run_ad
1856
1857    SUBROUTINE wrf_run_ad_standalone( )
1858 !<DESCRIPTION>
1859 !     WRF adjoint code standalone run
1860 !</DESCRIPTION>
1861
1862 !</DESCRIPTION>
1863
1864    INTEGER :: rc, time_step, id, ierr
1865    CHARACTER*256 :: timestr
1866    CHARACTER (LEN=80) :: bdyname
1867
1868    ! Return immediately if not dyn_em_ad
1869    IF ( config_flags%dyn_opt .NE. dyn_em_ad ) RETURN
1870
1871    ! Force to turn off history output in this case
1872    CALL WRFU_AlarmRingerOff( head_grid%alarms( HISTORY_ALARM ), rc=rc )
1873
1874    IF ( head_grid%trajectory_io ) THEN
1875
1876    ! Force to read the lateral boundary condition here.
1877    !CALL med_before_solve_io ( head_grid, config_flags )
1878
1879    ! Close boundary file, as it will be read again from start later
1880    !CALL close_dataset ( head_grid%lbc_fid , config_flags , "DATASET=BOUNDARY" )
1881
1882    CALL init_domain_size ( head_grid, config_flags )
1883
1884    ! Release linked list for trajectory
1885    call xtraj_io_initialize
1886
1887    CALL       wrf_message ( "wrf: calling nonlinear integrate" )
1888
1889    ! Set up basic states output
1890    CALL nl_get_time_step ( head_grid%id, time_step )
1891    CALL nl_set_auxhist6_interval_s ( head_grid%id, time_step )
1892    CALL nl_set_io_form_auxhist6 ( head_grid%id, 2 )
1893    CALL nl_set_frames_per_auxhist6 ( head_grid%id, 1 )
1894    IF ( ASSOCIATED( head_grid%domain_clock ) ) THEN
1895      IF ( head_grid%domain_clock_created ) THEN
1896        CALL WRFU_ClockDestroy( head_grid%domain_clock )
1897        head_grid%domain_clock_created = .FALSE.
1898      ENDIF
1899    ENDIF
1900    IF ( ASSOCIATED( head_grid%alarms ) .AND. &
1901         ASSOCIATED( head_grid%alarms_created ) ) THEN
1902      DO alarmid = 1, MAX_WRF_ALARMS
1903        IF ( head_grid%alarms_created( alarmid ) ) THEN
1904          CALL WRFU_AlarmDestroy( head_grid%alarms( alarmid ) )
1905          head_grid%alarms_created( alarmid ) = .FALSE.
1906        ENDIF
1907      ENDDO
1908    ENDIF
1909    CALL Setup_Timekeeping ( head_grid )
1910
1911    CALL wrf_run
1912
1913    ! Turn off basic states output
1914    CALL nl_set_io_form_auxhist6 ( head_grid%id, 0 )
1915    CALL nl_set_auxhist6_interval_s ( head_grid%id, 0 )
1916    IF ( ASSOCIATED( head_grid%domain_clock ) ) THEN
1917      IF ( head_grid%domain_clock_created ) THEN
1918        CALL WRFU_ClockDestroy( head_grid%domain_clock )
1919        head_grid%domain_clock_created = .FALSE.
1920      ENDIF
1921    ENDIF
1922    IF ( ASSOCIATED( head_grid%alarms ) .AND. &
1923         ASSOCIATED( head_grid%alarms_created ) ) THEN
1924      DO alarmid = 1, MAX_WRF_ALARMS
1925        IF ( head_grid%alarms_created( alarmid ) ) THEN
1926          CALL WRFU_AlarmDestroy( head_grid%alarms( alarmid ) )
1927          head_grid%alarms_created( alarmid ) = .FALSE.
1928        ENDIF
1929      ENDDO
1930    ENDIF
1931    CALL Setup_Timekeeping ( head_grid )
1932
1933    CALL       wrf_message ( "wrf: back from nonlinear integrate" )
1934
1935    ENDIF
1936
1937    ! Set the file names and interval for reading basic states.
1938    model_config_rec%auxinput6_inname = "auxhist6_d<domain>_<date>"
1939    CALL nl_get_time_step ( head_grid%id, time_step )
1940    CALL nl_set_auxinput6_interval_s (head_grid%id, time_step )
1941    CALL nl_set_io_form_auxinput6 ( head_grid%id, 2 )
1942    CALL nl_set_frames_per_auxinput6 ( head_grid%id, 1 )
1943
1944    IF ( config_flags%mp_physics .NE. 0 .and. config_flags%mp_physics .NE. 99 .and. &
1945         config_flags%mp_physics .NE. 98 ) CALL nl_set_mp_physics (head_grid%id, config_flags%mp_physics_ad)
1946    IF ( config_flags%bl_pbl_physics .GT. 0 ) CALL nl_set_bl_pbl_physics (head_grid%id, 98)
1947    IF ( config_flags%cu_physics .GT. 0 ) THEN
1948       CALL nl_set_cu_physics (head_grid%id, 98)
1949       head_grid%cudt = 0
1950    ENDIF
1951    CALL nl_set_ra_lw_physics (head_grid%id, 0)
1952    CALL nl_set_ra_sw_physics (head_grid%id, 0)
1953    CALL nl_set_sf_sfclay_physics (head_grid%id, 0)
1954
1955    CALL zero_out_ad ( head_grid )
1956
1957    head_grid%g_u_1 = 0.0
1958    head_grid%g_v_1 = 0.0
1959    head_grid%g_w_1 = 0.0
1960    head_grid%g_t_1 = 0.0
1961    head_grid%g_ph_1 = 0.0
1962    head_grid%g_mu_1 = 0.0
1963
1964    head_grid%g_u_2 = 0.0
1965    head_grid%g_v_2 = 0.0
1966    head_grid%g_w_2 = 0.0
1967    head_grid%g_t_2 = 0.0
1968    head_grid%g_ph_2 = 0.0
1969    head_grid%g_mu_2 = 0.0
1970
1971    head_grid%g_p = 0.0
1972
1973    head_grid%g_moist = 0.0
1974    head_grid%g_tracer = 0.0
1975    head_grid%g_scalar = 0.0
1976    head_grid%g_rainnc  = 0.0
1977    head_grid%g_rainncv = 0.0
1978    head_grid%g_rainc  = 0.0
1979    head_grid%g_raincv = 0.0
1980
1981    CALL domain_clock_get( head_grid, stop_timestr=timestr )
1982    CALL domain_clock_set( head_grid, current_timestr=timestr )
1983    config_flags%auxinput7_inname = "final_sens_d01"
1984    config_flags%io_form_auxinput7 = 2
1985    CALL nl_set_io_form_auxinput7 ( head_grid%id, 2 )
1986    config_flags%frames_per_auxinput7 = 1
1987    CALL med_auxinput_in ( head_grid, auxinput7_alarm, config_flags )
1988    CALL  wrf_debug ( 0 , 'Read in final sensitivuty' )
1989    config_flags%io_form_auxinput7 = 0
1990    CALL domain_clock_get( head_grid, start_timestr=timestr )
1991    CALL domain_clock_set( head_grid, current_timestr=timestr )
1992
1993    gradient_out = .TRUE.
1994
1995    CALL wrf_run_ad
1996
1997 !-- Output gradient in boundary
1998    CALL construct_filename1( bdyname , 'gradient_wrfbdy' , head_grid%id , 2 )
1999    CALL open_w_dataset ( id, TRIM(bdyname) , head_grid , config_flags , output_boundary , "DATASET=BOUNDARY", ierr )
2000    IF ( ierr .NE. 0 ) THEN
2001       CALL wrf_error_fatal( 'real: error opening wrfbdy for writing' )
2002    END IF
2003    print *,'Output LBC gradient valid between these times ',current_date, ' ',start_date
2004    CALL output_boundary ( id, head_grid , config_flags , ierr )
2005
2006    ! Release linked list for trajectory
2007    call xtraj_io_initialize
2008
2009    END SUBROUTINE wrf_run_ad_standalone
2010 !------------------AD/TL code end--------------------
2011 #endif
2012
2013    SUBROUTINE wrf_run( )
2014 !<DESCRIPTION>
2015 !     WRF run routine.
2016 !</DESCRIPTION>
2017
2018 #if (WRFPLUS == 1)
2019    integer :: io_auxh7, io_auxh8
2020    CHARACTER (LEN=80)      :: bdyname
2021    INTEGER                 :: open_status
2022 #endif
2023 !<DESCRIPTION>
2024 ! Once the top-level domain has been allocated, configured, and
2025 ! initialized, the model time integration is ready to proceed.  The start
2026 ! and stop times for the domain are set to the start and stop time of the
2027 ! model run, and then <a href=integrate.html>integrate</a> is called to
2028 ! advance the domain forward through that specified time interval.  On
2029 ! return, the simulation is completed.
2030 !
2031 !</DESCRIPTION>
2032
2033    !  The forecast integration for the most coarse grid is now started.  The
2034    !  integration is from the first step (1) to the last step of the simulation.
2035
2036    CALL       wrf_debug ( 100 , 'wrf: calling integrate' )
2037 #if (WRFPLUS == 1)
2038    model_config_rec%dyn_opt = dyn_em
2039    IF ( model_config_rec%bl_pbl_physics(head_grid%id) .EQ. SURFDRAGSCHEME ) THEN
2040       head_grid%rublten = 0.0
2041       head_grid%rvblten = 0.0
2042       head_grid%rthblten = 0.0
2043       head_grid%rqvblten = 0.0
2044    ENDIF
2045    IF ( model_config_rec%mp_physics(head_grid%id) .EQ. LSCONDSCHEME ) THEN
2046       head_grid%h_diabatic = 0.0
2047       head_grid%qv_diabatic = 0.0
2048       head_grid%qc_diabatic = 0.0
2049       head_grid%rainnc = 0.0
2050       head_grid%rainncv = 0.0
2051    ENDIF
2052    IF ( model_config_rec%mp_physics(head_grid%id) .EQ. MKESSLERSCHEME ) THEN
2053       head_grid%h_diabatic = 0.0
2054       head_grid%qv_diabatic = 0.0
2055       head_grid%qc_diabatic = 0.0
2056       head_grid%rainnc = 0.0
2057       head_grid%rainncv = 0.0
2058    ENDIF
2059    IF ( model_config_rec%cu_physics(head_grid%id) .EQ. DUCUSCHEME ) THEN
2060       head_grid%rthcuten = 0.0
2061       head_grid%rqvcuten = 0.0
2062       head_grid%rainc = 0.0
2063       head_grid%raincv = 0.0
2064    ENDIF
2065
2066    CALL nl_get_io_form_auxhist7( head_grid%id, io_auxh7 )
2067    CALL nl_get_io_form_auxhist8( head_grid%id, io_auxh8 )
2068    CALL nl_set_io_form_auxhist7( head_grid%id, 0 )
2069    CALL nl_set_io_form_auxhist8( head_grid%id, 0 )
2070
2071    head_grid%itimestep = 0
2072    CALL start_domain ( head_grid , .TRUE. )
2073 #endif
2074    CALL integrate ( head_grid )
2075
2076 #if (WRFPLUS == 1)
2077    CALL nl_set_io_form_auxhist7( head_grid%id, io_auxh7 )
2078    CALL nl_set_io_form_auxhist8( head_grid%id, io_auxh8 )
2079
2080    ! Close boundary file, as it will be read again from start
2081    CALL construct_filename2a ( bdyname , config_flags%bdy_inname , head_grid%id , 2 , 'dummydate' )
2082    CALL wrf_inquire_opened(head_grid%lbc_fid , TRIM(bdyname) , open_status , ierr )
2083    IF ( open_status .EQ. WRF_FILE_OPENED_FOR_READ ) THEN
2084       CALL close_dataset ( head_grid%lbc_fid , config_flags , "DATASET=BOUNDARY" )
2085    ENDIF
2086 #endif
2087
2088    CALL       wrf_debug ( 100 , 'wrf: back from integrate' )
2089
2090    END SUBROUTINE wrf_run
2091
2092
2093
2094    SUBROUTINE wrf_finalize( no_shutdown )
2095 !<DESCRIPTION>
2096 !     WRF finalize routine.
2097 !</DESCRIPTION>
2098
2099 !<DESCRIPTION>
2100 ! A Mediation Layer-provided
2101 ! subroutine, <a href=med_shutdown_io.html>med_shutdown_io</a> is called
2102 ! to allow the the model to do any I/O specific cleanup and shutdown, and
2103 ! then the WRF Driver Layer routine <a
2104 ! href=wrf_shutdown.html>wrf_shutdown</a> (quilt servers would be
2105 ! directed to shut down here) is called to properly end the run,
2106 ! including shutting down the communications (for example, most comm
2107 ! layers would call MPI_FINALIZE at this point if they're using MPI).
2108 !
2109 !</DESCRIPTION>
2110      LOGICAL, OPTIONAL, INTENT(IN) :: no_shutdown
2111
2112    ! shut down I/O
2113    CALL med_shutdown_io ( head_grid , config_flags )
2114    CALL       wrf_debug ( 100 , 'wrf: back from med_shutdown_io' )
2115
2116    CALL       wrf_debug (   0 , 'wrf: SUCCESS COMPLETE WRF' )
2117
2118    ! Call wrf_shutdown() (which calls MPI_FINALIZE()
2119    ! for DM parallel runs).
2120    IF ( .NOT. PRESENT( no_shutdown ) ) THEN
2121      ! Finalize time manager
2122       IF (coupler_on) THEN
2123          CALL cpl_finalize()
2124       ELSE
2125          CALL WRFU_Finalize
2126          CALL wrf_shutdown
2127       ENDIF
2128    ENDIF
2129
2130    END SUBROUTINE wrf_finalize
2131
2132
2133    SUBROUTINE wrf_dfi()
2134 !<DESCRIPTION>
2135 ! Runs a digital filter initialization procedure.
2136 !</DESCRIPTION>
2137       IMPLICIT NONE
2138
2139 ! #if (EM_CORE == 1)
2140       ! Initialization procedure
2141       IF ( config_flags%dfi_opt .NE. DFI_NODFI ) THEN
2142
2143          SELECT CASE ( config_flags%dfi_opt )
2144
2145             CASE (DFI_DFL)
2146                wrf_err_message = 'Initializing with DFL'
2147                CALL wrf_message(TRIM(wrf_err_message))
2148
2149                wrf_err_message = '   Filtering forward in time'
2150                CALL wrf_message(TRIM(wrf_err_message))
2151
2152                CALL wrf_dfi_fwd_init()
2153                CALL wrf_run()
2154
2155                CALL wrf_dfi_array_reset()
2156
2157                CALL wrf_dfi_fst_init()
2158
2159                IF ( config_flags%dfi_write_filtered_input ) THEN
2160                   CALL wrf_dfi_write_initialized_state()
2161                END IF
2162
2163             CASE (DFI_DDFI)
2164                wrf_err_message = 'Initializing with DDFI'
2165                CALL wrf_message(TRIM(wrf_err_message))
2166
2167                wrf_err_message = '   Integrating backward in time'
2168                CALL wrf_message(TRIM(wrf_err_message))
2169
2170                CALL wrf_dfi_bck_init()
2171                CALL wrf_run()
2172
2173                wrf_err_message = '   Filtering forward in time'
2174                CALL wrf_message(TRIM(wrf_err_message))
2175
2176                CALL wrf_dfi_fwd_init()
2177                CALL wrf_run()
2178
2179                CALL wrf_dfi_array_reset()
2180
2181                CALL wrf_dfi_fst_init()
2182
2183                IF ( config_flags%dfi_write_filtered_input ) THEN
2184                   CALL wrf_dfi_write_initialized_state()
2185                END IF
2186
2187             CASE (DFI_TDFI)
2188                wrf_err_message = 'Initializing with TDFI'
2189                CALL wrf_message(TRIM(wrf_err_message))
2190
2191                wrf_err_message = '   Integrating backward in time'
2192                CALL wrf_message(TRIM(wrf_err_message))
2193
2194                CALL wrf_dfi_bck_init()
2195                CALL wrf_run()
2196
2197                CALL wrf_dfi_array_reset()
2198
2199                wrf_err_message = '   Filtering forward in time'
2200                CALL wrf_message(TRIM(wrf_err_message))
2201
2202                CALL wrf_dfi_fwd_init()
2203                CALL wrf_run()
2204
2205                CALL wrf_dfi_array_reset()
2206
2207                CALL wrf_dfi_fst_init()
2208
2209                IF ( config_flags%dfi_write_filtered_input ) THEN
2210                   CALL wrf_dfi_write_initialized_state()
2211                END IF
2212
2213             CASE DEFAULT
2214                wrf_err_message = 'Unrecognized DFI_OPT in namelist'
2215                CALL wrf_error_fatal(TRIM(wrf_err_message))
2216
2217          END SELECT
2218
2219       END IF
2220 ! #endif
2221
2222    END SUBROUTINE wrf_dfi
2223
2224    SUBROUTINE set_derived_rconfigs
2225 !<DESCRIPTION>
2226 ! Some derived rconfig entries need to be set based on the value of other,
2227 ! non-derived entries before package-dependent memory allocation takes place.
2228 ! This might be employed when, for example, we want to allocate arrays in
2229 ! a package that depends on the setting of two or more namelist variables.
2230 ! In this subroutine, we do just that.
2231 !</DESCRIPTION>
2232
2233       IMPLICIT NONE
2234
2235       INTEGER :: i
2236
2237
2238 ! #if (EM_CORE == 1)
2239       IF ( model_config_rec % dfi_opt .EQ. DFI_NODFI ) THEN
2240         DO i = 1, model_config_rec % max_dom
2241            model_config_rec % mp_physics_dfi(i) = -1
2242         ENDDO
2243       ELSE
2244         DO i = 1, model_config_rec % max_dom
2245            model_config_rec % mp_physics_dfi(i) = model_config_rec % mp_physics(i)
2246         ENDDO
2247       END IF
2248 ! #endif
2249
2250 #if (EM_CORE == 1)
2251       IF ( model_config_rec % dfi_opt .EQ. DFI_NODFI ) THEN
2252         DO i = 1, model_config_rec % max_dom
2253            model_config_rec % bl_pbl_physics_dfi(i) = -1
2254         ENDDO
2255       ELSE
2256         DO i = 1, model_config_rec % max_dom
2257            model_config_rec % bl_pbl_physics_dfi(i) = model_config_rec % bl_pbl_physics(i)
2258         ENDDO
2259       END IF
2260 #endif
2261
2262 #if (DA_CORE == 1)
2263       IF ( model_config_rec % dyn_opt .EQ. 2 ) THEN
2264         DO i = 1, model_config_rec % max_dom
2265            model_config_rec % mp_physics_4dvar(i) = -1
2266         ENDDO
2267       ELSE
2268         DO i = 1, model_config_rec % max_dom
2269            model_config_rec % mp_physics_4dvar(i) = model_config_rec % mp_physics(i)
2270         ENDDO
2271       END IF
2272 #endif
2273
2274    END SUBROUTINE set_derived_rconfigs
2275
2276    RECURSIVE SUBROUTINE alloc_doms_for_dfi ( grid )
2277
2278       !  Input variables.
2279
2280       TYPE (domain) , pointer :: grid
2281
2282       !  Local variables.
2283
2284       TYPE (domain) , pointer :: new_nest_loc
2285       TYPE (grid_config_rec_type) :: parent_config_flags
2286       INTEGER :: nestid_loc , kid_loc
2287
2288          !  Are there any subdomains from this level.  The output is the nestid (the domain
2289          !  ID of the nest), and kid (an index to which of the parent's children this new nested
2290          !  domain represents).
2291
2292          DO WHILE ( nests_to_open( grid , nestid_loc , kid_loc ) )
2293
2294             !  If we found another child domain, we continue on: allocate, set up time keeping,
2295             !  initialize.
2296
2297             CALL alloc_and_configure_domain ( domain_id  = nestid_loc   , &
2298                                               grid       = new_nest_loc , &
2299                                               parent     = grid         , &
2300                                               kid        = kid_loc        )
2301
2302 print *,'for parent domain id #',grid%id,', found child domain #',nestid_loc
2303             !  Since this is a DFI run, set the DFI switches to the same for all domains.
2304
2305             new_nest_loc%dfi_opt = head_grid%dfi_opt
2306             new_nest_loc%dfi_stage = DFI_SETUP
2307
2308             !  Set up time keeping for the fine grid space that was just allocated.
2309
2310             CALL Setup_Timekeeping (new_nest_loc)
2311
2312             !  With space allocated, and timers set, the fine grid can be initialized with data.
2313
2314             CALL model_to_grid_config_rec ( grid%id , model_config_rec , parent_config_flags )
2315             CALL med_nest_initial ( grid , new_nest_loc , config_flags )
2316
2317             !  Here's the recursive part.  For each of these child domains, we call this same routine.
2318             !  This will find all of "new_nest_loc" first generation progeny.
2319
2320             CALL alloc_doms_for_dfi ( new_nest_loc )
2321
2322          END DO
2323
2324    END SUBROUTINE alloc_doms_for_dfi
2325
2326 END MODULE module_wrf_top