var/convertor/nmm_interface_convertor/gridmod.F

   1 !-------------------------------------------------------------------------
   2 !    NOAA/NCEP, National Centers for Environmental Prediction GSI        !
   3 !-------------------------------------------------------------------------
   4 !BOP
   5 !
   6 ! !MODULE:  gridmod --- GSI grid related variable declarations
   7 !
   8 ! !INTERFACE:
   9 !
  10 module gridmod
  11
  12 ! !USES:
  13
  14   use kinds, only: i_byte,r_kind,r_single,i_kind
  15   implicit none
  16
  17 ! !DESCRIPTION: module containing grid related variable declarations
  18 !
  19 ! !REVISION HISTORY:
  20 !   2003-09-25  kleist
  21 !   2003-xx-xx  parrish,wu  regional components added
  22 !   2004-05-13  kleist, documentation
  23 !   2004-07-15  todling, protex-compliant prologue
  24 !   2004-12-23  treadon - add routines get_ij and get_ijk
  25 !   2005-01-20  okamoto - add nsig5p1
  26 !   2005-03-04  derber  - add nsig3p3,nsig3p2
  27 !   2005-03-07  dee     - add gmao_intfc option for gmao interface
  28 !   2005-05-24  pondeca - regional surface component added
  29 !   2005-06-01  treadon - add variables msig and array nlayers
  30 !   2005-09-28  derber  - put grid calculations into get_ij and get_ijk
  31 !   2006-01-09  derber  - add sigsum
  32 !   2006-02-01  parrish - correct error to dx_an, dy_an when using filled_grid
  33 !   2006-04-14  treadon - remove global sigi,sigl; add ntracter,ncloud,ck5
  34 !   2006-04-17  treadon - remove regional sigi_ll,sigl_ll
  35 !   2006-10-17  kleist  - add lat dependent coriolis parameter
  36 !   2007-05-07  treadon - add ncep_sigio, ncepgfs_head(v)
  37 !   2007-05-08  kleist  - add variables for fully generalized vertical coordinate
  38 !   2007-10-24  parrish - fix error in wind rotation reference angle field
  39 !
  40 ! !AUTHOR:
  41 !   kleist           org: np20                date: 2003-09-25
  42 !
  43 !EOP
  44 !-------------------------------------------------------------------------
  45
  46   logical regional          ! .t. for regional background/analysis
  47   logical diagnostic_reg    ! .t. to activate regional analysis diagnostics
  48
  49   logical gmao_intfc        ! .t. if model is GMAO gcm
  50
  51   logical ncep_sigio        ! .t. if using ncep sgio format file
  52
  53   logical eta_regional      !
  54   logical wrf_nmm_regional  !
  55   logical wrf_mass_regional !
  56   logical twodvar_regional  ! .t. to run code in regional 2D-var mode
  57   logical netcdf            ! .t. for regional netcdf i/o
  58
  59   logical hybrid            ! .t. to set hybrid vertical coordinates
  60   logical filled_grid       !
  61   logical half_grid         !
  62   logical update_regsfc     !
  63
  64   integer(i_kind) nsig1o            ! no. of levels distributed on each processor
  65   integer(i_kind) nlat              ! no. of latitudes
  66   integer(i_kind) nlon              ! no. of longitudes
  67   integer(i_kind) nsig              ! no. of levels
  68   integer(i_kind) idvc5             ! vertical coordinate identifier
  69 !                                        1: sigma
  70 !                                        2: sigma-pressure
  71 !                                        3: sigma-pressure-theta
  72   integer(i_kind) idvm5
  73   integer(i_kind) idpsfc5           ! surface pressure identifier
  74 !                                      0/1: ln(ps)
  75 !                                        2: ps
  76   integer(i_kind) idthrm5           ! thermodynamic variable identifier
  77 !                                      0/1: virtual temperature
  78 !                                        2: sensible temperature
  79 !                                        3: enthalpy (CpT)
  80   integer(i_kind) idsl5             ! midlayer pressure definition
  81 !                                        1: Philips
  82 !                                        2: average
  83   integer(i_kind) nsig2             ! 2 times number of levels
  84   integer(i_kind) nsig3             ! 3 times number of levels
  85   integer(i_kind) nsig3p1           ! 3 times number of levels plus 1
  86   integer(i_kind) nsig3p2           ! 3 times number of levels plus 2
  87   integer(i_kind) nsig3p3           ! 3 times number of levels plus 3
  88   integer(i_kind) nsig4             ! 4 times number of levels
  89   integer(i_kind) nsig5             ! 5 times number of levels
  90   integer(i_kind) nsig5p1           ! 5 times number of levels plus 1
  91   integer(i_kind) nsig_hlf          ! half number of levels
  92
  93   integer(i_kind) ntracer           ! number of tracers
  94   integer(i_kind) ncloud            ! number of cloud types
  95
  96   integer(i_kind) ns1               ! 2 times number of levels plus 1
  97   integer(i_kind) n1                ! no. of levels plus 1
  98   integer(i_kind) lat1              ! no. of lats on subdomain (no buffer)
  99   integer(i_kind) lon1              ! no. of lons on subdomain (no buffer)
 100   integer(i_kind) lat2              ! no. of lats on subdomain (buffer points on ends)
 101   integer(i_kind) lon2              ! no. of lons on subdomain (buffer points on ends)
 102   integer(i_kind) latlon11          ! horizontal points in subdomain (with buffer)
 103   integer(i_kind) latlon1n          ! no. of points in subdomain (with buffer)
 104   integer(i_kind) iglobal           ! number of horizontal points on global grid
 105   integer(i_kind) itotsub           ! number of horizontal points of all subdomains combined
 106   integer(i_kind) msig              ! number of profile layers to use when calling RTM
 107
 108   logical periodic                              ! logical flag for periodic e/w domains
 109   logical,allocatable,dimension(:):: periodic_s ! logical flag for periodic e/w subdomain (all tasks)
 110
 111   integer(i_kind),allocatable,dimension(:):: jstart  ! start lon of the whole array on each pe
 112   integer(i_kind),allocatable,dimension(:):: istart  ! start lat of the whole array on each pe
 113   integer(i_kind),allocatable,dimension(:):: ilat1   ! no. of lats for each subdomain (no buffer)
 114   integer(i_kind),allocatable,dimension(:):: jlon1   ! no. of lons for each subdomain (no buffer)
 115   integer(i_kind),allocatable,dimension(:):: ijn_s   ! no. of horiz. points for each subdomain (with buffer)
 116   integer(i_kind),allocatable,dimension(:):: ijn     ! no. of horiz. points for each subdomain (no buffer)
 117   integer(i_kind),allocatable,dimension(:):: ijn_s3  ! 3 times no. of horiz. pnts for each subdomain (with buffer)
 118   integer(i_kind),allocatable,dimension(:):: isc_g   ! no. array, count for send to global; size of subdomain
 119   integer(i_kind),allocatable,dimension(:):: ijn3    ! 3 times no. of horiz. pnts for each subdomain (no buffer)
 120
 121                                                ! comm. array ...
 122   integer(i_kind),allocatable,dimension(:):: irc_s     !   count for receive on subdomain
 123   integer(i_kind),allocatable,dimension(:):: irc_s3    !   count for receive on subdomain
 124   integer(i_kind),allocatable,dimension(:):: isc_g3    !   count for send to global; 3*size of subdomain
 125   integer(i_kind),allocatable,dimension(:):: ird_s     !   displacement for receive on subdomain
 126   integer(i_kind),allocatable,dimension(:):: isd_g     !   displacement for send to global
 127   integer(i_kind),allocatable,dimension(:):: ird_s3    !   displacement for receive on subdomain
 128   integer(i_kind),allocatable,dimension(:):: isd_g3    !   displacement for send to global
 129   integer(i_kind),allocatable,dimension(:):: displs_s  !   displacement for send from subdomain
 130   integer(i_kind),allocatable,dimension(:):: displs_s3 !   displacement for send from subdomain
 131   integer(i_kind),allocatable,dimension(:):: displs_g  !   displacement for receive on global grid
 132   integer(i_kind),allocatable,dimension(:):: displs_g3 !   displacement for receive on global grid
 133
 134                                              ! array element indices for location of ...
 135   integer(i_kind),allocatable,dimension(:):: ltosi   !   lats in iglobal array excluding buffer
 136   integer(i_kind),allocatable,dimension(:):: ltosj   !   lons in iglobal array excluding buffer
 137   integer(i_kind),allocatable,dimension(:):: ltosi_s !   lats in itotsub array including buffer
 138   integer(i_kind),allocatable,dimension(:):: ltosj_s !   lons in itotsub array including buffer
 139
 140   integer(i_kind),dimension(100):: nlayers        ! number of RTM layers per model layer
 141                                                   ! (k=1 is near surface layer), default is 1
 142   real(r_kind) dlm0,dph0
 143   real(r_kind) gencode
 144
 145   real(r_kind),allocatable,dimension(:):: sigsum  ! used to estimate 3d pressure increment for qoption = 2
 146   real(r_kind),allocatable,dimension(:):: rlats   ! grid latitudes (radians)
 147   real(r_kind),allocatable,dimension(:):: rlons   ! grid longitudes (radians)
 148   real(r_kind),allocatable,dimension(:):: ak5,bk5,ck5,tref5 ! coefficients for generalized vertical coordinate
 149   real(r_kind),allocatable,dimension(:):: cp5     ! specific heat for tracers
 150   real(r_kind),allocatable,dimension(:):: coslon  ! cos(grid longitudes (radians))
 151   real(r_kind),allocatable,dimension(:):: sinlon  ! sin(grid longitudes (radians))
 152   real(r_kind),allocatable,dimension(:):: wgtlats !  gaussian integration weights
 153   real(r_kind),allocatable,dimension(:):: corlats ! coriolis parameter by latitude
 154   real(r_kind),allocatable,dimension(:):: rbs2    ! 1./sin(grid latitudes))**2
 155
 156 ! additional variables for regional mode
 157   real(r_kind),allocatable:: deta1_ll(:)          !
 158   real(r_kind),allocatable::  eta1_ll(:)          !
 159   real(r_kind),allocatable:: aeta1_ll(:)          !
 160   real(r_kind),allocatable:: deta2_ll(:)          !
 161   real(r_kind),allocatable::  eta2_ll(:)          !
 162   real(r_kind),allocatable:: aeta2_ll(:)          !
 163   real(r_kind),allocatable::region_lon(:,:)       !
 164   real(r_kind),allocatable::region_lat(:,:)       !
 165   real(r_kind),allocatable::region_dx(:,:)        !
 166   real(r_kind),allocatable::region_dy(:,:)        !
 167   real(r_kind),allocatable::coeffx(:,:)           !
 168   real(r_kind),allocatable::coeffy(:,:)           !
 169
 170   real(r_kind) dlon_ll,dlat_ll,dlon_regional,dlat_regional
 171   real(r_kind) rlon_min_ll,rlon_max_ll,rlat_min_ll,rlat_max_ll
 172   real(r_kind) rlon_min_dd,rlon_max_dd,rlat_min_dd,rlat_max_dd
 173   real(r_kind) rlon_min_regional,rlon_max_regional
 174   real(r_kind) rlat_min_regional,rlat_max_regional
 175   real(r_kind) rlon0_origin_ll,rlat0_origin_ll,clat0_origin_ll,slat0_origin_ll,pi_reg_glob_ll
 176   real(r_kind) dt_ll,pdtop_ll,pt_ll
 177
 178   integer(i_kind) nlon_regional,nlat_regional,nsoil_regional
 179   integer(i_kind) itb_regional,jtb_regional
 180   integer(i_kind) order_a2e,order_e2a
 181   integer(i_kind) nhr_assimilation
 182   real(r_kind) regional_fhr
 183   integer(i_kind) regional_time(6)
 184
 185 ! The following is for the generalized transform
 186   real(r_kind) pihalf,sign_pole,rlambda0
 187   real(r_kind) atilde_x,btilde_x,atilde_y,btilde_y
 188   real(r_kind) btilde_xinv,btilde_yinv
 189   integer(i_kind) nxtilde,nytilde
 190   real(r_kind),allocatable::xtilde0(:,:),ytilde0(:,:)
 191   real(r_kind),allocatable::beta_ref(:,:),cos_beta_ref(:,:),sin_beta_ref(:,:)
 192   integer(i_kind),allocatable::i0_tilde(:,:),j0_tilde(:,:)
 193   integer(i_byte),allocatable::ip_tilde(:,:),jp_tilde(:,:)
 194 !----------temporary variables to keep track of number of observations falling in beta_ref jump zone
 195   real(r_kind):: count_beta_diff,count_beta_diff_gt_20
 196   real(r_kind) beta_diff_max,beta_diff_min,beta_diff_rms
 197   real(r_kind) beta_diff_max_gt_20
 198
 199 ! Define structure to hold NCEP sigio/gfsio header information
 200   type:: ncepgfs_head
 201      integer(i_kind):: ivs
 202      integer(i_kind):: version
 203      real(r_single) :: fhour
 204      integer(i_kind):: idate(4)
 205      integer(i_kind):: nrec
 206      integer(i_kind):: latb
 207      integer(i_kind):: lonb
 208      integer(i_kind):: levs
 209      integer(i_kind):: jcap
 210      integer(i_kind):: itrun
 211      integer(i_kind):: iorder
 212      integer(i_kind):: irealf
 213      integer(i_kind):: igen
 214      integer(i_kind):: latf
 215      integer(i_kind):: lonf
 216      integer(i_kind):: latr
 217      integer(i_kind):: lonr
 218      integer(i_kind):: ntrac
 219      integer(i_kind):: icen2
 220      integer(i_kind):: iens(2)
 221      integer(i_kind):: idpp
 222      integer(i_kind):: idsl
 223      integer(i_kind):: idvc
 224      integer(i_kind):: idvm
 225      integer(i_kind):: idvt
 226      integer(i_kind):: idrun
 227      integer(i_kind):: idusr
 228      real(r_single) :: pdryini
 229      integer(i_kind):: ncldt
 230      integer(i_kind):: ixgr
 231      integer(i_kind):: nvcoord
 232      integer(i_kind):: idrt
 233   end type ncepgfs_head
 234
 235   type:: ncepgfs_headv
 236      real(r_single),allocatable:: vcoord(:,:)
 237      real(r_single),allocatable:: cpi(:)
 238   end type ncepgfs_headv
 239
 240 contains
 241
 242 !-------------------------------------------------------------------------
 243 !    NOAA/NCEP, National Centers for Environmental Prediction GSI        !
 244 !-------------------------------------------------------------------------
 245 !BOP
 246 !
 247 ! !IROUTINE:  init_grid --- Initialize defaults for grid related variables
 248 !
 249 ! !INTERFACE:
 250 !
 251   subroutine init_grid
 252
 253 ! !DESCRIPTION: initialize defaults for grid related variables
 254 !
 255 ! !REVISION HISTORY:
 256 !   2003-09-25  kleist
 257 !   2004-05-13  kleist, documentation
 258 !   2004-07-15  todling, protex-compliant prologue
 259 !   2005-03-03  treadon - add implicit none
 260 !   2005-06-01  treadon - add initialization of msig and nlayers
 261 !
 262 ! !REMARKS:
 263 !   language: f90
 264 !   machine:  ibm rs/6000 sp; SGI Origin 2000; Compaq/HP
 265 !
 266 ! !AUTHOR:
 267 !   kleist           org: np20                date: 2003-09-25
 268 !
 269 !EOP
 270 !-------------------------------------------------------------------------
 271     implicit none
 272
 273     integer k
 274
 275     nsig = 42
 276     nsig1o = 7
 277     nlat = 96
 278     nlon = 384
 279     idvc5 = 1
 280     idvm5=0
 281     idpsfc5=1
 282     idthrm5=1
 283     idsl5 = 1
 284     ntracer=1
 285     ncloud=0
 286     gencode = 80
 287     regional = .false.
 288     gmao_intfc = .false.
 289     ncep_sigio = .true.
 290     periodic = .false.
 291     wrf_nmm_regional = .false.
 292     wrf_mass_regional = .false.
 293     twodvar_regional = .false.
 294     netcdf = .false.
 295     hybrid = .false.
 296     filled_grid = .false.
 297     half_grid = .false.
 298     lat1 = nlat
 299     lon1 = nlon
 300     lat2 = lat1+2
 301     lon2 = lon1+2
 302
 303     diagnostic_reg=.false.
 304     update_regsfc=.false.
 305     nlon_regional=0
 306     nlat_regional=0
 307     nhr_assimilation=6
 308
 309     msig = nsig
 310     do k=1,100
 311        nlayers(k) = 1
 312     end do
 313
 314     return
 315   end subroutine init_grid
 316
 317 !-------------------------------------------------------------------------
 318 !    NOAA/NCEP, National Centers for Environmental Prediction GSI        !
 319 !-------------------------------------------------------------------------
 320 !BOP
 321 !
 322 ! !IROUTINE:  init_grid_vars --- Set grid related variables
 323 !
 324 ! !INTERFACE:
 325 !
 326   subroutine init_grid_vars(jcap,npe)
 327
 328 ! !USES:
 329
 330     implicit none
 331
 332 ! !INPUT PARAMETERS:
 333
 334    integer(i_kind),intent(in)::jcap   ! spectral truncation
 335    integer(i_kind),intent(in)::npe    ! number of mpi tasks
 336
 337 ! !DESCRIPTION: set grid related variables (post namelist read)
 338 !
 339 ! !REVISION HISTORY:
 340 !   2003-09-25  kleist
 341 !   2004-05-13  kleist, documentation
 342 !   2004-07-15  todling, protex-compliant prologue
 343 !   2005-06-01  treadon - add computation of msig
 344 !
 345 !   input argument list:
 346 !
 347 !   output argument list:
 348 !
 349 ! !REMARKS:
 350 !   language: f90
 351 !   machine:  ibm rs/6000 sp; SGI Origin 2000; Compaq/HP
 352 !
 353 ! !AUTHOR:
 354 !   kleist           org: np20                date: 2003-09-25
 355 !
 356 !EOP
 357 !-------------------------------------------------------------------------
 358     integer(i_kind) vlevs,k
 359
 360     if(jcap==62) gencode=80.0
 361     ns1=2*nsig+1
 362     nsig2=2*nsig
 363     nsig3=3*nsig
 364     nsig3p1=3*nsig+1
 365     nsig3p2=3*nsig+2
 366     nsig3p3=3*nsig+3
 367     nsig4=4*nsig
 368     nsig5=5*nsig
 369     nsig5p1=5*nsig+1
 370     nsig_hlf=nsig/2
 371     iglobal=nlat*nlon
 372
 373 ! Initialize nsig1o to distribute levs/variables
 374 ! as evenly as possible over the tasks
 375     vlevs=(6*nsig)+4
 376     nsig1o=vlevs/npe
 377     if(mod(vlevs,npe)/=0) nsig1o=nsig1o+1
 378
 379 ! Sum total number of vertical layers for RTM call
 380     msig = 0
 381     do k=1,nsig
 382        msig = msig + nlayers(k)
 383     end do
 384
 385     return
 386   end subroutine init_grid_vars
 387
 388 !-------------------------------------------------------------------------
 389 !    NOAA/NCEP, National Centers for Environmental Prediction GSI        !
 390 !-------------------------------------------------------------------------
 391 !BOP
 392 !
 393 ! !IROUTINE:  init_subdomain_vars --- Initialize variables related to subdomains
 394 !
 395 ! !INTERFACE:
 396 !
 397   subroutine init_subdomain_vars
 398
 399 ! !DESCRIPTION: initialize variables related to subdomains
 400 !
 401 ! !REVISION HISTORY:
 402 !   2003-09-25  kleist
 403 !   2004-05-13  kleist, documentation
 404 !   2004-07-15  todling, protex-compliant prologue
 405 !   2005-03-03  treadon - add implicit none
 406 !
 407 ! !REMARKS:
 408 !   language: f90
 409 !   machine:  ibm rs/6000 sp; SGI Origin 2000; Compaq/HP
 410 !
 411 ! !AUTHOR:
 412 !   kleist           org: np20                date: 2003-09-25
 413 !
 414 !EOP
 415 !-------------------------------------------------------------------------
 416     implicit none
 417
 418     lat2 = lat1+2
 419     lon2 = lon1+2
 420     latlon11 = lat2*lon2
 421     latlon1n = latlon11*nsig
 422
 423     return
 424   end subroutine init_subdomain_vars
 425
 426 !-------------------------------------------------------------------------
 427 !    NOAA/NCEP, National Centers for Environmental Prediction GSI        !
 428 !-------------------------------------------------------------------------
 429 !BOP
 430 !
 431 ! !IROUTINE:  create_grid_vars --- Allocate memory for grid related variables
 432 !
 433 ! !INTERFACE:
 434 !
 435   subroutine create_grid_vars
 436
 437 ! !DESCRIPTION: allocate memory for grid related variables
 438 !
 439 ! !REVISION HISTORY:
 440 !   2003-09-25  kleist
 441 !   2004-05-13  kleist, documentation
 442 !   2004-07-15  todling, protex-compliant prologue
 443 !   2005-03-03  treadon - add implicit none
 444 !
 445 ! !REMARKS:
 446 !   language: f90
 447 !   machine:  ibm rs/6000 sp; SGI Origin 2000; Compaq/HP
 448 !
 449 ! !AUTHOR:
 450 !   kleist           org: np20                date: 2003-09-25
 451 !
 452 !EOP
 453 !-------------------------------------------------------------------------
 454     implicit none
 455
 456     allocate(sigsum(nsig))
 457     allocate(rlats(nlat),rlons(nlon),coslon(nlon),sinlon(nlon),&
 458              wgtlats(nlat),rbs2(nlat),corlats(nlat))
 459     allocate(ak5(nsig+1),bk5(nsig+1),ck5(nsig+1),tref5(nsig))
 460     return
 461   end subroutine create_grid_vars
 462
 463 !-------------------------------------------------------------------------
 464 !    NOAA/NCEP, National Centers for Environmental Prediction GSI        !
 465 !-------------------------------------------------------------------------
 466 !BOP
 467 !
 468 ! !IROUTINE:  destroy_grid_vars --- Deallocate memory for grid related variables
 469 !
 470 ! !INTERFACE:
 471 !
 472   subroutine destroy_grid_vars
 473
 474 ! !DESCRIPTION: deallocate memory for grid related variables
 475 !
 476 ! !REVISION HISTORY:
 477 !   2003-09-25  kleist
 478 !   2004-05-13  kleist, documentation
 479 !   2004-07-15  todling, protex-compliant prologue
 480 !   2005-03-03  treadon - add implicit none
 481 !
 482 ! !REMARKS:
 483 !   language: f90
 484 !   machine:  ibm rs/6000 sp; SGI Origin 2000; Compaq/HP
 485 !
 486 ! !AUTHOR:
 487 !   kleist           org: np20                date: 2003-09-25
 488 !
 489 !EOP
 490 !-------------------------------------------------------------------------
 491     implicit none
 492     deallocate(sigsum)
 493     deallocate(rlats,rlons,corlats,coslon,sinlon,wgtlats,rbs2)
 494     deallocate(ak5,bk5,ck5,tref5)
 495     if (allocated(cp5)) deallocate(cp5)
 496     return
 497   end subroutine destroy_grid_vars
 498
 499 !-------------------------------------------------------------------------
 500 !    NOAA/NCEP, National Centers for Environmental Prediction GSI        !
 501 !-------------------------------------------------------------------------
 502 !BOP
 503 !
 504 ! !IROUTINE:  create_mapping --- Init vars mapping between global domain/subd.
 505 !
 506 ! !INTERFACE:
 507 !
 508   subroutine create_mapping(nlat,nlon,npe)
 509
 510 ! !USES:
 511
 512     use constants, only: izero
 513     implicit none
 514
 515 ! !INPUT PARAMETERS:
 516
 517     integer(i_kind),intent(in):: nlat  ! number of latitudes
 518     integer(i_kind),intent(in):: nlon  ! number of longitudes
 519     integer(i_kind),intent(in):: npe   ! number of mpi tasks
 520
 521 ! !DESCRIPTION: allocate and initialize variables that create mapping
 522 !           between global domain and subdomains
 523 !
 524 ! !REVISION HISTORY:
 525 !   2003-09-25  kleist
 526 !   2004-05-13  kleist, documentation
 527 !   2004-07-15  todling, protex-compliant prologue
 528 !
 529 ! !REMARKS:
 530 !   language: f90
 531 !   machine:  ibm rs/6000 sp; SGI Origin 2000; Compaq/HP
 532 !
 533 ! !AUTHOR:
 534 !   kleist           org: np20                date: 2003-09-25
 535 !
 536 !EOP
 537 !-------------------------------------------------------------------------
 538     integer(i_kind) i
 539
 540     allocate(ltosi(nlat*nlon),ltosj(nlat*nlon),&
 541          ltosi_s(2*nlat*nlon),ltosj_s(2*nlat*nlon))
 542
 543     allocate(periodic_s(npe),jstart(npe),istart(npe),&
 544          ilat1(npe),jlon1(npe),&
 545        ijn_s(npe),irc_s(npe),ird_s(npe),displs_s(npe),&
 546        ijn_s3(npe),irc_s3(npe),ird_s3(npe),displs_s3(npe),&
 547        ijn(npe),isc_g(npe),isd_g(npe),displs_g(npe),&
 548        ijn3(npe),isc_g3(npe),isd_g3(npe),displs_g3(npe))
 549
 550     do i=1,nlat*nlon
 551       ltosi(i)=izero
 552       ltosj(i)=izero
 553     end do
 554
 555     do i=1,2*nlat*nlon
 556       ltosi_s(i)=izero
 557       ltosj_s(i)=izero
 558     end do
 559
 560     do i=1,npe
 561       periodic_s(i)= .false.
 562       jstart(i)    = izero
 563       istart(i)    = izero
 564       ilat1(i)     = izero
 565       jlon1(i)     = izero
 566       ijn_s(i)     = izero
 567       irc_s(i)     = izero
 568       ird_s(i)     = izero
 569       displs_s(i)  = izero
 570       ijn_s3(i)    = izero
 571       irc_s3(i)    = izero
 572       ird_s3(i)    = izero
 573       displs_s3(i) = izero
 574       ijn(i)       = izero
 575       isc_g(i)     = izero
 576       isd_g(i)     = izero
 577       displs_g(i)  = izero
 578       ijn3(i)      = izero
 579       isc_g3(i)    = izero
 580       isd_g3(i)    = izero
 581       displs_g3(i) = izero
 582     end do
 583
 584     return
 585   end subroutine create_mapping
 586
 587 !-------------------------------------------------------------------------
 588 !    NOAA/NCEP, National Centers for Environmental Prediction GSI        !
 589 !-------------------------------------------------------------------------
 590 !BOP
 591 !
 592 ! !IROUTINE:  destroy_mapping --- Dealloc global/subdomain mapping arrays
 593 !
 594 ! !INTERFACE:
 595 !
 596   subroutine destroy_mapping
 597
 598 ! !DESCRIPTION: deallocate memory for global/subdomain mapping variables
 599 !
 600 ! !REVISION HISTORY:
 601 !   2003-09-25  kleist
 602 !   2004-05-13  kleist, documentation
 603 !   2004-07-15  todling, protex-compliant prologue
 604 !   2005-03-03  treadon - add implicit none
 605 !
 606 ! !REMARKS:
 607 !   language: f90
 608 !   machine:  ibm rs/6000 sp; SGI Origin 2000; Compaq/HP
 609 !
 610 ! !AUTHOR:
 611 !   kleist           org: np20                date: 2003-09-25
 612 !
 613 !EOP
 614 !-------------------------------------------------------------------------
 615     implicit none
 616     deallocate(ltosi,ltosj,ltosi_s,ltosj_s)
 617     deallocate(periodic_s,jstart,istart,ilat1,jlon1,&
 618        ijn_s,displs_s,&
 619        ijn_s3,irc_s3,ird_s3,displs_s3,&
 620        ijn,isc_g,isd_g,displs_g,&
 621        ijn3,isc_g3,isd_g3,displs_g3)
 622
 623     return
 624   end subroutine destroy_mapping
 625
 626 !-------------------------------------------------------------------------
 627 !    NOAA/NCEP, National Centers for Environmental Prediction GSI        !
 628 !-------------------------------------------------------------------------
 629 !BOP
 630 !
 631 ! !IROUTINE:  init_reg_glob_ll --- In case regional, initialize setting
 632 !
 633 ! !INTERFACE:
 634 !
 635   subroutine init_reg_glob_ll(mype,lendian_in)
 636
 637 ! !USES:
 638
 639     use kinds, only: r_kind,r_single,i_kind
 640     use constants, only: zero, one, three, deg2rad,pi,half, two
 641     implicit none
 642
 643 ! !INPUT PARAMETERS:
 644
 645     integer(i_kind), intent(in) :: mype          ! mpi task id
 646     integer(i_kind), intent(in) :: lendian_in    ! unit number reserved for
 647                                                  !  little endian input
 648
 649 ! !DESCRIPTION: decide if regional run or not, and initialize constants
 650 !           required for rotation transformation
 651 !
 652 !
 653 !   output argument list:
 654 !
 655 !   Notes about grid definition:
 656 !   \begin{enumerate}
 657 !   \item  The origin of the analysis coordinate system is always $rlon=180.$, $rlat=0.$,
 658 !          whether this is a global or regional run.  The point $rlon=180$, $rlat=0$ in
 659 !          the rotated coordinate coincides with the point rlon0\_origin, rlat0\_origin
 660 !          in earth coordinates.  This is why $rlon0_origin=180$. in the global case.
 661 !
 662 !   \item  For regional runs, the rotated coordinate origin and extent of the domain are read
 663 !          in from the NMM restart file.
 664 !
 665 !   \item  The reason for having the longitude of the origin = 180 is because there are
 666 !          places in the global analysis that depend on $0 < lon < 360$.  So to minimize changes
 667 !          to the global code, this approach has been adopted.
 668 !
 669 !   \item  The regional analysis domain is larger than the corresponding NMM grid.  A halo is included
 670 !          whose width is a function of the interpolation order for transfers between grids.
 671 !          This is so the analysis increment is always being interpolated and added on to the
 672 !          full model domain.
 673 !   \end{enumerate}
 674 !
 675 ! !REVISION HISTORY:
 676 !   2003-08-28  parrish
 677 !   2004-05-13  kleist, documentation
 678 !   2004-07-15  todling, protex-compliant prologue
 679 !   2004-12-15  treadon - explicity set value for inges
 680 !   2005-05-24  pondeca - add the surface analysis option
 681 !   2006-04-06  middlecoff - changed inges from 21 to lendian_in so it can be set to little endian.
 682 !
 683 ! !REMARKS:
 684 !   language: f90
 685 !   machine:  ibm rs/6000 sp; SGI Origin 2000; Compaq/HP
 686 !
 687 ! !AUTHOR:
 688 !   parrish          org: np22                date: 2003-08-28
 689 !
 690 !EOP
 691 !-------------------------------------------------------------------------
 692
 693     real(r_kind) del,pthis,rlongl
 694     logical fexist
 695     real(r_single) dt
 696     integer(i_kind) nfcst,nbc,list
 697     integer(i_kind) im,jm,lm
 698     integer(i_kind) i,j,k,kk
 699     real(r_single)dy,cpgfv,en,ent,r,pt,tddamp,f4d,f4q,ef4t,pdtop
 700     real(r_single),allocatable:: deta1(:),aeta1(:),eta1(:),deta2(:),aeta2(:),eta2(:)
 701     real(r_single),allocatable:: ptbl(:,:),ttbl(:,:)
 702     real(r_single) r1,pt1,tsph,wbd,sbd,tlm0d,tph0d,dlmd,dphd
 703     real(r_single),allocatable:: glat(:,:),glon(:,:)
 704     real(r_single),allocatable:: dx_nmm(:,:),dy_nmm(:,:)
 705     real(r_single),allocatable:: dx_mc(:,:),dy_mc(:,:)
 706
 707     real(r_kind),parameter:: r0_01=0.01_r_kind
 708     real(r_kind),parameter:: r1_5=1.5_r_kind
 709     real(r_kind),parameter:: six=6.0_r_kind
 710     real(r_kind),parameter:: r90=90.0_r_kind
 711     real(r_kind),parameter:: r360=360.0_r_kind
 712     real(r_kind),parameter:: r1013=1013.0_r_kind
 713
 714     real(r_kind),allocatable::glat_an(:,:),glon_an(:,:)
 715     real(r_kind),allocatable:: dx_an(:,:),dy_an(:,:)
 716     integer(i_kind) ierror
 717     character(6) filename
 718     real(r_kind) rlong_sw,rlatg_sw,rlatgl,errmax
 719     integer(i_kind) ihr
 720     logical outside
 721     real(r_kind),allocatable::gxtemp(:,:),gytemp(:,:)
 722     real(r_kind),allocatable::gxtemp_an(:,:),gytemp_an(:,:)
 723     real(r_kind) rtemp
 724     real(r_kind) rlon_min_ll,rlon_max_ll,rlat_min_ll,rlat_max_ll
 725
 726     if(.not.regional) then
 727 ! This is global run
 728       rlat_min_ll=-r90*deg2rad
 729       rlat_max_ll=r90*deg2rad
 730       rlon_min_ll=zero*deg2rad
 731       rlon_max_ll=r360*deg2rad
 732       rlon_min_dd=rlon_min_ll-deg2rad
 733       rlon_max_dd=rlon_max_ll+deg2rad
 734       rlat_min_dd=rlat_min_ll-deg2rad
 735       rlat_max_dd=rlat_max_ll+deg2rad
 736       dt_ll=zero
 737     end if
 738
 739     if(wrf_nmm_regional) then     ! begin wrf_nmm section
 740 ! This is a wrf_nmm regional run.
 741       if(diagnostic_reg.and.mype.eq.0)  &
 742         write(6,*)' in init_reg_glob_ll, initializing for wrf nmm regional run'
 743
 744 ! Get regional constants
 745       ihr=-999
 746       do i=0,12
 747         write(filename,'("sigf",i2.2)')i
 748         inquire(file=filename,exist=fexist)
 749         if(fexist) then
 750           ihr=i
 751           exit
 752         end if
 753       end do
 754       if(ihr.lt.0) then
 755         stop ' NO INPUT FILE AVAILABLE FOR REGIONAL (WRFNMM) ANALYSIS.  PROGRAM STOPS'
 756       end if
 757
 758       if(diagnostic_reg.and.mype.eq.0) write(6,*)' in init_reg_glob_ll, lendian_in=',lendian_in
 759       open(lendian_in,file=filename,form='unformatted')
 760       rewind lendian_in
 761       read(lendian_in) regional_time,nlon_regional,nlat_regional,nsig, &
 762                   dlmd,dphd,pt,pdtop
 763       regional_fhr=zero  !  with wrf nmm fcst hr is not currently available.
 764
 765       if(diagnostic_reg.and.mype.eq.0) write(6,'(" in init_reg_glob_ll, yr,mn,dy,h,m,s=",6i6)') &
 766                regional_time
 767       if(diagnostic_reg.and.mype.eq.0) write(6,'(" in init_reg_glob_ll, nlon_regional=",i6)') &
 768                nlon_regional
 769       if(diagnostic_reg.and.mype.eq.0) write(6,'(" in init_reg_glob_ll, nlat_regional=",i6)') &
 770                nlat_regional
 771       if(diagnostic_reg.and.mype.eq.0) write(6,'(" in init_reg_glob_ll, nsig=",i6)') nsig
 772
 773 ! Get vertical info for hybrid coordinate and sigma coordinate we will interpolate to
 774       allocate(aeta1_ll(nsig),eta1_ll(nsig+1),aeta2_ll(nsig),eta2_ll(nsig+1))
 775       allocate(deta1(nsig),aeta1(nsig),eta1(nsig+1),deta2(nsig),aeta2(nsig),eta2(nsig+1))
 776       allocate(glat(nlon_regional,nlat_regional),glon(nlon_regional,nlat_regional))
 777       allocate(dx_nmm(nlon_regional,nlat_regional),dy_nmm(nlon_regional,nlat_regional))
 778       read(lendian_in) deta1
 779       read(lendian_in) aeta1
 780       read(lendian_in) eta1
 781       read(lendian_in) deta2
 782       read(lendian_in) aeta2
 783       read(lendian_in) eta2
 784
 785       if(diagnostic_reg.and.mype.eq.0) write(6,*)' in init_reg_glob_ll, pdtop,pt=',pdtop,pt
 786       if(diagnostic_reg.and.mype.eq.0) then
 787         write(6,*)' in init_reg_glob_ll, aeta1 aeta2 follow:'
 788         do k=1,nsig
 789           write(6,'(" k,aeta1,aeta2=",i3,2f10.4)') k,aeta1(k),aeta2(k)
 790         end do
 791         write(6,*)' in init_reg_glob_ll, deta1 deta2 follow:'
 792         do k=1,nsig
 793           write(6,'(" k,deta1,deta2=",i3,2f10.4)') k,deta1(k),deta2(k)
 794         end do
 795         write(6,*)' in init_reg_glob_ll, deta1 deta2 follow:'
 796         do k=1,nsig+1
 797           write(6,'(" k,eta1,eta2=",i3,2f10.4)') k,eta1(k),eta2(k)
 798         end do
 799       end if
 800
 801       pdtop_ll=r0_01*pdtop                    !  check units--this converts to mb
 802       pt_ll=r0_01*pt                          !  same here
 803
 804       if(diagnostic_reg.and.mype.eq.0) write(6,*)' in init_reg_glob_ll, pdtop_ll,pt_ll=',pdtop_ll,pt_ll
 805       eta1_ll=eta1
 806       aeta1_ll=aeta1
 807       eta2_ll=eta2
 808       aeta2_ll=aeta2
 809       read(lendian_in) glat,dx_nmm
 810       read(lendian_in) glon,dy_nmm
 811       close(lendian_in)
 812
 813       rlon_min_ll=1
 814       rlat_min_ll=1
 815       if(filled_grid) then
 816         nlon=2*nlon_regional-1
 817         nlat=nlat_regional
 818         rlon_max_ll=nlon
 819         rlat_max_ll=nlat
 820         rlat_min_dd=rlat_min_ll+three
 821         rlat_max_dd=rlat_max_ll-three
 822         rlon_min_dd=rlon_min_ll+six
 823         rlon_max_dd=rlon_max_ll-six
 824       end if
 825       if(half_grid) then
 826         nlon=nlon_regional
 827         nlat=1+nlat_regional/2
 828         rlon_max_ll=nlon
 829         rlat_max_ll=nlat
 830         rlat_min_dd=rlat_min_ll+r1_5
 831         rlat_max_dd=rlat_max_ll-r1_5
 832         rlon_min_dd=rlon_min_ll+three
 833         rlon_max_dd=rlon_max_ll-three
 834       end if
 835
 836       if(diagnostic_reg.and.mype.eq.0) then
 837         write(6,*)' in init_reg_glob_ll, rlat_min_dd=',rlat_min_dd
 838         write(6,*)' in init_reg_glob_ll, rlat_max_dd=',rlat_max_dd
 839         write(6,*)' in init_reg_glob_ll, rlon_min_dd=',rlon_min_dd
 840         write(6,*)' in init_reg_glob_ll, rlon_max_dd=',rlon_max_dd
 841         write(6,*)' in init_reg_glob_ll, rlat_min_ll=',rlat_min_ll
 842         write(6,*)' in init_reg_glob_ll, rlat_max_ll=',rlat_max_ll
 843         write(6,*)' in init_reg_glob_ll, rlon_min_ll=',rlon_min_ll
 844         write(6,*)' in init_reg_glob_ll, rlon_max_ll=',rlon_max_ll
 845         write(6,*)' in init_reg_glob_ll, filled_grid,half_grid=',filled_grid,half_grid
 846         write(6,*)' in init_reg_glob_ll, nlon,nlat=',nlon,nlat
 847       end if
 848
 849       allocate(region_lat(nlat,nlon),region_lon(nlat,nlon))
 850       allocate(region_dy(nlat,nlon),region_dx(nlat,nlon))
 851       allocate(coeffy(nlat,nlon),coeffx(nlat,nlon))
 852
 853 !   generate earth lats and lons on analysis grid
 854
 855       allocate(glat_an(nlon,nlat),glon_an(nlon,nlat))
 856       allocate(dx_an(nlon,nlat),dy_an(nlon,nlat))
 857
 858       if(half_grid) then
 859        call half_nmm_grid2a(glon,nlon_regional,nlat_regional,glon_an,1)
 860        call half_nmm_grid2a(glat,nlon_regional,nlat_regional,glat_an,1)
 861        call half_nmm_grid2a(dx_nmm,nlon_regional,nlat_regional,dx_an,1)
 862        call half_nmm_grid2a(dy_nmm,nlon_regional,nlat_regional,dy_an,1)
 863        dx_an=two*dx_an
 864        dy_an=two*dy_an
 865       end if
 866
 867       if(filled_grid) then
 868        allocate(gxtemp(nlon_regional,nlat_regional))
 869        allocate(gytemp(nlon_regional,nlat_regional))
 870        allocate(gxtemp_an(nlon,nlat))
 871        allocate(gytemp_an(nlon,nlat))
 872        sign_pole=-one
 873        pihalf=half*pi
 874        if(maxval(glat)/deg2rad.lt.zero) sign_pole=one
 875        do j=1,nlat_regional
 876         do i=1,nlon_regional
 877          rtemp=pihalf-sign_pole*glat(i,j)
 878          gxtemp(i,j)=rtemp*cos(one*glon(i,j))
 879          gytemp(i,j)=rtemp*sin(one*glon(i,j))
 880         end do
 881        end do
 882        call fill_nmm_grid2a3(gxtemp,nlon_regional,nlat_regional,gxtemp_an)
 883        call fill_nmm_grid2a3(gytemp,nlon_regional,nlat_regional,gytemp_an)
 884        do j=1,nlat
 885         do i=1,nlon
 886          rtemp=sqrt(gxtemp_an(i,j)**2+gytemp_an(i,j)**2)
 887          glat_an(i,j)=sign_pole*(pihalf-rtemp)
 888          glon_an(i,j)=atan2(gytemp_an(i,j),gxtemp_an(i,j))
 889         end do
 890        end do
 891        gxtemp=dx_nmm
 892        gytemp=dy_nmm
 893        call fill_nmm_grid2a3(gxtemp,nlon_regional,nlat_regional,dx_an)
 894        call fill_nmm_grid2a3(gytemp,nlon_regional,nlat_regional,dy_an)
 895        deallocate(gxtemp,gytemp,gxtemp_an,gytemp_an)
 896       end if
 897
 898       do k=1,nlon
 899         do i=1,nlat
 900           region_lat(i,k)=glat_an(k,i)
 901           region_lon(i,k)=glon_an(k,i)
 902           region_dy(i,k)=dy_an(k,i)
 903           region_dx(i,k)=dx_an(k,i)
 904           coeffy(i,k)=half/dy_an(k,i)
 905           coeffx(i,k)=half/dx_an(k,i)
 906         end do
 907       end do
 908
 909 ! ???????  later change glat_an,glon_an to region_lat,region_lon, with dimensions flipped
 910       call init_general_transform(glat_an,glon_an,mype)
 911
 912       deallocate(deta1,aeta1,eta1,deta2,aeta2,eta2,glat,glon,glat_an,glon_an)
 913       deallocate(dx_nmm,dy_nmm,dx_an,dy_an)
 914
 915     end if   ! end if wrf_nmm section
 916
 917     if(wrf_mass_regional) then     ! begin wrf mass core section
 918 ! This is a wrf_mass regional run.
 919       if(diagnostic_reg.and.mype.eq.0) &
 920           write(6,*)' in init_reg_glob_ll, initializing for wrf mass core regional run'
 921
 922 ! Get regional constants
 923       ihr=-999
 924       do i=0,12
 925         write(filename,'("sigf",i2.2)')i
 926         inquire(file=filename,exist=fexist)
 927         if(fexist) then
 928           ihr=i
 929           exit
 930         end if
 931       end do
 932       if(ihr.lt.0) then
 933         stop ' NO INPUT FILE AVAILABLE FOR REGIONAL (WRF MASS CORE) ANALYSIS.  PROGRAM STOPS'
 934       end if
 935       if(diagnostic_reg.and.mype.eq.0) write(6,*)' in init_reg_glob_ll, lendian_in=',lendian_in
 936       open(lendian_in,file=filename,form='unformatted')
 937       rewind lendian_in
 938       read(lendian_in) regional_time,nlon_regional,nlat_regional,nsig,pt
 939       regional_fhr=zero  !  with wrf mass core fcst hr is not currently available.
 940
 941       if(diagnostic_reg.and.mype.eq.0) write(6,'(" in init_reg_glob_ll, yr,mn,dy,h,m,s=",6i6)') &
 942                regional_time
 943       if(diagnostic_reg.and.mype.eq.0) write(6,'(" in init_reg_glob_ll, nlon_regional=",i6)') &
 944                nlon_regional
 945       if(diagnostic_reg.and.mype.eq.0) write(6,'(" in init_reg_glob_ll, nlat_regional=",i6)') &
 946                nlat_regional
 947       if(diagnostic_reg.and.mype.eq.0) write(6,'(" in init_reg_glob_ll, nsig=",i6)') nsig
 948
 949 ! Get vertical info for wrf mass core
 950       allocate(aeta1_ll(nsig),eta1_ll(nsig+1))
 951       allocate(aeta1(nsig),eta1(nsig+1))
 952       allocate(glat(nlon_regional,nlat_regional),glon(nlon_regional,nlat_regional))
 953       allocate(dx_mc(nlon_regional,nlat_regional),dy_mc(nlon_regional,nlat_regional))
 954       read(lendian_in) aeta1
 955       read(lendian_in) eta1
 956
 957       if(diagnostic_reg.and.mype.eq.0) write(6,*)' in init_reg_glob_ll, pt=',pt
 958       if(diagnostic_reg.and.mype.eq.0) then
 959         write(6,*)' in init_reg_glob_ll, aeta1 follows:'
 960         do k=1,nsig
 961           write(6,'(" k,aeta1=",i3,f10.4)') k,aeta1(k)
 962         end do
 963         write(6,*)' in init_reg_glob_ll, eta1 follows:'
 964         do k=1,nsig+1
 965           write(6,'(" k,eta1=",i3,f10.4)') k,eta1(k)
 966         end do
 967       end if
 968
 969       pt_ll=r0_01*pt                    !  check units--this converts to mb
 970
 971       if(diagnostic_reg.and.mype.eq.0) write(6,*)' in init_reg_glob_ll, pt_ll=',pt_ll
 972       eta1_ll=eta1
 973       aeta1_ll=aeta1
 974       read(lendian_in) glat,dx_mc
 975       read(lendian_in) glon,dy_mc
 976       close(lendian_in)
 977
 978       rlon_min_ll=1
 979       rlat_min_ll=1
 980       nlon=nlon_regional
 981       nlat=nlat_regional
 982       rlon_max_ll=nlon
 983       rlat_max_ll=nlat
 984       rlat_min_dd=rlat_min_ll+r1_5
 985       rlat_max_dd=rlat_max_ll-r1_5
 986       rlon_min_dd=rlon_min_ll+r1_5
 987       rlon_max_dd=rlon_max_ll-r1_5
 988
 989       if(diagnostic_reg.and.mype.eq.0) then
 990         write(6,*)' in init_reg_glob_ll, rlat_min_dd=',rlat_min_dd
 991         write(6,*)' in init_reg_glob_ll, rlat_max_dd=',rlat_max_dd
 992         write(6,*)' in init_reg_glob_ll, rlon_min_dd=',rlon_min_dd
 993         write(6,*)' in init_reg_glob_ll, rlon_max_dd=',rlon_max_dd
 994         write(6,*)' in init_reg_glob_ll, rlat_min_ll=',rlat_min_ll
 995         write(6,*)' in init_reg_glob_ll, rlat_max_ll=',rlat_max_ll
 996         write(6,*)' in init_reg_glob_ll, rlon_min_ll=',rlon_min_ll
 997         write(6,*)' in init_reg_glob_ll, rlon_max_ll=',rlon_max_ll
 998         write(6,*)' in init_reg_glob_ll, nlon,nlat=',nlon,nlat
 999       end if
