| blob | 9b78583a6816f518d741e6eb673ae30876191b2f |
1 ! Create an initial data set for the WRF model based on real data. This
2 ! program is specifically set up for the Eulerian, mass-based coordinate.
3 PROGRAM tc_data
4 USE module_machine
5 USE module_domain, ONLY : domain, alloc_and_configure_domain, &
6 domain_clock_set, head_grid, program_name, domain_clockprint
7 USE module_io_domain
10 USE module_driver_constants
11 USE module_configure, ONLY : grid_config_rec_type, model_config_rec, &
12 initial_config, get_config_as_buffer, set_config_as_buffer
13 USE module_timing
14 USE module_state_description, ONLY : realonly
15 USE module_symbols_util, ONLY: wrfu_cal_gregorian
16 USE module_utility, ONLY : WRFU_finalize
18 IMPLICIT NONE
21 REAL :: time , bdyfrq
23 INTEGER :: loop , levels_to_process , debug_level
26 TYPE(domain) , POINTER :: null_domain
27 TYPE(domain) , POINTER :: grid , another_grid
28 TYPE(domain) , POINTER :: grid_ptr , grid_ptr2
29 TYPE (grid_config_rec_type) :: config_flags
30 INTEGER :: number_at_same_level
32 INTEGER :: max_dom, domain_id , grid_id , parent_id , parent_id1 , id
33 INTEGER :: e_we , e_sn , i_parent_start , j_parent_start
34 INTEGER :: idum1, idum2
35 #ifdef DM_PARALLEL
36 INTEGER :: nbytes
37 INTEGER, PARAMETER :: configbuflen = 4* CONFIG_BUF_LEN
38 INTEGER :: configbuf( configbuflen )
39 LOGICAL , EXTERNAL :: wrf_dm_on_monitor
40 #endif
41 LOGICAL found_the_id
43 INTEGER :: ids , ide , jds , jde , kds , kde
44 INTEGER :: ims , ime , jms , jme , kms , kme
45 INTEGER :: ips , ipe , jps , jpe , kps , kpe
46 INTEGER :: ijds , ijde , spec_bdy_width
47 INTEGER :: i , j , k , idts, rc
48 INTEGER :: sibling_count , parent_id_hold , dom_loop
50 CHARACTER (LEN=80) :: message
52 INTEGER :: start_year , start_month , start_day , start_hour , start_minute , start_second
53 INTEGER :: end_year , end_month , end_day , end_hour , end_minute , end_second
54 INTEGER :: interval_seconds , real_data_init_type
55 INTEGER :: time_loop_max , time_loop
56 real::t1,t2
57 INTERFACE
58 SUBROUTINE Setup_Timekeeping( grid )
59 USE module_domain, ONLY : domain
60 TYPE(domain), POINTER :: grid
61 END SUBROUTINE Setup_Timekeeping
62 END INTERFACE
64 #include "version_decl"
66 ! Define the name of this program (program_name defined in module_domain)
68 program_name = "TC_EM " // TRIM(release_version) // " PREPROCESSOR"
70 ! The TC bogus algorithm assumes that the user defines a central point, and then
71 ! allows the program to remove a typhoon based on a distance in km. This is
72 ! implemented on a single processor only.
74 #ifdef DM_PARALLEL
75 IF ( .NOT. wrf_dm_on_monitor() ) THEN
76 CALL wrf_error_fatal( 'TC bogus must run with a single processor only, re-run with num procs set to 1' )
77 END IF
78 #endif
80 #ifdef DM_PARALLEL
81 CALL disable_quilting
82 #endif
84 ! Initialize the modules used by the WRF system. Many of the CALLs made from the
85 ! init_modules routine are NO-OPs. Typical initializations are: the size of a
86 ! REAL, setting the file handles to a pre-use value, defining moisture and
87 ! chemistry indices, etc.
89 CALL wrf_debug ( 100 , 'real_em: calling init_modules ' )
90 CALL init_modules(1) ! Phase 1 returns after MPI_INIT() (if it is called)
91 CALL WRFU_Initialize( defaultCalendar=WRFU_CAL_GREGORIAN, rc=rc )
92 CALL init_modules(2) ! Phase 2 resumes after MPI_INIT() (if it is called)
94 ! The configuration switches mostly come from the NAMELIST input.
96 #ifdef DM_PARALLEL
97 IF ( wrf_dm_on_monitor() ) THEN
98 CALL initial_config
99 END IF
100 CALL get_config_as_buffer( configbuf, configbuflen, nbytes )
101 CALL wrf_dm_bcast_bytes( configbuf, nbytes )
102 CALL set_config_as_buffer( configbuf, configbuflen )
103 ! CALL wrf_dm_initialize
104 #else
105 CALL initial_config
106 #endif
109 CALL nl_get_debug_level ( 1, debug_level )
110 CALL set_wrf_debug_level ( debug_level )
112 CALL wrf_message ( program_name )
114 ! There are variables in the Registry that are only required for the real
115 ! program, fields that come from the WPS package. We define the run-time
116 ! flag that says to allocate space for these input-from-WPS-only arrays.
118 CALL nl_set_use_wps_input ( 1 , REALONLY )
120 ! Allocate the space for the mother of all domains.
122 NULLIFY( null_domain )
123 CALL wrf_debug ( 100 , 'real_em: calling alloc_and_configure_domain ' )
124 CALL alloc_and_configure_domain ( domain_id = 1 , &
125 grid = head_grid , &
126 parent = null_domain , &
127 kid = -1 )
129 grid => head_grid
130 CALL nl_get_max_dom ( 1 , max_dom )
132 IF ( model_config_rec%interval_seconds .LE. 0 ) THEN
133 CALL wrf_error_fatal( 'namelist value for interval_seconds must be > 0')
134 END IF
136 all_domains : DO domain_id = 1 , max_dom
138 IF ( ( model_config_rec%input_from_file(domain_id) ) .OR. &
139 ( domain_id .EQ. 1 ) ) THEN
141 CALL Setup_Timekeeping ( grid )
142 CALL set_current_grid_ptr( grid )
143 CALL domain_clockprint ( 150, grid, &
144 'DEBUG real: clock after Setup_Timekeeping,' )
145 CALL domain_clock_set( grid, &
146 time_step_seconds=model_config_rec%interval_seconds )
147 CALL domain_clockprint ( 150, grid, &
148 'DEBUG real: clock after timeStep set,' )
151 CALL wrf_debug ( 100 , 'real_em: calling set_scalar_indices_from_config ' )
152 CALL set_scalar_indices_from_config ( grid%id , idum1, idum2 )
154 CALL wrf_debug ( 100 , 'real_em: calling model_to_grid_config_rec ' )
155 CALL model_to_grid_config_rec ( grid%id , model_config_rec , config_flags )
157 ! Initialize the WRF IO: open files, init file handles, etc.
159 CALL wrf_debug ( 100 , 'real_em: calling init_wrfio' )
160 CALL init_wrfio
162 ! Some of the configuration values may have been modified from the initial READ
163 ! of the NAMELIST, so we re-broadcast the configuration records.
165 #ifdef DM_PARALLEL
166 CALL wrf_debug ( 100 , 'real_em: re-broadcast the configuration records' )
167 CALL get_config_as_buffer( configbuf, configbuflen, nbytes )
168 CALL wrf_dm_bcast_bytes( configbuf, nbytes )
169 CALL set_config_as_buffer( configbuf, configbuflen )
170 #endif
172 ! No looping in this layer.
174 CALL wrf_debug ( 100 , 'calling tc_med_sidata_input' )
175 CALL tc_med_sidata_input ( grid , config_flags )
176 CALL wrf_debug ( 100 , 'backfrom tc_med_sidata_input' )
178 ELSE
179 CYCLE all_domains
180 END IF
182 END DO all_domains
184 CALL set_current_grid_ptr( head_grid )
186 ! We are done.
188 CALL wrf_debug ( 0 , 'real_em: SUCCESS COMPLETE REAL_EM INIT' )
190 CALL wrf_shutdown
192 CALL WRFU_Finalize( rc=rc )
195 END PROGRAM tc_data
198 !-----------------------------------------------------------------
199 SUBROUTINE tc_med_sidata_input ( grid , config_flags )
200 ! Driver layer
201 USE module_domain
202 USE module_io_domain
203 ! Model layer
204 USE module_configure
205 USE module_bc_time_utilities
206 USE module_optional_input
208 USE module_date_time
209 USE module_utility
211 IMPLICIT NONE
214 ! Interface
215 INTERFACE
216 SUBROUTINE start_domain ( grid , allowed_to_read ) ! comes from module_start in appropriate dyn_ directory
217 USE module_domain
218 TYPE (domain) grid
219 LOGICAL, INTENT(IN) :: allowed_to_read
220 END SUBROUTINE start_domain
221 END INTERFACE
223 ! Arguments
224 TYPE(domain) :: grid
225 TYPE (grid_config_rec_type) :: config_flags
226 ! Local
227 INTEGER :: time_step_begin_restart
228 INTEGER :: idsi , ierr , myproc, internal_time_loop
229 CHARACTER (LEN=80) :: si_inpname
230 CHARACTER (LEN=80) :: message
232 CHARACTER(LEN=19) :: start_date_char , end_date_char , current_date_char , next_date_char
234 INTEGER :: time_loop_max , loop, rc
235 INTEGER :: julyr , julday
236 REAL :: gmt
237 real::t1,t2,t3,t4
239 grid%input_from_file = .true.
240 grid%input_from_file = .false.
242 CALL tc_compute_si_start ( model_config_rec%start_year (grid%id) , &
243 model_config_rec%start_month (grid%id) , &
244 model_config_rec%start_day (grid%id) , &
245 model_config_rec%start_hour (grid%id) , &
246 model_config_rec%start_minute(grid%id) , &
247 model_config_rec%start_second(grid%id) , &
248 model_config_rec%interval_seconds , &
249 model_config_rec%real_data_init_type , &
250 start_date_char)
252 end_date_char = start_date_char
253 IF ( end_date_char .LT. start_date_char ) THEN
254 CALL wrf_error_fatal( 'Ending date in namelist ' // end_date_char // ' prior to beginning date ' // start_date_char )
255 END IF
256 print *,"the start date char ",start_date_char
257 print *,"the end date char ",end_date_char
259 time_loop_max = 1
260 ! Override stop time with value computed above.
261 CALL domain_clock_set( grid, stop_timestr=end_date_char )
263 ! TBH: for now, turn off stop time and let it run data-driven
264 CALL WRFU_ClockStopTimeDisable( grid%domain_clock, rc=rc )
265 CALL wrf_check_error( WRFU_SUCCESS, rc, &
266 'WRFU_ClockStopTimeDisable(grid%domain_clock) FAILED', &
267 __FILE__ , &
268 __LINE__ )
269 CALL domain_clockprint ( 150, grid, &
270 'DEBUG med_sidata_input: clock after stopTime set,' )
272 ! Here we define the initial time to process, for later use by the code.
274 current_date_char = start_date_char
275 start_date = start_date_char // '.0000'
276 current_date = start_date
278 CALL nl_set_bdyfrq ( grid%id , REAL(model_config_rec%interval_seconds) )
280 !!!!!!! Loop over each time period to process.
282 CALL cpu_time ( t1 )
283 DO loop = 1 , time_loop_max
285 internal_time_loop = loop
286 IF ( ( grid%id .GT. 1 ) .AND. ( loop .GT. 1 ) .AND. &
287 ( model_config_rec%grid_fdda(grid%id) .EQ. 0 ) .AND. &
288 ( model_config_rec%sst_update .EQ. 0 ) ) EXIT
290 print *,' '
291 print *,'-----------------------------------------------------------------------------'
292 print *,' '
293 print '(A,I2,A,A,A,I4,A,I4)' , &
294 ' Domain ',grid%id,': Current date being processed: ',current_date, ', which is loop #',loop,' out of ',time_loop_max
296 ! After current_date has been set, fill in the julgmt stuff.
298 CALL geth_julgmt ( config_flags%julyr , config_flags%julday , config_flags%gmt )
300 print *,'configflags%julyr, %julday, %gmt:',config_flags%julyr, config_flags%julday, config_flags%gmt
301 ! Now that the specific Julian info is available, save these in the model config record.
303 CALL nl_set_gmt (grid%id, config_flags%gmt)
304 CALL nl_set_julyr (grid%id, config_flags%julyr)
305 CALL nl_set_julday (grid%id, config_flags%julday)
307 ! Open the input file for tc stuff. Either the "new" one or the "old" one. The "new" one could have
308 ! a suffix for the type of the data format. Check to see if either is around.
310 CALL cpu_time ( t3 )
311 WRITE ( wrf_err_message , FMT='(A,A)' )'med_sidata_input: calling open_r_dataset for ', &
312 TRIM(config_flags%auxinput1_inname)
313 CALL wrf_debug ( 100 , wrf_err_message )
314 IF ( config_flags%auxinput1_inname(1:8) .NE. 'wrf_real' ) THEN
315 CALL construct_filename4a( si_inpname , config_flags%auxinput1_inname , grid%id , 2 , &
316 current_date_char , config_flags%io_form_auxinput1 )
317 ELSE
318 CALL construct_filename2a( si_inpname , config_flags%auxinput1_inname , grid%id , 2 , &
319 current_date_char )
320 END IF
321 CALL open_r_dataset ( idsi, TRIM(si_inpname) , grid , config_flags , "DATASET=AUXINPUT1", ierr )
322 IF ( ierr .NE. 0 ) THEN
323 CALL wrf_error_fatal( 'error opening ' // TRIM(si_inpname) // &
324 ' for input; bad date in namelist or file not in directory' )
325 END IF
327 ! Input data.
329 CALL wrf_debug ( 100 , 'med_sidata_input: calling input_aux_model_input1' )
330 CALL input_aux_model_input1 ( idsi , grid , config_flags , ierr )
331 WRITE ( wrf_err_message , FMT='(A,I10,A)' ) 'Timing for input ',NINT(t4-t3) ,' s.'
332 CALL wrf_debug( 0, wrf_err_message )
333 CALL close_dataset ( idsi , config_flags , "DATASET=AUXINPUT1" )
334 CALL cpu_time ( t4 )
335 print *,"the out name ",config_flags%auxinput1_outname
337 ! Possible optional SI input. This sets flags used by init_domain.
338 CALL cpu_time ( t3 )
339 IF ( loop .EQ. 1 ) THEN
340 already_been_here = .FALSE.
341 CALL wrf_debug ( 100 , 'med_sidata_input: calling init_module_optional_input' )
342 CALL init_module_optional_input ( grid , config_flags )
343 END IF
344 CALL wrf_debug ( 100 , 'med_sidata_input: calling optional_input' )
345 CALL optional_input ( grid , idsi , config_flags )
347 CALL model_to_grid_config_rec ( grid%id , model_config_rec , config_flags )
350 CALL cpu_time ( t3 )
351 print *,"before assemble output"
353 CALL assemble_output ( grid , config_flags , loop , time_loop_max, current_date_char )
354 CALL cpu_time ( t4 )
355 WRITE ( wrf_err_message , FMT='(A,I10,A)' ) 'Timing for output ',NINT(t4-t3) ,' s.'
356 CALL wrf_debug( 0, wrf_err_message )
357 CALL cpu_time ( t2 )
358 WRITE ( wrf_err_message , FMT='(A,I4,A,I10,A)' ) 'Timing for loop # ',loop,' = ',NINT(t2-t1) ,' s.'
359 CALL wrf_debug( 0, wrf_err_message )
361 CALL cpu_time ( t1 )
362 END DO
364 END SUBROUTINE tc_med_sidata_input
367 !-------------------------------------------------------------------------------------
368 SUBROUTINE tc_compute_si_start( &
369 start_year , start_month , start_day , start_hour , start_minute , start_second , &
370 interval_seconds , real_data_init_type , &
371 start_date_char)
373 USE module_date_time
375 IMPLICIT NONE
377 INTEGER :: start_year , start_month , start_day , start_hour , start_minute , start_second
378 INTEGER :: end_year , end_month , end_day , end_hour , end_minute , end_second
379 INTEGER :: interval_seconds , real_data_init_type
380 INTEGER :: time_loop_max , time_loop
382 CHARACTER(LEN=19) :: current_date_char , start_date_char , end_date_char , next_date_char
384 #ifdef PLANET
385 WRITE ( start_date_char , FMT = '(I4.4,"-",I5.5,"_",I2.2,":",I2.2,":",I2.2)' ) &
386 start_year,start_day,start_hour,start_minute,start_second
387 #else
388 WRITE ( start_date_char , FMT = '(I4.4,"-",I2.2,"-",I2.2,"_",I2.2,":",I2.2,":",I2.2)' ) &
389 start_year,start_month,start_day,start_hour,start_minute,start_second
390 #endif
393 END SUBROUTINE tc_compute_si_start
395 !-----------------------------------------------------------------------
396 SUBROUTINE assemble_output ( grid , config_flags , loop , time_loop_max,current_date_char )
398 USE module_big_step_utilities_em
399 USE module_domain
400 USE module_io_domain
401 USE module_configure
402 USE module_date_time
403 USE module_bc
404 IMPLICIT NONE
406 TYPE(domain) :: grid
407 TYPE (grid_config_rec_type) :: config_flags
409 INTEGER , INTENT(IN) :: loop , time_loop_max
411 !These values are in the name list and are avaiable from
412 !from the config_flags.
413 real :: latc_loc,lonc_loc,vmax,vmax_ratio,rmax
416 INTEGER :: ids , ide , jds , jde , kds , kde
417 INTEGER :: ims , ime , jms , jme , kms , kme
418 INTEGER :: ips , ipe , jps , jpe , kps , kpe
419 INTEGER :: ijds , ijde , spec_bdy_width
420 INTEGER :: i , j , k , idts
422 INTEGER :: id1 , interval_seconds , ierr, rc, sst_update, grid_fdda
423 INTEGER , SAVE :: id, id2, id4
424 CHARACTER (LEN=80) :: tcoutname , bdyname
425 CHARACTER(LEN= 4) :: loop_char
426 CHARACTER(LEN=19) :: current_date_char
428 character *19 :: temp19
429 character *24 :: temp24 , temp24b
431 real::t1,t2
433 ! Various sizes that we need to be concerned about.
435 ids = grid%sd31
436 ide = grid%ed31
437 kds = grid%sd32
438 kde = grid%ed32
439 jds = grid%sd33
440 jde = grid%ed33
442 ims = grid%sm31
443 ime = grid%em31
444 kms = grid%sm32
445 kme = grid%em32
446 jms = grid%sm33
447 jme = grid%em33
449 ips = grid%sp31
450 ipe = grid%ep31
451 kps = grid%sp32
452 kpe = grid%ep32
453 jps = grid%sp33
454 jpe = grid%ep33
456 ijds = MIN ( ids , jds )
457 ijde = MAX ( ide , jde )
459 ! Boundary width, scalar value.
461 spec_bdy_width = model_config_rec%spec_bdy_width
462 interval_seconds = model_config_rec%interval_seconds
463 sst_update = model_config_rec%sst_update
464 grid_fdda = model_config_rec%grid_fdda(grid%id)
467 ! here we are going to call the temp subroutine.
468 print *,"here before temp sub"
469 grid%stand_lon = 141.
470 grid%cen_lat = 26.99999
471 grid%map_proj = 3
474 latc_loc = config_flags%latc_loc
475 lonc_loc = config_flags%lonc_loc
476 vmax = config_flags%vmax_meters_per_second
477 vmax_ratio = config_flags%vmax_ratio
478 rmax = config_flags%rmax
479 call tc_bogus(grid%cen_lat,grid%stand_lon,grid%map_proj,grid%truelat1,grid%truelat2, &
480 grid%dx,ipe,jpe,grid%num_metgrid_levels,latc_loc,lonc_loc,vmax,vmax_ratio,rmax,grid)
481 stop
482 !Notes jps is the starting y. Usally 1
483 ! jpe is the ending y position on the staggered grid north south or ny
484 ! ipe is the ending x postiion on the staggered grid east west of nx
486 ! Open the tc bogused output file. cd
487 CALL construct_filename4a( tcoutname , config_flags%auxinput1_outname , grid%id , 2 , &
488 current_date_char , config_flags%io_form_auxinput1 )
490 print *,"outfile name from construct filename ",tcoutname
491 CALL open_w_dataset ( id1, TRIM(tcoutname) , grid , config_flags ,output_aux_model_input1,"DATASET=AUXINPUT1",ierr )
492 IF ( ierr .NE. 0 ) THEN
493 CALL wrf_error_fatal( 'tc_em: error opening tc bogus file for writing' )
494 END IF
495 CALL output_aux_model_input1( id1, grid , config_flags , ierr )
496 CALL close_dataset ( id1 , config_flags , "DATASET=AUXINPUT1" )
499 END SUBROUTINE assemble_output
501 !----------------------------------------------------------------------------------------------
503 SUBROUTINE tc_bogus(centerlat,stdlon,nproj,truelat1,truelat2,dsm,ew,ns,nz, &
504 latc_loc,lonc_loc,vmax,vmax_ratio,rmax,grid)
506 !center lat is the center latitude from the global attributes
507 !stdlon is the center longitude from the global attributes
508 !nproj is the map projection from the global attributes
509 !dsm is the spacing in meters from the global attributes
510 !ew is the west_east_stag from the dimensions. We want the full domain.
511 !ns is the south_north_stag from the dimensions. We want the full domain.
512 !nz is the number of metgrid levels from the dimensions
513 !latc_loc is the center latitude of the bogus strorm. It comes from the namelist entry
514 !lonc_loc is the center longitude of the bogus strorm. It comes from the namelist entry
515 !vmax is the max vortex it comes from the namelist entry
516 !vmax_ratio
517 !rmax
519 !module_llxy resides in the share directory.
520 USE module_llxy
521 !This is for the large structure (grid)
522 USE module_domain
526 IMPLICIT NONE
527 TYPE(domain) :: grid
528 integer ew,ns,nz
529 integer nproj
530 integer num_storm,nstrm
531 real :: centerlat,stdlon,conef,truelat1,truelat2,dsm,dx
532 real :: latc_loc,lonc_loc,vmax,vmax_ratio,rmax
534 real :: press(ns,ew,nz),p(nz), rhmx(nz), vwgt(nz)
535 real, dimension(:,:,:) , allocatable :: u11,v11,t11,rh11,phi11
536 real, dimension(:,:,:) , allocatable :: umass,vmass
537 real, dimension(:,:,:) , allocatable :: u1 , v1 , t1 , rh1 , phi1
538 real, dimension(ns,ew) :: mfx,mfd,lond,terrain,cor,xlat,pslx
541 CHARACTER*2 jproj
542 LOGICAL :: l_tcbogus
544 !These values are hard coded in the procedure proc_namelist.F
545 !Since we don't include this module we will declare them here
546 !and hard code them.
547 real :: r_search,r_vor,beta,devps,humidity_max
548 real :: devpc,const,r_vor2,cnst,alphar,epsilon,vormx , rad , sum_q
549 real :: avg_q ,q_old,ror,q_new,sum_rh,avg_rh,rh_min,rhbkg,rhbog,dph
550 real :: rh_max,min_RH_value,r_ratio,ps
551 integer :: vert_variation
552 integer :: i,k,j,jx,ix,kx
553 integer :: k00,kfrm ,kto ,k85,n_iter,i_mvc,j_mvc,nct,itr
554 integer :: strmci(nz), strmcj(nz)
555 real :: disx,disy,alpha,degran,pie,rovcp,cp
556 REAL :: rho,pprm,phip0,x0,y0,vmx,xico,xjco,xicn,xjcn,p85,xlo,rconst
557 real :: ptop_in_pa
558 real :: latinc,loninc
559 real :: rtemp,colat0,colat
560 REAL :: q1(ns,ew,nz), psi1(ns,ew,nz)
562 ! This is the entire map projection enchilada.
563 TYPE(proj_info) :: proj
567 REAL :: lat1 , lon1
568 ! These values are read in from the data set.
569 real :: knowni,knownj
571 ! TC bogus
572 REAL utcr(ns,ew,nz), vtcr(ns,ew,nz)
573 REAL utcp(ns,ew,nz), vtcp(ns,ew,nz)
574 REAL psitc(ns,ew,nz), psiv(nz)
575 REAL vortc(ns,ew,nz), vorv(nz)
576 REAL tptc(ns,ew,nz), rhtc(ns,ew,nz)
577 REAL phiptc(ns,ew,nz)
579 ! Work arrays
580 REAL uuwork(nz), vvwork(nz)
581 REAL vort(ns,ew,nz), div(ns,ew,nz)
582 REAL vortsv(ns,ew,nz)
583 REAL theta(ns,ew,nz), t_reduce(ns,ew,nz)
584 REAL ug(ns,ew,nz), vg(ns,ew,nz), vorg(ns,ew,nz)
585 REAL uchi(ns,ew,nz), vchi(ns,ew,nz)
586 REAL delpx(ns,ew)
588 !subroutines for relaxation
589 REAL outold(ns,ew), outnew(ns,ew)
590 REAL rd(ns,ew), ff(ns,ew)
591 REAL tmp1(ns,ew), tmp2(ns,ew)
593 ! Background fields.
594 REAL , DIMENSION (ns,ew,nz) :: u0, v0, t0, t00, rh0, q0, &
595 phi0, psi0, chi
598 ! Perturbations
599 REAL psipos(ns,ew,nz), upos(ns,ew,nz), vpos(ns,ew,nz), &
600 tpos(ns,ew,nz), psi(ns,ew,nz), &
601 phipos(ns,ew,nz), phip(ns,ew,nz)
603 ! Final fields.
604 REAL u2(ns,ew,nz), v2(ns,ew,nz), &
605 t2(ns,ew,nz), &
606 z2(ns,ew,nz), phi2(ns,ew,nz), &
607 rh2(ns,ew,nz), q2(ns,ew,nz)
610 print *,"the dimensions: north-south = ",ns," east-west =",ew
611 IF (nproj .EQ. 1) THEN
612 jproj = 'LC'
613 print *,"Lambert Conformal projection"
614 ELSE IF (nproj .EQ. 2) THEN
615 jproj = 'ST'
616 ELSE IF (nproj .EQ. 3) THEN
617 jproj = 'ME'
618 print *,"A mercator projection"
619 END IF
620 num_storm = 1
622 !NOTE: the ptop in the namelist is not in pa
623 ptop_in_pa = 10000
624 knowni = 1.
625 knownj = 1.
626 pie = 3.141592653589793
627 degran = pie/180.
628 rconst = 287.04
629 min_RH_value = 5.0
630 cp = 1004.0
631 rovcp = rconst/cp
633 r_search = 400000.0
634 r_vor = 300000.0
635 r_vor2 = r_vor * 4
636 beta = 0.5
637 devpc= 40.0
638 vert_variation = 1
639 humidity_max = 95.0
640 alphar = 1.8
641 latinc = -999.
642 loninc = -999.
644 ! Set up initializations for map projection so that the lat/lon
645 ! of the tropical storm can be put into model (i,j) space. This needs to be done once per
646 ! map projection definition. Since this is the domain that we are "GOING TO", it is a once
647 ! per regridder requirement. If the user somehow ends up calling this routine for several
648 ! time periods, there is no problemos, just a bit of overhead with redundant calls.
650 dx = dsm
651 lat1 = grid%xlat_gc(1,1)
652 lon1 = grid%xlong_gc(1,1)
653 IF( jproj .EQ. 'ME' )THEN
654 IF ( lon1 .LT. -180. ) lon1 = lon1 + 360.
655 IF ( lon1 .GT. 180. ) lon1 = lon1 - 360.
656 IF ( stdlon .LT. -180. ) stdlon = stdlon + 360.
657 IF ( stdlon .GT. 180. ) stdlon = stdlon - 360.
658 CALL map_set ( proj_merc, proj, lat1, lon1, lat1, lon1, knowni, knownj, dx, &
659 latinc,loninc,stdlon , truelat1 , truelat2)
660 conef = 0.
661 ELSE IF ( jproj .EQ. 'LC' ) THEN
662 CALL map_set (proj_lc,proj, lat1, lon1, lat1, lon1, knowni, knownj, dx, &
663 latinc,loninc,stdlon , truelat1 , truelat2)
664 conef = proj%cone
665 ELSE IF ( jproj .EQ. 'ST' ) THEN
666 conef = 1.
667 CALL map_set ( proj_ps,proj,lat1, lon1, lat1, lon1, knowni, knownj, dx, &
668 latinc,loninc,stdlon , truelat1 , truelat2)
669 END IF
672 ! Put the vertical column of pressures into hPa from Pa. The first level is defined as the surface.
673 ! The second level is usually 1000 mb. The last level, kx, is defined as ptop. These come in from
674 ! the namelist
675 kx = nz
676 DO k=1,kx
677 p(k) = grid%p_gc(1,k,1)*0.01
678 END DO
680 ! Initialize the vertical profiles for humidity and weighting.
681 !The ptop variable will be read in from the namelist
682 IF ( ( ptop_in_pa .EQ. 40000. ) .OR. ( ptop_in_pa .EQ. 60000. ) ) THEN
683 PRINT '(A)','Hold on pardner, your value for PTOP is gonna cause problems for the TC bogus option.'
684 PRINT '(A)','Make it higher up than 400 mb.'
685 STOP 'ptop_woes_for_tc_bogus'
686 END IF
688 IF ( vert_variation .EQ. 1 ) THEN
689 DO k=1,kx
690 IF ( p(k) .GT. 400. ) THEN
691 rhmx(k) = humidity_max
692 ELSE
693 rhmx(k) = humidity_max * MAX( 0.1 , (p(k) - ptop_in_pa/100.)/(400.-ptop_in_pa/100.) )
694 END IF
696 IF ( p(k) .GT. 600. ) THEN
697 vwgt(k) = 1.0
698 ELSE IF ( p(k) .LE. 100. ) THEN
699 vwgt(k) = 0.0001
700 ELSE
701 vwgt(k) = MAX ( 0.0001 , (p(k)-ptop_in_pa/100.)/(600.-ptop_in_pa/100.) )
702 END IF
703 END DO
705 ELSE IF ( vert_variation .EQ. 2 ) THEN
706 IF ( kx .eq. 24 ) THEN
707 rhmx = (/ 95., 95., 95., 95., 95., 95., 95., 95., &
708 95., 95., 95., 95., 95., 90., 85., 80., 75., &
709 70., 66., 60., 39., 10., 10., 10./)
710 vwgt = (/ 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 0.9850, &
711 0.9680, 0.9500, 0.9290, 0.9060, 0.8810, 0.8500, 0.7580, 0.6500, 0.5100, &
712 0.3500, 0.2120, 0.0500, 0.0270, 0.0001, 0.0001, 0.0001/)
713 ELSE
714 PRINT '(A)','Number of vertical levels assumed to be 24 for AFWA TC bogus option'
715 STOP 'AFWA_TC_BOGUS_LEVEL_ERROR'
716 END IF
717 END IF
719 ns = ns
720 ew = ew
721 !Put the U and V into the new arrays.
722 allocate(u11 (1:ns, 1:ew, 1:nz))
723 allocate(v11 (1:ns, 1:ew, 1:nz))
724 allocate(umass(1:ew, 1:nz, 1:ns))
725 allocate(vmass(1:ew, 1:nz, 1:ns))
727 !first U
728 do i = 1,ns-1
729 do j = 1,ew
730 do k = 1,nz
731 u11(i,j,k) = grid%u_gc(j,k,i)
732 end do
733 end do
734 end do
738 !now V
739 do i = 1,ns
740 do j = 1,ew-1
741 do k = 1,nz
742 v11(i,j,k) = grid%v_gc(j,k,i)
743 end do
744 end do
745 end do
747 do i = 1,ns
748 print *,v11(i,1,1)
749 end do
750 stop
753 ! do i = 1,ew
754 ! do k = 1,nz
755 ! do j = 1,ns-1
756 ! umass(i,k,j) = 0.5 * (grid%u_gc(i,k,j) + grid%u_gc(i + 1,k,j))
757 ! end do
758 ! end do
759 ! end do
762 ! do i = 1,ew-1
763 ! do k = 1,nz
764 ! do j = 1,ns
765 ! vmass(i,k,j) = 0.5 * (grid%v_gc(i,k,j) + grid%v_gc(i,k,j+1))
766 ! end do
767 ! end do
768 ! end do
771 !here we need to transpose the level(k) and the ns direction
772 !So the ordering will be j,i,k - ns,ew,level
773 !we want MM5 ordering
774 allocate(u11 (1:ns, 1:ew, 1:nz))
775 allocate(v11 (1:ns, 1:ew, 1:nz))
776 allocate(t11 (1:ns, 1:ew, 1:nz))
777 allocate(rh11 (1:ns, 1:ew, 1:nz))
778 allocate(phi11(1:ns, 1:ew, 1:nz))
780 do i = 1,ns-1
781 do j = 1,ew
782 do k = 1,nz
783 u11(i,j,k) = grid%u_gc(j,k,i)
784 end do
785 end do
786 end do
788 do i = 1,ns-1
789 print *,i,u11(i,1,1)
790 end do
792 stop
793 do i = 1,ns
794 do j = 1,ew
795 do k = 1,nz
796 u11(i,j,k) = umass(j,k,i)
797 v11(i,j,k) = vmass(j,k,i)
798 t11(i,j,k) = grid%t_gc(j,k,i)
799 rh11(i,j,k) = grid%rh_gc(j,k,i)
800 phi11(i,j,k) = grid%ght_gc(j,k,i)
801 end do
802 end do
803 end do
804 deallocate(umass)
805 deallocate(vmass)
809 ! Save initial fields.
810 allocate(u1 (1:ns, 1:ew, 1:nz))
811 allocate(v1 (1:ns, 1:ew, 1:nz))
812 allocate(t1 (1:ns, 1:ew, 1:nz))
813 allocate(rh1 (1:ns, 1:ew, 1:nz))
814 allocate(phi1(1:ns, 1:ew, 1:nz))
815 !for wrf we need this
816 jx = ew
817 ix = ns
818 u1 = u11
819 v1 = v11
820 t1 = t11
821 rh1 = rh11
822 phi1 = phi11 * 9.81
824 !The two d fields
825 do i = 1,ns
826 do j = 1,ew
827 ! pslx(i,j) = grid%pslv_gc(j,i) * 0.01
828 ! cor(i,j) = grid%f(j,i) !coreolous
829 ! lond(i,j) = grid%xlong_gc(j,i)
830 terrain(i,j) = grid%ht_gc(j,i)
831 xlat(i,j) = grid%xlat_gc(1,1)
832 end do
833 end do
834 cor = transpose(grid%f)
835 pslx = transpose(grid%pslv_gc)
836 pslx = pslx * 0.01
837 lond = transpose(grid%xlong_gc)
838 call crs2dot(lond,ns,ew)
841 DO k=1,kx
842 CALL expand ( rh1(1,1,k) , ns-1 , ew-1 , ns , ew )
843 CALL expand ( t1 (1,1,k) , ns-1 , ew-1 , ns , ew )
844 END DO
847 ! Loop over the number of storms to process.
849 l_tcbogus = .FALSE.
851 k00 = 2
852 kfrm = k00
853 p85 = 850.
855 kto = kfrm
856 DO k=kfrm+1,kx
857 IF ( p(k) .GE. p85 ) THEN
858 kto = kto + 1
859 END IF
860 END DO
861 k85 = kto
864 ! Parameters for max wind
865 rho = 1.2
866 pprm = devpc*100.
867 phip0= pprm/rho
869 !latc_loc and lonc_loc will come in from the namelist. Hard code
870 !them for now.
872 CALL latlon_to_ij ( proj , latc_loc , lonc_loc , x0 , y0 )
873 IF ( ( x0 .LT. 1. ) .OR. ( x0 .GT. REAL(jx-1) ) .OR. &
874 ( y0 .LT. 1. ) .OR. ( y0 .GT. REAL(ix-1) ) ) THEN
875 PRINT '(A,I3,A,A,A)',' Storm position is outside the computational domain.'
876 PRINT '(A,2F6.2,A)' ,' Storm postion: (x,y) = ',x0,y0,'.'
877 stop
878 END IF
880 l_tcbogus = .TRUE.
881 ! Bogus vortex specifications, vmax (m/s); rmax (m);
882 vmx = vmax * vmax_ratio
883 print *,"the vmx ",vmax,vmax_ratio,vmx
885 IF ( latc_loc .LT. 0. ) THEN
886 vmx = -vmx
887 END IF
889 IF ( vmax .LE. 0. ) THEN
890 vmx = SQRT( 2.*(1-beta)*ABS(phip0) )
891 END IF
893 xico = y0
894 xjco = x0
895 xicn = xico
896 xjcn = xjco
897 n_iter = 1
900 ! Start computing.
902 PRINT '(/,A,I3,A,A,A)' ,'---> TC: Processing storm number= ',nstrm,' for date= simlaku'
903 PRINT '(A,F6.2,A,F7.2,A)' ,' Storm center lat= ',latc_loc,' lon= ',lonc_loc,'.'
904 PRINT '(A,2F6.2,A)' ,' Storm center grid position (x,y)= ',xjcn,xicn,'.'
905 PRINT '(A,F5.2,F9.2,A)' ,' Storm max wind (m/s) and max radius (m)= ',vmx,rmax,'.'
906 PRINT '(A,F5.2,A)' ,' Estmated central press dev (mb)= ',devpc,'.'
909 ! Initialize storm center to (1,1)
911 DO k=1,kx
912 strmci(k) = 1
913 strmcj(k) = 1
914 END DO
916 ! Define complete field of bogus storm
917 !Note dx is spacing in meters
918 DO i=1,ns
919 DO j=1,ew
920 disx = REAL(j) - xjcn
921 disy = REAL(i) - xicn
922 CALL rankine(disx,disy,dx,kx,vwgt,rmax,vmx,uuwork,vvwork,psiv,vorv)
923 DO k=1,kx
924 utcp(i,j,k) = uuwork(k)
925 vtcp(i,j,k) = vvwork(k)
926 psitc(i,j,k) = psiv(k)
927 vortc(i,j,k) = vorv(k)
928 END DO
929 END DO
930 END DO
933 ! Rotate wind to map proj
935 DO i=1,ns-1
936 DO j=1,ew-1
937 !i don't think we want xlong_gc here check with dave.
938 xlo = stdlon-lond(i,j)
939 ! print *,xlo
940 IF ( xlo .GT. 180.)xlo = xlo-360.
941 IF ( xlo .LT.-180.)xlo = xlo+360.
943 alpha = xlo*conef*degran*SIGN(1.,centerlat)
944 DO k=1,kx
945 utcr(i,j,k) = vtcp(i,j,k)*SIN(alpha)+utcp(i,j,k)*COS(alpha)
946 vtcr(i,j,k) = vtcp(i,j,k)*COS(alpha)-utcp(i,j,k)*SIN(alpha)
947 if(vtcr(i,j,k) .gt. 500.) then
948 print *,i,j,k,"a very bad value of vtcr"
949 stop
950 end if
951 if(utcr(i,j,k) .gt. 500.) then
952 print *,i,j,k,"a very bad value of utcr"
953 stop
954 end if
955 END DO
956 END DO
957 END DO
960 DO k=1,kx
961 DO i=1,ns-1
962 utcr(i,ew,k) = utcr(i,ew-1,k)
963 vtcr(i,ew,k) = vtcr(i,ew-1,k)
964 END DO
965 DO j=1,ew-1
966 utcr(ns,j,k) = utcr(ns-1,j,k)
967 vtcr(ns,j,k) = vtcr(ns-1,j,k)
968 END DO
970 utcr( 1,ew,k) = utcr( 2,ew-1,k)
971 utcr(ns, 1,k) = utcr(ns-1, 2,k)
972 utcr(ns,ew,k) = utcr(ns-1,ew-1,k)
974 vtcr( 1,ew,k) = vtcr( 2,ew-1,k)
975 vtcr(ns, 1,k) = vtcr(ns-1, 2,k)
976 vtcr(ns,ew,k) = vtcr(ns-1,ew-1,k)
977 END DO
979 !This is the map scale factor on cross points.
980 do i = 1,ns
981 do j = 1,ew
982 mfx(i,j) = grid%msft(j,i)
983 mfd(i,j) = grid%msft(j,i)
984 end do
985 end do
986 call crs2dot(mfd,ns,ew)
989 ! Compute vorticity of FG
990 CALL vor(u1,v1,mfd,mfx,ns,ew,kx,dx,vort)
993 ! Compute divergence of FG
994 CALL diverg(u1,v1,mfd,mfx,ns,ew,kx,dx,div)
997 ! Compute mixing ratio of FG
998 CALL mxratprs(rh1,t1,p*100.,ns,ew,kx,q1,min_RH_value)
999 q1(:,:,1) = q1(:,:,2)
1001 DO k=1,kx
1002 CALL expand ( q1(1,1,k) , ns-1 , ew-1 , ns , ew )
1003 END DO
1005 ! Compute initial streamfunction - PSI1
1006 vortsv = vort
1007 q0 = q1
1010 ! Solve for streamfunction.
1011 DO k=1,kx
1012 DO i=1,ns
1013 DO j=1,ew
1014 ff(i,j) = vort(i,j,k)
1015 tmp1(i,j)= 0.0
1016 END DO
1017 END DO
1018 epsilon = 1.E-2
1019 CALL relax(tmp1,ff,rd,ns,ew,dx,epsilon,alphar)
1020 DO i=1,ns
1021 DO j=1,ew
1022 psi1(i,j,k) = tmp1(i,j)
1023 END DO
1024 END DO
1025 END DO
1026 print *,"after relax one"
1028 DO k=1,kx !start of the k loop
1029 IF ( latc_loc .GE. 0. ) THEN
1030 vormx = -1.e10
1031 ELSE
1032 vormx = 1.e10
1033 END IF
1035 i_mvc = 1
1036 j_mvc = 1
1038 DO i=1,ns
1039 DO j=1,ew
1040 rad = SQRT((REAL(i)-xico)**2.+(REAL(j)-xjco)**2.)*dx
1041 IF ( rad .LE. r_search ) THEN
1042 IF ( latc_loc .GE. 0. ) THEN
1043 IF ( vortsv(i,j,k) .GT. vormx ) THEN
1044 vormx = vortsv(i,j,k)
1045 i_mvc = i
1046 j_mvc = j
1047 END IF
1048 ELSE IF (latc_loc .LT. 0. ) THEN
1049 IF ( vortsv(i,j,k) .LT. vormx ) THEN
1050 vormx = vortsv(i,j,k)
1051 i_mvc = i
1052 j_mvc = j
1053 END IF
1054 END IF
1055 END IF
1056 END DO
1057 END DO
1058 strmci(k) = i_mvc
1059 strmcj(k) = j_mvc
1061 DO i=1,ns
1062 DO j=1,ew
1063 rad = SQRT(REAL((i-i_mvc)**2.+(j-j_mvc)**2.))*dx
1064 IF ( rad .GT. r_vor ) THEN
1065 vort(i,j,k) = 0.
1066 div(i,j,k) = 0.
1067 END IF
1068 END DO
1069 END DO
1071 DO itr=1,n_iter
1072 sum_q = 0.
1073 nct = 0
1074 DO i=1,ns
1075 DO j=1,ew
1076 rad = SQRT(REAL(i-i_mvc)**2.+REAL(j-j_mvc)**2.)*dx
1077 IF ( (rad .LT. r_vor2).AND.(rad .GE. 0.8*r_vor2) ) THEN
1078 sum_q = sum_q + q0(i,j,k)
1079 nct = nct + 1
1080 END IF
1081 END DO
1082 END DO
1083 avg_q = sum_q/MAX(REAL(nct),1.)
1085 DO i=1,ns
1086 DO j=1,ew
1087 q_old = q0(i,j,k)
1088 rad = SQRT(REAL(i-i_mvc)**2.+REAL(j-j_mvc)**2.)*dx
1089 IF ( rad .LT. r_vor2 ) THEN
1090 ror = rad/r_vor2
1091 q_new = ((1.-ror)*avg_q) + (ror*q_old)
1092 q0(i,j,k) = q_new
1093 END IF
1094 END DO
1095 END DO
1096 END DO !end of itr loop
1097 END DO !of the k loop
1100 ! Compute divergent wind
1101 DO k=1,kx
1102 DO i=1,ns
1103 DO j=1,ew
1104 ff(i,j) = div(i,j,k)
1105 tmp1(i,j)= 0.0
1106 END DO
1107 END DO
1108 epsilon = 1.e-2
1109 ! epsilon = 1.E-5
1110 CALL relax(tmp1,ff,rd,ns,ew,dx,epsilon,alphar)
1111 DO i=1,ns
1112 DO j=1,ew
1113 chi(i,j,k) = tmp1(i,j)
1114 END DO
1115 END DO
1116 END DO !of the k loop for divergent winds
1118 DO k=1,kx !start of k loop
1119 DO i=2,ns-1
1120 DO j=2,ew-1
1121 uchi(i,j,k) = ((chi(i ,j ,k)+chi(i-1,j ,k))- (chi(i ,j-1,k)+chi(i-1,j-1,k)))/(2.*dx)
1122 vchi(i,j,k) = ((chi(i ,j ,k)+chi(i ,j-1,k))- (chi(i-1,j ,k)+chi(i-1,j-1,k)))/(2.*dx)
1123 END DO
1124 END DO
1126 DO i=2,ns-1
1127 uchi(i,1,k) = (chi(i , 2,k)-chi(i , 1,k))/(2.*dx)
1128 uchi(i,ew,k) = (chi(i ,ew ,k)-chi(i ,ew-1,k))/(2.*dx)
1129 vchi(i,1,k) = (chi(i , 2,k)-chi(i -1, 2,k))/(2.*dx)
1130 vchi(i,ew,k) = (chi(i ,ew-1,k)-chi(i -1,ew-1,k))/(2.*dx)
1131 END DO
1132 DO j=2,ew-1
1133 uchi(1,j,k) = (chi( 2,j ,k)-chi( 2,j -1,k))/(2.*dx)
1134 uchi(ns,j,k) = (chi(ns-1,j ,k)-chi(ns-1,j -1,k))/(2.*dx)
1135 vchi(1,j,k) = (chi( 2,j -1,k)-chi( 1,j -1,k))/(2.*dx)
1136 vchi(ns,j,k) = (chi(ns ,j -1,k)-chi(ns-1,j -1,k))/(2.*dx)
1137 END DO
1139 uchi( 1, 1,k) = uchi( 2, 2,k)
1140 uchi( 1,ew,k) = uchi( 2,ew-1,k)
1141 uchi(ns, 1,k) = uchi(ns-1, 2,k)
1142 uchi(ns,ns,k) = uchi(ns-1,ew-1,k)
1144 vchi( 1, 1,k) = vchi( 2, 2,k)
1145 vchi( 1,ew,k) = vchi( 2,ew-1,k)
1146 vchi(ns, 1,k) = vchi(ns-1, 2,k)
1147 vchi(ns,ew,k) = vchi(ns-1,ew-1,k)
1148 END DO !end of k loop
1151 ! Compute background streamfunction (PSI0) and perturbation field (PSI)
1152 DO k=1,kx
1153 DO i=1,ns
1154 DO j=1,ew
1155 ff(i,j)=vort(i,j,k)
1156 tmp1(i,j)=0.0
1157 END DO
1158 END DO
1159 epsilon = 1.e-2
1160 ! epsilon = 1.E-5
1161 CALL relax(tmp1,ff,rd,ns,ew,dx,epsilon,alphar)
1162 DO i=1,ns
1163 DO j=1,ew
1164 psi(i,j,k)=tmp1(i,j)
1165 END DO
1166 END DO
1167 END DO
1169 DO k=1,kx
1170 DO i=2,ns-1
1171 DO j=2,ew-1
1172 psi0(i,j,k) = psi1(i,j,k)-psi(i,j,k)
1173 END DO
1174 END DO
1175 END DO
1177 DO k=k00,kx
1178 DO i=1,ns
1179 DO j=1,ew
1180 psipos(i,j,k)=psi(i,j,k)
1181 END DO
1182 END DO
1183 END DO
1186 DO k=1,kx
1187 DO i=2,ns-1
1188 DO j=2,ew-1
1189 upos(i,j,k) = -((psi(i ,j ,k)+psi(i ,j-1,k))-(psi(i-1,j-1,k)+psi(i-1,j ,k)))/(2.*dx)
1190 vpos(i,j,k) = ((psi(i ,j ,k)+psi(i-1,j ,k))-(psi(i-1,j-1,k)+psi(i ,j-1,k)))/(2.*dx)
1191 END DO
1192 END DO
1194 DO i=2,ns-1
1195 upos(i,1,k) = -(psi(i , 2,k)-psi(i -1, 2,k))/(2.*dx)
1196 upos(i,ew,k) = -(psi(i ,ew-1,k)-psi(i -1,ew-1,k))/(2.*dx)
1197 vpos(i,1,k) = (psi(i , 2,k)-psi(i , 1,k))/(2.*dx)
1198 vpos(i,ew,k) = (psi(i ,ew ,k)-psi(i ,ew-1,k))/(2.*dx)
1199 END DO
1200 DO j=2,ew-1
1201 upos(1,j,k) = -(psi( 2,j ,k)-psi( 1,j ,k))/(2.*dx)
1202 upos(ns,j,k) = -(psi(ns ,j ,k)-psi(ns-1,j ,k))/(2.*dx)
1203 vpos(1,j,k) = (psi( 2,j ,k)-psi( 2,j -1,k))/(2.*dx)
1204 vpos(ns,j,k) = (psi(ns-1,j ,k)-psi(ns-1,j -1,k))/(2.*dx)
1205 END DO
1207 upos( 1, 1,k) = upos( 2, 2,k)
1208 upos( 1,ew,k) = upos( 2,ew-1,k)
1209 upos(ns, 1,k) = upos(ns-1, 2,k)
1210 upos(ns,ew,k) = upos(ns-1,ew-1,k)
1212 vpos( 1, 1,k) = vpos( 2, 2,k)
1213 vpos( 1,ew,k) = vpos( 2,ew-1,k)
1214 vpos(ns, 1,k) = vpos(ns-1, 2,k)
1215 vpos(ns,ew,k) = vpos(ns-1,ew-1,k)
1217 END DO
1219 ! Background u, v fields.
1220 DO k=1,kx
1221 DO i=1,ns
1222 DO j=1,ew
1223 u0(i,j,k) = u1(i,j,k)-(upos(i,j,k)+uchi(i,j,k))
1224 v0(i,j,k) = v1(i,j,k)-(vpos(i,j,k)+vchi(i,j,k))
1225 END DO
1226 END DO
1227 END DO
1230 ! Geostrophic vorticity.
1231 CALL geowind(phi1,mfx,mfd,cor,ns,ew,kx,dx,ug,vg)
1232 CALL vor(ug,vg,mfd,mfx,ns,ew,kx,dx,vorg)
1235 DO k=1,kx
1236 i_mvc = strmci(k)
1237 j_mvc = strmcj(k)
1239 DO i=1,ns
1240 DO j=1,ew
1241 rad = SQRT(REAL(i-i_mvc)**2.+REAL(j-j_mvc)**2.)*dx
1242 IF ( rad .GT. r_vor ) THEN
1243 vorg(i,j,k) = 0.
1244 END IF
1245 END DO
1246 END DO
1247 END DO
1249 DO k=k00,kx
1250 DO i=1,ns
1251 DO j=1,ew
1252 ff(i,j) = vorg(i,j,k)
1253 tmp1(i,j)= 0.0
1254 END DO
1255 END DO
1256 epsilon = 1.e-3
1257 CALL relax(tmp1,ff,rd,ns,ew,dx,epsilon,alphar)
1258 DO i=1,ns
1259 DO j=1,ew
1260 phip(i,j,k) = tmp1(i,j)
1261 END DO
1262 END DO
1263 END DO
1265 print *,"the background geopotential"
1266 ! Background geopotential.
1267 DO k=k00,kx
1268 DO i=1,ns
1269 DO j=1,ew
1270 phi0(i,j,k) = phi1(i,j,k) - phip(i,j,k)
1271 ! print *,i,j,k,phi0(i,j,k)
1272 END DO
1273 END DO
1274 END DO
1276 print *,"the background temperature"
1277 ! Background temperature
1278 DO k=k00,kx
1279 DO i=1,ns
1280 DO j=1,ew
1281 IF( k .EQ. 2 ) THEN
1282 tpos(i,j,k) = (-1./rconst)*(phip(i,j,k+1)-phip(i,j,k ))/LOG(p(k+1)/p(k ))
1283 ELSE IF ( k .EQ. kx ) THEN
1284 tpos(i,j,k) = (-1./rconst)*(phip(i,j,k )-phip(i,j,k-1))/LOG(p(k )/p(k-1))
1285 ELSE
1286 tpos(i,j,k) = (-1./rconst)*(phip(i,j,k+1)-phip(i,j,k-1))/LOG(p(k+1)/p(k-1))
1287 END IF
1289 t0(i,j,k) = t1(i,j,k)-tpos(i,j,k)
1290 t00(i,j,k) = t0(i,j,k)
1291 END DO
1292 END DO
1293 END DO
1295 print *,"the background rh"
1296 ! Background RH.
1297 CALL qvtorh (q0,t0,p*100.,k00,ix,jx,kx,rh0,min_RH_value)
1299 DO k=k00,kx
1300 CALL expand ( rh0(1,1,k) , ix-1 , jx-1 , ix , jx )
1301 END DO
1303 DO k=k00,kx
1304 rh_max= rhmx(k)
1305 i_mvc = strmci(k)
1306 j_mvc = strmcj(k)
1308 sum_rh = 0.
1309 nct = 0
1310 DO i=1,ns
1311 DO j=1,ew
1312 rad = SQRT(REAL(i-i_mvc)**2.+REAL(j-j_mvc)**2.)*dx
1313 IF ( (rad .LT. r_vor2).AND.(rad .GE. 0.8*r_vor2) ) THEN
1314 sum_rh = sum_rh + rh0(i,j,k)
1315 nct = nct + 1
1316 END IF
1317 END DO
1318 END DO
1319 avg_rh = sum_rh/MAX(REAL(nct),1.)
1321 DO i=1,ns
1322 DO j=1,ew
1323 rh_min = avg_rh
1324 rad = SQRT((REAL(i)-xico)**2.+(REAL(j)-xjco)**2.)*dx
1325 IF ( rad .LE. rmax ) THEN
1326 rhtc(i,j,k) = rh_max
1327 ELSE
1328 rhtc(i,j,k) = (rmax/rad)*rh_max+(1.-(rmax/rad))*rh_min
1329 END IF
1330 END DO
1331 END DO
1332 END DO
1334 ! adjust T0
1335 DO k=k00,kx
1336 DO i=1,ns
1337 DO j=1,ew
1338 theta(i,j,k) = t1(i,j,k)*(1000./p(k))**rovcp
1339 END DO
1340 END DO
1341 END DO
1343 i_mvc = strmci(k00)
1344 j_mvc = strmcj(k00)
1347 DO k=kfrm,kto
1348 DO i=1,ns
1349 DO j=1,ew
1350 rad = SQRT(REAL(i-i_mvc)**2.+REAL(j-j_mvc)**2.)*dx
1351 IF ( rad .LT. r_vor2 ) THEN
1352 t_reduce(i,j,k) = theta(i,j,k85)-0.03*(p(k)-p(k85))
1353 t0(i,j,k) = t00(i,j,k)*(rad/r_vor2) + (((p(k)/1000.)**rovcp)*t_reduce(i,j,k))*(1.-(rad/r_vor2))
1354 END IF
1355 END DO
1356 END DO
1357 END DO
1359 ! New RH.
1360 DO k=k00,kx
1361 DO i=1,ns
1362 DO j=1,ew
1363 rhbkg = rh0(i,j,k)
1364 rhbog = rhtc(i,j,k)
1365 rad = SQRT((REAL(i)-xico)**2.+(REAL(j)-xjco)**2.)*dx
1366 IF ( (rad.GT.rmax) .AND. (rad.LE.r_vor2) ) THEN
1367 r_ratio = (rad-rmax)/(r_vor2-rmax)
1368 rh2(i,j,k) = ((1.-r_ratio)*rhbog) + (r_ratio*rhbkg)
1369 ELSE IF (rad .LE. rmax ) THEN
1370 rh2(i,j,k) = rhbog
1371 ELSE
1372 rh2(i,j,k) = rhbkg
1373 END IF
1374 END DO
1375 END DO
1376 END DO
1379 ! New wind field.
1380 DO k=1,kx
1381 DO i=1,ns
1382 DO j=1,ew
1383 u2(i,j,k) = u0(i,j,k)+utcr(i,j,k)
1384 v2(i,j,k) = v0(i,j,k)+vtcr(i,j,k)
1385 END DO
1386 END DO
1387 END DO
1389 ! Geopotential perturbation
1390 DO k=1,kx
1391 DO i=1,ns
1392 DO j=1,ew
1393 tmp1(i,j)=psitc(i,j,k)
1394 END DO
1395 END DO
1396 ! print *,"before balance ",dx
1397 CALL balance(cor,tmp1,ns,ew,dx,outold)
1398 DO i=1,ns
1399 DO j=1,ew
1400 ff(i,j)=outold(i,j)
1401 print *,"ff ",i,j,k,ff(i,j)
1402 tmp1(i,j)=0.0
1403 END DO
1404 END DO
1405 epsilon = 1.e-3
1406 CALL relax (tmp1,ff,rd,ns,ew,dx,epsilon,alphar)
1407 print *,"after the geop relax",k
1408 DO i=1,ns
1409 DO j=1,ew
1410 phiptc(i,j,k) = tmp1(i,j)
1411 print *,i,j,phiptc(i,j,k)
1412 END DO
1413 END DO
1414 END DO
1415 print *,"after the geopotential perturbation"
1416 stop
1417 ! New geopotential field.
1418 DO k=1,kx
1419 DO i=1,ns
1420 DO j=1,ew
1421 phi2(i,j,k) = phi0(i,j,k) + phiptc(i,j,k)
1422 END DO
1423 END DO
1424 END DO
1426 print *,"new temperature field"
1427 ! New temperature field.
1428 DO k=k00,kx
1429 DO i=1,ns
1430 DO j=1,ew
1431 IF( k .EQ. 2 ) THEN
1432 tptc(i,j,k)=(-1./rconst)*(phiptc(i,j,k+1)-phiptc(i,j,k ))/LOG(p(k+1)/p(k ))
1433 ELSE IF ( k .EQ. kx ) THEN
1434 tptc(i,j,k)=(-1./rconst)*(phiptc(i,j,k )-phiptc(i,j,k-1))/LOG(p(k )/p(k-1))
1435 ELSE
1436 tptc(i,j,k)=(-1./rconst)*(phiptc(i,j,k+1)-phiptc(i,j,k-1))/LOG(p(k+1)/p(k-1))
1437 END IF
1438 t2(i,j,k) = t0(i,j,k) + tptc(i,j,k)
1439 print *,i,j,k,t2(i,j,k)
1440 END DO
1441 END DO
1442 END DO
1445 ! Surface pressure change.
1446 DO i=1,ns
1447 DO j=1,ew
1448 dph = phi2(i,j,k00)-phi1(i,j,k00)
1449 delpx(i,j) = rho*dph*0.01
1450 END DO
1451 END DO
1453 ! New SLP.
1454 DO i=1,ns
1455 DO j=1,ew
1456 pslx(i,j) = pslx(i,j)+delpx(i,j)
1457 END DO
1458 END DO
1460 ! tmp1=delpx
1461 ! CALL crs2dot(tmp1,ix,jx)
1462 ! DO i=1,ix
1463 ! DO j=1,jx
1464 ! psld(i,j) = psld(i,j)+tmp1(i,j)
1465 ! END DO
1466 ! END DO
1469 ! Set new geopotential at surface to terrain elevation.
1471 DO i=1,ns-1
1472 DO j=1,ew-1
1473 z2(i,j,1) = terrain(i,j)
1474 END DO
1475 END DO
1476 CALL expand(z2(1,1,1),ix-1,jx-1,ix,jx)
1478 ! Geopotential back to height.
1480 DO k=k00,kx
1481 DO i=1,ns
1482 DO j=1,ew
1483 z2(i,j,k) = phi2(i,j,k)/9.81
1484 END DO
1485 END DO
1486 CALL expand(z2(1,1,k),ix-1,jx-1,ix,jx)
1487 END DO
1490 ! New surface temperature, assuming same theta as from 1000 mb.
1491 print *,"surface temperature"
1492 DO i=1,ns
1493 DO j=1,ew
1494 ps = pslx(i,j)
1495 t2(i,j,1) = t2(i,j,k00)*((ps/1000.)**rovcp)
1496 print *,i,j,t2(i,j,1),ps
1497 END DO
1498 END DO
1499 stop
1501 ! Set surface RH to the value from 1000 mb.
1502 DO i=1,ns
1503 DO j=1,ew
1504 rh2(i,j,1) = rh2(i,j,k00)
1505 END DO
1506 END DO
1508 ! Modification of tropical storm complete.
1509 PRINT '(A,I3,A)' ,' Bogus storm number ',nstrm,' completed.'
1512 deallocate(u11)
1513 deallocate(v11)
1514 deallocate(t11)
1515 deallocate(rh11)
1516 deallocate(phi11)
1517 deallocate(u1)
1518 deallocate(v1)
1519 deallocate(t1)
1520 deallocate(rh1)
1521 deallocate(phi1)
1524 do i = 1,ns
1525 do j = 1,ew
1526 do k = 1,nz
1527 grid%u_gc(j,k,i) = u2(i,j,k)
1528 grid%v_gc(j,k,i) = v2(i,j,k)
1529 grid%t_gc(j,k,i) = t2(i,j,k)
1530 print *,i,j,k,t2(i,j,k)
1531 grid%rh_gc(j,k,i) = rh2(i,j,k)
1532 grid%ght_gc(j,k,i) = z2(i,j,k)
1533 end do
1534 end do
1535 end do
1537 END SUBROUTINE tc_bogus
1539 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1540 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1542 SUBROUTINE expand(slab,istart,jstart,itot,jtot)
1544 ! Fill the nearest data to the empty rows or columns of a slab
1546 IMPLICIT NONE
1548 INTEGER :: istart, jstart, itot, jtot
1549 REAL , DIMENSION(itot,jtot) :: slab
1551 INTEGER :: i1, j1, i, j
1553 DO j1=jstart,jtot-1
1554 DO i=1,istart
1555 slab(i,j1+1)=slab(i,j1)
1556 END DO
1557 END DO
1559 DO i1=istart,itot-1
1560 DO j=1,jtot
1561 slab(i1+1,j)=slab(i1,j)
1562 END DO
1563 END DO
1565 END SUBROUTINE expand
1566 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1567 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1569 SUBROUTINE rankine(dx,dy,ds,nlvl,vwgt,rmax,vmax,uu,vv,psi,vor)
1571 ! Define analytical bogus vortex
1573 IMPLICIT NONE
1575 INTEGER nlvl
1576 REAL , DIMENSION(nlvl) :: uu, vv, psi, vor
1577 REAL , DIMENSION(nlvl) :: vwgt
1578 REAL :: dx,dy,ds,rmax,vmax
1580 REAL , PARAMETER :: alpha1= 1.
1581 REAL , PARAMETER :: alpha2= -0.75
1582 real :: pi
1585 INTEGER :: k
1586 REAL :: vr , ang , rr , term1 , bb , term2 , alpha
1589 pi = 3.141592653589793
1590 ! Wind component
1592 DO k=1,nlvl
1593 rr = SQRT(dx**2+dy**2)*ds
1594 IF ( rr .LT. rmax ) THEN
1595 alpha = 1.
1596 ELSE IF ( rr .GE. rmax ) THEN
1597 alpha = alpha2
1598 END IF
1599 vr = vmax * (rr/rmax)**(alpha)
1600 IF ( dx.GE.0. ) THEN
1601 ang = (pi/2.) - ATAN2(dy,MAX(dx,1.e-6))
1602 uu(k) = vwgt(k)*(-vr*COS(ang))
1603 vv(k) = vwgt(k)*( vr*SIN(ang))
1604 ELSE IF ( dx.LT.0. ) THEN
1605 ang = ((3.*pi)/2.) + ATAN2(dy,dx)
1606 uu(k) = vwgt(k)*(-vr*COS(ang))
1607 vv(k) = vwgt(k)*(-vr*SIN(ang))
1608 END IF
1609 END DO
1611 ! psi
1613 DO k=1,nlvl
1614 rr = SQRT(dx**2+dy**2)*ds
1615 IF ( rr .LT. rmax ) THEN
1616 psi(k) = vwgt(k) * (vmax*rr*rr)/(2.*rmax)
1617 ELSE IF ( rr .GE. rmax ) THEN
1618 IF (alpha1.EQ.1.0 .AND. alpha2.eq.-1.0) THEN
1619 psi(k) = vwgt(k) * vmax*rmax*(0.5+LOG(rr/rmax))
1620 ELSE IF (alpha1.EQ.1.0 .AND. alpha2.NE.-1.0) THEN
1621 term1 = vmax/(rmax**alpha1)*(rmax**(alpha1+1)/(alpha1+1))
1622 bb = (rr**(alpha2+1)/(alpha2+1))-(rmax**(alpha2+1))/(alpha2+1)
1623 term2 = vmax/(rmax**alpha2)*bb
1624 psi(k) = vwgt(k) * (term1 + term2)
1625 END IF
1626 END IF
1627 END DO
1629 ! vort
1631 DO k=1,nlvl
1632 rr = SQRT(dx**2+dy**2)*ds
1633 IF ( rr .LT. rmax ) THEN
1634 vor(k) = vwgt(k) * (2.*vmax)/rmax
1635 ELSE IF ( rr .GE. rmax ) THEN
1636 vor(k) = vwgt(k) * ( (vmax/rmax**alpha2)*(rr**(alpha2-1.))*(1.+alpha2) )
1637 END IF
1638 END DO
1640 END SUBROUTINE rankine
1641 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1642 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1644 SUBROUTINE vor(u,v,dmf,xmf,i1,j1,k1,ds,vort)
1646 ! Compute k component of del cross velocity
1647 ! vort = m*m (dv/dx - du/dy), where u and v are coupled
1648 ! with map factors (dot point) and m is the map factors
1649 ! on cross points
1651 IMPLICIT NONE
1653 INTEGER :: i1 , j1 , k1
1655 REAL , DIMENSION(i1,j1,k1) :: u, v, vort
1656 REAL , DIMENSION(i1,j1 ) :: xmf, dmf
1658 REAL :: ds
1660 REAL :: ds2r , u1 , u2 , u3 , u4 , v1 , v2 , v3 , v4
1661 INTEGER :: i , j , k
1663 ds2r=1./(2.*ds)
1665 DO k=1,k1
1666 DO j=1,j1-1
1667 DO i=1,i1-1
1668 u1=u(i ,j ,k)/dmf(i ,j )
1669 u2=u(i+1,j ,k)/dmf(i+1,j )
1670 u3=u(i ,j+1,k)/dmf(i ,j+1)
1671 u4=u(i+1,j+1,k)/dmf(i+1,j+1)
1672 v1=v(i ,j ,k)/dmf(i ,j )
1673 v2=v(i+1,j ,k)/dmf(i+1,j )
1674 v3=v(i ,j+1,k)/dmf(i ,j+1)
1675 v4=v(i+1,j+1,k)/dmf(i+1,j+1)
1676 vort(i,j,k)=xmf(i,j)*xmf(i,j)*ds2r*((v4-v2+v3-v1)-(u2-u1+u4-u3))
1677 print *,i,j,k,vort(i,j,k)
1678 END DO
1679 END DO
1680 END DO
1682 CALL fillit(vort,i1,j1,k1,i1,j1,1,i1-1,1,j1-1)
1684 END SUBROUTINE vor
1686 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1687 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1689 SUBROUTINE diverg(u,v,dmf,xmf,i1,j1,k1,ds,div)
1691 ! Computes divergence
1692 ! div = m*m (du/dx + dv/dy)
1694 IMPLICIT NONE
1696 INTEGER :: i1, j1 , k1
1698 REAL , DIMENSION(i1,j1,k1) :: u, v, div
1699 REAL , DIMENSION(i1,j1 ) :: xmf, dmf
1700 REAL :: ds
1702 REAL :: ds2r , u1 , u2 , u3 , u4 , v1 , v2 , v3 , v4
1703 INTEGER :: i , j , k
1705 ds2r = 1./(2.*ds)
1707 DO k = 1, k1
1708 DO j = 1, j1-1
1709 DO i = 1, i1-1
1710 u1=u(i ,j ,k)/dmf(i ,j )
1711 u2=u(i+1,j ,k)/dmf(i+1,j )
1712 u3=u(i ,j+1,k)/dmf(i ,j+1)
1713 u4=u(i+1,j+1,k)/dmf(i+1,j+1)
1714 v1=v(i ,j ,k)/dmf(i ,j )
1715 v2=v(i+1,j ,k)/dmf(i+1,j )
1716 v2=v(i+1,j ,k)/dmf(i+1,j )
1717 v3=v(i ,j+1,k)/dmf(i ,j+1)
1718 v4=v(i+1,j+1,k)/dmf(i+1,j+1)
1719 div(i,j,k) = xmf(i,j)*xmf(i,j)*ds2r*((u3-u1+u4-u2)+(v2-v1+v4-v3))
1720 END DO
1721 END DO
1722 END DO
1723 END SUBROUTINE diverg
1725 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1726 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1728 SUBROUTINE mxratprs (rh, t, ppa, ix, jx, kx, q, min_RH_value)
1731 IMPLICIT NONE
1733 INTEGER :: i , ix , j , jx , k , kx
1736 REAL :: min_RH_value
1737 REAL :: ppa(kx)
1738 REAL :: p( KX )
1739 REAL :: q (ix,jx,kx),rh(ix,jx,kx),t(ix,jx,kx)
1741 REAL :: es
1742 REAL :: qs
1743 REAL :: cp = 1004.0
1744 REAL :: svp1,svp2,svp3
1745 REAL :: celkel
1746 REAL :: eps
1749 ! This function is designed to compute (q) from basic variables
1750 ! p (mb), t(K) and rh(0-100%) to give (q) in (kg/kg).
1753 p = ppa * 0.01
1755 DO k = 1, kx
1756 DO j = 1, jx - 1
1757 DO i = 1, ix - 1
1758 rh(i,j,k) = MIN ( MAX ( rh(i,j,k) ,min_RH_value ) , 100. )
1759 END DO
1760 END DO
1761 END DO
1763 rh(ix,:,:) = rh(ix-1,:,:)
1764 rh(:,jx,:) = rh(:,jx-1,:)
1765 t (ix,:,:) = t (ix-1,:,:)
1766 t (:,jx,:) = t (:,jx-1,:)
1768 svp3 = 29.65
1769 svp1 = 0.6112
1770 svp2 = 17.67
1771 celkel = 273.15
1772 eps = 0.622
1773 DO k = 1, kx
1774 DO j = 1, jx
1775 DO i = 1, ix
1776 es = svp1 * 10. * EXP(svp2 * (t(i,j,k) - celkel ) / (t(i,j,k) - svp3 ))
1777 qs = eps * es / (p(k) - es)
1778 q(i,j,k) = MAX(0.01 * rh(i,j,k) * qs,0.0)
1779 END DO
1780 END DO
1781 END DO
1783 END SUBROUTINE mxratprs
1785 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1786 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1788 SUBROUTINE fillit (f, ix, jx, kx, imx, jmx, ifirst, ilast, jfirst, jlast)
1790 IMPLICIT NONE
1792 INTEGER :: i
1793 INTEGER :: ifirst
1794 INTEGER :: ilast
1795 INTEGER :: imx
1796 INTEGER :: ix
1797 INTEGER :: j
1798 INTEGER :: jfirst
1799 INTEGER :: jlast
1800 INTEGER :: jmx
1801 INTEGER :: jx
1802 INTEGER :: k
1803 INTEGER :: kx
1805 REAL , dimension(ix,jx,kx) :: f
1807 DO k = 1 , kx
1808 DO j = jfirst, jlast
1809 DO i = 1, ifirst - 1
1810 f(i,j,k) = f(ifirst,j,k)
1811 END DO
1812 DO i = ilast + 1, imx
1813 f(i,j,k) = f(ilast,j,k)
1814 END DO
1815 END DO
1817 DO j = 1, jfirst - 1
1818 f(:,j,k) = f(:,jfirst,k)
1819 END DO
1820 DO j = jlast + 1, jmx
1821 f(:,j,k) = f(:,jlast,k)
1822 END DO
1823 END DO
1825 END SUBROUTINE fillit
1826 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1827 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1828 SUBROUTINE crs2dot(field,dim1,dim2)
1830 IMPLICIT NONE
1832 INTEGER :: dim1 , dim2
1833 REAL , DIMENSION(dim1,dim2) :: field,dummy
1834 INTEGER :: i , j
1836 dummy(2:dim1-1,2:dim2-1) = ( field(1:dim1-2,1:dim2-2) + &
1837 field(1:dim1-2,2:dim2-1) + &
1838 field(2:dim1-1,1:dim2-2) + &
1839 field(2:dim1-1,2:dim2-1) ) * 0.25
1841 dummy(2:dim1-1,1:dim2:dim2-1) = ( field(1:dim1-2,1:dim2-1:dim2-2) + &
1842 field(2:dim1-1,1:dim2-1:dim2-2) ) * 0.5
1844 dummy(1:dim1:dim1-1,2:dim2-1) = ( field(1:dim1-1:dim1-2,1:dim2-2) + &
1845 field(1:dim1-1:dim1-2,2:dim2-1) ) * 0.5
1847 dummy(1:dim1:dim1-1,1:dim2:dim2-1) = field(1:dim1-1:dim1-2,1:dim2-1:dim2-2)
1849 field = dummy
1851 END SUBROUTINE crs2dot
1852 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1853 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1855 SUBROUTINE relax (chi, ff, rd, imx, jmx, ds, smallres, alpha)
1857 IMPLICIT NONE
1859 INTEGER, PARAMETER :: mm = 20000
1861 INTEGER :: i
1862 INTEGER :: ie
1863 INTEGER :: imx
1864 INTEGER :: iter
1865 INTEGER :: j
1866 INTEGER :: je
1867 INTEGER :: jm
1868 INTEGER :: jmx
1869 INTEGER :: mi
1871 REAL :: alpha
1872 REAL :: alphaov4
1873 REAL :: chi(imx,jmx)
1874 REAL :: chimx( jmx )
1875 REAL :: ds
1876 REAL :: epx
1877 REAL :: fac
1878 REAL :: ff(imx,jmx)
1879 REAL :: rd(imx,jmx)
1880 REAL :: rdmax( jmx )
1881 REAL :: smallres
1883 LOGICAL :: converged = .FALSE.
1885 fac = ds * ds
1886 alphaov4 = alpha * 0.25
1888 ie=imx-2
1889 je=jmx-2
1891 DO j = 1, jmx
1892 DO i = 1, imx
1893 ff(i,j) = fac * ff(i,j)
1894 rd(i,j) = 0.0
1895 END DO
1896 END DO
1898 iter_loop : DO iter = 1, mm
1899 mi = iter
1900 chimx = 0.0
1903 DO j = 2, je
1904 DO i = 2, ie
1905 chimx(j) = MAX(ABS(chi(i,j)),chimx(j))
1906 END DO
1907 END DO
1909 epx = MAXVAL(chimx) * SMALLRES * 4.0 / alpha
1911 DO j = 2, je
1912 DO i = 2, ie
1913 rd(i,j) = chi(i,j+1) + chi(i,j-1) + chi(i+1,j) + chi(i-1,j) - 4.0 * chi(i,j) - ff(i,j)
1914 chi(i,j) = chi(i,j) + rd(i,j) * alphaov4
1915 END DO
1916 END DO
1918 rdmax = 0.0
1920 DO j = 2, je
1921 DO i = 2, ie
1922 rdmax(j) = MAX(ABS(rd(i,j)),rdmax(j))
1923 END DO
1924 END DO
1926 IF (MAXVAL(rdmax) .lt. epx) THEN
1927 converged = .TRUE.
1928 EXIT iter_loop
1929 END IF
1931 END DO iter_loop
1933 IF (converged ) THEN
1934 ! PRINT '(A,I5,A)','Relaxation converged in ',mi,' iterations.'
1935 ELSE
1936 PRINT '(A,I5,A)','Relaxation did not converge in',mm,' iterations.'
1937 STOP 'no_converge'
1938 END IF
1940 END SUBROUTINE relax
1941 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1942 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1943 SUBROUTINE geowind(height,xmf,dmf,cor,imx,jmx,kx,ds,ug,vg)
1945 ! Computes the geostrophic wind components from the height gradient.
1946 ! There is no Coriolis parameter used - this is the tropics.
1948 IMPLICIT NONE
1950 ! input height geopotential cross 3d
1951 ! xmf map factors cross 2d
1952 ! dmf map factors dot 2d
1953 ! imx dot point dimension n-s
1954 ! jmx dot point dimension e-w
1955 ! kx number of vertical levels
1956 !
1957 ! output ug u component of geo wind cross 3d
1958 ! vg v component of geo wind cross 3D
1960 INTEGER :: imx , jmx , kx
1961 REAL :: ds
1962 REAL , DIMENSION(imx,jmx,kx) :: height
1963 REAL , DIMENSION(imx,jmx ) :: xmf , dmf , cor
1965 REAL , DIMENSION(imx,jmx,kx) :: ug , vg
1967 REAL :: ds2r , h1 , h2 , h3 , h4
1968 INTEGER :: i , j , k
1970 ds2r=1./(2.*ds)
1972 DO k=1,kx
1973 DO j=2,jmx-1
1974 DO i=2,imx-1
1975 h1=height(i-1,j-1,k)
1976 h2=height(i ,j-1,k)
1977 h3=height(i-1,j ,k)
1978 h4=height(i ,j ,k)
1979 ! ug(i,j,k)=-1.*g*dmf(i,j)/cor(i,j)*ds2r*(h4-h3+h2-h1)
1980 ! vg(i,j,k)= g*dmf(i,j)/cor(i,j)*ds2r*(h4-h2+h3-h1)
1981 ug(i,j,k)=-1.*dmf(i,j)*ds2r*(h4-h3+h2-h1)
1982 vg(i,j,k)= dmf(i,j)*ds2r*(h4-h2+h3-h1)
1983 END DO
1984 END DO
1985 END DO
1987 CALL fillit(ug,imx,jmx,kx,imx,jmx,2,imx-1,2,jmx-1)
1988 CALL fillit(vg,imx,jmx,kx,imx,jmx,2,imx-1,2,jmx-1)
1990 END SUBROUTINE geowind
1991 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1992 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1993 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1994 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1996 SUBROUTINE balance (f,psi,ix,jx,ds,out)
1998 ! Calculates the forcing terms in balance equation
2000 IMPLICIT NONE
2002 ! f coriolis force
2003 ! psi stream function
2004 ! ix, jx grid points in east west, north south direction, respectively
2005 ! ds grid distance
2006 ! out output array
2008 INTEGER :: ix , jx
2009 REAL , DIMENSION(ix,jx) :: f,psi,out
2010 REAL :: ds
2012 REAL :: psixx , psiyy , psiy , psixy
2013 REAL :: dssq , ds2 , dssq4
2015 INTEGER :: i , j
2017 dssq = ds * ds
2018 ds2 = ds * 2.
2019 dssq4 = ds * ds * 4.
2021 ! print *,"in balance",dssq,ds2,dssq4,ds
2022 DO i=2,ix-2
2023 DO j=2,jx-2
2024 psixx = ( psi(i,j+1) + psi(i,j-1) - 2.*psi(i,j) ) / dssq
2025 psiyy = ( psi(i+1,j) + psi(i-1,j) - 2.*psi(i,j) ) / dssq
2026 psiy = ( psi(i+1,j) - psi(i-1,j) ) / ds2
2027 psixy = ( psi(i+1,j+1)+psi(i-1,j-1)-psi(i-1,j+1)-psi(i+1,j-1)) / dssq4
2028 ! print *,f(i,j),f(i+1,j),f(i+1,j+1),f(i,j+1),psixx,psiyy,psiy,psixy
2029 out(i,j)=0.25*(f(i,j)+f(i+1,j)+f(i+1,j+1)+f(i,j+1))*(psixx+psiyy) &
2030 +psiy*(f(i+1,j+1)+f(i+1,j)-f(i,j)-f(i,j+1))/ ds2 &
2031 -2.*(psixy*psixy-psixx*psiyy)
2032 END DO
2033 END DO
2035 CALL fill(out,ix,jx,ix,jx,2,ix-2,2,jx-2)
2037 END SUBROUTINE balance
2038 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2039 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2041 SUBROUTINE fill (f, ix, jx, imx, jmx, ifirst, ilast, jfirst, jlast)
2043 IMPLICIT NONE
2045 INTEGER :: I
2046 INTEGER :: IFIRST
2047 INTEGER :: ILAST
2048 INTEGER :: IMX
2049 INTEGER :: IX
2050 INTEGER :: J
2051 INTEGER :: JFIRST
2052 INTEGER :: JLAST
2053 INTEGER :: JMX
2054 INTEGER :: JX
2056 REAL :: F(ix,jx)
2058 DO j = jfirst, jlast
2059 DO i = 1, ifirst - 1
2060 f(i,j) = f(ifirst,j)
2061 END DO
2062 DO i = ilast + 1, imx
2063 f(i,j) = f(ilast,j)
2064 END DO
2065 END DO
2067 DO j = 1, jfirst - 1
2068 f(:,j) = f(:,jfirst)
2069 END DO
2070 DO j = jlast + 1, jmx
2071 f(:,j) = f(:,jlast)
2072 END DO
2074 END SUBROUTINE fill
2075 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2076 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2078 SUBROUTINE qvtorh ( q , t , p , k00, imx , jmx , kxs , rh, min_RH_value )
2080 IMPLICIT NONE
2082 INTEGER , INTENT(IN) :: imx , jmx , kxs , k00
2083 REAL , INTENT(IN) , DIMENSION(imx,jmx,kxs) :: q , t
2084 REAL , INTENT(IN) , DIMENSION(kxs) :: p
2086 REAL , INTENT(OUT) , DIMENSION(imx,jmx,kxs) :: rh
2088 real min_RH_value
2090 ! Local variables.
2092 INTEGER :: i , j , k
2093 REAL :: es
2094 REAL :: qs
2095 REAL :: cp = 1004.0
2096 REAL :: svp1,svp2,svp3
2097 REAL :: celkel
2098 REAL :: eps
2099 svp3 = 29.65
2100 svp1 = 0.6112
2101 svp2 = 17.67
2102 celkel = 273.15
2103 eps = 0.622
2105 DO k = k00 , kxs
2106 DO j = 1 , jmx - 1
2107 DO i = 1 , imx - 1
2108 es = svp1 * 10. * EXP(svp2 * (t(i,j,k) - celkel ) / (t(i,j,k) - svp3 ))
2109 qs = eps*es/(0.01*p(k) - es)
2110 rh(i,j,k) = MIN ( 100. , MAX ( 100.*q(i,j,k)/qs , min_RH_value ) )
2111 END DO
2112 END DO
2113 END DO
2115 END SUBROUTINE qvtorh