1000
1001       allocate(region_lat(nlat,nlon),region_lon(nlat,nlon))
1002       allocate(region_dy(nlat,nlon),region_dx(nlat,nlon))
1003       allocate(coeffy(nlat,nlon),coeffx(nlat,nlon))
1004
1005 !   trasfer earth lats and lons to arrays region_lat, region_lon
1006
1007       allocate(glat_an(nlon,nlat),glon_an(nlon,nlat))
1008       do k=1,nlon
1009         do i=1,nlat
1010           glat_an(k,i)=glat(k,i)
1011           glon_an(k,i)=glon(k,i)
1012           region_lat(i,k)=glat(k,i)
1013           region_lon(i,k)=glon(k,i)
1014           region_dx(i,k)=dx_mc(k,i)
1015           region_dy(i,k)=dy_mc(k,i)
1016           coeffx(i,k)=half/dx_mc(k,i)
1017           coeffy(i,k)=half/dy_mc(k,i)
1018         end do
1019       end do
1020
1021 ! ???????  later change glat_an,glon_an to region_lat,region_lon, with dimensions flipped
1022       call init_general_transform(glat_an,glon_an,mype)
1023
1024       deallocate(aeta1,eta1,glat,glon,glat_an,glon_an)
1025       deallocate(dx_mc,dy_mc)
1026
1027     end if   ! end if wrf_nmm section
1028
1029 !   Begin surface analysis section (regional 2D-var)
1030     if(twodvar_regional) then
1031
1032 ! This is a surface analysis regional run.
1033       if(diagnostic_reg.and.mype.eq.0) &
1034           write(6,*)' in init_reg_glob_ll, initializing for surface analysis regional run'
1035
1036 ! Get regional constants
1037       ihr=-999
1038       do i=0,12
1039         write(filename,'("sigf",i2.2)')i
1040         inquire(file=filename,exist=fexist)
1041         if(fexist) then
1042           ihr=i
1043           exit
1044         end if
1045       end do
1046       if(ihr.lt.0) then
1047         stop ' NO INPUT FILE AVAILABLE FOR REGIONAL (SURFACE) ANALYSIS.  PROGRAM STOPS'
1048       end if
1049       if(diagnostic_reg.and.mype.eq.0) write(6,*)' in init_reg_glob_ll, lendian_in=',lendian_in
1050       open(lendian_in,file=filename,form='unformatted')
1051       rewind lendian_in
1052       read(lendian_in) regional_time,nlon_regional,nlat_regional,nsig
1053       regional_fhr=zero  !  with twodvar analysis fcst hr is not currently available.
1054
1055       if(diagnostic_reg.and.mype.eq.0) write(6,'(" in init_reg_glob_ll, yr,mn,dy,h,m,s=",6i6)') &
1056                regional_time
1057       if(diagnostic_reg.and.mype.eq.0) write(6,'(" in init_reg_glob_ll, nlon_regional=",i6)') &
1058                nlon_regional
1059       if(diagnostic_reg.and.mype.eq.0) write(6,'(" in init_reg_glob_ll, nlat_regional=",i6)') &
1060                nlat_regional
1061       if(diagnostic_reg.and.mype.eq.0) write(6,'(" in init_reg_glob_ll, nsig=",i6)') nsig
1062
1063 ! Get vertical info
1064       allocate(aeta1_ll(nsig),eta1_ll(nsig+1))
1065       allocate(aeta1(nsig),eta1(nsig+1))
1066       allocate(glat(nlon_regional,nlat_regional),glon(nlon_regional,nlat_regional))
1067       allocate(dx_mc(nlon_regional,nlat_regional),dy_mc(nlon_regional,nlat_regional))
1068
1069       aeta1=one                ! set to this value for convenience
1070       eta1=one                 ! set to this value for convenience
1071       pt=0._r_single           ! set to this value for convenience
1072
1073       pt_ll=r0_01*pt
1074       eta1_ll=eta1
1075       aeta1_ll=aeta1
1076
1077       read(lendian_in) glat,dx_mc
1078       read(lendian_in) glon,dy_mc
1079       close(lendian_in)
1080
1081       rlon_min_ll=1
1082       rlat_min_ll=1
1083       nlon=nlon_regional
1084       nlat=nlat_regional
1085       rlon_max_ll=nlon
1086       rlat_max_ll=nlat
1087       rlat_min_dd=rlat_min_ll+r1_5
1088       rlat_max_dd=rlat_max_ll-r1_5
1089       rlon_min_dd=rlon_min_ll+r1_5
1090       rlon_max_dd=rlon_max_ll-r1_5
1091
1092       if(diagnostic_reg.and.mype.eq.0) then
1093         write(6,*)' in init_reg_glob_ll, rlat_min_dd=',rlat_min_dd
1094         write(6,*)' in init_reg_glob_ll, rlat_max_dd=',rlat_max_dd
1095         write(6,*)' in init_reg_glob_ll, rlon_min_dd=',rlon_min_dd
1096         write(6,*)' in init_reg_glob_ll, rlon_max_dd=',rlon_max_dd
1097         write(6,*)' in init_reg_glob_ll, rlat_min_ll=',rlat_min_ll
1098         write(6,*)' in init_reg_glob_ll, rlat_max_ll=',rlat_max_ll
1099         write(6,*)' in init_reg_glob_ll, rlon_min_ll=',rlon_min_ll
1100         write(6,*)' in init_reg_glob_ll, rlon_max_ll=',rlon_max_ll
1101         write(6,*)' in init_reg_glob_ll, nlon,nlat=',nlon,nlat
1102       end if
1103
1104       allocate(region_lat(nlat,nlon),region_lon(nlat,nlon))
1105       allocate(region_dy(nlat,nlon),region_dx(nlat,nlon))
1106       allocate(coeffy(nlat,nlon),coeffx(nlat,nlon))
1107
1108 !   transfer earth lats and lons to arrays region_lat, region_lon
1109
1110       allocate(glat_an(nlon,nlat),glon_an(nlon,nlat))
1111       do k=1,nlon
1112         do i=1,nlat
1113           glat_an(k,i)=glat(k,i)
1114           glon_an(k,i)=glon(k,i)
1115           region_lat(i,k)=glat(k,i)
1116           region_lon(i,k)=glon(k,i)
1117           region_dx(i,k)=dx_mc(k,i)
1118           region_dy(i,k)=dy_mc(k,i)
1119           coeffx(i,k)=half/dx_mc(k,i)
1120           coeffy(i,k)=half/dy_mc(k,i)
1121         end do
1122       end do
1123
1124 ! ???????  later change glat_an,glon_an to region_lat,region_lon, with dimensions flipped
1125       call init_general_transform(glat_an,glon_an,mype)
1126
1127       deallocate(aeta1,eta1,glat,glon,glat_an,glon_an)
1128       deallocate(dx_mc,dy_mc)
1129
1130     end if   ! end if twodvar analysis section
1131
1132     return
1133   end subroutine init_reg_glob_ll
1134
1135  subroutine init_general_transform(glats,glons,mype)
1136
1137   use constants, only: zero,one,half,pi,deg2rad,two
1138   implicit none
1139   real(r_kind) glats(nlon,nlat),glons(nlon,nlat)
1140   integer(i_kind) mype
1141
1142   real(r_kind),parameter:: r0_01=0.01_r_kind
1143   real(r_kind) xbar_min,xbar_max,ybar_min,ybar_max
1144   real(r_kind) clon,slon,r_of_lat,xbar,ybar
1145   integer(i_kind) i,j,istart0,iend,iinc,itemp,ilast,jlast
1146   real(r_kind) cosalpha,sinalpha,denom,epslon,r0,r1,x0,x1,x2,y0,y1,y2
1147   integer(i_kind) ip
1148
1149   pihalf=half*pi
1150
1151 !  define xtilde, ytilde grid, transform
1152
1153 !      glons,glats are lons, lats of input grid points of dimension nlon,nlat
1154        if(mype.eq.0) write(6,*)' at  1 in init_general_transform'
1155   call get_xytilde_domain(nlon,nlat,glons,glats,nxtilde,nytilde, &
1156                    xbar_min,xbar_max,ybar_min,ybar_max)
1157   allocate(i0_tilde(nxtilde,nytilde),j0_tilde(nxtilde,nytilde))
1158   allocate(ip_tilde(nxtilde,nytilde),jp_tilde(nxtilde,nytilde))
1159   allocate(xtilde0(nlon,nlat),ytilde0(nlon,nlat))
1160
1161 ! define atilde_x, btilde_x, atilde_y, btilde_y
1162
1163   btilde_x=(nxtilde-one)/(xbar_max-xbar_min)
1164   btilde_xinv=(xbar_max-xbar_min)/(nxtilde-one)
1165   atilde_x=one-btilde_x*xbar_min
1166   btilde_y=(nytilde-one)/(ybar_max-ybar_min)
1167   btilde_yinv=(ybar_max-ybar_min)/(nytilde-one)
1168   atilde_y=one-btilde_y*ybar_min
1169
1170 ! define xtilde0,ytilde0
1171   do j=1,nlat
1172    do i=1,nlon
1173     r_of_lat=pihalf+sign_pole*glats(i,j)
1174     clon=cos(glons(i,j)+rlambda0)
1175     slon=sin(glons(i,j)+rlambda0)
1176     xbar=r_of_lat*clon
1177     ybar=r_of_lat*slon
1178     xtilde0(i,j)=atilde_x+btilde_x*xbar
1179     ytilde0(i,j)=atilde_y+btilde_y*ybar
1180    end do
1181   end do
1182
1183 !  now get i0_tilde, j0_tilde, ip_tilde,jp_tilde
1184   ilast=1 ; jlast=1
1185   istart0=nxtilde
1186   iend=1
1187   iinc=-1
1188   do j=1,nytilde
1189    itemp=istart0
1190    istart0=iend
1191    iend=itemp
1192    iinc=-iinc
1193    ybar=j
1194    do i=istart0,iend,iinc
1195     xbar=i
1196     call nearest_3(ilast,jlast,i0_tilde(i,j),j0_tilde(i,j), &
1197                    ip_tilde(i,j),jp_tilde(i,j),xbar,ybar,nlon,nlat,xtilde0,ytilde0)
1198    end do
1199   end do
1200
1201 !  now compute beta_ref, used in alpha = beta_ref + sign_pole*earth_lon, and alpha is
1202 !   angle between earth positive east and grid positive x.  This is needed
1203 !   for rotation of u,v from earth to grid coordinate.
1204   allocate(beta_ref(nlon,nlat),cos_beta_ref(nlon,nlat),sin_beta_ref(nlon,nlat))
1205   epslon=r0_01*deg2rad
1206   do j=1,nlat
1207    do i=1,nlon-1
1208     ip=i+1
1209     r0=two*cos(glats(i,j))/(one-sign_pole*sin(glats(i,j)))
1210     x0=r0*cos(glons(i,j))
1211     y0=r0*sin(glons(i,j))
1212     r1=two*cos(glats(ip,j))/(one-sign_pole*sin(glats(ip,j)))
1213     x1=r1*cos(glons(ip,j))
1214     y1=r1*sin(glons(ip,j))
1215     x2=r0*cos(glons(i,j)+epslon)
1216     y2=r0*sin(glons(i,j)+epslon)
1217     denom=one/sqrt(((x1-x0)**2+(y1-y0)**2)*((x2-x0)**2+(y2-y0)**2))
1218     cosalpha=((x2-x0)*(x1-x0)+(y2-y0)*(y1-y0))*denom
1219     sinalpha=((x2-x0)*(y1-y0)-(y2-y0)*(x1-x0))*denom
1220     beta_ref(i,j)=atan2(sinalpha,cosalpha)-sign_pole*glons(i,j)
1221     cos_beta_ref(i,j)=cos(beta_ref(i,j))
1222     sin_beta_ref(i,j)=sin(beta_ref(i,j))
1223    end do
1224    i=nlon
1225    ip=nlon-1
1226    r0=two*cos(glats(i,j))/(one-sign_pole*sin(glats(i,j)))
1227    x0=r0*cos(glons(i,j))
1228    y0=r0*sin(glons(i,j))
1229    r1=two*cos(glats(ip,j))/(one-sign_pole*sin(glats(ip,j)))
1230    x1=r1*cos(glons(ip,j))
1231    y1=r1*sin(glons(ip,j))
1232    x2=r0*cos(glons(i,j)-epslon)
1233    y2=r0*sin(glons(i,j)-epslon)
1234    denom=one/sqrt(((x1-x0)**2+(y1-y0)**2)*((x2-x0)**2+(y2-y0)**2))
1235    cosalpha=((x2-x0)*(x1-x0)+(y2-y0)*(y1-y0))*denom
1236    sinalpha=((x2-x0)*(y1-y0)-(y2-y0)*(x1-x0))*denom
1237    beta_ref(i,j)=atan2(sinalpha,cosalpha)-sign_pole*glons(i,j)
1238    cos_beta_ref(i,j)=cos(beta_ref(i,j))
1239    sin_beta_ref(i,j)=sin(beta_ref(i,j))
1240   end do
1241   beta_diff_max=-huge(beta_diff_max)
1242   beta_diff_max_gt_20=-huge(beta_diff_max_gt_20)
1243   beta_diff_min= huge(beta_diff_min)
1244   beta_diff_rms=zero
1245   count_beta_diff=zero
1246   count_beta_diff_gt_20=zero
1247
1248  end subroutine init_general_transform
1249
1250 !-------------------------------------------------------------------------
1251 !    NOAA/NCEP, National Centers for Environmental Prediction GSI        !
1252 !-------------------------------------------------------------------------
1253 !BOP
1254 !
1255 ! !IROUTINE:  tll2xy --- convert earth lon-lat to x-y grid coordinates
1256 !
1257 ! !INTERFACE:
1258 !
1259   subroutine tll2xy(rlon,rlat,x,y,outside)
1260
1261 ! !USES:
1262
1263     use kinds, only: r_kind,i_kind
1264     use constants, only: one
1265     implicit none
1266
1267     real(r_kind),intent(in)::rlon  ! earth longitude (radians)
1268     real(r_kind),intent(in)::rlat  ! earth latitude  (radians)
1269
1270 ! !OUTPUT PARAMETERS:
1271
1272     real(r_kind),intent(out)::x  ! x-grid coordinate (grid units)
1273     real(r_kind),intent(out)::y  ! y-grid coordinate (grid units)
1274     logical,intent(out)::outside     ! .false., then point is inside x-y domain
1275                                      ! .true.,  then point is outside x-y domain
1276
1277 ! !DESCRIPTION: to convert earth lon-lat to x-y grid units of a
1278 !           general regional rectangular domain.  Also, decide if
1279 !           point is inside this domain.  As a result, there is
1280 !           no restriction on type of horizontal coordinate for
1281 !           a regional run, other than that it not have periodicity
1282 !           or polar singularities.
1283 !           This is done by first converting rlon, rlat to an
1284 !           intermediate coordinate xtilde,ytilde, which has
1285 !           precomputed pointers and constants for final conversion
1286 !           to the desired x,y via 3 point inverse interpolation.
1287 !           All of the information needed is derived from arrays
1288 !           specifying earth latitude and longitude of every point
1289 !           on the input grid.  Currently, the input x-y grid that
1290 !           this is based on must be non-staggered.  This restriction
1291 !           will eventually be lifted so we can run directly from
1292 !           model grids that are staggered without first resorting
1293 !           to interpolation of the guess to a non-staggered grid.
1294 !
1295 ! !REVISION HISTORY:
1296 !   2003-08-28  parrish
1297 !   2004-05-13  kleist, documentation
1298 !   2004-07-15  todling, protex-compliant prologue
1299 !   2004-07-23  parrish - new routine
1300 !
1301 ! !REMARKS:
1302 !   language: f90
1303 !   machine:  ibm rs/6000 sp; SGI Origin 2000; Compaq/HP
1304 !
1305 ! !AUTHOR:
1306 !   parrish          org: np22                date: 2003-08-28
1307 !
1308 !EOP
1309 !-------------------------------------------------------------------------
1310
1311     real(r_kind) clon,slon,r_of_lat,xtilde,ytilde
1312     real(r_kind) dtilde,etilde
1313     real(r_kind) d1tilde,d2tilde,e1tilde,e2tilde,detinv
1314     integer(i_kind) itilde,jtilde
1315     integer(i_kind) i0,j0,ip,jp
1316
1317 !   first compute xtilde, ytilde
1318
1319     clon=cos(rlon+rlambda0)
1320     slon=sin(rlon+rlambda0)
1321     r_of_lat=pihalf+sign_pole*rlat
1322
1323     xtilde=atilde_x+btilde_x*r_of_lat*clon
1324     ytilde=atilde_y+btilde_y*r_of_lat*slon
1325
1326 !  next get interpolation information
1327
1328     itilde=max(1,min(nint(xtilde),nxtilde))
1329     jtilde=max(1,min(nint(ytilde),nytilde))
1330
1331     i0=i0_tilde(itilde,jtilde)
1332     j0=j0_tilde(itilde,jtilde)
1333     ip=i0+ip_tilde(itilde,jtilde)
1334     jp=j0+jp_tilde(itilde,jtilde)
1335     dtilde=xtilde-xtilde0(i0,j0)
1336     etilde=ytilde-ytilde0(i0,j0)
1337     d1tilde=(xtilde0(ip,j0)-xtilde0(i0,j0))*(ip-i0)
1338     d2tilde=(xtilde0(i0,jp)-xtilde0(i0,j0))*(jp-j0)
1339     e1tilde=(ytilde0(ip,j0)-ytilde0(i0,j0))*(ip-i0)
1340     e2tilde=(ytilde0(i0,jp)-ytilde0(i0,j0))*(jp-j0)
1341     detinv=one/(d1tilde*e2tilde-d2tilde*e1tilde)
1342     x = i0+detinv*(e2tilde*dtilde-d2tilde*etilde)
1343     y = j0+detinv*(d1tilde*etilde-e1tilde*dtilde)
1344
1345     outside=x < rlon_min_dd .or. x > rlon_max_dd .or. &
1346             y < rlat_min_dd .or. y > rlat_max_dd
1347
1348  end subroutine tll2xy
1349
1350 !-------------------------------------------------------------------------
1351 !    NOAA/NCEP, National Centers for Environmental Prediction GSI        !
1352 !-------------------------------------------------------------------------
1353 !BOP
1354 !
1355 ! !IROUTINE:  txy2ll ---  convert x-y grid units to earth lat-lon coordinates
1356 !
1357 ! !INTERFACE:
1358 !
1359   subroutine txy2ll(x,y,rlon,rlat)
1360
1361 ! !USES:
1362
1363     use kinds, only: r_kind,i_kind
1364     use constants, only: one
1365     implicit none
1366
1367 ! !INPUT PARAMETERS:
1368
1369     real(r_kind),intent(in):: x      ! x-grid coordinate (grid units)
1370     real(r_kind),intent(in):: y      ! y_grid coordinate (grid units)
1371
1372 ! !OUTPUT PARAMETERS:
1373
1374     real(r_kind),intent(out)::rlon   ! earth longitude (radians)
1375     real(r_kind),intent(out)::rlat   ! earth latitude  (radians)
1376
1377 ! !DESCRIPTION: to convert earth lon-lat to x-y grid units of a
1378 !           general regional rectangular domain.  Also, decide if
1379 !           point is inside this domain.  As a result, there is
1380 !           no restriction on type of horizontal coordinate for
1381 !           a regional run, other than that it not have periodicity
1382 !           or polar singularities.
1383 !           This is done by first converting rlon, rlat to an
1384 !           intermediate coordinate xtilde,ytilde, which has
1385 !           precomputed pointers and constants for final conversion
1386 !           to the desired x,y via 3 point inverse interpolation.
1387 !           All of the information needed is derived from arrays
1388 !           specifying earth latitude and longitude of every point
1389 !           on the input grid.  Currently, the input x-y grid that
1390 !           this is based on must be non-staggered.  This restriction
1391 !           will eventually be lifted so we can run directly from
1392 !           model grids that are staggered without first resorting
1393 !           to interpolation of the guess to a non-staggered grid.
1394 !
1395 ! !REVISION HISTORY:
1396 !   2003-08-28  parrish
1397 !   2004-05-13  kleist, documentation
1398 !   2004-07-15  todling, protex-compliant prologue
1399 !   2004-07-20  todling, fixed description
1400 !   2004-07-23  parrish - new routine
1401 !
1402 ! !REMARKS:
1403 !   language: f90
1404 !   machine:  ibm rs/6000 sp; SGI Origin 2000; Compaq/HP
1405 !
1406 ! !AUTHOR:
1407 !   parrish          org: np22                date: 2003-08-28
1408 !
1409 !EOP
1410 !-------------------------------------------------------------------------
1411
1412     real(r_kind) clon,slon,r_of_lat,xtilde,ytilde
1413     real(r_kind) dtilde,etilde,xbar,ybar
1414     real(r_kind) d1tilde,d2tilde,e1tilde,e2tilde
1415     integer(i_kind) itilde,jtilde
1416     integer(i_kind) i0,j0,ip,jp
1417
1418     i0=nint(x)
1419     j0=nint(y)
1420     i0=max(1,min(i0,nlon))
1421     j0=max(1,min(j0,nlat))
1422     ip=i0+nint(sign(one,x-i0))
1423     jp=j0+nint(sign(one,y-j0))
1424     if(ip.lt.1) then
1425      i0=2
1426      ip=1
1427     end if
1428     if(jp.lt.1) then
1429      j0=2
1430      jp=1
1431     end if
1432     if(ip.gt.nlon) then
1433      i0=nlon-1
1434      ip=nlon
1435     end if
1436     if(jp.gt.nlat) then
1437      j0=nlat-1
1438      jp=nlat
1439     end if
1440     d1tilde=(xtilde0(ip,j0)-xtilde0(i0,j0))*(ip-i0)
1441     d2tilde=(xtilde0(i0,jp)-xtilde0(i0,j0))*(jp-j0)
1442     e1tilde=(ytilde0(ip,j0)-ytilde0(i0,j0))*(ip-i0)
1443     e2tilde=(ytilde0(i0,jp)-ytilde0(i0,j0))*(jp-j0)
1444     dtilde=d1tilde*(x-i0)+d2tilde*(y-j0)
1445     etilde=e1tilde*(x-i0)+e2tilde*(y-j0)
1446     xtilde=dtilde+xtilde0(i0,j0)
1447     ytilde=etilde+ytilde0(i0,j0)
1448
1449     xbar=(xtilde-atilde_x)*btilde_xinv
1450     ybar=(ytilde-atilde_y)*btilde_yinv
1451     r_of_lat=sqrt(xbar**2+ybar**2)
1452     rlat=(r_of_lat-pihalf)*sign_pole
1453     rlon=atan2(ybar,xbar)-rlambda0
1454
1455  end subroutine txy2ll
1456
1457  subroutine nearest_3(ilast,jlast,i0,j0,ip,jp,x,y,nx0,ny0,x0,y0)
1458
1459 !   find closest 3 points to (x,y) on grid defined by x0,y0
1460
1461   implicit none
1462   integer(i_kind),intent(inout)::ilast,jlast
1463   integer(i_kind),intent(out)::i0,j0
1464   integer(i_byte),intent(out)::ip,jp
1465   integer(i_kind),intent(in)::nx0,ny0
1466   real(r_kind),intent(in)::x,y
1467   real(r_kind),intent(in)::x0(nx0,ny0),y0(nx0,ny0)
1468
1469   real(r_kind) dista,distb,dist2,dist2min
1470   integer(i_kind) i,i1mirror,inext,j,j1mirror,jnext
1471
1472   do
1473    i0=ilast
1474    j0=jlast
1475    dist2min=huge(dist2min)
1476    inext=0
1477    jnext=0
1478    do j=max(j0-1,1),min(j0+1,ny0)
1479     do i=max(i0-1,1),min(i0+1,nx0)
1480      dist2=(x-x0(i,j))**2+(y-y0(i,j))**2
1481      if(dist2.lt.dist2min) then
1482       dist2min=dist2
1483       inext=i
1484       jnext=j
1485      end if
1486     end do
1487    end do
1488    if(inext.eq.i0.and.jnext.eq.j0) exit
1489    ilast=inext
1490    jlast=jnext
1491   end do
1492
1493 !  now find which way to go in x for second point
1494
1495   ip=0
1496   if(i0.eq.nx0) ip=-1
1497   if(i0.eq.1) ip=1
1498   if(ip.eq.0) then
1499    dista=(x-x0(i0-1,j0))**2+(y-y0(i0-1,j0))**2
1500    distb=(x-x0(i0+1,j0))**2+(y-y0(i0+1,j0))**2
1501    if(distb.lt.dista) then
1502     ip=1
1503    else
1504     ip=-1
1505    end if
1506   end if
1507
1508 !  repeat for y for 3rd point
1509
1510   jp=0
1511   if(j0.eq.ny0) jp=-1
1512   if(j0.eq.1) jp=1
1513   if(jp.eq.0) then
1514    dista=(x-x0(i0,j0-1))**2+(y-y0(i0,j0-1))**2
1515    distb=(x-x0(i0,j0+1))**2+(y-y0(i0,j0+1))**2
1516    if(distb.lt.dista) then
1517     jp=1
1518    else
1519     jp=-1
1520    end if
1521   end if
1522
1523   ilast=i0
1524   jlast=j0
1525
1526  end subroutine nearest_3
1527
1528  subroutine get_xytilde_domain(nx0,ny0,rlons0,rlats0, &
1529                                   nx,ny,xminout,xmaxout,yminout,ymaxout)
1530
1531    use constants, only: one, deg2rad,half,zero
1532 !  define parameters for xy domain which optimally overlays input grid
1533
1534   implicit none
1535   integer(i_kind),intent(in)::nx0,ny0
1536   real(r_kind),intent(in)::rlons0(nx0,ny0),rlats0(nx0,ny0)
1537
1538   integer(i_kind),intent(out)::nx,ny
1539   real(r_kind),intent(out)::xminout,xmaxout,yminout,ymaxout
1540
1541   real(r_kind),parameter:: r10=10.0_r_kind
1542   real(r_kind),parameter:: r37=37.0_r_kind
1543
1544   real(r_kind) area,areamax,areamin,extra,rlats0max,rlats0min,testlambda
1545   real(r_kind) xthis,ythis
1546   integer(i_kind) i,ip1,j,jp1,m
1547
1548   real(r_kind) coslon0(nx0,ny0),sinlon0(nx0,ny0)
1549   real(r_kind) coslat0(nx0,ny0),sinlat0(nx0,ny0)
1550   real(r_kind) count,delbar
1551   real(r_kind) dx,dy,disti,distj,distmin,distmax
1552   real(r_kind) xmin,xmax,ymin,ymax
1553
1554 !  get range of lats for input grid
1555
1556   rlats0max=maxval(rlats0) ; rlats0min=minval(rlats0)
1557
1558 !   assign hemisphere ( parameter sign_pole )
1559
1560   if(rlats0min.gt.-r37*deg2rad) sign_pole=-one   !  northern hemisphere xy domain
1561   if(rlats0max.lt. r37*deg2rad) sign_pole= one   !  southern hemisphere xy domain
1562
1563 !   get optimum rotation angle rlambda0
1564
1565   areamin=huge(areamin)
1566   areamax=-huge(areamax)
1567   do m=0,359
1568    testlambda=m*deg2rad
1569    xmax=-huge(xmax)
1570    xmin=huge(xmin)
1571    ymax=-huge(ymax)
1572    ymin=huge(ymin)
1573    do j=1,ny0,ny0-1
1574     do i=1,nx0
1575      xthis=(pihalf+sign_pole*rlats0(i,j))*cos(rlons0(i,j)+testlambda)
1576      ythis=(pihalf+sign_pole*rlats0(i,j))*sin(rlons0(i,j)+testlambda)
1577      xmax=max(xmax,xthis)
1578      ymax=max(ymax,ythis)
1579      xmin=min(xmin,xthis)
1580      ymin=min(ymin,ythis)
1581     end do
1582    end do
1583    do j=1,ny0
1584     do i=1,nx0,nx0-1
1585      xthis=(pihalf+sign_pole*rlats0(i,j))*cos(rlons0(i,j)+testlambda)
1586      ythis=(pihalf+sign_pole*rlats0(i,j))*sin(rlons0(i,j)+testlambda)
1587      xmax=max(xmax,xthis)
1588      ymax=max(ymax,ythis)
1589      xmin=min(xmin,xthis)
1590      ymin=min(ymin,ythis)
1591     end do
1592    end do
1593    area=(xmax-xmin)*(ymax-ymin)
1594    areamax=max(area,areamax)
1595    if(area.lt.areamin) then
1596     areamin=area
1597     rlambda0=testlambda
1598     xmaxout=xmax
1599     xminout=xmin
1600     ymaxout=ymax
1601     yminout=ymin
1602    end if
1603   end do
1604
1605 !   now determine resolution of input grid and choose nx,ny of xy grid accordingly
1606 !                 (currently hard-wired at 1/2 the average input grid increment)
1607
1608   do j=1,ny0
1609    do i=1,nx0
1610     coslon0(i,j)=cos(one*rlons0(i,j)) ; sinlon0(i,j)=sin(one*rlons0(i,j))
1611     coslat0(i,j)=cos(one*rlats0(i,j)) ; sinlat0(i,j)=sin(one*rlats0(i,j))
1612    end do
1613   end do
1614
1615   delbar=zero
1616   count=zero
1617   do j=1,ny0-1
1618    jp1=j+1
1619    do i=1,nx0-1
1620     ip1=i+1
1621     disti=acos(sinlat0(i,j)*sinlat0(ip1,j)+coslat0(i,j)*coslat0(ip1,j)* &
1622                (sinlon0(i,j)*sinlon0(ip1,j)+coslon0(i,j)*coslon0(ip1,j)))
1623     distj=acos(sinlat0(i,j)*sinlat0(i,jp1)+coslat0(i,j)*coslat0(i,jp1)* &
1624                (sinlon0(i,j)*sinlon0(i,jp1)+coslon0(i,j)*coslon0(i,jp1)))
1625     distmax=max(disti,distj)
1626     distmin=min(disti,distj)
1627     delbar=delbar+distmax
1628     count=count+one
1629    end do
1630   end do
1631   delbar=delbar/count
1632   dx=half*delbar
1633   dy=dx
1634
1635 !   add extra space to computational grid to push any boundary problems away from
1636 !     area of interest
1637
1638   extra=r10*dx
1639   xmaxout=xmaxout+extra
1640   xminout=xminout-extra
1641   ymaxout=ymaxout+extra
1642   yminout=yminout-extra
1643   nx=1+(xmaxout-xminout)/dx
1644   ny=1+(ymaxout-yminout)/dy
1645
1646  end subroutine get_xytilde_domain
1647
1648  subroutine half_nmm_grid2a(gin,nx,ny,gout,igtype)
1649
1650 !$$$  subprogram documentation block
1651 !                .      .    .                                       .
1652 ! subprogram:    half_nmm_grid2a same as half_nmm_grid2, but output not reorganized
1653 !   prgmmr: parrish         org: w/emc2               date: 2004-06-22
1654 !
1655 ! abstract: creates an unstaggered A grid from the staggered E grid used by the wrf nmm.
1656 !           This is done by keeping every other row of the original E grid.  If this
1657 !           is a mass variable (igtype=1), then no interpolation is required.  If this
1658 !           is a wind variable (igtype=2), then interpolation is necessary.  This procedure
1659 !           is necessary because the gsi is not yet able to work with anything other than
1660 !           unstaggered grids.  This solution introduces greater interpolation error
1661 !           compared to the option fill_nmm_grid2, but has the advantage of 4 times fewer
1662 !           grid points compared to the output of fill_nmm__grid2.  This routine will be
1663 !           eliminated when the gsi has the capability to work directly with staggered grids.
1664 !
1665 ! program history log:
1666 !   2004-06-22  parrish, document
1667 !   2005-03-03  treadon - add implicit none
1668 !
1669 !   input argument list:
1670 !     gin      - input staggered E grid field over entire horizontal domain
1671 !     nx,ny    - input grid dimensions
1672 !     igtype   - =1, then (1,1) on staggered grid is at corner of grid (mass point for nmm)
1673 !              - =2, then (1,1) is staggered (wind point for nmm, see illustration below)
1674
1675 !                   igtype=1:
1676 !
1677 !
1678 !
1679 !       ^   3             x     x     x     x
1680 !       |
1681 !       y   2                x     x     x     x
1682 !
1683 !           1             x     x     x     x
1684 !
1685 !                         1     2     3
1686 !
1687 !                           x -->
1688
1689 !                   igtype=2:
1690 !
1691 !
1692 !
1693 !       ^   3             x     x     x     x
1694 !       |
1695 !       y   2          x     x     x     x
1696 !
1697 !           1             x     x     x     x
1698 !
1699 !                         1     2     3
1700 !
1701 !                           x -->
1702 !
1703 !   output argument list
1704 !     gout     - output unstaggered half grid  (reorganized for distibution to local domains)
1705 !
1706 ! attributes:
1707 !   language: f90
1708 !   machine:  ibm RS/6000 SP
1709 !
1710 !$$$
1711   use constants, only: quarter
1712   implicit none
1713   integer(i_kind) nx,ny,igtype
1714   real(r_single) gin(nx,ny)
1715   real(r_kind) gout(nx,*)
1716
1717   integer(i_kind) i,i0,im,j,jj,jm,jp
1718
1719   if(igtype.eq.1) then
1720    jj=0
1721    do j=1,ny,2
1722     jj=jj+1
1723     do i=1,nx
1724      gout(i,jj)=gin(i,j)
1725     end do
1726    end do
1727   else
1728    jj=0
1729    do j=1,ny,2
1730     jj=jj+1
1731     jp=j+1 ; if(jp.gt.ny) jp=j-1
1732     jm=j-1 ; if(jm.lt.1) jm=j+1
1733     do i=1,nx
1734      im=i-1 ; if(im.lt.1) im=i
1735      i0=i ; if(i.eq.nx) i0=im
1736      gout(i,jj)=quarter*(gin(im,j)+gin(i0,j)+gin(i,jp)+gin(i,jm))
1737     end do
1738    end do
1739   end if
1740
1741  end subroutine half_nmm_grid2a
1742
1743  subroutine fill_nmm_grid2a3(gin,nx,ny,gout)
1744
1745   implicit none
1746   integer(i_kind) nx,ny
1747   real(r_kind) gin(nx,ny)
1748   real(r_kind) gout(2*nx-1,ny)
1749
1750   integer(i_kind) i,j
1751   integer(i_kind) i1a(2*nx-1),i2a(2*nx-1)
1752   integer(i_kind) i3a(2*nx-1),i4a(2*nx-1)
1753   real(r_kind) r1a(2*nx-1),r2a(2*nx-1)
1754   real(r_kind) r3a(2*nx-1),r4a(2*nx-1)
1755   real(r_kind) x,x1,x2,x3,x4
1756
1757 !  first transfer all staggered points to appropriate
1758 !   points on filled output grid
1759
1760   do j=1,ny,2
1761    do i=1,nx
1762     gout(2*i-1,j)=gin(i,j)
1763    end do
1764   end do
1765   do j=2,ny,2
1766    do i=1,nx-1
1767     gout(2*i,j)=gin(i,j)
1768    end do
1769   end do
1770
1771 !   compute all interpolation constants for even x points on odd y rows
1772
1773   i=2
1774    i1a(i)=i-1 ; i2a(i)=i+1 ; i3a(i)=i+3 ; i4a(i)=i+5
1775    x=i ; x1=i1a(i) ; x2=i2a(i) ; x3=i3a(i) ; x4=i4a(i)
1776    r1a(i)=       (x-x2)*(x-x3)*(x-x4)/(        (x1-x2)*(x1-x3)*(x1-x4))
1777    r2a(i)=(x-x1)       *(x-x3)*(x-x4)/((x2-x1)        *(x2-x3)*(x2-x4))
1778    r3a(i)=(x-x1)*(x-x2)       *(x-x4)/((x3-x1)*(x3-x2)        *(x3-x4))
1779    r4a(i)=(x-x1)*(x-x2)*(x-x3)       /((x4-x1)*(x4-x2)*(x4-x3)        )
1780
1781   do i=4,2*nx-4,2
1782    i1a(i)=i-3 ; i2a(i)=i-1 ; i3a(i)=i+1 ; i4a(i)=i+3
1783    x=i ; x1=i1a(i) ; x2=i2a(i) ; x3=i3a(i) ; x4=i4a(i)
1784    r1a(i)=       (x-x2)*(x-x3)*(x-x4)/(        (x1-x2)*(x1-x3)*(x1-x4))
1785    r2a(i)=(x-x1)       *(x-x3)*(x-x4)/((x2-x1)        *(x2-x3)*(x2-x4))
1786    r3a(i)=(x-x1)*(x-x2)       *(x-x4)/((x3-x1)*(x3-x2)        *(x3-x4))
1787    r4a(i)=(x-x1)*(x-x2)*(x-x3)       /((x4-x1)*(x4-x2)*(x4-x3)        )
1788   end do
1789
1790   i=2*nx-2
1791    i1a(i)=i-5 ; i2a(i)=i-3 ; i3a(i)=i-1 ; i4a(i)=i+1
1792    x=i ; x1=i1a(i) ; x2=i2a(i) ; x3=i3a(i) ; x4=i4a(i)
1793    r1a(i)=       (x-x2)*(x-x3)*(x-x4)/(        (x1-x2)*(x1-x3)*(x1-x4))
1794    r2a(i)=(x-x1)       *(x-x3)*(x-x4)/((x2-x1)        *(x2-x3)*(x2-x4))
1795    r3a(i)=(x-x1)*(x-x2)       *(x-x4)/((x3-x1)*(x3-x2)        *(x3-x4))
1796    r4a(i)=(x-x1)*(x-x2)*(x-x3)       /((x4-x1)*(x4-x2)*(x4-x3)        )
1797
1798 !   now get all interpolation constants for odd x points on even y rows
1799
1800   i=1
1801    i1a(i)=i+1 ; i2a(i)=i+3 ; i3a(i)=i+5 ; i4a(i)=i+7
1802    x=i ; x1=i1a(i) ; x2=i2a(i) ; x3=i3a(i) ; x4=i4a(i)
1803    r1a(i)=       (x-x2)*(x-x3)*(x-x4)/(        (x1-x2)*(x1-x3)*(x1-x4))
1804    r2a(i)=(x-x1)       *(x-x3)*(x-x4)/((x2-x1)        *(x2-x3)*(x2-x4))
1805    r3a(i)=(x-x1)*(x-x2)       *(x-x4)/((x3-x1)*(x3-x2)        *(x3-x4))
1806    r4a(i)=(x-x1)*(x-x2)*(x-x3)       /((x4-x1)*(x4-x2)*(x4-x3)        )
1807
1808   i=3
1809    i1a(i)=i-1 ; i2a(i)=i+1 ; i3a(i)=i+3 ; i4a(i)=i+5
1810    x=i ; x1=i1a(i) ; x2=i2a(i) ; x3=i3a(i) ; x4=i4a(i)
1811    r1a(i)=       (x-x2)*(x-x3)*(x-x4)/(        (x1-x2)*(x1-x3)*(x1-x4))
1812    r2a(i)=(x-x1)       *(x-x3)*(x-x4)/((x2-x1)        *(x2-x3)*(x2-x4))
1813    r3a(i)=(x-x1)*(x-x2)       *(x-x4)/((x3-x1)*(x3-x2)        *(x3-x4))
1814    r4a(i)=(x-x1)*(x-x2)*(x-x3)       /((x4-x1)*(x4-x2)*(x4-x3)        )
1815
1816   do i=5,2*nx-5,2
1817    i1a(i)=i-3 ; i2a(i)=i-1 ; i3a(i)=i+1 ; i4a(i)=i+3
1818    x=i ; x1=i1a(i) ; x2=i2a(i) ; x3=i3a(i) ; x4=i4a(i)
1819    r1a(i)=       (x-x2)*(x-x3)*(x-x4)/(        (x1-x2)*(x1-x3)*(x1-x4))
1820    r2a(i)=(x-x1)       *(x-x3)*(x-x4)/((x2-x1)        *(x2-x3)*(x2-x4))
1821    r3a(i)=(x-x1)*(x-x2)       *(x-x4)/((x3-x1)*(x3-x2)        *(x3-x4))
1822    r4a(i)=(x-x1)*(x-x2)*(x-x3)       /((x4-x1)*(x4-x2)*(x4-x3)        )
1823   end do
1824
1825   i=2*nx-3
1826    i1a(i)=i-5 ; i2a(i)=i-3 ; i3a(i)=i-1 ; i4a(i)=i+1
1827    x=i ; x1=i1a(i) ; x2=i2a(i) ; x3=i3a(i) ; x4=i4a(i)
1828    r1a(i)=       (x-x2)*(x-x3)*(x-x4)/(        (x1-x2)*(x1-x3)*(x1-x4))
1829    r2a(i)=(x-x1)       *(x-x3)*(x-x4)/((x2-x1)        *(x2-x3)*(x2-x4))
1830    r3a(i)=(x-x1)*(x-x2)       *(x-x4)/((x3-x1)*(x3-x2)        *(x3-x4))
1831    r4a(i)=(x-x1)*(x-x2)*(x-x3)       /((x4-x1)*(x4-x2)*(x4-x3)        )
1832
1833   i=2*nx-1
1834    i1a(i)=i-7 ; i2a(i)=i-5 ; i3a(i)=i-3 ; i4a(i)=i-1
1835    x=i ; x1=i1a(i) ; x2=i2a(i) ; x3=i3a(i) ; x4=i4a(i)
1836    r1a(i)=       (x-x2)*(x-x3)*(x-x4)/(        (x1-x2)*(x1-x3)*(x1-x4))
1837    r2a(i)=(x-x1)       *(x-x3)*(x-x4)/((x2-x1)        *(x2-x3)*(x2-x4))
1838    r3a(i)=(x-x1)*(x-x2)       *(x-x4)/((x3-x1)*(x3-x2)        *(x3-x4))
1839    r4a(i)=(x-x1)*(x-x2)*(x-x3)       /((x4-x1)*(x4-x2)*(x4-x3)        )
1840
1841   do j=1,ny,2
1842    do i=2,2*nx-2,2
1843     gout(i,j)=r1a(i)*gout(i1a(i),j)+r2a(i)*gout(i2a(i),j)+ &
1844               r3a(i)*gout(i3a(i),j)+r4a(i)*gout(i4a(i),j)
1845    end do
1846   end do
1847   do j=2,ny,2
1848    do i=1,2*nx-1,2
1849     gout(i,j)=r1a(i)*gout(i1a(i),j)+r2a(i)*gout(i2a(i),j)+ &
1850               r3a(i)*gout(i3a(i),j)+r4a(i)*gout(i4a(i),j)
1851    end do
1852   end do
1853
1854  end subroutine fill_nmm_grid2a3
1855
1856 !-------------------------------------------------------------------------
1857 !    NOAA/NCEP, National Centers for Environmental Prediction GSI        !
1858 !-------------------------------------------------------------------------
1859 !BOP
1860 !
1861 ! !IROUTINE:  rotate_wind_ll2xy ---  Rotate earth vector wind
1862 !
1863 ! !INTERFACE:
1864 !
1865   subroutine rotate_wind_ll2xy(u0,v0,u,v,rlon0,rlat0,x,y)
1866
1867 ! !USES:
1868
1869     use kinds, only: r_kind,i_kind
1870     use constants, only: one,two,pi,rad2deg
1871     implicit none
1872
1873 ! !INPUT PARAMETERS:
1874
1875     real(r_kind),intent(in)::u0,v0        ! earth wind component
1876     real(r_kind),intent(in)::rlon0,rlat0  ! earth   lon/lat (radians)
1877     real(r_kind),intent(in)::x,y          ! local x,y coordinate (grid units)
1878
1879 ! !OUTPUT PARAMETERS:
1880
1881     real(r_kind),intent(out)::u,v         ! rotated coordinate of winds
1882
1883 ! !DESCRIPTION: to convert earth vector wind components to corresponding
1884 !           local x,y coordinate
1885 !
1886 ! !REVISION HISTORY:
1887 !   2003-09-30  parrish
1888 !   2004-05-13  kleist, documentation
1889 !   2004-07-15  todling, protex-compliant prologue
1890 !
1891 ! !REMARKS:
1892 !   language: f90
1893 !   machine:  ibm rs/6000 sp; SGI Origin 2000; Compaq/HP
1894 !
1895 ! !AUTHOR:
1896 !   parrish          org: np22                date: 2003-09-30
1897 !
1898 !EOP
1899 !-------------------------------------------------------------------------
1900
1901   real(r_kind) beta,delx,delxp,dely,delyp
1902   real(r_kind) beta_old,two_pi,sin_beta,cos_beta,thisdiff
1903   integer(i_kind) ix,iy,k
1904
1905 !  interpolate departure from longitude part of angle between earth positive east and local positive x
1906
1907   ix=x
1908   iy=y
1909   ix=max(1,min(ix,nlon-1))
1910   iy=max(1,min(iy,nlat-1))
1911   delx=x-ix
1912   dely=y-iy
1913   delxp=one-delx
1914   delyp=one-dely
1915   beta_old=beta_ref(ix  ,iy  )*delxp*delyp+beta_ref(ix+1,iy  )*delx *delyp+ &
1916        beta_ref(ix  ,iy+1)*delxp*dely +beta_ref(ix+1,iy+1)*delx *dely
1917   cos_beta=cos_beta_ref(ix  ,iy  )*delxp*delyp+cos_beta_ref(ix+1,iy  )*delx *delyp+ &
1918        cos_beta_ref(ix  ,iy+1)*delxp*dely +cos_beta_ref(ix+1,iy+1)*delx *dely
1919   sin_beta=sin_beta_ref(ix  ,iy  )*delxp*delyp+sin_beta_ref(ix+1,iy  )*delx *delyp+ &
1920        sin_beta_ref(ix  ,iy+1)*delxp*dely +sin_beta_ref(ix+1,iy+1)*delx *dely
1921   beta=atan2(sin_beta,cos_beta)
1922   thisdiff=huge(thisdiff)
1923   two_pi=two*pi
1924   do k=-6,6
1925     thisdiff=min(abs(beta-beta_old+k*two_pi),thisdiff)
1926   end do
1927   if(thisdiff*rad2deg.gt.0.1_r_kind) then
1928     count_beta_diff_gt_20=count_beta_diff_gt_20 + one
1929     beta_diff_max_gt_20=max(beta_diff_max_gt_20,thisdiff)
1930   else
1931     beta_diff_max=max(beta_diff_max,thisdiff)
1932     beta_diff_min=min(beta_diff_min,thisdiff)
1933     beta_diff_rms=beta_diff_rms+thisdiff**2
1934     count_beta_diff=count_beta_diff+one
1935   end if
1936
1937 !  now rotate;
1938
1939   u= u0*cos(beta+rlon0*sign_pole)+v0*sin(beta+rlon0*sign_pole)
1940   v=-u0*sin(beta+rlon0*sign_pole)+v0*cos(beta+rlon0*sign_pole)
1941
1942  end subroutine rotate_wind_ll2xy
1943
1944 !-------------------------------------------------------------------------
1945 !    NOAA/NCEP, National Centers for Environmental Prediction GSI        !
1946 !-------------------------------------------------------------------------
1947 !BOP
1948 !
1949 ! !IROUTINE:  rotate_wind_xy2ll ---  Unrotate earth vector wind
1950 !
1951 ! !INTERFACE:
1952 !
1953   subroutine rotate_wind_xy2ll(u,v,u0,v0,rlon0,rlat0,x,y)
1954
1955 ! !USES:
1956
1957     use kinds, only: r_kind,i_kind
1958     use constants, only: one
1959     implicit none
1960
1961 ! !INPUT PARAMETERS:
1962
1963     real(r_kind),intent(in)::u,v         ! rotated coordinate winds
1964     real(r_kind),intent(in)::rlon0,rlat0 ! earth   lon/lat (radians)
1965     real(r_kind),intent(in)::x,y         ! rotated lon/lat (radians)
1966
1967 ! !OUTPUT PARAMETERS:
1968
1969     real(r_kind),intent(out)::u0,v0      ! earth winds
1970
1971 ! !DESCRIPTION: rotate u,v in local x,y coordinate to u0,v0 in earth
1972 !           lat, lon coordinate
1973 !
1974 ! !REVISION HISTORY:
1975 !   2003-09-30  parrish
1976 !   2004-05-13  kleist, documentation
1977 !   2004-07-15  todling, protex-compliant prologue
1978 !   2004-07-20  todling, fixed description
1979 !
1980 ! !REMARKS:
1981 !   language: f90
1982 !   machine:  ibm rs/6000 sp; SGI Origin 2000; Compaq/HP
1983 !
1984 ! !AUTHOR:
1985 !   parrish          org: np22                date: 2003-09-30
1986 !
1987 !EOP
1988 !-------------------------------------------------------------------------
1989   real(r_kind) beta,delx,delxp,dely,delyp
1990   real(r_kind) sin_beta,cos_beta
1991   integer(i_kind) ix,iy
1992
1993 !  interpolate departure from longitude part of angle between earth
1994 !  positive east and local positive x
1995
1996   ix=x
1997   iy=y
1998   ix=max(1,min(ix,nlon-1))
1999   iy=max(1,min(iy,nlat-1))
2000   delx=x-ix
2001   dely=y-iy
2002   delxp=one-delx
2003   delyp=one-dely
2004   cos_beta=cos_beta_ref(ix  ,iy  )*delxp*delyp+cos_beta_ref(ix+1,iy  )*delx *delyp+ &
2005        cos_beta_ref(ix  ,iy+1)*delxp*dely +cos_beta_ref(ix+1,iy+1)*delx *dely
2006   sin_beta=sin_beta_ref(ix  ,iy  )*delxp*delyp+sin_beta_ref(ix+1,iy  )*delx *delyp+ &
2007        sin_beta_ref(ix  ,iy+1)*delxp*dely +sin_beta_ref(ix+1,iy+1)*delx *dely
2008   beta=atan2(sin_beta,cos_beta)
2009
2010 !  now rotate;
2011
2012   u0= u*cos(beta+rlon0*sign_pole)-v*sin(beta+rlon0*sign_pole)
2013   v0= u*sin(beta+rlon0*sign_pole)+v*cos(beta+rlon0*sign_pole)
2014
2015  end subroutine rotate_wind_xy2ll
2016
2017 !-------------------------------------------------------------------------
2018 !    NOAA/NCEP, National Centers for Environmental Prediction GSI        !
2019 !-------------------------------------------------------------------------
2020 !BOP
2021 !
2022 ! !IROUTINE:  load_grid --- strip off south/north latitude rows
2023 !
2024 ! !INTERFACE:
2025 !
2026  subroutine load_grid(grid_in,grid_out)
2027
2028 ! !USES:
2029
2030    use kinds, only: r_kind,i_kind
2031    implicit none
2032
2033 ! !INPUT PARAMETERS:
2034
2035    real(r_kind),dimension(max(iglobal,itotsub)),intent(in):: grid_in  ! input grid
2036    real(r_kind),dimension(nlon,nlat-2),intent(out):: grid_out         ! output grid
2037
2038 ! !DESCRIPTION: This routine prepares grids for use in splib
2039 !               grid to spectral tranforms.  This preparation
2040 !               entails to two steps
2041 !                  1) reorder indexing of the latitude direction.
2042 !                     The GSI ordering is south to north.  The
2043 !                     ordering assumed in splib routines is north
2044 !                     to south.
2045 !                  2) The global GSI adds two latitude rows, one
2046 !                     for each pole.  These latitude rows are not
2047 !                     needed in the grid to spectral transforms of
2048 !                     splib.  The code below strips off these
2049 !                     "pole rows"
2050 !
2051 ! !REVISION HISTORY:
2052 !   2004-08-27  treadon
2053 !
2054 ! !REMARKS:
2055 !   language: f90
2056 !   machine:  ibm rs/6000
2057 !
2058 ! !AUTHOR:
2059 !   treadon          org: np23                date: 2004-08-27
2060 !
2061 !EOP
2062 !-------------------------------------------------------------------------
2063    integer(i_kind) i,j,k,nlatm1,jj
2064    real(r_kind),dimension(nlon,nlat):: grid
2065
2066 !  Transfer input grid from 1d to 2d local array.  As loading
2067 !  local array, reverse direction of latitude index.  Coming
2068 !  into the routine the order is south --> north.  On exit
2069 !  the order is north --> south
2070    do k=1,iglobal
2071       i=nlat-ltosi(k)+1
2072       j=ltosj(k)
2073       grid(j,i)=grid_in(k)
2074    end do
2075
2076 !  Transfer contents of local array to output array.
2077    nlatm1=nlat-1
2078    do j=2,nlatm1
2079       jj=j-1
2080       do i=1,nlon
2081          grid_out(i,jj)=grid(i,j)
2082       end do
2083    end do
2084
2085    return
2086  end subroutine load_grid
2087
2088 !-------------------------------------------------------------------------
2089 !    NOAA/NCEP, National Centers for Environmental Prediction GSI        !
2090 !-------------------------------------------------------------------------
2091 !BOP
2092 !
2093 ! !IROUTINE:  fill_ns --- add southern/northern latitude rows
2094 !
2095 ! !INTERFACE:
2096 !
2097  subroutine fill_ns(grid_in,grid_out)
2098
2099 ! !USES:
2100
2101    use kinds, only: r_kind,i_kind
2102    use constants, only: zero,one
2103    implicit none
2104
2105 ! !INPUT PARAMETERS:
2106
2107    real(r_kind),dimension(nlon,nlat-2),intent(in):: grid_in  ! input grid
2108    real(r_kind),dimension(itotsub),intent(out):: grid_out    ! output grid
2109
2110 ! !DESCRIPTION: This routine adds a southern and northern latitude
2111 !               row to the input grid.  The southern row contains
2112 !               the longitudinal mean of the adjacent latitude row.
2113 !               The northern row contains the longitudinal mean of
2114 !               the adjacent northern row.
2115 !
2116 !               The added rows correpsond to the south and north poles.
2117 !
2118 !               In addition to adding latitude rows corresponding to the
2119 !               south and north poles, the routine reorder the output
2120 !               array so that it is a one-dimensional array read in
2121 !               an order consisten with that assumed for total domain
2122 !               gsi grids.
2123 !
2124 !               The assumed order for the input grid is longitude as
2125 !               the first dimension with array index increasing from
2126 !               east to west.  The second dimension is latitude with
2127 !               the index increasing from north to south.  This ordering
2128 !               differs from that used in the GSI.
2129 !
2130 !               The GSI ordering is latitude first with the index
2131 !               increasing from south to north.  The second dimension is
2132 !               longitude with the index increasing from east to west.
2133 !
2134 !               Thus, the code below also rearranges the indexing and
2135 !               order of the dimensions to make the output grid
2136 !               consistent with that which is expected in the rest of
2137 !               gsi.
2138 !
2139 !
2140 ! !REVISION HISTORY:
2141 !   2004-08-27  treadon
2142 !
2143 ! !REMARKS:
2144 !   language: f90
2145 !   machine:  ibm rs/6000
2146 !
2147 ! !AUTHOR:
2148 !   treadon          org: np23                date: 2004-08-27
2149 !
2150 !EOP
2151 !-------------------------------------------------------------------------
2152 !  Declare local variables
2153    integer(i_kind) i,j,k,jj,nlatm2
2154    real(r_kind) rnlon,sumn,sums
2155    real(r_kind),dimension(nlon,nlat):: grid
2156
2157 !  Transfer contents of input grid to local work array
2158 !  Reverse ordering in j direction from n-->s to s-->n
2159    do j=2,nlat-1
2160       jj=nlat-j
2161       do i=1,nlon
2162          grid(i,j)=grid_in(i,jj)
2163       end do
2164    end do
2165
2166 !  Compute mean along southern and northern latitudes
2167    sumn=zero
2168    sums=zero
2169    nlatm2=nlat-2
2170    do i=1,nlon
2171       sumn=sumn+grid_in(i,1)
2172       sums=sums+grid_in(i,nlatm2)
2173    end do
2174    rnlon=one/float(nlon)
2175    sumn=sumn*rnlon
2176    sums=sums*rnlon
2177
2178 !  Load means into local work array
2179    do i=1,nlon
2180       grid(i,1)   =sums
2181       grid(i,nlat)=sumn
2182    end do
2183
2184 !  Transfer local work array to output grid
2185    do k=1,itotsub
2186       i=ltosi_s(k)
2187       j=ltosj_s(k)
2188       grid_out(k)=grid(j,i)
2189    end do
2190
2191    return
2192  end subroutine fill_ns
2193
2194 !-------------------------------------------------------------------------
2195 !    NOAA/NCEP, National Centers for Environmental Prediction GSI        !
2196 !-------------------------------------------------------------------------
2197 !BOP
2198 !
2199 ! !IROUTINE:  filluv_ns --- add southern/northern latitude rows
2200 !
2201 ! !INTERFACE:
2202 !
2203  subroutine filluv_ns(gridu_in,gridv_in,gridu_out,gridv_out)
2204
2205 ! !USES:
2206
2207    use kinds, only: r_kind,i_kind
2208    use constants, only: zero,one
2209    implicit none
2210
2211 ! !INPUT PARAMETERS:
2212
2213    real(r_kind),dimension(nlon,nlat-2),intent(in):: gridu_in,gridv_in  ! input grid
2214    real(r_kind),dimension(itotsub),intent(out):: gridu_out,gridv_out    ! output grid
2215
2216 ! !DESCRIPTION: This routine adds a southern and northern latitude
2217 !               row to the input grid.  The southern row contains
2218 !               the longitudinal mean of the adjacent latitude row.
2219 !               The northern row contains the longitudinal mean of
2220 !               the adjacent northern row.
2221 !
2222 !               The added rows correpsond to the south and north poles.
2223 !
2224 !               In addition to adding latitude rows corresponding to the
2225 !               south and north poles, the routine reorder the output
2226 !               array so that it is a one-dimensional array read in
2227 !               an order consisten with that assumed for total domain
2228 !               gsi grids.
2229 !
2230 !               The assumed order for the input grid is longitude as
2231 !               the first dimension with array index increasing from
2232 !               east to west.  The second dimension is latitude with
2233 !               the index increasing from north to south.  This ordering
2234 !               differs from that used in the GSI.
2235 !
2236 !               The GSI ordering is latitude first with the index
2237 !               increasing from south to north.  The second dimension is
2238 !               longitude with the index increasing from east to west.
2239 !
2240 !               Thus, the code below also rearranges the indexing and
2241 !               order of the dimensions to make the output grid
2242 !               consistent with that which is expected in the rest of
2243 !               gsi.
2244 !
2245 !
2246 ! !REVISION HISTORY:
2247 !   2004-08-27  treadon
2248 !
2249 ! !REMARKS:
2250 !   language: f90
2251 !   machine:  ibm rs/6000
2252 !
2253 ! !AUTHOR:
2254 !   treadon          org: np23                date: 2004-08-27
2255 !
2256 !EOP
2257 !-------------------------------------------------------------------------
2258 !  Declare local variables
2259    integer(i_kind) i,j,k,jj,nlatm2
2260    real(r_kind) rnlon,polnu,polnv,polsu,polsv
2261    real(r_kind),dimension(nlon,nlat):: grid,grid2
2262
2263 !  Transfer contents of input grid to local work array
2264 !  Reverse ordering in j direction from n-->s to s-->n
2265    do j=2,nlat-1
2266       jj=nlat-j
2267       do i=1,nlon
2268          grid(i,j)=gridu_in(i,jj)
2269          grid2(i,j)=gridv_in(i,jj)
2270       end do
2271    end do
2272
2273 !  Compute mean along southern and northern latitudes
2274    polnu=zero
2275    polnv=zero
2276    polsu=zero
2277    polsv=zero
2278    do i=1,nlon
2279       polnu=polnu+grid(i,nlat-1)*coslon(i)-grid2(i,nlat-1)*sinlon(i)
2280       polnv=polnv+grid(i,nlat-1)*sinlon(i)+grid2(i,nlat-1)*coslon(i)
2281       polsu=polsu+grid(i,2)*coslon(i)+grid2(i,2)*sinlon(i)
2282       polsv=polsv+grid(i,2)*sinlon(i)-grid2(i,2)*coslon(i)
2283    end do
2284    polnu=polnu/float(nlon)
2285    polnv=polnv/float(nlon)
2286    polsu=polsu/float(nlon)
2287    polsv=polsv/float(nlon)
2288    do i=1,nlon
2289       grid(i,nlat)= polnu*coslon(i)+polnv*sinlon(i)
2290       grid2(i,nlat)=-polnu*sinlon(i)+polnv*coslon(i)
2291       grid(i,1)= polsu*coslon(i)+polsv*sinlon(i)
2292       grid2(i,1)= polsu*sinlon(i)-polsv*coslon(i)
2293    end do
2294
2295 !  Transfer local work array to output grid
2296    do k=1,itotsub
2297       i=ltosi_s(k)
2298       j=ltosj_s(k)
2299       gridu_out(k)=grid(j,i)
2300       gridv_out(k)=grid2(j,i)
2301    end do
2302
2303    return
2304  end subroutine filluv_ns
2305
2306 !-------------------------------------------------------------------------
2307 !    NOAA/NCEP, National Centers for Environmental Prediction GSI        !
2308 !-------------------------------------------------------------------------
2309 !BOP
2310 !
2311 ! !IROUTINE:  get_ij --- get (i,j) grid indices and interpolation weights
2312 !
2313 ! !INTERFACE:
2314 !
2315  subroutine get_ij(mm1,obs_lat,obs_lon,jgrd,wgrd,jjlat,jjlon)
2316
2317 ! !USES:
2318
2319    use kinds, only: r_kind,i_kind
2320    use constants, only: one
2321    implicit none
2322
2323 ! !INPUT PARAMETERS:
2324
2325    integer(i_kind),intent(in)::  mm1
2326    integer(i_kind),dimension(4),intent(out):: jgrd
2327    integer(i_kind),intent(out),optional:: jjlat,jjlon
2328
2329    real(r_kind),intent(in):: obs_lat,obs_lon
2330    real(r_kind),dimension(4),intent(out):: wgrd
2331    real(r_kind):: dx,dy,dx1,dy1
2332
2333    integer(i_kind):: jlat,jlon
2334
2335 ! !DESCRIPTION: This routine returns the sub-domain grid relative
2336 !               i,j index of a given observation (lat,lon).  The
2337 !               routine also returns weights needed for bilinear
2338 !               from the four surrounding analysis grid points to
2339 !               the observation location.
2340 !
2341 ! !REVISION HISTORY:
2342 !   2004-12-23  treadon
2343 !   2006-01-06  treadon - add optional arguments jjlat,jjlon
2344 !
2345 ! !REMARKS:
2346 !   language: f90
2347 !   machine:  ibm rs/6000
2348 !
2349 ! !AUTHOR:
2350 !   treadon          org: np23                date: 2004-08-27
2351 !
2352 !EOP
2353 !-------------------------------------------------------------------------
2354
2355 !  Set (i,j) indices of guess gridpoint that bound obs location
2356    jlat = obs_lat
2357    jlon = obs_lon
2358
2359 !  Compute weights for bilinear interpolation
2360    dy  = obs_lat-jlat
2361    dx  = obs_lon-jlon
2362    dx1 = one-dx
2363    dy1 = one-dy
2364
2365 !  Bound lat and lon indices to fall within analysis grid limits
2366    jlat = min(max(1,jlat),nlat)
2367    jlon = min(max(0,jlon),nlon)
2368
2369 !  Handle special case of e/w periodicity
2370    if (jstart(mm1)==1 .and. jlon==nlon) jlon=0
2371    if (jstart(mm1)+jlon1(mm1)==nlon+1 .and. jlon==0) jlon=nlon
2372
2373 !  Convert global (i,j) indices to sub-domain specific (i,j) indices
2374    jlat=jlat-istart(mm1)+2
2375    jlon=jlon-jstart(mm1)+2
2376
2377    jgrd(1)=jlat+(jlon-1)*lat2
2378    jgrd(2)=jgrd(1)+1
2379    jgrd(3)=jgrd(1)+lat2
2380    jgrd(4)=jgrd(3)+1
2381
2382    wgrd(1)=dx1*dy1
2383    wgrd(2)=dx1*dy
2384    wgrd(3)=dx*dy1
2385    wgrd(4)=dx*dy
2386
2387    if (present(jjlat)) jjlat=jlat
2388    if (present(jjlon)) jjlon=jlon
2389
2390    return
2391  end subroutine get_ij
2392
2393 !-------------------------------------------------------------------------
2394 !    NOAA/NCEP, National Centers for Environmental Prediction GSI        !
2395 !-------------------------------------------------------------------------
2396 !BOP
2397 !
2398 ! !IROUTINE:  get_ijk --- get (i,j,k) grid indices and interpolation weights
2399 !
2400 ! !INTERFACE:
2401 !
2402  subroutine get_ijk(mm1,obs_lat,obs_lon,obs_sig,jgrd,wgrd)
2403
2404 ! !USES:
2405
2406    use kinds, only: r_kind,i_kind
2407    use constants, only: one
2408    implicit none
2409
2410 ! !INPUT PARAMETERS:
2411
2412    integer(i_kind),intent(in)::  mm1
2413    integer(i_kind),dimension(8),intent(out):: jgrd
2414
2415    real(r_kind),intent(in):: obs_lat,obs_lon,obs_sig
2416    real(r_kind),dimension(8),intent(out):: wgrd
2417    real(r_kind) :: dx,dy,dx1,dy1,ds,ds1
2418
2419    integer(i_kind):: jlat,jlon,jsig,latlon11_l
2420
2421 ! !DESCRIPTION: This routine returns the sub-domain grid relative
2422 !               i,j,k index of a given observation (lat,lon,sig).
2423 !               The routine also returns weights needed for bilinear
2424 !               from the eight surrounding analysis grid points to
2425 !               the observation location
2426 !
2427 ! !REVISION HISTORY:
2428 !   2004-12-23  treadon
2429 !
2430 ! !REMARKS:
2431 !   language: f90
2432 !   machine:  ibm rs/6000
2433 !
2434 ! !AUTHOR:
2435 !   treadon          org: np23                date: 2004-08-27
2436 !
2437 !EOP
2438 !-------------------------------------------------------------------------
2439 !  Declare local variables
2440    real(r_kind) obs_s
2441
2442 !  Special handling for vertical coordinate
2443    obs_s = obs_sig
2444    if (obs_s < one) obs_s = one
2445
2446 !  Set (i,j,k) indices of guess gridpoint that bound obs location
2447    jlat = obs_lat
2448    jlon = obs_lon
2449    jsig = obs_s
2450
2451 !  Compute weights for bilinear interpolation
2452    dy  = obs_lat-jlat
2453    dx  = obs_lon-jlon
2454    ds  = obs_s-jsig
2455
2456    dx1 = one-dx
2457    dy1 = one-dy
2458    ds1 = one-ds
2459
2460 !  Bound lat and lon indices to fall within analysis grid limits
2461    jlat = min(max(1,jlat),nlat)
2462    jlon = min(max(0,jlon),nlon)
2463
2464 !  Handle special case of e/w periodicity
2465    if (jstart(mm1)==1 .and. jlon==nlon) jlon=0
2466    if (jstart(mm1)+jlon1(mm1)==nlon+1 .and. jlon==0) jlon=nlon
2467
2468 !  Convert global (i,j) indices to sub-domain specific (i,j) indices
2469    jlat=jlat-istart(mm1)+2
2470    jlon=jlon-jstart(mm1)+2
2471
2472 !  Set number of points on horizontal layer
2473    latlon11_l = latlon11
2474    if(jsig==nsig) latlon11_l=0
2475    jgrd(1)=jlat+(jlon-1)*lat2+(jsig-1)*latlon11
2476    jgrd(2)=jgrd(1)+1
2477    jgrd(3)=jgrd(1)+lat2
2478    jgrd(4)=jgrd(3)+1
2479    jgrd(5)=jgrd(1)+latlon11_l
2480    jgrd(6)=jgrd(5)+1
2481    jgrd(7)=jgrd(5)+lat2
2482    jgrd(8)=jgrd(7)+1
2483
2484    wgrd(1)=dx1*dy1*ds1
2485    wgrd(2)=dx1*dy*ds1
2486    wgrd(3)=dx*dy1*ds1
2487    wgrd(4)=dx*dy*ds1
2488    wgrd(5)=dx1*dy1*ds
2489    wgrd(6)=dx1*dy*ds
2490    wgrd(7)=dx*dy1*ds
2491    wgrd(8)=dx*dy*ds
2492
2493    return
2494  end subroutine get_ijk
2495
2496  subroutine check_rotate_wind(string)
2497    use constants, only: zero,one,rad2deg
2498    implicit none
2499    character(len=*),intent(in):: string
2500
2501    if(count_beta_diff.gt.zero.or.count_beta_diff_gt_20.gt.zero) then
2502       beta_diff_rms=sqrt(beta_diff_rms/(max(one,count_beta_diff)))
2503       write(6,*)'CHECK_ROTATE_WIND:  called from routine ',trim(string)
2504       write(6,100) beta_diff_max*rad2deg, beta_diff_min*rad2deg, beta_diff_rms*rad2deg
2505       write(6,110) count_beta_diff, count_beta_diff_gt_20
2506       write(6,115) beta_diff_max_gt_20*rad2deg
2507 100   format('   beta_diff_mass,min,rms(deg)          = ',3(g18.12,2x))
2508 110   format('   count_beta_diff,count_beta_diff_gt_20= ',2(g18.12,2x))
2509 115   format('   beta_diff_max_gt_20(deg)             = ',g18.12)
2510    end if
2511  end subroutine check_rotate_wind
2512
2513 end module gridmod