!REAL:MODEL_LAYER:INITIALIZATION ! This MODULE holds the routines which are used to perform various initializations ! for individual domains utilizing the NMM dynamical core. !----------------------------------------------------------------------- MODULE module_initialize_ideal USE module_bc USE module_configure USE module_domain USE module_io_domain USE module_model_constants ! USE module_si_io_nmm USE module_state_description USE module_timing USE module_soil_pre #ifdef DM_PARALLEL USE module_dm, ONLY : LOCAL_COMMUNICATOR & ,MYTASK,NTASKS,NTASKS_X & ,NTASKS_Y USE module_comm_dm USE module_ext_internal #endif INTEGER :: internal_time_loop INTEGER:: MPI_COMM_COMP,MYPE INTEGER:: loopinc, iloopinc CONTAINS !------------------------------------------------------------------- SUBROUTINE init_domain ( grid ) IMPLICIT NONE ! Input space and data. No gridded meteorological data has been stored, though. ! TYPE (domain), POINTER :: grid TYPE (domain) :: grid ! Local data. INTEGER :: idum1, idum2 CALL set_scalar_indices_from_config ( head_grid%id , idum1, idum2 ) CALL init_domain_nmm (grid & ! #include "actual_new_args.inc" ! ) END SUBROUTINE init_domain !------------------------------------------------------------------- !--------------------------------------------------------------------- SUBROUTINE init_domain_nmm ( grid & ! # include "dummy_new_args.inc" ! ) USE module_optional_input IMPLICIT NONE ! Input space and data. No gridded meteorological data has been stored, though. ! TYPE (domain), POINTER :: grid TYPE (domain) :: grid # include "dummy_new_decl.inc" real ,DIMENSION(grid%sm31:grid%em31,grid%sm32:grid%em32,1:grid%num_metgrid_levels) :: ght_out real ,DIMENSION(grid%sm31:grid%em31,grid%sm32:grid%em32,1:grid%num_metgrid_levels) :: rh_out real ,DIMENSION(grid%sm31:grid%em31,grid%sm32:grid%em32,1:grid%num_metgrid_levels) :: t_out real ,DIMENSION(grid%sm31:grid%em31,grid%sm32:grid%em32,1:grid%num_metgrid_levels) :: u_out real ,DIMENSION(grid%sm31:grid%em31,grid%sm32:grid%em32,1:grid%num_metgrid_levels) :: v_out real ,DIMENSION(grid%sm31:grid%em31,grid%sm32:grid%em32,1:grid%num_land_cat) :: landusef_out real ,DIMENSION(grid%sm31:grid%em31,grid%sm32:grid%em32,1:grid%num_soil_cat) :: soilcbot_out real ,DIMENSION(grid%sm31:grid%em31,grid%sm32:grid%em32,1:grid%num_soil_cat) :: soilctop_out TYPE (grid_config_rec_type) :: config_flags ! Local domain indices and counters. INTEGER :: num_veg_cat , num_soil_top_cat , num_soil_bot_cat INTEGER :: num_veg_gc , num_soil_top_gc , num_soil_bot_gc INTEGER :: & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte, & ips, ipe, jps, jpe, kps, kpe, & i, j, k, NNXP, NNYP ! Local data CHARACTER(LEN=19):: start_date #ifdef DM_PARALLEL LOGICAL,EXTERNAL :: WRF_DM_ON_MONITOR logical :: test ! INTEGER :: DOMDESC REAL,ALLOCATABLE :: SICE_G(:,:), SM_G(:,:) INTEGER, ALLOCATABLE:: IHE_G(:),IHW_G(:) INTEGER, ALLOCATABLE:: ITEMP(:,:) INTEGER :: my_e,my_n,my_s,my_w,my_ne,my_nw,my_se,my_sw,myi,myj,npe INTEGER :: istat,inpes,jnpes #else integer, allocatable:: ihw(:),ihe(:) #endif CHARACTER (LEN=255) :: message INTEGER :: error REAL :: p_surf, p_level REAL :: cof1, cof2 REAL :: qvf , qvf1 , qvf2 , pd_surf REAL :: p00 , t00 , a REAL :: hold_znw, rmin,rmax REAL :: p_top_requested , ptsgm INTEGER :: num_metgrid_levels, ICOUNT REAL , DIMENSION(max_eta) :: eta_levels LOGICAL :: stretch_grid, dry_sounding, debug, log_flag_sst, hyb_coor REAL, ALLOCATABLE,DIMENSION(:,:):: ADUM2D,SNOWC,HT,TG_ALT, & PDVP,PSFC_OUTV REAL, ALLOCATABLE,DIMENSION(:,:,:):: P3D_OUT,P3DV_OUT,P3DV_IN, & QTMP,QTMP2 INTEGER, ALLOCATABLE, DIMENSION(:):: KHL2,KVL2,KHH2,KVH2, & KHLA,KHHA,KVLA,KVHA ! INTEGER, ALLOCATABLE, DIMENSION(:,:):: grid%lu_index REAL, ALLOCATABLE, DIMENSION(:):: DXJ,WPDARJ,CPGFUJ,CURVJ, & FCPJ,FDIVJ,EMJ,EMTJ,FADJ, & HDACJ,DDMPUJ,DDMPVJ REAL, ALLOCATABLE,DIMENSION(:),SAVE:: SG1,SG2,DSG1,DSG2, & SGML1,SGML2 !-- Carsel and Parrish [1988] REAL , DIMENSION(100) :: lqmi integer iicount REAL:: TPH0D,TLM0D REAL:: TPH0,WB,SB,TDLM,TDPH REAL:: WBI,SBI,EBI,ANBI,STPH0,CTPH0 REAL:: TSPH,DTAD,DTCF REAL:: ACDT,CDDAMP,DXP,FP REAL:: WBD,SBD REAL:: RSNOW,SNOFAC REAL, PARAMETER:: SALP=2.60 REAL, PARAMETER:: SNUP=0.040 REAL:: SMCSUM,STCSUM,SEAICESUM,FISX REAL:: cur_smc, aposs_smc REAL:: COAC, CODAMP INTEGER,PARAMETER:: DOUBLE=SELECTED_REAL_KIND(15,300) REAL(KIND=DOUBLE):: TERM1,APH,TLM,TPH,DLM,DPH,STPH,CTPH INTEGER:: KHH,KVH,JAM,JA, IHL, IHH, L INTEGER:: II,JJ,ISRCH,ISUM,ITER,Ilook,Jlook REAL, PARAMETER:: DTR=0.01745329 REAL, PARAMETER:: W_NMM=0.08 #if ( HWRF == 1 ) REAL, PARAMETER:: DDFC=1.0 #else REAL, PARAMETER:: DDFC=8.0 #endif REAL, PARAMETER:: TWOM=.00014584 REAL, PARAMETER:: CP=1004.6 REAL, PARAMETER:: DFC=1.0 REAL, PARAMETER:: ROI=916.6 REAL, PARAMETER:: R=287.04 REAL, PARAMETER:: CI=2060.0 REAL, PARAMETER:: ROS=1500. REAL, PARAMETER:: CS=1339.2 REAL, PARAMETER:: DS=0.050 REAL, PARAMETER:: AKS=.0000005 REAL, PARAMETER:: DZG=2.85 REAL, PARAMETER:: DI=.1000 REAL, PARAMETER:: AKI=0.000001075 REAL, PARAMETER:: DZI=2.0 REAL, PARAMETER:: THL=210. REAL, PARAMETER:: PLQ=70000. REAL, PARAMETER:: ERAD=6371200. REAL, PARAMETER:: TG0=258.16 REAL, PARAMETER:: TGA=30.0 integer :: numzero,numexamined #if ( HWRF == 1 ) !============================================================================ ! gopal's doing for ocean coupling !============================================================================ REAL, DIMENSION(:,:), ALLOCATABLE :: NHLAT,NHLON,NVLAT,NVLON,HRES_SM REAL :: NDLMD,NDPHD,NWBD,NSBD INTEGER :: NIDE,NJDE,ILOC,JLOC INTEGER fid, ierr, nprocs CHARACTER*255 f65name, SysString !============================================================================ ! end gopal's doing for ocean coupling !============================================================================ #endif if (ALLOCATED(ADUM2D)) DEALLOCATE(ADUM2D) if (ALLOCATED(TG_ALT)) DEALLOCATE(TG_ALT) !#define COPY_IN !#include "scalar_derefs.inc" #ifdef DM_PARALLEL # include "data_calls.inc" #endif SELECT CASE ( model_data_order ) CASE ( DATA_ORDER_ZXY ) kds = grid%sd31 ; kde = grid%ed31 ; ids = grid%sd32 ; ide = grid%ed32 ; jds = grid%sd33 ; jde = grid%ed33 ; kms = grid%sm31 ; kme = grid%em31 ; ims = grid%sm32 ; ime = grid%em32 ; jms = grid%sm33 ; jme = grid%em33 ; kts = grid%sp31 ; kte = grid%ep31 ; ! tile is entire patch its = grid%sp32 ; ite = grid%ep32 ; ! tile is entire patch jts = grid%sp33 ; jte = grid%ep33 ; ! tile is entire patch CASE ( DATA_ORDER_XYZ ) ids = grid%sd31 ; ide = grid%ed31 ; jds = grid%sd32 ; jde = grid%ed32 ; kds = grid%sd33 ; kde = grid%ed33 ; ims = grid%sm31 ; ime = grid%em31 ; jms = grid%sm32 ; jme = grid%em32 ; kms = grid%sm33 ; kme = grid%em33 ; its = grid%sp31 ; ite = grid%ep31 ; ! tile is entire patch jts = grid%sp32 ; jte = grid%ep32 ; ! tile is entire patch kts = grid%sp33 ; kte = grid%ep33 ; ! tile is entire patch CASE ( DATA_ORDER_XZY ) ids = grid%sd31 ; ide = grid%ed31 ; kds = grid%sd32 ; kde = grid%ed32 ; jds = grid%sd33 ; jde = grid%ed33 ; ims = grid%sm31 ; ime = grid%em31 ; kms = grid%sm32 ; kme = grid%em32 ; jms = grid%sm33 ; jme = grid%em33 ; its = grid%sp31 ; ite = grid%ep31 ; ! tile is entire patch kts = grid%sp32 ; kte = grid%ep32 ; ! tile is entire patch jts = grid%sp33 ; jte = grid%ep33 ; ! tile is entire patch END SELECT #ifdef DM_PARALLEL CALL WRF_GET_MYPROC(MYPE) CALL WRF_GET_DM_COMMUNICATOR(MPI_COMM_COMP) call wrf_get_nprocx(inpes) call wrf_get_nprocy(jnpes) ! allocate(itemp(inpes,jnpes),stat=istat) npe=0 ! do j=1,jnpes do i=1,inpes itemp(i,j)=npe if(npe==mype)then myi=i myj=j endif npe=npe+1 enddo enddo ! my_n=-1 if(myj+1<=jnpes)my_n=itemp(myi,myj+1) ! my_e=-1 if(myi+1<=inpes)my_e=itemp(myi+1,myj) ! my_s=-1 if(myj-1>=1)my_s=itemp(myi,myj-1) ! my_w=-1 if(myi-1>=1)my_w=itemp(myi-1,myj) ! my_ne=-1 if((myi+1<=inpes).and.(myj+1<=jnpes)) & my_ne=itemp(myi+1,myj+1) ! my_se=-1 if((myi+1<=inpes).and.(myj-1>=1)) & my_se=itemp(myi+1,myj-1) ! my_sw=-1 if((myi-1>=1).and.(myj-1>=1)) & my_sw=itemp(myi-1,myj-1) ! my_nw=-1 if((myi-1>=1).and.(myj+1<=jnpes)) & my_nw=itemp(myi-1,myj+1) ! ! my_neb(1)=my_n ! my_neb(2)=my_e ! my_neb(3)=my_s ! my_neb(4)=my_w ! my_neb(5)=my_ne ! my_neb(6)=my_se ! my_neb(7)=my_sw ! my_neb(8)=my_nw ! deallocate(itemp) #endif grid%DT=float(grid%TIME_STEP) NNXP=min(ITE,IDE-1) NNYP=min(JTE,JDE-1) write(message,*) 'IDE, JDE: ', IDE,JDE write(message,*) 'NNXP, NNYP: ', NNXP,NNYP CALL wrf_message(message) JAM=6+2*(JDE-JDS-10) if (internal_time_loop .eq. 1) then ALLOCATE(ADUM2D(grid%sm31:grid%em31,jms:jme)) ALLOCATE(KHL2(JTS:NNYP),KVL2(JTS:NNYP),KHH2(JTS:NNYP),KVH2(JTS:NNYP)) ALLOCATE(DXJ(JTS:NNYP),WPDARJ(JTS:NNYP),CPGFUJ(JTS:NNYP),CURVJ(JTS:NNYP)) ALLOCATE(FCPJ(JTS:NNYP),FDIVJ(JTS:NNYP),& FADJ(JTS:NNYP)) ALLOCATE(HDACJ(JTS:NNYP),DDMPUJ(JTS:NNYP),DDMPVJ(JTS:NNYP)) ALLOCATE(KHLA(JAM),KHHA(JAM)) ALLOCATE(KVLA(JAM),KVHA(JAM)) endif CALL model_to_grid_config_rec ( grid%id , model_config_rec , config_flags ) IF ( CONFIG_FLAGS%FRACTIONAL_SEAICE == 1 ) THEN CALL WRF_ERROR_FATAL("NMM cannot use FRACTIONAL_SEAICE = 1") ENDIF if ( config_flags%bl_pbl_physics == BOULACSCHEME ) then call wrf_error_fatal("Cannot use BOULAC PBL with NMM") endif write(message,*) 'cen_lat: ', config_flags%cen_lat CALL wrf_debug(100,message) write(message,*) 'cen_lon: ', config_flags%cen_lon CALL wrf_debug(100,message) write(message,*) 'dx: ', config_flags%dx CALL wrf_debug(100,message) write(message,*) 'dy: ', config_flags%dy CALL wrf_debug(100,message) write(message,*) 'config_flags%start_year: ', config_flags%start_year CALL wrf_debug(100,message) write(message,*) 'config_flags%start_month: ', config_flags%start_month CALL wrf_debug(100,message) write(message,*) 'config_flags%start_day: ', config_flags%start_day CALL wrf_debug(100,message) write(message,*) 'config_flags%start_hour: ', config_flags%start_hour CALL wrf_debug(100,message) write(start_date,435) config_flags%start_year, config_flags%start_month, & config_flags%start_day, config_flags%start_hour 435 format(I4,'-',I2.2,'-',I2.2,'_',I2.2,':00:00') grid%dlmd=config_flags%dx grid%dphd=config_flags%dy tph0d=config_flags%cen_lat tlm0d=config_flags%cen_lon !========================================================================== !! ! Check to see if the boundary conditions are set ! properly in the namelist file. ! This checks for sufficiency and redundancy. CALL boundary_condition_check( config_flags, bdyzone, error, grid%id ) ! Some sort of "this is the first time" initialization. Who knows. grid%itimestep=0 ! Pull in the info in the namelist to compare it to the input data. grid%real_data_init_type = model_config_rec%real_data_init_type ! write(message,*) 'what is flag_metgrid: ', flag_metgrid ! CALL wrf_message(message) IF ( flag_metgrid .EQ. 1 ) THEN ! <----- START OF VERTICAL INTERPOLATION PART ----> num_metgrid_levels = grid%num_metgrid_levels !--------------------------------------------------------------------- ! bug fix for albedo and emissivity DO j = jts, MIN(jte,jde-1) DO i = its, MIN(ite,ide-1) grid%landmask(I,J) = 0. END DO END DO ! end bug fix for albedo and emissivity ! ! gopal's doing for ideal cases ! num_veg_gc = SIZE ( grid%landusef_gc , DIM=3 ) num_soil_top_gc = SIZE ( grid%soilctop_gc , DIM=3 ) num_soil_bot_gc = SIZE ( grid%soilcbot_gc , DIM=3 ) #ifdef DM_PARALLEL ips=its ; ipe=ite ; jps=jts ; jpe=jte ; kps=kts ; kpe=kte !JWB # include "HALO_NMM_IDEAL_1.inc" #endif ! WRITE(message,*)'--------------- ght_gc before calling vortex --------------------------' ! call wrf_debug(1,message) ! WRITE(message,*)grid%ght_gc(100,100,:) CALL vortex ( grid%ght_gc,grid%rh_gc,grid%t_gc,grid%u_gc,grid%v_gc,grid%p_gc & &, ght_out,rh_out,t_out,u_out,v_out & &, grid%ht_gc,grid%tsk_gc,grid%xice_gc & &, grid%hlat_gc,grid%hlon_gc,grid%vlat_gc,grid%vlon_gc & &, grid%greenfrac_gc,grid%albedo12m_gc,grid%landusef_gc & &, grid%soilctop_gc,grid%soilcbot_gc & &, landusef_out,soilctop_out,soilcbot_out & &, num_veg_gc,num_soil_top_gc,num_soil_bot_gc & &, config_flags%dx,internal_time_loop & &, 1,grid%num_metgrid_levels,config_flags%sf_surface_physics & &, ids,ide,jds,jde,kds,kde & &, ims,ime,jms,jme,kms,kme & &, its,ite,jts,jte,kts,kte ) !---------------------------------------------------------------------- IF (grid%ght_gc(its,jts,num_metgrid_levels/2) .lt. grid%ght_gc(its,jts,num_metgrid_levels/2+1)) THEN write(message,*) 'normal ground up file order' hyb_coor=.false. CALL wrf_message(message) ELSE hyb_coor=.true. write(message,*) 'reverse the order of coordinate' CALL wrf_message(message) CALL reverse_vert_coord(grid%ght_gc, 2, num_metgrid_levels & &, IDS,IDE,JDS,JDE,KDS,KDE & &, IMS,IME,JMS,JME,KMS,KME & &, ITS,ITE,JTS,JTE,KTS,KTE ) #if ( HWRF == 1 ) if(.not. grid%use_prep_hybrid) then #endif CALL reverse_vert_coord(grid%p_gc, 2, num_metgrid_levels & &, IDS,IDE,JDS,JDE,KDS,KDE & &, IMS,IME,JMS,JME,KMS,KME & &, ITS,ITE,JTS,JTE,KTS,KTE ) CALL reverse_vert_coord(grid%t_gc, 2, num_metgrid_levels & &, IDS,IDE,JDS,JDE,KDS,KDE & &, IMS,IME,JMS,JME,KMS,KME & &, ITS,ITE,JTS,JTE,KTS,KTE ) CALL reverse_vert_coord(grid%u_gc, 2, num_metgrid_levels & &, IDS,IDE,JDS,JDE,KDS,KDE & &, IMS,IME,JMS,JME,KMS,KME & &, ITS,ITE,JTS,JTE,KTS,KTE ) CALL reverse_vert_coord(grid%v_gc, 2, num_metgrid_levels & &, IDS,IDE,JDS,JDE,KDS,KDE & &, IMS,IME,JMS,JME,KMS,KME & &, ITS,ITE,JTS,JTE,KTS,KTE ) CALL reverse_vert_coord(grid%rh_gc, 2, num_metgrid_levels & &, IDS,IDE,JDS,JDE,KDS,KDE & &, IMS,IME,JMS,JME,KMS,KME & &, ITS,ITE,JTS,JTE,KTS,KTE ) #if ( HWRF == 1 ) endif #endif endif IF (hyb_coor) THEN ! limit extreme deviations from source model topography ! due to potential for nasty extrapolation/interpolation issues ! write(message,*) 'min, max of grid%ht_gc before adjust: ', minval(grid%ht_gc), maxval(grid%ht_gc) CALL wrf_debug(100,message) ICOUNT=0 DO J=JTS,min(JTE,JDE-1) DO I=ITS,min(ITE,IDE-1) IF ((grid%ht_gc(I,J) - grid%ght_gc(I,J,2)) .LT. -150.) THEN grid%ht_gc(I,J)=grid%ght_gc(I,J,2)-150. IF (ICOUNT .LT. 20) THEN write(message,*) 'increasing NMM topo toward RUC ', I,J CALL wrf_debug(100,message) ICOUNT=ICOUNT+1 ENDIF ELSEIF ((grid%ht_gc(I,J) - grid%ght_gc(I,J,2)) .GT. 150.) THEN grid%ht_gc(I,J)=grid%ght_gc(I,J,2)+150. IF (ICOUNT .LT. 20) THEN write(message,*) 'decreasing NMM topo toward RUC ', I,J CALL wrf_debug(100,message) ICOUNT=ICOUNT+1 ENDIF ENDIF END DO END DO write(message,*) 'min, max of ht_gc after correct: ', minval(grid%ht_gc), maxval(grid%ht_gc) CALL wrf_debug(100,message) ENDIF CALL boundary_smooth(grid%ht_gc,grid%landmask, grid, 12 , 12 & &, IDS,IDE,JDS,JDE,KDS,KDE & &, IMS,IME,JMS,JME,KMS,KME & &, ITS,ITE,JTS,JTE,KTS,KTE ) DO j = jts, MIN(jte,jde-1) DO i = its, MIN(ite,ide-1) if (grid%landmask(I,J) .gt. 0.5) grid%sm(I,J)=0. if (grid%landmask(I,J) .le. 0.5) grid%sm(I,J)=1. if (grid%tsk_gc(I,J) .gt. 0.) then grid%nmm_tsk(I,J)=grid%tsk_gc(I,J) else #if ( HWRF == 1 ) if(grid%use_prep_hybrid) then if(grid%t(I,J,1)<100) then write(message,*) 'NO VALID SURFACE TEMPERATURE: I,J,TSK_GC(I,J),T(I,J,1) = ', & I,J,grid%TSK_GC(I,J),grid%T(I,J,1) call wrf_debug(1,message) else grid%nmm_tsk(I,J)=grid%t(I,J,1) ! stopgap measure end if else #endif grid%nmm_tsk(I,J)=grid%t_gc(I,J,1) ! stopgap measure #if ( HWRF == 1 ) endif #endif endif ! grid%glat(I,J)=grid%hlat_gc(I,J)*DEGRAD grid%glon(I,J)=grid%hlon_gc(I,J)*DEGRAD grid%weasd(I,J)=grid%snow(I,J) grid%xice(I,J)=grid%xice_gc(I,J) ENDDO ENDDO ! First item is to define the target vertical coordinate num_metgrid_levels = grid%num_metgrid_levels eta_levels(1:kde) = model_config_rec%eta_levels(1:kde) ptsgm = model_config_rec%ptsgm p_top_requested = grid%p_top_requested grid%pt=p_top_requested if (internal_time_loop .eq. 1) then if (eta_levels(1) .ne. 1.0) then #if ( HWRF == 1 ) if(grid%use_prep_hybrid) then call wrf_error_fatal('PREP_HYBRID ERROR: eta_levels is not specified, but use_prep_hybrid=.true.') end if #endif write(message,*) '********************************************************************* ' CALL wrf_message(message) write(message,*) '** eta_levels appears not to be specified in the namelist' CALL wrf_message(message) write(message,*) '** We will call compute_nmm_levels to define layer thicknesses.' CALL wrf_message(message) write(message,*) '** These levels should be reasonable for running the model, ' CALL wrf_message(message) write(message,*) '** but may not be ideal for the simulation being made. Consider ' CALL wrf_message(message) write(message,*) '** defining your own levels by specifying eta_levels in the model ' CALL wrf_message(message) write(message,*) '** namelist. ' CALL wrf_message(message) write(message,*) '********************************************************************** ' CALL wrf_message(message) CALL compute_nmm_levels(KDE,p_top_requested,eta_levels) DO L=1,KDE write(message,*) 'L, eta_levels(L) returned :: ', L,eta_levels(L) CALL wrf_message(message) ENDDO endif write(message,*) 'KDE-1: ', KDE-1 CALL wrf_debug(1,message) allocate(SG1(1:KDE-1)) allocate(SG2(1:KDE-1)) allocate(DSG1(1:KDE-1)) allocate(DSG2(1:KDE-1)) allocate(SGML1(1:KDE)) allocate(SGML2(1:KDE)) CALL define_nmm_vertical_coord (kde-1, ptsgm, grid%pt,grid%pdtop, eta_levels, & grid%eta1,grid%deta1,grid%aeta1, & grid%eta2,grid%deta2,grid%aeta2, grid%dfl, grid%dfrlg ) DO L=KDS,KDE-1 grid%deta(L)=eta_levels(L)-eta_levels(L+1) ENDDO endif write(message,*) 'num_metgrid_levels: ', num_metgrid_levels CALL wrf_message(message) DO j = jts, MIN(jte,jde-1) DO i = its, MIN(ite,ide-1) grid%fis(I,J)=grid%ht_gc(I,J)*g ! ! IF ( grid%p_gc(I,J,1) .ne. 200100. .AND. (grid%ht_gc(I,J) .eq. grid%ght_gc(I,J,1)) .AND. grid%ht_gc(I,J) .ne. 0) THEN IF ( grid%p_gc(I,J,1) .ne. 200100. .AND. (abs(grid%ht_gc(I,J)-grid%ght_gc(I,J,1)) .lt. 0.01) .AND. grid%ht_gc(I,J) .ne. 0) THEN IF (mod(I,10) .eq. 0 .and. mod(J,10) .eq. 0) THEN write(message,*) 'grid%ht_gc and grid%toposoil to swap, flag_soilhgt ::: ', & I,J, grid%ht_gc(I,J),grid%toposoil(I,J),flag_soilhgt CALL wrf_debug(10,message) ENDIF IF ( ( flag_soilhgt.EQ. 1 ) ) THEN grid%ght_gc(I,J,1)=grid%toposoil(I,J) ENDIF ENDIF ENDDO ENDDO numzero=0 numexamined=0 DO j = jts, MIN(jte,jde-1) DO i = its, MIN(ite,ide-1) numexamined=numexamined+1 if(grid%fis(i,j)<1e-5 .and. grid%fis(i,j)>-1e5 ) then numzero=numzero+1 end if enddo enddo write(message,*) 'TOTAL NEAR-ZERO FIS POINTS: ',numzero,' OF ',numexamined call wrf_debug(10,message) #if ( HWRF == 1 ) interp_notph: if(.not. grid%use_prep_hybrid) then #endif if (.NOT. allocated(PDVP)) allocate(PDVP(IMS:IME,JMS:JME)) if (.NOT. allocated(P3D_OUT)) allocate(P3D_OUT(IMS:IME,JMS:JME,KDS:KDE-1)) if (.NOT. allocated(PSFC_OUTV)) allocate(PSFC_OUTV(IMS:IME,JMS:JME)) if (.NOT. allocated(P3DV_OUT)) allocate(P3DV_OUT(IMS:IME,JMS:JME,KDS:KDE-1)) if (.NOT. allocated(P3DV_IN)) allocate(P3DV_IN(IMS:IME,JMS:JME,num_metgrid_levels)) CALL compute_nmm_surfacep (grid%ht_gc, grid%ght_gc, grid%p_gc , grid%t_gc & &, grid%psfc_out, num_metgrid_levels & &, IDS,IDE,JDS,JDE,KDS,KDE & &, IMS,IME,JMS,JME,KMS,KME & &, ITS,ITE,JTS,JTE,KTS,KTE ) ! H points if (internal_time_loop .eq. 1) then write(message,*) 'psfc points (final combined)' loopinc=max( (JTE-JTS)/20,1) iloopinc=max( (ITE-ITS)/10,1) CALL wrf_message(message) DO J=min(JTE,JDE-1),JTS,-loopinc write(message,633) (grid%psfc_out(I,J)/100.,I=its,min(ite,IDE-1),iloopinc) CALL wrf_message(message) ENDDO endif 633 format(35(f5.0,1x)) CALL compute_3d_pressure (grid%psfc_out,grid%aeta1,grid%aeta2 & &, grid%pdtop,grid%pt,grid%pd,p3d_out & &, IDS,IDE,JDS,JDE,KDS,KDE & &, IMS,IME,JMS,JME,KMS,KME & &, ITS,ITE,JTS,JTE,KTS,KTE ) #ifdef DM_PARALLEL ips=its ; ipe=ite ; jps=jts ; jpe=jte ; kps=kts ; kpe=kte # include "HALO_NMM_MG2.inc" #endif #ifdef DM_PARALLEL # include "HALO_NMM_MG3.inc" #endif do K=1,num_metgrid_levels do J=JTS,min(JTE,JDE-1) do I=ITS,min(ITE,IDE-1) IF (K .eq. KTS) THEN IF (J .eq. JDS .and. I .lt. IDE-1) THEN ! S boundary PDVP(I,J)=0.5*(grid%pd(I,J)+grid%pd(I+1,J)) PSFC_OUTV(I,J)=0.5*(grid%psfc_out(I,J)+grid%psfc_out(I+1,J)) ELSEIF (J .eq. JDE-1 .and. I .lt. IDE-1) THEN ! N boundary PDVP(I,J)=0.5*(grid%pd(I,J)+grid%pd(I+1,J)) PSFC_OUTV(I,J)=0.5*(grid%psfc_out(I,J)+grid%psfc_out(I+1,J)) ELSEIF (I .eq. IDS .and. mod(J,2) .eq. 0) THEN ! W boundary PDVP(I,J)=0.5*(grid%pd(I,J-1)+grid%pd(I,J+1)) PSFC_OUTV(I,J)=0.5*(grid%psfc_out(I,J-1)+grid%psfc_out(I,J+1)) ELSEIF (I .eq. IDE-1 .and. mod(J,2) .eq. 0) THEN ! E boundary PDVP(I,J)=0.5*(grid%pd(I,J-1)+grid%pd(I,J+1)) PSFC_OUTV(I,J)=0.5*(grid%psfc_out(I,J-1)+grid%psfc_out(I,J+1)) ELSEIF (I .eq. IDE-1 .and. mod(J,2) .eq. 1) THEN ! phantom E boundary PDVP(I,J)=grid%pd(I,J) PSFC_OUTV(I,J)=grid%psfc_out(I,J) ELSEIF (mod(J,2) .eq. 0) THEN ! interior even row PDVP(I,J)=0.25*(grid%pd(I,J)+grid%pd(I-1,J)+grid%pd(I,J+1)+grid%pd(I,J-1)) PSFC_OUTV(I,J)=0.25*(grid%psfc_out(I,J)+grid%psfc_out(I-1,J)+ & grid%psfc_out(I,J+1)+grid%psfc_out(I,J-1)) ELSE ! interior odd row PDVP(I,J)=0.25*(grid%pd(I,J)+grid%pd(I+1,J)+grid%pd(I,J+1)+grid%pd(I,J-1)) PSFC_OUTV(I,J)=0.25*(grid%psfc_out(I,J)+grid%psfc_out(I+1,J)+ & grid%psfc_out(I,J+1)+grid%psfc_out(I,J-1)) ENDIF ENDIF IF (J .eq. JDS .and. I .lt. IDE-1) THEN ! S boundary P3DV_IN(I,J,K)=0.5*(grid%p_gc(I,J,K)+grid%p_gc(I+1,J,K)) ELSEIF (J .eq. JDE-1 .and. I .lt. IDE-1) THEN ! N boundary P3DV_IN(I,J,K)=0.5*(grid%p_gc(I,J,K)+grid%p_gc(I+1,J,K)) ELSEIF (I .eq. IDS .and. mod(J,2) .eq. 0) THEN ! W boundary P3DV_IN(I,J,K)=0.5*(grid%p_gc(I,J-1,K)+grid%p_gc(I,J+1,K)) ELSEIF (I .eq. IDE-1 .and. mod(J,2) .eq. 0) THEN ! E boundary P3DV_IN(I,J,K)=0.5*(grid%p_gc(I,J-1,K)+grid%p_gc(I,J+1,K)) ELSEIF (I .eq. IDE-1 .and. mod(J,2) .eq. 1) THEN ! phantom E boundary P3DV_IN(I,J,K)=grid%p_gc(I,J,K) ELSEIF (mod(J,2) .eq. 0) THEN ! interior even row P3DV_IN(I,J,K)=0.25*(grid%p_gc(I,J,K)+grid%p_gc(I-1,J,K) + & grid%p_gc(I,J+1,K)+grid%p_gc(I,J-1,K)) ELSE ! interior odd row P3DV_IN(I,J,K)=0.25*(grid%p_gc(I,J,K)+grid%p_gc(I+1,J,K) + & grid%p_gc(I,J+1,K)+grid%p_gc(I,J-1,K)) ENDIF enddo enddo enddo CALL compute_3d_pressure (psfc_outv,grid%aeta1,grid%aeta2 & &, grid%pdtop,grid%pt,pdvp,p3dv_out & &, IDS,IDE,JDS,JDE,KDS,KDE & &, IMS,IME,JMS,JME,KMS,KME & &, ITS,ITE,JTS,JTE,KTS,KTE ) CALL interp_press2press_lin(grid%p_gc, p3d_out & &, grid%t_gc, grid%t,num_metgrid_levels & &, .TRUE.,.TRUE.,.TRUE. & ! extrap, ignore_lowest, t_field &, IDS,IDE,JDS,JDE,KDS,KDE & &, IMS,IME,JMS,JME,KMS,KME & &, ITS,ITE,JTS,JTE,KTS,KTE, internal_time_loop ) CALL interp_press2press_lin(p3dv_in, p3dv_out & &, grid%u_gc, grid%u,num_metgrid_levels & &, .FALSE.,.TRUE.,.FALSE. & &, IDS,IDE,JDS,JDE,KDS,KDE & &, IMS,IME,JMS,JME,KMS,KME & &, ITS,ITE,JTS,JTE,KTS,KTE, internal_time_loop ) CALL interp_press2press_lin(p3dv_in, p3dv_out & &, grid%v_gc, grid%v,num_metgrid_levels & &, .FALSE.,.TRUE.,.FALSE. & &, IDS,IDE,JDS,JDE,KDS,KDE & &, IMS,IME,JMS,JME,KMS,KME & &, ITS,ITE,JTS,JTE,KTS,KTE, internal_time_loop ) IF (hyb_coor) THEN CALL wind_adjust(p3dv_in,p3dv_out,grid%u_gc,grid%v_gc,grid%u,grid%v & &, num_metgrid_levels,5000. & &, IDS,IDE,JDS,JDE,KDS,KDE & &, IMS,IME,JMS,JME,KMS,KME & &, ITS,ITE,JTS,JTE,KTS,KTE ) ENDIF ALLOCATE(qtmp(IMS:IME,JMS:JME,num_metgrid_levels)) ALLOCATE(qtmp2(IMS:IME,JMS:JME,num_metgrid_levels)) CALL rh_to_mxrat (grid%rh_gc, grid%t_gc, grid%p_gc, qtmp , .TRUE. , & ids , ide , jds , jde , 1 , num_metgrid_levels , & ims , ime , jms , jme , 1 , num_metgrid_levels , & its , ite , jts , jte , 1 , num_metgrid_levels ) do K=1,num_metgrid_levels do J=JTS,min(JTE,JDE-1) do I=ITS,min(ITE,IDE-1) QTMP2(I,J,K)=QTMP(I,J,K)/(1.0+QTMP(I,J,K)) end do end do end do CALL interp_press2press_log(grid%p_gc, p3d_out & &, QTMP2, grid%q,num_metgrid_levels & &, .FALSE.,.TRUE. & &, IDS,IDE,JDS,JDE,KDS,KDE & &, IMS,IME,JMS,JME,KMS,KME & &, ITS,ITE,JTS,JTE,KTS,KTE, internal_time_loop ) IF (ALLOCATED(QTMP)) DEALLOCATE(QTMP) IF (ALLOCATED(QTMP)) DEALLOCATE(QTMP2) #if ( HWRF == 1 ) else ! we are using prep_hybrid ! Compute surface pressure: grid%psfc_out=grid%pdtop+grid%pd end if interp_notph #endif ! Get the monthly values interpolated to the current date ! for the traditional monthly ! fields of green-ness fraction and background grid%albedo. if (internal_time_loop .eq. 1 .or. config_flags%sst_update .eq. 1) then CALL monthly_interp_to_date ( grid%greenfrac_gc , current_date , grid%vegfra , & ids , ide , jds , jde , kds , kde , & ims , ime , jms , jme , kms , kme , & its , ite , jts , jte , kts , kte ) CALL monthly_interp_to_date ( grid%albedo12m_gc , current_date , grid%albbck , & ids , ide , jds , jde , kds , kde , & ims , ime , jms , jme , kms , kme , & its , ite , jts , jte , kts , kte ) ! Get the min/max of each i,j for the monthly green-ness fraction. CALL monthly_min_max ( grid%greenfrac_gc , grid%shdmin , grid%shdmax , & ids , ide , jds , jde , kds , kde , & ims , ime , jms , jme , kms , kme , & its , ite , jts , jte , kts , kte ) ! The model expects the green-ness values in percent, not fraction. DO j = jts, MIN(jte,jde-1) DO i = its, MIN(ite,ide-1) !! grid%vegfra(i,j) = grid%vegfra(i,j) * 100. grid%shdmax(i,j) = grid%shdmax(i,j) * 100. grid%shdmin(i,j) = grid%shdmin(i,j) * 100. grid%vegfrc(I,J)=grid%vegfra(I,J) END DO END DO ! The model expects the albedo fields as ! a fraction, not a percent. Set the water values to 8%. DO j = jts, MIN(jte,jde-1) DO i = its, MIN(ite,ide-1) if (grid%albbck(i,j) .lt. 5.) then write(message,*) 'reset albedo to 8%... I,J,albbck:: ', I,J,grid%albbck(I,J) CALL wrf_debug(10,message) grid%albbck(I,J)=8. endif grid%albbck(i,j) = grid%albbck(i,j) / 100. grid%snoalb(i,j) = grid%snoalb(i,j) / 100. IF ( grid%landmask(i,j) .LT. 0.5 ) THEN grid%albbck(i,j) = 0.08 grid%snoalb(i,j) = 0.08 END IF grid%albase(i,j)=grid%albbck(i,j) grid%mxsnal(i,j)=grid%snoalb(i,j) END DO END DO endif #if ( HWRF == 1 ) if(.not.grid%use_prep_hybrid) then #endif ! new deallocs DEALLOCATE(p3d_out,p3dv_out,p3dv_in) #if ( HWRF == 1 ) end if #endif END IF ! <----- END OF VERTICAL INTERPOLATION PART ----> !! compute SST at each time if updating SST if (config_flags%sst_update == 1) then DO j = jts, MIN(jde-1,jte) DO i = its, MIN(ide-1,ite) if (grid%SM(I,J) .lt. 0.5) then grid%SST(I,J)=0. endif if (grid%SM(I,J) .gt. 0.5) then grid%SST(I,J)=grid%NMM_TSK(I,J) grid%NMM_TSK(I,J)=0. endif IF ( (grid%NMM_TSK(I,J)+grid%SST(I,J)) .lt. 200. .or. & (grid%NMM_TSK(I,J)+grid%SST(I,J)) .gt. 350. ) THEN write(message,*) 'TSK, SST trouble at : ', I,J CALL wrf_message(message) write(message,*) 'SM, NMM_TSK,SST ', grid%SM(I,J),grid%NMM_TSK(I,J),grid%SST(I,J) CALL wrf_message(message) ENDIF ENDDO ENDDO endif ! sst_update test if (internal_time_loop .eq. 1) then !!! weasd has "snow water equivalent" in mm DO j = jts, MIN(jte,jde-1) DO i = its, MIN(ite,ide-1) IF(grid%sm(I,J).GT.0.9) THEN IF (grid%xice(I,J) .gt. 0) then grid%si(I,J)=1.0 ENDIF ! SEA grid%epsr(I,J)=.97 grid%embck(I,J)=.97 grid%gffc(I,J)=0. grid%albedo(I,J)=.06 grid%albase(I,J)=.06 IF(grid%si (I,J).GT.0. ) THEN ! SEA-ICE grid%sm(I,J)=0. grid%si(I,J)=0. grid%sice(I,J)=1. grid%gffc(I,J)=0. ! just leave zero as irrelevant grid%albedo(I,J)=.60 grid%albase(I,J)=.60 ENDIF ELSE grid%si(I,J)=5.0*grid%weasd(I,J)/1000. ! LAND grid%epsr(I,J)=1.0 grid%embck(I,J)=1.0 grid%gffc(I,J)=0.0 ! just leave zero as irrelevant grid%sice(I,J)=0. grid%sno(I,J)=grid%si(I,J)*.20 ENDIF ENDDO ENDDO ! DETERMINE grid%albedo OVER LAND DO j = jts, MIN(jte,jde-1) DO i = its, MIN(ite,ide-1) IF(grid%sm(I,J).LT.0.9.AND.grid%sice(I,J).LT.0.9) THEN ! SNOWFREE albedo IF ( (grid%sno(I,J) .EQ. 0.0) .OR. & (grid%albase(I,J) .GE. grid%mxsnal(I,J) ) ) THEN grid%albedo(I,J) = grid%albase(I,J) ELSE ! MODIFY albedo IF SNOWCOVER: ! BELOW SNOWDEPTH THRESHOLD... IF (grid%sno(I,J) .LT. SNUP) THEN RSNOW = grid%sno(I,J)/SNUP SNOFAC = 1. - ( EXP(-SALP*RSNOW) - RSNOW*EXP(-SALP)) ! ABOVE SNOWDEPTH THRESHOLD... ELSE SNOFAC = 1.0 ENDIF ! CALCULATE grid%albedo ACCOUNTING FOR SNOWDEPTH AND VGFRCK grid%albedo(I,J) = grid%albase(I,J) & + (1.0-grid%vegfra(I,J))*SNOFAC*(grid%mxsnal(I,J)-grid%albase(I,J)) ENDIF END IF grid%si(I,J)=5.0*grid%weasd(I,J) grid%sno(I,J)=grid%weasd(I,J) !! convert vegfra grid%vegfra(I,J)=grid%vegfra(I,J)*100. ! ENDDO ENDDO #ifdef DM_PARALLEL ALLOCATE(SM_G(IDS:IDE,JDS:JDE),SICE_G(IDS:IDE,JDS:JDE)) CALL WRF_PATCH_TO_GLOBAL_REAL( grid%sice(IMS,JMS) & &, SICE_G,grid%DOMDESC & &, 'z','xy' & &, IDS,IDE-1,JDS,JDE-1,1,1 & &, IMS,IME,JMS,JME,1,1 & &, ITS,ITE,JTS,JTE,1,1 ) CALL WRF_PATCH_TO_GLOBAL_REAL( grid%sm(IMS,JMS) & &, SM_G,grid%DOMDESC & &, 'z','xy' & &, IDS,IDE-1,JDS,JDE-1,1,1 & &, IMS,IME,JMS,JME,1,1 & &, ITS,ITE,JTS,JTE,1,1 ) IF (WRF_DM_ON_MONITOR()) THEN 637 format(40(f3.0,1x)) allocate(IHE_G(JDS:JDE-1),IHW_G(JDS:JDE-1)) DO j = JDS, JDE-1 IHE_G(J)=MOD(J+1,2) IHW_G(J)=IHE_G(J)-1 ENDDO DO ITER=1,10 DO j = jds+1, (jde-1)-1 DO i = ids+1, (ide-1)-1 ! any sea ice around point in question? IF (SM_G(I,J) .ge. 0.9) THEN SEAICESUM=SICE_G(I+IHE_G(J),J+1)+SICE_G(I+IHW_G(J),J+1)+ & SICE_G(I+IHE_G(J),J-1)+SICE_G(I+IHW_G(J),J-1) IF (SEAICESUM .ge. 1. .and. SEAICESUM .lt. 3.) THEN IF ((SICE_G(I+IHE_G(J),J+1).lt.0.1 .and. SM_G(I+IHE_G(J),J+1).lt.0.1) .OR. & (SICE_G(I+IHW_G(J),J+1).lt.0.1 .and. SM_G(I+IHW_G(J),J+1).lt.0.1) .OR. & (SICE_G(I+IHE_G(J),J-1).lt.0.1 .and. SM_G(I+IHE_G(J),J-1).lt.0.1) .OR. & (SICE_G(I+IHW_G(J),J-1).lt.0.1 .and. SM_G(I+IHW_G(J),J-1).lt.0.1)) THEN ! HAVE SEA ICE AND A SURROUNDING LAND POINT - CONVERT TO SEA ICE write(message,*) 'making seaice (1): ', I,J CALL wrf_debug(100,message) SICE_G(I,J)=1.0 SM_G(I,J)=0. ENDIF ELSEIF (SEAICESUM .ge. 3) THEN ! WATER POINT SURROUNDED BY ICE - CONVERT TO SEA ICE write(message,*) 'making seaice (2): ', I,J CALL wrf_debug(100,message) SICE_G(I,J)=1.0 SM_G(I,J)=0. ENDIF ENDIF ENDDO ENDDO ENDDO ENDIF CALL WRF_GLOBAL_TO_PATCH_REAL( SICE_G, grid%sice & &, grid%DOMDESC & &, 'z','xy' & &, IDS,IDE-1,JDS,JDE-1,1,1 & &, IMS,IME,JMS,JME,1,1 & &, ITS,ITE,JTS,JTE,1,1 ) CALL WRF_GLOBAL_TO_PATCH_REAL( SM_G,grid%sm & &, grid%DOMDESC & &, 'z','xy' & &, IDS,IDE-1,JDS,JDE-1,1,1 & &, IMS,IME,JMS,JME,1,1 & &, ITS,ITE,JTS,JTE,1,1 ) IF (WRF_DM_ON_MONITOR()) THEN #if ( HWRF == 1 ) ! SM_G is still needed for the high-res grid #else DEALLOCATE(SM_G) #endif deallocate(SICE_G) DEALLOCATE(IHE_G,IHW_G) ENDIF ! write(message,*) 'revised sea ice on patch' ! CALL wrf_debug(100,message) ! DO J=JTE,JTS,-(((JTE-JTS)/25)+1) ! write(message,637) (grid%sice(I,J),I=ITS,ITE,ITE/20) ! CALL wrf_debug(100,message) ! END DO #else ! serial sea ice reprocessing allocate(IHE(JDS:JDE-1),IHW(JDS:JDE-1)) DO j = jts, MIN(jte,jde-1) IHE(J)=MOD(J+1,2) IHW(J)=IHE(J)-1 ENDDO DO ITER=1,10 DO j = jts+1, MIN(jte,jde-1)-1 DO i = its+1, MIN(ite,ide-1)-1 ! any sea ice around point in question? IF (grid%sm(I,J) .gt. 0.9) THEN SEAICESUM=grid%sice(I+IHE(J),J+1)+grid%sice(I+IHW(J),J+1)+ & grid%sice(I+IHE(J),J-1)+grid%sice(I+IHW(J),J-1) IF (SEAICESUM .ge. 1. .and. SEAICESUM .lt. 3.) THEN IF ((grid%sice(I+IHE(J),J+1).lt.0.1 .and. grid%sm(I+IHE(J),J+1).lt.0.1) .OR. & (grid%sice(I+IHW(J),J+1).lt.0.1 .and. grid%sm(I+IHW(J),J+1).lt.0.1) .OR. & (grid%sice(I+IHE(J),J-1).lt.0.1 .and. grid%sm(I+IHE(J),J-1).lt.0.1) .OR. & (grid%sice(I+IHW(J),J-1).lt.0.1 .and. grid%sm(I+IHW(J),J-1).lt.0.1)) THEN ! HAVE SEA ICE AND A SURROUNDING LAND POINT - CONVERT TO SEA ICE grid%sice(I,J)=1.0 grid%sm(I,J)=0. ENDIF ELSEIF (SEAICESUM .ge. 3) THEN ! WATER POINT SURROUNDED BY ICE - CONVERT TO SEA ICE grid%sice(I,J)=1.0 grid%sm(I,J)=0. ENDIF ENDIF ENDDO ENDDO ENDDO DEALLOCATE(IHE,IHW) #endif ! this block meant to guarantee land/sea agreement between sm and landmask DO j = jts, MIN(jte,jde-1) DO i = its, MIN(ite,ide-1) IF (grid%sm(I,J) .gt. 0.5) THEN grid%landmask(I,J)=0.0 ELSEIF (grid%sm(I,J) .lt. 0.5 .and. grid%sice(I,J) .gt. 0.9) then grid%landmask(I,J)=0.0 ELSEIF (grid%sm(I,J) .lt. 0.5 .and. grid%sice(I,J) .lt. 0.1) then grid%landmask(I,J)=1.0 ELSE write(message,*) 'missed point in grid%landmask definition ' , I,J CALL wrf_message(message) grid%landmask(I,J)=0.0 ENDIF ! IF (grid%sice(I,J) .gt. 0.5 .and. grid%nmm_tsk(I,J) .lt. 0.1 .and. grid%sst(I,J) .gt. 0.) THEN write(message,*) 'set grid%nmm_tsk to: ', grid%sst(I,J) CALL wrf_message(message) grid%nmm_tsk(I,J)=grid%sst(I,J) grid%sst(I,J)=0. endif ENDDO ENDDO ! For sf_surface_physics = 1, we want to use close to a 10 cm value ! for the bottom level of the soil temps. IF ( ( model_config_rec%sf_surface_physics(grid%id) .EQ. 1 ) .AND. & ( flag_st000010 .EQ. 1 ) ) THEN DO j = jts , MIN(jde-1,jte) DO i = its , MIN(ide-1,ite) grid%soiltb(i,j) = grid%st000010(i,j) END DO END DO END IF ! Adjust the various soil temperature values depending on the difference in ! in elevation between the current model's elevation and the incoming data's ! orography. IF ( ( flag_toposoil .EQ. 1 ) ) THEN ALLOCATE(HT(ims:ime,jms:jme)) DO J=jms,jme DO I=ims,ime HT(I,J)=grid%fis(I,J)/9.81 END DO END DO ! if (maxval(grid%toposoil) .gt. 100.) then ! ! Being avoided. Something to revisit eventually. ! !1219 might be simply a matter of including toposoil ! ! CODE NOT TESTED AT NCEP USING THIS FUNCTIONALITY, ! SO TO BE SAFE WILL AVOID FOR RETRO RUNS. ! ! CALL adjust_soil_temp_new ( grid%soiltb , 2 , & ! grid%nmm_tsk , ht , grid%toposoil , grid%landmask, flag_toposoil , & ! grid%st000010 , st010040 , st040100 , st100200 , st010200 , & ! flag_st000010 , flag_st010040 , flag_st040100 , & ! flag_st100200 , flag_st010200 , & ! soilt000 , soilt005 , soilt020 , soilt040 , & ! soilt160 , soilt300 , & ! flag_soilt000 , flag_soilt005 , flag_soilt020 , & ! flag_soilt040 , flag_soilt160 , flag_soilt300 , & ! ids , ide , jds , jde , kds , kde , & ! ims , ime , jms , jme , kms , kme , & ! its , ite , jts , jte , kts , kte ) ! endif END IF ! Process the LSM data. ! surface_input_source=1 => use data from static file ! (fractional category as input) ! surface_input_source=2 => use data from grib file ! (dominant category as input) IF ( config_flags%surface_input_source .EQ. 1 ) THEN grid%vegcat (its,jts) = 0 grid%soilcat(its,jts) = 0 END IF ! Generate the vegetation and soil category information ! from the fractional input ! data, or use the existing dominant category fields if they exist. IF ((grid%soilcat(its,jts) .LT. 0.5) .AND. (grid%vegcat(its,jts) .LT. 0.5)) THEN num_veg_cat = SIZE ( grid%landusef_gc , DIM=3 ) num_soil_top_cat = SIZE ( grid%soilctop_gc , DIM=3 ) num_soil_bot_cat = SIZE ( grid%soilcbot_gc , DIM=3 ) do J=JMS,JME do K=1,num_veg_cat do I=IMS,IME grid%landusef(I,K,J)=grid%landusef_gc(I,J,K) enddo enddo enddo do J=JMS,JME do K=1,num_soil_top_cat do I=IMS,IME grid%soilctop(I,K,J)=grid%soilctop_gc(I,J,K) enddo enddo enddo do J=JMS,JME do K=1,num_soil_bot_cat do I=IMS,IME grid%soilcbot(I,K,J)=grid%soilcbot_gc(I,J,K) enddo enddo enddo ! grid%sm (1=water, 0=land) ! grid%landmask(0=water, 1=land) write(message,*) 'landmask into process_percent_cat_new' CALL wrf_debug(1,message) do J=JTE,JTS,-(((JTE-JTS)/20)+1) write(message,641) (grid%landmask(I,J),I=ITS,min(ITE,IDE-1),((ITE-ITS)/15)+1) CALL wrf_debug(1,message) enddo 641 format(25(f3.0,1x)) CALL process_percent_cat_new ( grid%landmask , & grid%landusef , grid%soilctop , grid%soilcbot , & grid%isltyp , grid%ivgtyp , & num_veg_cat , num_soil_top_cat , num_soil_bot_cat , & ids , ide , jds , jde , kds , kde , & ims , ime , jms , jme , kms , kme , & its , ite , jts , jte , kts , kte , & model_config_rec%iswater(grid%id) ) DO j = jts , MIN(jde-1,jte) DO i = its , MIN(ide-1,ite) grid%vegcat(i,j) = grid%ivgtyp(i,j) grid%soilcat(i,j) = grid%isltyp(i,j) END DO END DO ELSE ! Do we have dominant soil and veg data from the input already? IF ( grid%soilcat(its,jts) .GT. 0.5 ) THEN DO j = jts, MIN(jde-1,jte) DO i = its, MIN(ide-1,ite) grid%isltyp(i,j) = NINT( grid%soilcat(i,j) ) END DO END DO END IF IF ( grid%vegcat(its,jts) .GT. 0.5 ) THEN DO j = jts, MIN(jde-1,jte) DO i = its, MIN(ide-1,ite) grid%ivgtyp(i,j) = NINT( grid%vegcat(i,j) ) END DO END DO END IF ENDIF DO j = jts, MIN(jde-1,jte) DO i = its, MIN(ide-1,ite) IF (grid%sice(I,J) .lt. 0.1) THEN IF (grid%landmask(I,J) .gt. 0.5 .and. grid%sm(I,J) .gt. 0.5) THEN write(message,*) 'land mask and grid%sm both > 0.5: ', & I,J,grid%landmask(I,J),grid%sm(I,J) CALL wrf_message(message) grid%sm(I,J)=0. ELSEIF (grid%landmask(I,J) .lt. 0.5 .and. grid%sm(I,J) .lt. 0.5) THEN write(message,*) 'land mask and grid%sm both < 0.5: ', & I,J, grid%landmask(I,J),grid%sm(I,J) CALL wrf_message(message) grid%sm(I,J)=1. ENDIF ELSE IF (grid%landmask(I,J) .gt. 0.5 .and. grid%sm(I,J)+grid%sice(I,J) .gt. 0.9) then write(message,*) 'landmask says LAND, sm/sice say SEAICE: ', I,J ENDIF ENDIF ENDDO ENDDO DO j = jts, MIN(jde-1,jte) DO i = its, MIN(ide-1,ite) if (grid%sice(I,J) .gt. 0.9) then grid%isltyp(I,J)=16 grid%ivgtyp(I,J)=24 endif ENDDO ENDDO DO j = jts, MIN(jde-1,jte) DO i = its, MIN(ide-1,ite) if (grid%sm(I,J) .lt. 0.5) then grid%sst(I,J)=0. endif if (grid%sm(I,J) .gt. 0.5) then if (grid%sst(I,J) .lt. 0.1) then grid%sst(I,J)=grid%nmm_tsk(I,J) endif grid%nmm_tsk(I,J)=0. endif IF ( (grid%nmm_tsk(I,J)+grid%sst(I,J)) .lt. 200. .or. & (grid%nmm_tsk(I,J)+grid%sst(I,J)) .gt. 350. ) THEN write(message,*) 'TSK, sst trouble at : ', I,J CALL wrf_message(message) write(message,*) 'sm, nmm_tsk,sst ', grid%sm(I,J),grid%nmm_tsk(I,J),grid%sst(I,J) CALL wrf_message(message) ENDIF ENDDO ENDDO write(message,*) 'grid%sm' CALL wrf_message(message) DO J=min(jde-1,jte),jts,-((jte-jts)/15+1) write(message,635) (grid%sm(i,J),I=its,ite,((ite-its)/10)+1) CALL wrf_message(message) END DO write(message,*) 'sst/nmm_tsk' CALL wrf_debug(10,message) DO J=min(jde-1,jte),jts,-((jte-jts)/15+1) write(message,635) (grid%sst(I,J)+grid%nmm_tsk(I,J),I=ITS,min(ide-1,ite),((ite-its)/10)+1) CALL wrf_debug(10,message) END DO 635 format(20(f5.1,1x)) DO j = jts, MIN(jde-1,jte) DO i = its, MIN(ide-1,ite) IF ( ( grid%landmask(i,j) .LT. 0.5 ) .AND. ( flag_sst .EQ. 1 ) ) THEN grid%soiltb(i,j) = grid%sst(i,j) ELSE IF ( grid%landmask(i,j) .GT. 0.5 ) THEN grid%soiltb(i,j) = grid%nmm_tsk(i,j) END IF END DO END DO ! END IF ! Land use categories, dominant soil and vegetation types (if available). ! allocate(grid%lu_index(ims:ime,jms:jme)) DO j = jts, MIN(jde-1,jte) DO i = its, MIN(ide-1,ite) grid%lu_index(i,j) = grid%ivgtyp(i,j) END DO END DO if (flag_sst .eq. 1) log_flag_sst=.true. if (flag_sst .eq. 0) log_flag_sst=.false. write(message,*) 'st_input dimensions: ', size(st_input,dim=1), & size(st_input,dim=2),size(st_input,dim=3) CALL wrf_debug(100,message) ! write(message,*) 'maxval st_input(1): ', maxval(st_input(:,1,:)) ! CALL wrf_message(message) ! write(message,*) 'maxval st_input(2): ', maxval(st_input(:,2,:)) ! CALL wrf_message(message) ! write(message,*) 'maxval st_input(3): ', maxval(st_input(:,3,:)) ! CALL wrf_message(message) ! write(message,*) 'maxval st_input(4): ', maxval(st_input(:,4,:)) ! CALL wrf_message(message) ! ============================================================= IF (.NOT. ALLOCATED(TG_ALT))ALLOCATE(TG_ALT(grid%sm31:grid%em31,jms:jme)) TPH0=TPH0D*DTR WBD=-(((ide-1)-1)*grid%dlmd) WB= WBD*DTR SBD=-(((jde-1)/2)*grid%dphd) SB= SBD*DTR DLM=grid%dlmd*DTR DPH=grid%dphd*DTR TDLM=DLM+DLM TDPH=DPH+DPH WBI=WB+TDLM SBI=SB+TDPH EBI=WB+(ide-2)*TDLM ANBI=SB+(jde-2)*DPH STPH0=SIN(TPH0) CTPH0=COS(TPH0) TSPH=3600./GRID%DT DO J=JTS,min(JTE,JDE-1) TLM=WB-TDLM+MOD(J,2)*DLM !For velocity points on the E grid TPH=SB+float(J-1)*DPH STPH=SIN(TPH) CTPH=COS(TPH) DO I=ITS,MIN(ITE,IDE-1) if (I .eq. ITS) THEN TLM=TLM+TDLM*ITS else TLM=TLM+TDLM endif TERM1=(STPH0*CTPH*COS(TLM)+CTPH0*STPH) FP=TWOM*(TERM1) ! jbao orig grid%f(I,J)=0.5*GRID%DT*FP ! jbao Coriolis correction for idealized! grid%f(I,J)=3.15656e-5*0.5*GRID%DT ! jbao Coriolis correction for idealized! ENDDO ENDDO DO J=JTS,min(JTE,JDE-1) TLM=WB-TDLM+MOD(J+1,2)*DLM !For mass points on the E grid TPH=SB+float(J-1)*DPH STPH=SIN(TPH) CTPH=COS(TPH) DO I=ITS,MIN(ITE,IDE-1) if (I .eq. ITS) THEN TLM=TLM+TDLM*ITS else TLM=TLM+TDLM endif TERM1=(STPH0*CTPH*COS(TLM)+CTPH0*STPH) TERM1=MIN(TERM1,1.0D0) TERM1=MAX(TERM1,-1.0D0) APH=ASIN(TERM1) TG_ALT(I,J)=TG0+TGA*COS(APH)-grid%fis(I,J)/3333. ENDDO ENDDO DO j = jts, MIN(jde-1,jte) DO i = its, MIN(ide-1,ite) ! IF ( ( grid%landmask(i,j) .LT. 0.5 ) .AND. ( flag_sst .EQ. 1 ) .AND. & ! grid%sice(I,J) .eq. 0. ) THEN ! grid%tg(i,j) = grid%sst(i,j) ! ELSEIF (grid%sice(I,J) .eq. 1) THEN ! grid%tg(i,j) = 271.16 ! END IF if (grid%tg(I,J) .lt. 200.) then ! only use default TG_ALT definition if ! not getting TGROUND from grid%si grid%tg(I,J)=TG_ALT(I,J) endif if (grid%tg(I,J) .lt. 200. .or. grid%tg(I,J) .gt. 320.) then write(message,*) 'problematic grid%tg point at : ', I,J CALL wrf_message( message ) endif adum2d(i,j)=grid%nmm_tsk(I,J)+grid%sst(I,J) END DO END DO DEALLOCATE(TG_ALT) write(message,*) 'call process_soil_real with num_st_levels_input: ', num_st_levels_input CALL wrf_message( message ) ! ============================================================= CALL process_soil_real ( adum2d, grid%tg , & grid%landmask, grid%sst, & st_input, sm_input, sw_input, & st_levels_input , sm_levels_input , & sw_levels_input , & grid%sldpth , grid%dzsoil , grid%stc , grid%smc , grid%sh2o, & flag_sst , flag_soilt000, flag_soilm000, & ids , ide , jds , jde , kds , kde , & ims , ime , jms , jme , kms , kme , & its , ite , jts , jte , kts , kte , & model_config_rec%sf_surface_physics(grid%id) , & model_config_rec%num_soil_layers , & model_config_rec%real_data_init_type , & num_st_levels_input , num_sm_levels_input , & num_sw_levels_input , & num_st_levels_alloc , num_sm_levels_alloc , & num_sw_levels_alloc ) ! ============================================================= ! Minimum soil values, residual, from RUC LSM scheme. ! For input from Noah and using ! RUC LSM scheme, this must be subtracted from the input ! total soil moisture. For input RUC data and using the Noah LSM scheme, ! this value must be added to the soil moisture_input. lqmi(1:num_soil_top_cat) = & (/0.045, 0.057, 0.065, 0.067, 0.034, 0.078, 0.10, & 0.089, 0.095, 0.10, 0.070, 0.068, 0.078, 0.0, & 0.004, 0.065 /) !dusan , 0.020, 0.004, 0.008 /) ! At the initial time we care about values of soil moisture and temperature, ! other times are ignored by the model, so we ignore them, too. account_for_zero_soil_moisture : SELECT CASE ( model_config_rec%sf_surface_physics(grid%id) ) CASE ( LSMSCHEME ) iicount = 0 IF ( FLAG_SM000010 .EQ. 1 ) THEN DO j = jts, MIN(jde-1,jte) DO i = its, MIN(ide-1,ite) IF ((grid%landmask(i,j).gt.0.5) .and. (grid%stc(i,1,j) .gt. 200) .and. & (grid%stc(i,1,j) .lt. 400) .and. (grid%smc(i,1,j) .lt. 0.005)) then write(message,*) 'Noah > Noah: bad soil moisture at i,j = ',i,j,grid%smc(i,:,j) CALL wrf_message(message) iicount = iicount + 1 grid%smc(i,:,j) = 0.005 END IF END DO END DO IF ( iicount .GT. 0 ) THEN write(message,*) 'Noah -> Noah: total number of small soil moisture locations= ',& iicount CALL wrf_message(message) END IF ELSE IF ( FLAG_SOILM000 .EQ. 1 ) THEN DO j = jts, MIN(jde-1,jte) DO i = its, MIN(ide-1,ite) grid%smc(i,:,j) = grid%smc(i,:,j) + lqmi(grid%isltyp(i,j)) END DO END DO DO j = jts, MIN(jde-1,jte) DO i = its, MIN(ide-1,ite) IF ((grid%landmask(i,j).gt.0.5) .and. (grid%stc(i,1,j) .gt. 200) .and. & (grid%stc(i,1,j) .lt. 400) .and. (grid%smc(i,1,j) .lt. 0.004)) then write(message,*) 'RUC -> Noah: bad soil moisture at i,j = ' & ,i,j,grid%smc(i,:,j) CALL wrf_message(message) iicount = iicount + 1 grid%smc(i,:,j) = 0.004 END IF END DO END DO IF ( iicount .GT. 0 ) THEN write(message,*) 'RUC -> Noah: total number of small soil moisture locations = ',& iicount CALL wrf_message(message) END IF END IF CASE ( RUCLSMSCHEME ) iicount = 0 IF ( FLAG_SM000010 .EQ. 1 ) THEN DO j = jts, MIN(jde-1,jte) DO i = its, MIN(ide-1,ite) grid%smc(i,:,j) = MAX ( grid%smc(i,:,j) - lqmi(grid%isltyp(i,j)) , 0. ) END DO END DO ELSE IF ( FLAG_SOILM000 .EQ. 1 ) THEN ! no op END IF END SELECT account_for_zero_soil_moisture !!! zero out grid%nmm_tsk at water points again DO j = jts, MIN(jde-1,jte) DO i = its, MIN(ide-1,ite) if (grid%sm(I,J) .gt. 0.5) then grid%nmm_tsk(I,J)=0. endif END DO END DO !! check on grid%stc DO j = jts, MIN(jde-1,jte) DO i = its, MIN(ide-1,ite) IF (grid%sice(I,J) .gt. 0.9) then DO L = 1, grid%num_soil_layers grid%stc(I,L,J)=271.16 ! grid%tg value used by Eta/NMM END DO END IF IF (grid%sm(I,J) .gt. 0.9) then DO L = 1, grid%num_soil_layers grid%stc(I,L,J)=273.16 ! grid%tg value used by Eta/NMM END DO END IF END DO END DO DO j = jts, MIN(jde-1,jte) DO i = its, MIN(ide-1,ite) if (grid%sm(I,J) .lt. 0.1 .and. grid%stc(I,1,J) .lt. 0.1) THEN write(message,*) 'troublesome grid%sm,grid%stc,grid%smc value: ', I,J,grid%sm(I,J), grid%stc(I,1,J),grid%smc(I,1,J) CALL wrf_message(message) do JJ=J-1,J+1 do L=1, grid%num_soil_layers do II=I-1,I+1 if (II .ge. its .and. II .le. MIN(ide-1,ite) .and. & JJ .ge. jts .and. JJ .le. MIN(jde-1,jte)) then grid%stc(I,L,J)=amax1(grid%stc(I,L,J),grid%stc(II,L,JJ)) cur_smc=grid%smc(I,L,J) if ( grid%smc(II,L,JJ) .gt. 0.005 .and. grid%smc(II,L,JJ) .lt. 1.0) then aposs_smc=grid%smc(II,L,JJ) if ( cur_smc .eq. 0 ) then cur_smc=aposs_smc grid%smc(I,L,J)=cur_smc else cur_smc=amin1(cur_smc,aposs_smc) cur_smc=amin1(cur_smc,aposs_smc) grid%smc(I,L,J)=cur_smc endif endif endif ! bounds check enddo enddo enddo write(message,*) 'grid%stc, grid%smc(1) now: ', grid%stc(I,1,J),grid%smc(I,1,J) CALL wrf_message(message) endif if (grid%stc(I,1,J) .lt. 0.1) then write(message,*) 'QUITTING DUE TO STILL troublesome grid%stc value: ', I,J, grid%stc(I,1,J),grid%smc(I,1,J) call wrf_error_fatal(message) endif ENDDO ENDDO !hardwire soil stuff for time being ! RTDPTH=0. ! RTDPTH(1)=0.1 ! RTDPTH(2)=0.3 ! RTDPTH(3)=0.6 ! grid%sldpth=0. ! grid%sldpth(1)=0.1 ! grid%sldpth(2)=0.3 ! grid%sldpth(3)=0.6 ! grid%sldpth(4)=1.0 !!! main body of nmm_specific starts here ! do J=jts,min(jte,jde-1) do I=its,min(ite,ide-1) grid%res(I,J)=1. enddo enddo !! grid%hbm2 grid%hbm2=0. do J=jts,min(jte,jde-1) do I=its,min(ite,ide-1) IF ( (J .ge. 3 .and. J .le. (jde-1)-2) .AND. & (I .ge. 2 .and. I .le. (ide-1)-2+mod(J,2)) ) THEN grid%hbm2(I,J)=1. ENDIF enddo enddo !! grid%hbm3 grid%hbm3=0. !! LOOP OVER LOCAL DIMENSIONS do J=jts,min(jte,jde-1) grid%ihwg(J)=mod(J+1,2)-1 IF (J .ge. 4 .and. J .le. (jde-1)-3) THEN IHL=(ids+1)-grid%ihwg(J) IHH=(ide-1)-2 do I=its,min(ite,ide-1) IF (I .ge. IHL .and. I .le. IHH) grid%hbm3(I,J)=1. enddo ENDIF enddo !! grid%vbm2 grid%vbm2=0. do J=jts,min(jte,jde-1) do I=its,min(ite,ide-1) IF ( (J .ge. 3 .and. J .le. (jde-1)-2) .AND. & (I .ge. 2 .and. I .le. (ide-1)-1-mod(J,2)) ) THEN grid%vbm2(I,J)=1. ENDIF enddo enddo !! grid%vbm3 grid%vbm3=0. do J=jts,min(jte,jde-1) do I=its,min(ite,ide-1) IF ( (J .ge. 4 .and. J .le. (jde-1)-3) .AND. & (I .ge. 3-mod(J,2) .and. I .le. (ide-1)-2) ) THEN grid%vbm3(I,J)=1. ENDIF enddo enddo COAC=model_config_rec%coac(grid%id) CODAMP=model_config_rec%codamp(grid%id) DTAD=1.0 ! IDTCF=DTCF, IDTCF=4 DTCF=4.0 ! used? grid%dy_nmm=ERAD*DPH grid%cpgfv=-GRID%DT/(48.*grid%dy_nmm) grid%en= GRID%DT/( 4.*grid%dy_nmm)*DTAD grid%ent=GRID%DT/(16.*grid%dy_nmm)*DTAD DO J=jts,nnyp KHL2(J)=(IDE-1)*(J-1)-(J-1)/2+2 KVL2(J)=(IDE-1)*(J-1)-J/2+2 KHH2(J)=(IDE-1)*J-J/2-1 KVH2(J)=(IDE-1)*J-(J+1)/2-1 ENDDO TPH=SB-DPH DO J=jts,min(jte,jde-1) TPH=SB+float(J-1)*DPH DXP=ERAD*DLM*COS(TPH) DXJ(J)=DXP WPDARJ(J)=-W_NMM * & ((ERAD*DLM*AMIN1(COS(ANBI),COS(SBI)))**2+grid%dy_nmm**2)/ & (GRID%DT*32.*DXP*grid%dy_nmm) CPGFUJ(J)=-GRID%DT/(48.*DXP) CURVJ(J)=.5*GRID%DT*TAN(TPH)/ERAD FCPJ(J)=GRID%DT/(CP*192.*DXP*grid%dy_nmm) FDIVJ(J)=1./(12.*DXP*grid%dy_nmm) ! EMJ(J)= GRID%DT/( 4.*DXP)*DTAD ! EMTJ(J)=GRID%DT/(16.*DXP)*DTAD FADJ(J)=-GRID%DT/(48.*DXP*grid%dy_nmm)*DTAD ACDT=GRID%DT*SQRT((ERAD*DLM*AMIN1(COS(ANBI),COS(SBI)))**2+grid%dy_nmm**2) CDDAMP=CODAMP*ACDT HDACJ(J)=COAC*ACDT/(4.*DXP*grid%dy_nmm) DDMPUJ(J)=CDDAMP/DXP DDMPVJ(J)=CDDAMP/grid%dy_nmm ENDDO DO J=JTS,min(JTE,JDE-1) TLM=WB-TDLM+MOD(J,2)*DLM TPH=SB+float(J-1)*DPH STPH=SIN(TPH) CTPH=COS(TPH) DO I=ITS,MIN(ITE,IDE-1) if (I .eq. ITS) THEN TLM=TLM+TDLM*ITS else TLM=TLM+TDLM endif FP=TWOM*(CTPH0*STPH+STPH0*CTPH*COS(TLM)) !jbao orig grid%f(I,J)=0.5*GRID%DT*FP ! jbao Coriolis correction for idealized! grid%f(I,J)=3.15656e-5*0.5*GRID%DT print*,'real 2griddt is ',grid%dt ! jbao Coriolis correction for idealized! ENDDO ENDDO ! --------------DERIVED VERTICAL GRID CONSTANTS-------------------------- grid%ef4t=.5*GRID%DT/CP grid%f4q = -GRID%DT*DTAD grid%f4d =-.5*GRID%DT*DTAD DO L=KDS,KDE-1 grid%rdeta(L)=1./grid%deta(L) grid%f4q2(L)=-.25*GRID%DT*DTAD/grid%deta(L) ENDDO DO J=JTS,min(JTE,JDE-1) DO I=ITS,min(ITE,IDE-1) grid%dx_nmm(I,J)=DXJ(J) grid%wpdar(I,J)=WPDARJ(J)*grid%hbm2(I,J) grid%cpgfu(I,J)=CPGFUJ(J)*grid%vbm2(I,J) grid%curv(I,J)=CURVJ(J)*grid%vbm2(I,J) grid%fcp(I,J)=FCPJ(J)*grid%hbm2(I,J) grid%fdiv(I,J)=FDIVJ(J)*grid%hbm2(I,J) grid%fad(I,J)=FADJ(J) grid%hdacv(I,J)=HDACJ(J)*grid%vbm2(I,J) grid%hdac(I,J)=HDACJ(J)*1.25*grid%hbm2(I,J) ENDDO ENDDO DO J=JTS, MIN(JDE-1,JTE) IF (J.LE.5.OR.J.GE.(JDE-1)-4) THEN KHH=(IDE-1)-2+MOD(J,2) ! KHH is global...loop over I that have DO I=ITS,MIN(IDE-1,ITE) IF (I .ge. 2 .and. I .le. KHH) THEN grid%hdac(I,J)=grid%hdac(I,J)* DFC ENDIF ENDDO ELSE KHH=2+MOD(J,2) DO I=ITS,MIN(IDE-1,ITE) IF (I .ge. 2 .and. I .le. KHH) THEN grid%hdac(I,J)=grid%hdac(I,J)* DFC ENDIF ENDDO KHH=(IDE-1)-2+MOD(J,2) DO I=ITS,MIN(IDE-1,ITE) IF (I .ge. (IDE-1)-2 .and. I .le. KHH) THEN grid%hdac(I,J)=grid%hdac(I,J)* DFC ENDIF ENDDO ENDIF ENDDO DO J=JTS,min(JTE,JDE-1) DO I=ITS,min(ITE,IDE-1) grid%ddmpu(I,J)=DDMPUJ(J)*grid%vbm2(I,J) grid%ddmpv(I,J)=DDMPVJ(J)*grid%vbm2(I,J) grid%hdacv(I,J)=grid%hdacv(I,J)*grid%vbm2(I,J) ENDDO ENDDO ! --------------INCREASING DIFFUSION ALONG THE BOUNDARIES---------------- DO J=JTS,MIN(JDE-1,JTE) IF (J.LE.5.OR.J.GE.JDE-1-4) THEN KVH=(IDE-1)-1-MOD(J,2) DO I=ITS,min(IDE-1,ITE) IF (I .ge. 2 .and. I .le. KVH) THEN grid%ddmpu(I,J)=grid%ddmpu(I,J)*DDFC grid%ddmpv(I,J)=grid%ddmpv(I,J)*DDFC grid%hdacv(I,J)=grid%hdacv(I,J)* DFC ENDIF ENDDO ELSE KVH=3-MOD(J,2) DO I=ITS,min(IDE-1,ITE) IF (I .ge. 2 .and. I .le. KVH) THEN grid%ddmpu(I,J)=grid%ddmpu(I,J)*DDFC grid%ddmpv(I,J)=grid%ddmpv(I,J)*DDFC grid%hdacv(I,J)=grid%hdacv(I,J)* DFC ENDIF ENDDO KVH=(IDE-1)-1-MOD(J,2) DO I=ITS,min(IDE-1,ITE) IF (I .ge. IDE-1-2 .and. I .le. KVH) THEN grid%ddmpu(I,J)=grid%ddmpu(I,J)*DDFC grid%ddmpv(I,J)=grid%ddmpv(I,J)*DDFC grid%hdacv(I,J)=grid%hdacv(I,J)* DFC ENDIF ENDDO ENDIF ENDDO write(message,*) 'grid%stc(1)' CALL wrf_message(message) DO J=min(jde-1,jte),jts,-((jte-jts)/15+1) write(message,635) (grid%stc(I,1,J),I=its,min(ite,ide-1),(ite-its)/12+1) CALL wrf_message(message) ENDDO write(message,*) 'grid%smc(1)' CALL wrf_message(message) DO J=min(jde-1,jte),jts,-((jte-jts)/15+1) write(message,635) (grid%smc(I,1,J),I=its,min(ite,ide-1),(ite-its)/12+1) CALL wrf_message(message) ENDDO DO j = jts, MIN(jde-1,jte) DO i= ITS, MIN(IDE-1,ITE) if (grid%sm(I,J) .lt. 0.1 .and. grid%smc(I,1,J) .gt. 0.5 .and. grid%sice(I,J) .lt. 0.1) then write(message,*) 'very moist on land point: ', I,J,grid%smc(I,1,J) CALL wrf_debug(10,message) endif enddo enddo !!! compute grid%emt, grid%em on global domain, and only on task 0. #ifdef DM_PARALLEL IF (wrf_dm_on_monitor()) THEN !!!! NECESSARY TO LIMIT THIS TO TASK ZERO? #else IF (JDS .eq. JTS) THEN !! set unfailable condition for serial job #endif ALLOCATE(EMJ(JDS:JDE-1),EMTJ(JDS:JDE-1)) DO J=JDS,JDE-1 TPH=SB+float(J-1)*DPH DXP=ERAD*DLM*COS(TPH) EMJ(J)= GRID%DT/( 4.*DXP)*DTAD EMTJ(J)=GRID%DT/(16.*DXP)*DTAD ENDDO JA=0 DO 161 J=3,5 JA=JA+1 KHLA(JA)=2 KHHA(JA)=(IDE-1)-1-MOD(J+1,2) 161 grid%emt(JA)=EMTJ(J) DO 162 J=(JDE-1)-4,(JDE-1)-2 JA=JA+1 KHLA(JA)=2 KHHA(JA)=(IDE-1)-1-MOD(J+1,2) 162 grid%emt(JA)=EMTJ(J) DO 163 J=6,(JDE-1)-5 JA=JA+1 KHLA(JA)=2 KHHA(JA)=2+MOD(J,2) 163 grid%emt(JA)=EMTJ(J) DO 164 J=6,(JDE-1)-5 JA=JA+1 KHLA(JA)=(IDE-1)-2 KHHA(JA)=(IDE-1)-1-MOD(J+1,2) 164 grid%emt(JA)=EMTJ(J) ! --------------SPREADING OF UPSTREAM VELOCITY-POINT ADVECTION FACTOR---- JA=0 DO 171 J=3,5 JA=JA+1 KVLA(JA)=2 KVHA(JA)=(IDE-1)-1-MOD(J,2) 171 grid%em(JA)=EMJ(J) DO 172 J=(JDE-1)-4,(JDE-1)-2 JA=JA+1 KVLA(JA)=2 KVHA(JA)=(IDE-1)-1-MOD(J,2) 172 grid%em(JA)=EMJ(J) DO 173 J=6,(JDE-1)-5 JA=JA+1 KVLA(JA)=2 KVHA(JA)=2+MOD(J+1,2) 173 grid%em(JA)=EMJ(J) DO 174 J=6,(JDE-1)-5 JA=JA+1 KVLA(JA)=(IDE-1)-2 KVHA(JA)=(IDE-1)-1-MOD(J,2) 174 grid%em(JA)=EMJ(J) 696 continue ENDIF ! wrf_dm_on_monitor/serial job call NMM_SH2O(IMS,IME,JMS,JME,ITS,NNXP,JTS,NNYP,grid%num_soil_layers,grid%isltyp, & grid%sm,grid%sice,grid%stc,grid%smc,grid%sh2o) !! must be a better place to put this, but will eliminate "phantom" !! wind points here (no wind point on eastern boundary of odd numbered rows) IF ( abs(IDE-1-ITE) .eq. 1 ) THEN ! along eastern boundary write(message,*) 'zero phantom winds' CALL wrf_message(message) DO K=1,KDE-1 DO J=JDS,JDE-1,2 IF (J .ge. JTS .and. J .le. JTE) THEN grid%u(IDE-1,J,K)=0. grid%v(IDE-1,J,K)=0. ENDIF ENDDO ENDDO ENDIF 969 continue DO j = jms, jme DO i = ims, ime fisx=max(grid%fis(i,j),0.) grid%z0(I,J) =grid%sm(I,J)*Z0SEA+(1.-grid%sm(I,J))* & & (0.*Z0MAX+FISx *FCM+Z0LAND) ENDDO ENDDO write(message,*) 'grid%z0 over memory, leaving module_initialize_real' CALL wrf_message(message) DO J=JME,JMS,-((JME-JMS)/20+1) write(message,635) (grid%z0(I,J),I=IMS,IME,(IME-IMS)/14+1) CALL wrf_message(message) ENDDO endif ! on first_time check write(message,*) 'leaving init_domain_nmm' CALL wrf_message( TRIM(message) ) ! write(message,*)'STUFF MOVED TO REGISTRY:',grid%IDTAD, & & grid%NSOIL,grid%NRADL,grid%NRADS,grid%NPHS,grid%NCNVC,grid%sigma CALL wrf_message( TRIM(message) ) #if ( HWRF == 1 ) !========================================================================================= ! gopal's doing for ocean coupling. Produce a high resolution grid for the entire domain !========================================================================================= if(internal_time_loop.eq.1) then !Kwon's doing NDLMD=grid%dlmd/3. NDPHD=grid%dphd/3. NIDE=3*(IDE-1)-2 NJDE=3*(JDE-1)-2 ILOC=1 JLOC=1 NWBD= WBD ! + (ILOC -1)*2.*grid%dlmd + MOD(JLOC+1,2)*grid%dlmd NSBD= SBD ! + (JLOC -1)*grid%dphd ALLOCATE (NHLAT(NIDE,NJDE)) ALLOCATE (NHLON(NIDE,NJDE)) ALLOCATE (NVLAT(NIDE,NJDE)) ALLOCATE (NVLON(NIDE,NJDE)) ALLOCATE (HRES_SM(NIDE,NJDE)) #if defined(DM_PARALLEL) if(wrf_dm_on_monitor()) then ! Only the monitor process does the actual work (kinda ! stupid; should be parallelized, but it's better than ! writing garbage like it did before with >1 process) ! Get high-res lat & lon: CALL EARTH_LATLON_hwrf ( NHLAT,NHLON,NVLAT,NVLON, & ! rotated lat,lon at H and V points NDLMD,NDPHD,NWBD,NSBD, & tph0d,tlm0d, & 1,NIDE,1,NJDE,1,1, & 1,NIDE,1,NJDE,1,1, & 1,NIDE,1,NJDE,1,1 ) ! Interpolate landmask to high-res grid: CALL G2T2H_hwrf ( SM_G,HRES_SM, & ! output grid index and weights NHLAT,NHLON, & ! target (hres) input lat lon in degrees grid%DLMD,grid%DPHD,WBD,SBD, & ! parent res, west and south boundaries tph0d,tlm0d, & ! parent central lat,lon, all in degrees IDE,JDE,IDS,IDE,JDS,JDE, & ! parent imax and jmax, ime,jme 1,NIDE,1,NJDE,1,1, & 1,NIDE,1,NJDE,1,1, & 1,NIDE,1,NJDE,1,1 ) ! We're done with the low-res sm grid now: deallocate(SM_G) ! Write out high-res grid for coupler: WRITE(65)NHLAT(1:NIDE,1:NJDE) WRITE(65)NHLON(1:NIDE,1:NJDE) WRITE(65)NVLAT(1:NIDE,1:NJDE) WRITE(65)NVLON(1:NIDE,1:NJDE) WRITE(65)HRES_SM(1:NIDE,1:NJDE) endif #else ! This code is the same as above, but for the non-mpi version: CALL EARTH_LATLON_hwrf ( NHLAT,NHLON,NVLAT,NVLON, & ! rotated lat,lon at H and V points NDLMD,NDPHD,NWBD,NSBD, & tph0d,tlm0d, & 1,NIDE,1,NJDE,1,1, & 1,NIDE,1,NJDE,1,1, & 1,NIDE,1,NJDE,1,1 ) CALL G2T2H_hwrf ( grid%SM,HRES_SM, & ! output grid index and weights NHLAT,NHLON, & ! target (hres) input lat lon in degrees grid%DLMD,grid%DPHD,WBD,SBD, & ! parent res, west and south boundaries tph0d,tlm0d, & ! parent central lat,lon, all in degrees IDE,JDE,IMS,IME,JMS,JME, & ! parent imax and jmax, ime,jme 1,NIDE,1,NJDE,1,1, & 1,NIDE,1,NJDE,1,1, & 1,NIDE,1,NJDE,1,1 ) WRITE(65)NHLAT(1:NIDE,1:NJDE) WRITE(65)NHLON(1:NIDE,1:NJDE) WRITE(65)NVLAT(1:NIDE,1:NJDE) WRITE(65)NVLON(1:NIDE,1:NJDE) WRITE(65)HRES_SM(1:NIDE,1:NJDE) #endif DEALLOCATE (NHLAT) DEALLOCATE (NHLON) DEALLOCATE (NVLAT) DEALLOCATE (NVLON) DEALLOCATE (HRES_SM) endif !Kwon's doing !================================================================================== ! end gopal's doing for ocean coupling. !================================================================================== #endif !#define COPY_OUT !#include "scalar_derefs.inc" RETURN END SUBROUTINE init_domain_nmm !------------------------------------------------------ SUBROUTINE define_nmm_vertical_coord ( LM, PTSGM, pt, pdtop,HYBLEVS, & SG1,DSG1,SGML1, & SG2,DSG2,SGML2,dfl, dfrlg ) IMPLICIT NONE ! USE module_model_constants !!! certain physical parameters here probably don't need to be defined, as defined !!! elsewhere within WRF. Done for initial testing purposes. INTEGER :: LM, LPT2, L REAL :: PTSGM, pt, PL, PT2, pdtop REAL :: RGOG, PSIG,PHYB,PHYBM REAL, PARAMETER :: Rd = 287.04 ! J deg{-1} kg{-1} REAL, PARAMETER :: CP=1004.6,GAMMA=.0065,PRF0=101325.,T0=288. REAL, PARAMETER :: g=9.81 REAL, DIMENSION(LM) :: DSG,DSG1,DSG2 REAL, DIMENSION(LM) :: SGML1,SGML2 REAL, DIMENSION(LM+1) :: SG1,SG2,HYBLEVS,dfl,dfrlg CHARACTER(LEN=255) :: message LPT2=LM+1 write(message,*) 'pt= ', pt CALL wrf_message(message) DO L=LM+1,1,-1 pl=HYBLEVS(L)*(101325.-pt)+pt if(pl.lt.ptSGm) LPT2=l ENDDO IF(LPT2.lt.LM+1) THEN pt2=HYBLEVS(LPT2)*(101325.-pt)+pt ELSE pt2=pt ENDIF write(message,*) '*** Sigma system starts at ',pt2,' Pa, from level ',LPT2 CALL wrf_message(message) pdtop=pt2-pt write(message,*) 'allocating DSG,DSG1,DSG2 as ', LM CALL wrf_debug(10,message) DSG=-99. DO L=1,LM DSG(L)=HYBLEVS(L)- HYBLEVS(L+1) ENDDO DSG1=0. DSG2=0. DO L=LM,1,-1 IF(L.ge.LPT2) then DSG1(L)=DSG(L) ELSE DSG2(L)=DSG(L) ENDIF ENDDO SGML1=-99. SGML2=-99. IF(LPT2.le.LM+1) THEN DO L=LM+1,LPT2,-1 SG2(L)=0. ENDDO DO L=LPT2,2,-1 SG2(L-1)=SG2(L)+DSG2(L-1) ENDDO DO L=LPT2-1,1,-1 SG2(L)=SG2(L)/SG2(1) ENDDO SG2(1)=1. DO L=LPT2-1,1,-1 DSG2(L)=SG2(L)-SG2(L+1) SGML2(l)=(SG2(l)+SG2(l+1))*0.5 ENDDO ENDIF DO L=LM,LPT2,-1 DSG2(L)=0. SGML2(L)=0. ENDDO !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SG1(LM+1)=0. DO L=LM+1,LPT2,-1 SG1(L-1)=SG1(L)+DSG1(L-1) ENDDO DO L=LM,LPT2,-1 SG1(L)=SG1(L)/SG1(LPT2-1) ENDDO SG1(LPT2-1)=1. do l=LPT2-2,1,-1 SG1(l)=1. enddo DO L=LM,LPT2,-1 DSG1(L)=SG1(L)-SG1(L+1) SGML1(L)=(SG1(L)+SG1(L+1))*0.5 ENDDO DO L=LPT2-1,1,-1 DSG1(L)=0. SGML1(L)=1. ENDDO 1000 format('l,hyblevs,psig,SG1,SG2,phyb,phybm') 1100 format(' ',i4,f7.4,f10.2,2f7.4,2f10.2) write(message,1000) CALL wrf_debug(100,message) do l=1,LM+1 psig=HYBLEVS(L)*(101325.-pt)+pt phyb=SG1(l)*pdtop+SG2(l)*(101325.-pdtop-pt)+pt if(l.lt.LM+1) then phybm=SGML1(l)*pdtop+SGML2(l)*(101325.-pdtop-pt)+pt else phybm=-99. endif write(message,1100) l,HYBLEVS(L),psig & ,SG1(l),SG2(l),phyb,phybm CALL wrf_debug(100,message) enddo 632 format(f9.6) write(message,*) 'SG1' CALL wrf_debug(100,message) do L=LM+1,1,-1 write(message,632) SG1(L) CALL wrf_debug(100,message) enddo write(message,*) 'SG2' CALL wrf_debug(100,message) do L=LM+1,1,-1 write(message,632) SG2(L) CALL wrf_debug(100,message) enddo write(message,*) 'DSG1' CALL wrf_debug(100,message) do L=LM,1,-1 write(message,632) DSG1(L) CALL wrf_debug(100,message) enddo write(message,*) 'DSG2' CALL wrf_debug(100,message) do L=LM,1,-1 write(message,632) DSG2(L) CALL wrf_debug(100,message) enddo write(message,*) 'SGML1' CALL wrf_debug(100,message) do L=LM,1,-1 write(message,632) SGML1(L) CALL wrf_debug(100,message) enddo write(message,*) 'SGML2' CALL wrf_debug(100,message) do L=LM,1,-1 write(message,632) SGML2(L) CALL wrf_debug(100,message) enddo rgog=(rd*gamma)/g DO L=1,LM+1 dfl(L)=g*T0*(1.-((pt+SG1(L)*pdtop+SG2(L)*(101325.-pt2)) & /101325.)**rgog)/gamma dfrlg(L)=dfl(L)/g write(message,*) 'L, dfl(L): ', L, dfl(L) CALL wrf_debug(100,message) ENDDO END SUBROUTINE define_nmm_vertical_coord !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE compute_nmm_surfacep ( TERRAIN_HGT_T, Z3D_IN, PRESS3D_IN, T3D_IN & &, psfc_out,generic & &, IDS,IDE,JDS,JDE,KDS,KDE & &, IMS,IME,JMS,JME,KMS,KME & &, ITS,ITE,JTS,JTE,KTS,KTE ) IMPLICIT NONE real, allocatable:: dum2d(:,:),DUM2DB(:,:) integer :: IDS,IDE,JDS,JDE,KDS,KDE integer :: IMS,IME,JMS,JME,KMS,KME integer :: ITS,ITE,JTS,JTE,KTS,KTE,Ilook,Jlook integer :: I,J,II,generic,L,KINSERT,K,bot_lev,LL integer :: IHE(JMS:JME),IHW(JMS:JME) real :: TERRAIN_HGT_T(IMS:IME,JMS:JME) real :: Z3D_IN(IMS:IME,JMS:JME,generic) real :: T3D_IN(IMS:IME,JMS:JME,generic) real :: PRESS3D_IN(IMS:IME,JMS:JME,generic) real :: PSFC_IN(IMS:IME,JMS:JME),TOPO_IN(IMS:IME,JMS:JME) real :: psfc_out(IMS:IME,JMS:JME),rincr(IMS:IME,JMS:JME) real :: dif1,dif2,dif3,dif4,dlnpdz,BOT_INPUT_HGT,BOT_INPUT_PRESS,dpdz,rhs real :: zin(generic),pin(generic) character (len=255) :: message logical :: DEFINED_PSFC(IMS:IME,JMS:JME), DEFINED_PSFCB(IMS:IME,JMS:JME) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Ilook=25 Jlook=25 DO j = JMS, JME IHE(J)=MOD(J+1,2) IHW(J)=IHE(J)-1 ENDDO DO J=JMS,JME DO I=IMS,IME DEFINED_PSFC(I,J)=.FALSE. DEFINED_PSFCB(I,J)=.FALSE. IF (PRESS3D_IN(I,J,1) .ne. 200100.) THEN PSFC_IN(I,J)=PRESS3D_IN(I,J,1) TOPO_IN(I,J)=Z3D_IN(I,J,1) ELSE PSFC_IN(I,J)=PRESS3D_IN(I,J,2) TOPO_IN(I,J)=Z3D_IN(I,J,2) ENDIF ENDDO ENDDO ! input surface pressure smoothing over the ocean - still needed for NAM? II_loop: do II=1,8 CYCLE II_loop do J=JTS+1,min(JTE,JDE-1)-1 do I=ITS+1,min(ITE,IDE-1)-1 rincr(I,J)=0. if (PSFC_IN(I,J) .gt. 100000. .and. & PSFC_IN(I+IHE(J),J+1) .gt. 100000. .and. & PSFC_IN(I+IHE(J),J-1) .gt. 100000. .and. & PSFC_IN(I+IHW(J),J+1) .gt. 100000. .and. & PSFC_IN(I+IHW(J),J-1) .gt. 100000. ) then dif1=abs(PSFC_IN(I,J)-PSFC_IN(I+IHE(J),J+1)) dif2=abs(PSFC_IN(I,J)-PSFC_IN(I+IHE(J),J-1)) dif3=abs(PSFC_IN(I,J)-PSFC_IN(I+IHW(J),J+1)) dif4=abs(PSFC_IN(I,J)-PSFC_IN(I+IHW(J),J-1)) if (max(dif1,dif2,dif3,dif4) .lt. 200. .and. TOPO_IN(I,J).le. 0.5 .and. & TOPO_IN(I+IHE(J),J+1) .le. 0.5 .and. & TOPO_IN(I+IHW(J),J+1) .le. 0.5 .and. & TOPO_IN(I+IHE(J),J-1) .le. 0.5 .and. & TOPO_IN(I+IHW(J),J-1) .lt. 0.5) then rincr(I,J)=0.125*( 4.*PSFC_IN(I,J)+ & PSFC_IN(I+IHE(J),J+1)+PSFC_IN(I+IHE(J),J-1)+ & PSFC_IN(I+IHW(J),J+1)+PSFC_IN(I+IHW(J),J-1) ) & - PSFC_IN(I,J) ! if (rincr(I,J) .ne. 0 .and. abs(rincr(I,J)) .gt. 20.) then ! write(message,*) 'II, I,J,rincr: ', II, I,J,rincr(I,J) ! CALL wrf_message(message) ! endif endif endif ENDDO ENDDO DO J=JTS+1,min(JTE,JDE-1)-1 DO I=ITS+1,min(ITE,IDE-1)-1 PSFC_IN(I,J)=PSFC_IN(I,J) + rincr(I,J) ENDDO ENDDO ! write(message,*) ' -------------------------------------------------- ' ! CALL wrf_message(message) end do II_loop ALLOCATE(DUM2D(IMS:IME,JMS:JME)) DO J=JMS,JME DO I=IMS,IME DUM2D(I,J)=-9. END DO END DO DO J=JTS,min(JTE,JDE-1) I_loop: DO I=ITS,min(ITE,IDE-1) IF (PSFC_IN(I,J) .lt. 0.1) THEN write(message,*) 'QUITTING BECAUSE I,J, PSFC_IN: ', I,J,PSFC_IN(I,J) call wrf_error_fatal(message) ENDIF BOT_INPUT_PRESS=PSFC_IN(I,J) BOT_INPUT_HGT=TOPO_IN(I,J) IF (I .eq. Ilook .AND. J .eq. Jlook) THEN write(message,*) ' TERRAIN_HGT_T: ', I,J, TERRAIN_HGT_T(I,J) CALL wrf_message(message) write(message,*) ' PSFC_IN, TOPO_IN: ', & I, J, PSFC_IN(I,J),TOPO_IN(I,J) CALL wrf_message(message) DO L=1,generic write(message,*) ' L,PRESS3D_IN, Z3D_IN: ', & I,J,L, PRESS3D_IN(I,J,L),Z3D_IN(I,J,L) CALL wrf_debug(10,message) END DO ENDIF DO L=2,generic-1 IF ( PRESS3D_IN(i,j,L) .gt. PSFC_IN(I,J) .AND. & Z3D_IN(I,J,L) .lt. TERRAIN_HGT_T(I,J) .AND. & Z3D_IN(I,J,L+1) .gt. TERRAIN_HGT_T(I,J) ) THEN BOT_INPUT_PRESS=PRESS3D_IN(i,j,L) BOT_INPUT_HGT=Z3D_IN(I,J,L) ! IF (I .eq. Ilook .and. J .eq. Jlook) THEN ! write(message,*) 'BOT_INPUT_PRESS, BOT_INPUT_HGT NOW : ', & ! Ilook,Jlook, BOT_INPUT_PRESS, BOT_INPUT_HGT ! CALL wrf_message(message) ! ENDIF ENDIF END DO !!!!!!!!!!!!!!!!!!!!!! START HYDRO CHECK IF ( PRESS3D_IN(i,j,1) .ne. 200100. .AND. & (PSFC_IN(I,J) .gt. PRESS3D_IN(i,j,2) .OR. & TOPO_IN(I,J) .lt. Z3D_IN(I,J,2)) ) THEN ! extrapolate downward IF (J .eq. JTS .AND. I .eq. ITS) THEN write(message,*) 'hydro check - should only be for isobaric input' CALL wrf_message(message) ENDIF IF (Z3D_IN(I,J,2) .ne. TOPO_IN(I,J)) THEN dpdz=(PRESS3D_IN(i,j,2)-PSFC_IN(I,J))/(Z3D_IN(I,J,2)-TOPO_IN(I,J)) rhs=-9.81*((PRESS3D_IN(i,j,2)+ PSFC_IN(I,J))/2.)/(287.04* T3D_IN(I,J,2)) IF ( abs(PRESS3D_IN(i,j,2)-PSFC_IN(I,J)) .gt. 290.) THEN IF (dpdz .lt. 1.05*rhs .OR. dpdz .gt. 0.95*rhs) THEN write(message,*) 'I,J,P(2),Psfc,Z(2),Zsfc: ', & I,J,PRESS3D_IN(i,j,2),PSFC_IN(I,J),Z3D_IN(I,J,2),TOPO_IN(I,J) IF (mod(I,5).eq.0 .AND. mod(J,5).eq.0) CALL wrf_debug(50,message) CYCLE I_loop ENDIF ENDIF ELSE ! z(2) equals TOPO_IN IF (PRESS3D_IN(i,j,2) .eq. PSFC_IN(I,J)) THEN ! write(message,*) 'all equal at I,J: ', I,J ! CALL wrf_message(message) ELSE ! write(message,*) 'heights equal, pressures not: ', & ! PRESS3D_IN(i,j,2), PSFC_IN(I,J) ! CALL wrf_message(message) CYCLE I_loop ENDIF ENDIF IF ( abs(PRESS3D_IN(i,j,2)-PSFC_IN(I,J)) .gt. 290.) THEN IF (PRESS3D_IN(i,j,2) .lt. PSFC_IN(I,J) .and. & Z3D_IN(I,J,2) .lt. TOPO_IN(I,J)) THEN ! write(message,*) 'surface data mismatch(a) at I,J: ', I,J ! CALL wrf_message(message) CYCLE I_loop ELSEIF (PRESS3D_IN(i,j,2) .gt. PSFC_IN(I,J) .AND. & Z3D_IN(I,J,2) .gt. TOPO_IN(I,J)) THEN ! write(message,*) 'surface data mismatch(b) at I,J: ', I,J ! CALL wrf_message(message) CYCLE I_loop ENDIF ENDIF ENDIF !!!!!!! loop over a few more levels DO L=3,6 IF ( PRESS3D_IN(i,j,1) .ne. 200100. .AND. & (((PSFC_IN(I,J)-PRESS3D_IN(i,j,L)) .lt. 400.) .OR. & TOPO_IN(I,J) .lt. Z3D_IN(I,J,L))) then IF (Z3D_IN(I,J,L) .ne. TOPO_IN(I,J)) THEN dpdz=(PRESS3D_IN(i,j,L)-PSFC_IN(I,J))/ & (Z3D_IN(I,J,L)-TOPO_IN(I,J)) rhs=-9.81*((PRESS3D_IN(i,j,L)+ PSFC_IN(I,J))/2.)/ & (287.04*T3D_IN(I,J,L)) IF ( abs(PRESS3D_IN(i,j,L)-PSFC_IN(I,J)) .gt. 290.) THEN IF (dpdz .lt. 1.05*rhs .or. dpdz .gt. 0.95*rhs) THEN write(message,*) 'I,J,L,Piso,Psfc,Ziso,Zsfc: ', & I,J,L,PRESS3D_IN(i,j,L),PSFC_IN(I,J),& Z3D_IN(I,J,L),TOPO_IN(I,J) IF (mod(I,5).eq.0 .AND. mod(J,5).eq.0) & CALL wrf_debug(50,message) CYCLE I_loop ENDIF ENDIF ELSE IF (PRESS3D_IN(i,j,2) .eq. PSFC_IN(I,J)) THEN ! write(message,*) 'all equal at I,J: ', I,J ! CALL wrf_message(message) ELSE CYCLE I_loop ENDIF ENDIF ENDIF IF ( abs(PRESS3D_IN(i,j,L)-PSFC_IN(I,J)) .gt. 290.) THEN IF (PRESS3D_IN(i,j,L) .lt. PSFC_IN(I,J) .AND. & Z3D_IN(I,J,L) .lt. TOPO_IN(I,J)) THEN CYCLE I_loop ELSEIF (PRESS3D_IN(i,j,L) .gt. PSFC_IN(I,J) .AND. & Z3D_IN(I,J,L) .gt. TOPO_IN(I,J)) THEN CYCLE I_loop ENDIF ENDIF END DO !!!!!!!!!!!!!!!!!!!!!! END HYDRO CHECK IF (TERRAIN_HGT_T(I,J) .eq. BOT_INPUT_HGT ) THEN dum2d(I,J)=BOT_INPUT_PRESS IF (BOT_INPUT_HGT .ne. 0. .and. (BOT_INPUT_HGT-INT(BOT_INPUT_HGT) .ne. 0.) ) THEN write(message,*) 'with BOT_INPUT_HGT: ', BOT_INPUT_HGT, & 'set dum2d to bot_input_pres: ', I,J,dum2d(I,J) CALL wrf_message(message) ENDIF IF (dum2d(I,J) .lt. 50000. .OR. dum2d(I,J) .gt. 109000.) THEN write(message,*) 'bad dum2d(a): ', I,J,DUM2D(I,J) CALL wrf_message(message) ENDIF ELSEIF (TERRAIN_HGT_T(I,J) .lt. BOT_INPUT_HGT ) THEN ! target is below lowest possible input...extrapolate IF ( BOT_INPUT_PRESS-PRESS3D_IN(I,J,2) .gt. 500. ) THEN dlnpdz= (log(BOT_INPUT_PRESS)-log(PRESS3D_IN(i,j,2)) ) / & (BOT_INPUT_HGT-Z3D_IN(i,j,2)) IF (I .eq. Ilook .and. J .eq. Jlook) THEN write(message,*) 'I,J,dlnpdz(a): ', I,J,dlnpdz CALL wrf_message(message) ENDIF ELSE !! thin layer and/or just have lowest level - difference with 3rd level data IF ( abs(BOT_INPUT_PRESS - PRESS3D_IN(i,j,3)) .gt. 290. ) THEN dlnpdz= (log(BOT_INPUT_PRESS)-log(PRESS3D_IN(i,j,3)) ) / & (BOT_INPUT_HGT-Z3D_IN(i,j,3)) IF (I .eq. Ilook .and. J .eq. Jlook) then write(message,*) 'p diff: ', BOT_INPUT_PRESS, PRESS3D_IN(i,j,3) CALL wrf_message(message) write(message,*) 'z diff: ', BOT_INPUT_HGT, Z3D_IN(i,j,3) CALL wrf_message(message) ENDIF ELSE !! Loop up to level 7 looking for a sufficiently thick layer FIND_THICK: DO LL=4,7 IF( abs(BOT_INPUT_PRESS - PRESS3D_IN(i,j,LL)) .gt. 290.) THEN dlnpdz= (log(BOT_INPUT_PRESS)-log(PRESS3D_IN(i,j,LL)) ) / & (BOT_INPUT_HGT-Z3D_IN(i,j,LL)) EXIT FIND_THICK ENDIF END DO FIND_THICK ENDIF ENDIF dum2d(I,J)= exp(log(BOT_INPUT_PRESS) + dlnpdz * & (TERRAIN_HGT_T(I,J) - BOT_INPUT_HGT) ) IF (dum2d(I,J) .lt. 50000. .or. dum2d(I,J) .gt. 108000.) THEN write(message,*) 'bad dum2d(b): ', I,J,DUM2D(I,J) CALL wrf_message(message) write(message,*) 'BOT_INPUT_PRESS, dlnpdz, TERRAIN_HGT_T, BOT_INPUT_HGT: ', & BOT_INPUT_PRESS, dlnpdz, TERRAIN_HGT_T(I,J), BOT_INPUT_HGT CALL wrf_message(message) write(message,*) 'Z3D_IN: ', Z3D_IN(I,J,1:10) CALL wrf_message(message) write(message,*) 'PRESS3D_IN: ', PRESS3D_IN(I,J,1:10) CALL wrf_message(message) ENDIF ELSE ! target level bounded by input levels DO L=2,generic-1 IF (TERRAIN_HGT_T(I,J) .gt. Z3D_IN(i,j,L) .AND. & TERRAIN_HGT_T(I,J) .lt. Z3D_IN(i,j,L+1) ) THEN dlnpdz= (log(PRESS3D_IN(i,j,l))-log(PRESS3D_IN(i,j,L+1)) ) / & (Z3D_IN(i,j,l)-Z3D_IN(i,j,L+1)) dum2d(I,J)= log(PRESS3D_IN(i,j,l)) + & dlnpdz * (TERRAIN_HGT_T(I,J) - Z3D_IN(i,j,L) ) dum2d(i,j)=exp(dum2d(i,j)) IF (dum2d(I,J) .lt. 50000. .or. dum2d(I,J) .gt. 108000.) THEN write(message,*) 'bad dum2d(c): ', I,J,DUM2D(I,J) CALL wrf_message(message) ENDIF ENDIF ENDDO !!! account for situation where BOT_INPUT_HGT < TERRAIN_HGT_T < Z3D_IN(:,2,:) IF (dum2d(I,J) .eq. -9 .AND. BOT_INPUT_HGT .lt. TERRAIN_HGT_T(I,J) & .AND. TERRAIN_HGT_T(I,J) .lt. Z3D_IN(I,J,2)) then IF (mod(I,50) .eq. 0 .AND. mod(J,50) .eq. 0) THEN write(message,*) 'I,J,BOT_INPUT_HGT, bot_pres, TERRAIN_HGT_T: ', & I,J,BOT_INPUT_HGT, BOT_INPUT_PRESS, TERRAIN_HGT_T(I,J) CALL wrf_message(message) ENDIF dlnpdz= (log(PSFC_IN(i,j))-log(PRESS3D_IN(i,j,2)) ) / & (TOPO_IN(i,j)-Z3D_IN(i,j,2)) dum2d(I,J)= log(PSFC_IN(i,j)) + & dlnpdz * (TERRAIN_HGT_T(I,J) - TOPO_IN(i,j) ) dum2d(i,j)= exp(dum2d(i,j)) IF (dum2d(I,J) .lt. 50000. .or. dum2d(I,J) .gt. 108000.) THEN write(message,*) 'bad dum2d(d): ', I,J,DUM2D(I,J) CALL wrf_message(message) ENDIF ENDIF IF (dum2d(I,J) .eq. -9.) THEN write(message,*) 'must have flukey situation in new ', I,J CALL wrf_message(message) write(message,*) 'I,J,BOT_INPUT_HGT, bot_pres, TERRAIN_HGT_T: ', & I,J,BOT_INPUT_HGT, BOT_INPUT_PRESS, TERRAIN_HGT_T(I,J) CALL wrf_message(message) DO L=1,generic-1 IF ( TERRAIN_HGT_T(I,J) .eq. Z3D_IN(i,j,L) ) THEN ! problematic with HGT_M substitution for "input" surface height? dum2d(i,j)=PRESS3D_IN(I,J,L) IF (dum2d(I,J) .lt. 50000. .or. dum2d(I,J) .gt. 108000.) THEN write(message,*) 'bad dum2d(e): ', I,J,DUM2D(I,J) CALL wrf_message(message) ENDIF ENDIF ENDDO IF ( TERRAIN_HGT_T(I,J) .eq. TOPO_IN(I,J)) THEN dum2d(I,J)=PSFC_IN(I,J) IF (dum2d(I,J) .lt. 50000. .or. dum2d(I,J) .gt. 108000.) THEN write(message,*) 'bad dum2d(grid%f): ', I,J,DUM2D(I,J) CALL wrf_message(message) ENDIF write(message,*) 'matched input topo, psfc: ', I,J,TOPO_IN(I,J),PSFC_IN(I,J) CALL wrf_message(message) ENDIF IF (dum2d(I,J) .eq. -9.) THEN CALL wrf_error_fatal("quitting due to undefined surface pressure") ENDIF ENDIF DEFINED_PSFC(I,J)=.TRUE. IF (I .eq. Ilook .AND. J .eq. Jlook) THEN write(message,*) 'newstyle psfc: ', I,J,dum2d(I,J) CALL wrf_message(message) ENDIF ENDIF ENDDO I_loop ENDDO ! write(message,*) 'psfc points (new style)' ! CALL wrf_message(message) ! DO J=min(JTE,JDE-1),JTS,-loopinc ! write(message,633) (dum2d(I,J)/100.,I=ITS,min(ITE,IDE-1),iloopinc) ! CALL wrf_message(message) ! END DO 633 format(35(f5.0,1x)) write(message,*) 'PSFC extremes (new style)' CALL wrf_message(message) write(message,*) minval(dum2d,MASK=DEFINED_PSFC),maxval(dum2d,MASK=DEFINED_PSFC) CALL wrf_message(message) ! IF (minval(dum2d,MASK=DEFINED_PSFC) .lt. 50000. .or. maxval(dum2d,MASK=DEFINED_PSFC) .gt. 108000.) THEN DO J=JTS,min(JTE,JDE-1) DO I=ITS,min(ITE,IDE-1) IF (DEFINED_PSFC(I,J) .AND. dum2d(I,J) .lt. 50000. ) THEN IF (TERRAIN_HGT_T(I,J) .gt. 4500.) THEN WRITE(message,*) 'low surface pressure allowed because surface height is: ', TERRAIN_HGT_T(I,J) CALL wrf_debug(2,message) ELSE CALL wrf_error_fatal("quit due to unrealistic surface pressure") ENDIF ENDIF IF (DEFINED_PSFC(I,J) .AND. dum2d(I,J) .gt. 108000. ) THEN IF (TERRAIN_HGT_T(I,J) .lt. -10.) THEN WRITE(message,*) 'high surface pressure allowed because surface height is: ', TERRAIN_HGT_T(I,J) CALL wrf_debug(2,message) ELSE CALL wrf_error_fatal("quit due to unrealistic surface pressure") ENDIF ENDIF END DO END DO !! "traditional" isobaric only approach ------------------------------------------------ ALLOCATE (DUM2DB(IMS:IME,JMS:JME)) DO J=JMS,JME DO I=IMS,IME DUM2DB(I,J)=-9. END DO END DO DO J=JTS,min(JTE,JDE-1) DO I=ITS,min(ITE,IDE-1) IF (TERRAIN_HGT_T(I,J) .lt. Z3D_IN(i,j,2)) THEN ! targ below lowest IF ( abs(PRESS3D_IN(i,j,2)-PRESS3D_IN(i,j,3)) .gt. 290.) THEN dlnpdz= (log(PRESS3D_IN(i,j,2))-log(PRESS3D_IN(i,j,3)) ) / & (Z3D_IN(i,j,2)-Z3D_IN(i,j,3)) ELSE dlnpdz= (log(PRESS3D_IN(i,j,2))-log(PRESS3D_IN(i,j,4)) ) / & (Z3D_IN(i,j,2)-Z3D_IN(i,j,4)) ENDIF DUM2DB(I,J)= exp( log(PRESS3D_IN(i,j,2)) + dlnpdz * & (TERRAIN_HGT_T(I,J) - Z3D_IN(i,j,2)) ) IF (I .eq. Ilook .and. J .eq. Jlook) THEN write(message,*) 'I,K, trad: dlnpdz, press_in(2), terrain_t, Z3D_IN(2): ', I,J,dlnpdz, & PRESS3D_IN(i,j,2), TERRAIN_HGT_T(I,J), Z3D_IN(i,j,2) CALL wrf_message(message) ENDIF DEFINED_PSFCB(i,j)=.true. ELSEIF (TERRAIN_HGT_T(I,J) .gt. Z3D_IN(i,j,2)) THEN ! target level bounded by input levels DO L=2,generic-1 IF (TERRAIN_HGT_T(I,J) .gt. Z3D_IN(i,j,L) .AND. & TERRAIN_HGT_T(I,J) .lt. Z3D_IN(i,j,L+1) ) THEN dlnpdz= (log(PRESS3D_IN(i,j,l))-log(PRESS3D_IN(i,j,L+1)) ) / & (Z3D_IN(i,j,l)-Z3D_IN(i,j,L+1)) DUM2DB(I,J)= log(PRESS3D_IN(i,j,l)) + & dlnpdz * (TERRAIN_HGT_T(I,J) - Z3D_IN(i,j,L) ) DUM2DB(i,j)=exp(DUM2DB(i,j)) DEFINED_PSFCB(i,j)=.true. IF (DUM2DB(I,J) .lt. 13000.) THEN write(message,*) 'I,J,L,terrain,Z3d(L),z3d(L+1),p3d(L),p3d(l+1): ', I,J,L, & TERRAIN_HGT_T(I,J),Z3D_IN(I,J,L),Z3D_IN(I,J,L+1),PRESS3D_IN(I,J,L), & PRESS3D_IN(I,J,L+1) CALL wrf_error_fatal(message) ENDIF ENDIF ENDDO ELSEIF (TERRAIN_HGT_T(I,J) .eq. Z3D_IN(i,j,2)) THEN DUM2DB(i,j)=PRESS3D_IN(I,J,2) DEFINED_PSFCB(i,j)=.true. ENDIF IF (DUM2DB(I,J) .eq. -9.) THEN write(message,*) 'must have flukey situation in trad ', I,J CALL wrf_message(message) DO L=1,generic-1 IF ( TERRAIN_HGT_T(I,J) .eq. Z3D_IN(i,j,L) ) THEN DUM2DB(i,j)=PRESS3D_IN(I,J,L) DEFINED_PSFCB(i,j)=.true. ENDIF ENDDO ENDIF IF (DUM2DB(I,J) .eq. -9.) THEN write(message,*) 'HOPELESS PSFC, I QUIT' CALL wrf_error_fatal(message) ENDIF if (I .eq. Ilook .and. J .eq. Jlook) THEN write(message,*) ' traditional psfc: ', I,J,DUM2DB(I,J) CALL wrf_message(message) ENDIF ENDDO ENDDO ! write(message,*) 'psfc points (traditional)' ! CALL wrf_message(message) ! DO J=min(JTE,JDE-1),JTS,-loopinc ! write(message,633) (DUM2DB(I,J)/100.,I=its,min(ite,IDE-1),iloopinc) ! CALL wrf_message(message) ! ENDDO write(message,*) 'PSFC extremes (traditional)' CALL wrf_message(message) write(message,*) minval(DUM2DB,MASK=DEFINED_PSFCB),maxval(DUM2DB,MASK=DEFINED_PSFCB) CALL wrf_message(message) DO J=JTS,min(JTE,JDE-1) DO I=ITS,min(ITE,IDE-1) IF (DEFINED_PSFCB(I,J) .AND. dum2db(I,J) .lt. 50000. ) THEN IF (TERRAIN_HGT_T(I,J) .gt. 4500.) THEN WRITE(message,*) 'low surface pressure allowed because surface height is: ', TERRAIN_HGT_T(I,J) CALL wrf_debug(2,message) ELSE CALL wrf_error_fatal("quit due to unrealistic surface pressure") ENDIF ENDIF IF (DEFINED_PSFCB(I,J) .AND. dum2db(I,J) .gt. 108000. ) THEN IF (TERRAIN_HGT_T(I,J) .lt. -10.) THEN WRITE(message,*) 'high surface pressure allowed because surface height is: ', TERRAIN_HGT_T(I,J) CALL wrf_debug(2,message) ELSE CALL wrf_error_fatal("quit due to unrealistic surface pressure") ENDIF ENDIF ! IF (DEFINED_PSFCB(I,J) .AND. ( dum2db(I,J) .lt. 50000. .OR. dum2db(I,J) .gt. 108000. )) THEN ! IF (TERRAIN_HGT_T(I,J) .gt. -10. .and. TERRAIN_HGT_T(I,J) .lt. 5000.) THEN ! write(0,*) 'I,J, terrain_hgt_t, dum2db: ', I,J, terrain_hgt_t(I,J),dum2db(I,J) ! CALL wrf_error_fatal("quit due to unrealistic surface pressure") ! ELSE ! WRITE(message,*) 'surface pressure allowed because surface height is extreme value of: ', TERRAIN_HGT_T(I,J) ! CALL wrf_debug(2,message) ! ENDIF ! ENDIF ENDDO ENDDO !!!!! end traditional DO J=JTS,min(JTE,JDE-1) DO I=ITS,min(ITE,IDE-1) IF (DEFINED_PSFCB(I,J) .and. DEFINED_PSFC(I,J)) THEN IF ( abs(dum2d(I,J)-DUM2DB(I,J)) .gt. 400.) THEN write(message,*) 'BIG DIFF I,J, dum2d, DUM2DB: ', I,J,dum2d(I,J),DUM2DB(I,J) CALL wrf_message(message) ENDIF !! do we have enough confidence in new style to give it more than 50% weight? psfc_out(I,J)=0.5*(dum2d(I,J)+DUM2DB(I,J)) ELSEIF (DEFINED_PSFC(I,J)) THEN psfc_out(I,J)=dum2d(I,J) ELSEIF (DEFINED_PSFCB(I,J)) THEN psfc_out(I,J)=DUM2DB(I,J) ELSE write(message,*) 'I,J,dum2d,DUM2DB: ', I,J,dum2d(I,J),DUM2DB(I,J) CALL wrf_message(message) write(message,*) 'I,J,DEFINED_PSFC(I,J),DEFINED_PSFCB(I,J): ', I,J,DEFINED_PSFC(I,J),DEFINED_PSFCB(I,J) CALL wrf_message(message) call wrf_error_fatal("psfc_out completely undefined") ENDIF IF (I .eq. Ilook .AND. J .eq. Jlook) THEN write(message,*) ' combined psfc: ', I,J,psfc_out(I,J) CALL wrf_message(message) ENDIF IF (psfc_out(I,J) .lt. 50000. ) THEN IF (TERRAIN_HGT_T(I,J) .gt. 4500.) THEN WRITE(message,*) 'low surface pressure allowed because surface height is: ', TERRAIN_HGT_T(I,J) CALL wrf_debug(2,message) ELSE write(message,*) 'possibly bad combo on psfc_out: ', I,J, psfc_out(I,J) CALL wrf_debug(2,message) write(message,*) 'DEFINED_PSFC, dum2d: ', DEFINED_PSFC(I,J),dum2d(I,J) CALL wrf_debug(2,message) write(message,*) 'DEFINED_PSFCB, DUM2DB: ', DEFINED_PSFCB(I,J),DUM2DB(I,J) CALL wrf_debug(2,message) CALL wrf_error_fatal("quit due to unrealistic surface pressure") ENDIF ENDIF IF (psfc_out(I,J) .gt. 108000. ) THEN IF (TERRAIN_HGT_T(I,J) .lt. -10.) THEN WRITE(message,*) 'high surface pressure allowed because surface height is: ', TERRAIN_HGT_T(I,J) CALL wrf_debug(2,message) ELSE write(message,*) 'possibly bad combo on psfc_out: ', I,J, psfc_out(I,J) CALL wrf_debug(2,message) write(message,*) 'DEFINED_PSFC, dum2d: ', DEFINED_PSFC(I,J),dum2d(I,J) CALL wrf_debug(2,message) write(message,*) 'DEFINED_PSFCB, DUM2DB: ', DEFINED_PSFCB(I,J),DUM2DB(I,J) CALL wrf_debug(2,message) CALL wrf_error_fatal("quit due to unrealistic surface pressure") ENDIF ENDIF ENDDO ENDDO deallocate(dum2d,dum2db) END SUBROUTINE compute_nmm_surfacep !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE compute_3d_pressure(psfc_out,SGML1,SGML2,pdtop,pt & &, pd,p3d_out & &, IDS,IDE,JDS,JDE,KDS,KDE & &, IMS,IME,JMS,JME,KMS,KME & &, ITS,ITE,JTS,JTE,KTS,KTE ) INTEGER :: IDS,IDE,JDS,JDE,KDS,KDE INTEGER :: IMS,IME,JMS,JME,KMS,KME INTEGER :: ITS,ITE,JTS,JTE,KTS,KTE REAL, INTENT(IN) :: psfc_out(IMS:IME,JMS:JME) REAL, INTENT(IN) :: SGML1(KDE),SGML2(KDE),pdtop,pt REAL, INTENT(OUT):: p3d_out(IMS:IME,JMS:JME,KDS:KDE-1) REAL, INTENT(OUT):: pd(IMS:IME,JMS:JME) CHARACTER (len=255) :: message ! write(message,*) 'pdtop, pt, psfc_out(1,1): ', pdtop, pt, psfc_out(1,1) ! CALL wrf_message(message) DO J=JTS,min(JTE,JDE-1) DO I=ITS,min(ITE,IDE-1) pd(I,J)=psfc_out(I,J)-pdtop-pt ENDDO ENDDO DO J=JTS,min(JTE,JDE-1) DO K=KDS,KDE-1 DO I=ITS,min(ITE,IDE-1) p3d_out(I,J,K)=pd(I,J)*SGML2(K)+pdtop*SGML1(K)+pt IF (p3d_out(I,J,K) .ge. psfc_out(I,J) .or. p3d_out(I,J,K) .le. pt) THEN write(message,*) 'I,K,J,p3d_out: ', I,K,J,p3d_out(I,J,K) CALL wrf_error_fatal(message) ENDIF ENDDO ENDDO ENDDO END SUBROUTINE compute_3d_pressure !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE interp_press2press_lin(press_in,press_out, & data_in, data_out,generic & &, extrapolate,ignore_lowest,TFIELD & &, IDS,IDE,JDS,JDE,KDS,KDE & &, IMS,IME,JMS,JME,KMS,KME & &, ITS,ITE,JTS,JTE,KTS,KTE, internal_time ) ! Interpolates data from one set of pressure surfaces to ! another set of pressures INTEGER :: IDS,IDE,JDS,JDE,KDS,KDE INTEGER :: IMS,IME,JMS,JME,KMS,KME INTEGER :: ITS,ITE,JTS,JTE,KTS,KTE,generic INTEGER :: internal_time ! REAL, INTENT(IN) :: press_in(IMS:IME,generic,JMS:JME) REAL, INTENT(IN) :: press_in(IMS:IME,JMS:JME,generic) REAL, INTENT(IN) :: press_out(IMS:IME,JMS:JME,KDS:KDE-1) ! REAL, INTENT(IN) :: data_in(IMS:IME,generic,JMS:JME) REAL, INTENT(IN) :: data_in(IMS:IME,JMS:JME,generic) REAL, INTENT(OUT) :: data_out(IMS:IME,JMS:JME,KMS:KME) LOGICAL, INTENT(IN) :: extrapolate, ignore_lowest, TFIELD LOGICAL :: col_smooth INTEGER :: i,j INTEGER :: k,kk REAL :: desired_press REAL :: dvaldlnp,dlnp,tadiabat,tiso REAL, PARAMETER :: ADIAFAC=9.81/1004. REAL, PARAMETER :: TSTEXTRAPFAC=.0065 DO K=KMS,KME DO J=JMS,JME DO I=IMS,IME DATA_OUT(I,J,K)=-99999.9 ENDDO ENDDO ENDDO IF (ignore_lowest) then LMIN=2 ELSE LMIN=1 ENDIF DO j = JTS, min(JTE,JDE-1) test_i: DO i = ITS, min(ITE,IDE-1) IF (internal_time_loop .gt. 1) THEN IF (J .ne. JDS .and. J .ne. JDE-1 .and. & I .ne. IDS .and. I .ne. IDE-1 ) THEN !! not on boundary CYCLE test_i ENDIF ENDIF col_smooth=.false. output_loop: DO k = KDS,KDE-1 desired_press = press_out(i,j,k) if (K .gt. KDS) then if (TFIELD .and. col_smooth .and. desired_press .le. press_in(i,j,LMIN) & .and. press_out(i,j,k-1) .ge. press_in(i,j,LMIN)) then MAX_SMOOTH=K ! write(message,*) 'I,J, MAX_SMOOTH: ', I,J, MAX_SMOOTH ! CALL wrf_debug(100,message) endif endif ! keep track of where the extrapolation begins IF (desired_press .GT. press_in(i,j,LMIN)) THEN IF (TFIELD .and. K .eq. 1 .and. (desired_press - press_in(i,j,LMIN)) .gt. 3000.) then col_smooth=.TRUE. ! due to large extrapolation distance ENDIF IF ((desired_press - press_in(i,j,LMIN)).LT. 50.) THEN ! 0.5 mb data_out(i,j,k) = data_in(i,j,LMIN) ELSE IF (extrapolate) THEN ! Extrapolate downward because desired P level is below ! the lowest level in our input data. Extrapolate using simple ! 1st derivative of value with respect to ln P for the bottom 2 ! input layers. ! Add a check to make sure we are not using the gradient of ! a very thin layer if (TFIELD) then tiso=0.5*(data_in(i,j,1)+data_in(i,j,2)) endif IF ( (press_in(i,j,LMIN)-press_in(i,j,LMIN+1)) .GT. 500.) THEN ! likely isobaric data dlnp = log(press_in(i,j,LMIN))-log(press_in(i,j,LMIN+1)) dvaldlnp = (data_in(i,j,LMIN) - data_in(i,j,LMIN+1)) / dlnp ELSE ! assume terrain following dlnp = log(press_in(i,j,LMIN))-log(press_in(i,j,LMIN+5)) dvaldlnp = (data_in(i,j,LMIN) - data_in(i,j,LMIN+5)) / dlnp ENDIF data_out(i,j,k) = data_in(i,j,LMIN) + dvaldlnp * & ( log(desired_press)-log(press_in(i,j,LMIN)) ) if (TFIELD .and. data_out(i,j,k) .lt. tiso-0.2) then ! restrict slope to -1K/10 hPa dvaldlnp=max(dvaldlnp, -1.0/ & log( press_in(i,j,LMIN) / & ( press_in(i,j,LMIN)-1000.) )) data_out(I,J,K)= data_in(i,j,LMIN) + dvaldlnp * & ( log(desired_press)-log(press_in(i,j,LMIN)) ) elseif (TFIELD .and. data_out(i,j,k) .gt. tiso+0.2) then ! restrict slope to +0.8K/10 hPa dvaldlnp=min(dvaldlnp, 0.8/ & log( press_in(i,j,LMIN) / & ( press_in(i,j,LMIN)-1000.) )) data_out(I,J,K)= data_in(i,j,LMIN) + dvaldlnp * & ( log(desired_press)-log(press_in(i,j,LMIN)) ) endif ELSE data_out(i,j,k) = data_in(i,j,LMIN) ENDIF ENDIF ELSE IF (desired_press .LT. press_in(i,j,generic)) THEN IF ( (press_in(i,j,generic) - desired_press) .LT. 10.) THEN data_out(i,j,k) = data_in(i,j,generic) ELSE IF (extrapolate) THEN ! Extrapolate upward IF ((press_in(i,j,generic-1)-press_in(i,j,generic)).GT.50.) THEN dlnp =log(press_in(i,j,generic))-log(press_in(i,j,generic-1)) dvaldlnp=(data_in(i,j,generic)-data_in(i,j,generic-1))/dlnp ELSE dlnp =log(press_in(i,j,generic))-log(press_in(i,j,generic-2)) dvaldlnp=(data_in(i,j,generic)-data_in(i,j,generic-2))/dlnp ENDIF data_out(i,j,k) = data_in(i,j,generic) + & dvaldlnp * (log(desired_press)-log(press_in(i,j,generic))) ELSE data_out(i,j,k) = data_in(i,j,generic) ENDIF ENDIF ELSE ! We can trap between two levels and linearly interpolate input_loop: DO kk = LMIN, generic-1 IF (desired_press .EQ. press_in(i,j,kk) )THEN data_out(i,j,k) = data_in(i,j,kk) EXIT input_loop ELSE IF ( (desired_press .LT. press_in(i,j,kk)) .AND. & (desired_press .GT. press_in(i,j,kk+1)) ) THEN ! do trapped in lnp dlnp = log(press_in(i,j,kk)) - log(press_in(i,j,kk+1)) dvaldlnp = (data_in(i,j,kk)-data_in(i,j,kk+1))/dlnp data_out(i,j,k) = data_in(i,j,kk+1)+ & dvaldlnp*(log(desired_press)-log(press_in(i,j,kk+1))) EXIT input_loop ENDIF ENDDO input_loop ENDIF ENDDO output_loop if (col_smooth) then do K=max(KDS,MAX_SMOOTH-4),MAX_SMOOTH+4 data_out(I,J,K)=0.5*(data_out(I,J,K)+data_out(I,J,K+1)) enddo endif ENDDO test_i ENDDO END SUBROUTINE interp_press2press_lin SUBROUTINE wind_adjust(press_in,press_out, & U_in, V_in,U_out,V_out & &, generic,depth_replace & &, IDS,IDE,JDS,JDE,KDS,KDE & &, IMS,IME,JMS,JME,KMS,KME & &, ITS,ITE,JTS,JTE,KTS,KTE ) INTEGER :: IDS,IDE,JDS,JDE,KDS,KDE INTEGER :: IMS,IME,JMS,JME,KMS,KME INTEGER :: ITS,ITE,JTS,JTE,KTS,KTE,generic INTEGER :: MAXLIN,MAXLOUT REAL, INTENT(IN) :: press_in(IMS:IME,JMS:JME,generic) REAL, INTENT(IN) :: press_out(IMS:IME,JMS:JME,KDS:KDE-1) REAL, INTENT(IN) :: U_in(IMS:IME,JMS:JME,generic) REAL, INTENT(IN) :: V_in(IMS:IME,JMS:JME,generic) REAL, INTENT(INOUT) :: U_out(IMS:IME,KMS:KME,JMS:JME) REAL, INTENT(INOUT) :: V_out(IMS:IME,KMS:KME,JMS:JME) REAL :: p1d_in(generic) REAL :: p1d_out(KDS:KDE-1) DO j = JTS, min(JTE,JDE-1) DO i = ITS, min(ITE,IDE-1) ! IF (press_out(I,J,1) .lt. press_in(I,J,2)) then IF( (press_in(I,J,2)-press_out(I,J,1)) .gt. 200.) then U_out(I,1,J)=U_in(I,J,2) V_out(I,1,J)=V_in(I,J,2) INLOOP: DO L=2,generic p1d_in(L)=-9999. IF ( (press_in(I,J,2)-press_in(I,J,L)) .lt. depth_replace) THEN p1d_in(L)=(press_in(I,J,2)-press_in(I,J,L)) MAXLIN=L ELSE p1d_in(L)=(press_in(I,J,2)-press_in(I,J,L)) EXIT INLOOP ENDIF END DO INLOOP OUTLOOP: DO L=KDS,KDE-1 p1d_out(L)=-9999. IF ( (press_out(I,J,1)-press_out(I,J,L)) .lt. depth_replace) THEN p1d_out(L)=(press_out(I,J,1)-press_out(I,J,L)) MAXLOUT=L ELSE EXIT OUTLOOP ENDIF END DO OUTLOOP DO L=1,MAXLOUT ptarg=p1d_out(L) FINDLOOP: DO LL=2,MAXLIN if (p1d_in(LL) .lt. ptarg .and. p1d_in(LL+1) .gt. ptarg) then dlnp=log(p1d_in(LL))-log(p1d_in(LL+1)) dudlnp=(U_in(I,J,LL)-U_in(I,J,LL+1))/dlnp dvdlnp=(V_in(I,J,LL)-V_in(I,J,LL+1))/dlnp U_out(I,L,J)=U_in(I,J,LL)+dudlnp*(log(ptarg)-log(p1d_in(LL))) V_out(I,L,J)=V_in(I,J,LL)+dvdlnp*(log(ptarg)-log(p1d_in(LL))) EXIT FINDLOOP endif END DO FINDLOOP END DO ! MAXLOUT loop ENDIF ENDDO ENDDO END SUBROUTINE wind_adjust !-------------------------------------------------------------------- SUBROUTINE interp_press2press_log(press_in,press_out, & data_in, data_out, generic & &, extrapolate,ignore_lowest & &, IDS,IDE,JDS,JDE,KDS,KDE & &, IMS,IME,JMS,JME,KMS,KME & &, ITS,ITE,JTS,JTE,KTS,KTE, internal_time ) ! Interpolates ln(data) from one set of pressure surfaces to ! another set of pressures INTEGER :: IDS,IDE,JDS,JDE,KDS,KDE INTEGER :: IMS,IME,JMS,JME,KMS,KME INTEGER :: ITS,ITE,JTS,JTE,KTS,KTE,generic INTEGER :: internal_time ! REAL, INTENT(IN) :: press_in(IMS:IME,generic,JMS:JME) REAL, INTENT(IN) :: press_in(IMS:IME,JMS:JME,generic) REAL, INTENT(IN) :: press_out(IMS:IME,JMS:JME,KDS:KDE-1) ! REAL, INTENT(IN) :: data_in(IMS:IME,generic,JMS:JME) ! REAL, INTENT(IN) :: data_in(IMS:IME,JMS:JME,generic) REAL :: data_in(IMS:IME,JMS:JME,generic) REAL, INTENT(OUT) :: data_out(IMS:IME,JMS:JME,KMS:KME) LOGICAL, INTENT(IN) :: extrapolate, ignore_lowest INTEGER :: i,j INTEGER :: k,kk REAL :: desired_press REAL :: dlnvaldlnp,dlnp DO K=1,generic DO j = JTS, min(JTE,JDE-1) DO i = ITS, min(ITE,IDE-1) DATA_IN(I,J,K)=max(DATA_in(I,J,K),1.e-13) ENDDO ENDDO ENDDO DO K=KMS,KME DO J=JMS,JME DO I=IMS,IME DATA_OUT(I,J,K)=-99999.9 ENDDO ENDDO ENDDO IF (ignore_lowest) then LMIN=2 ELSE LMIN=1 ENDIF DO j = JTS, min(JTE,JDE-1) test_i: DO i = ITS, min(ITE,IDE-1) IF (internal_time .gt. 1) THEN IF (J .ne. JDS .and. J .ne. JDE-1 .and. & I .ne. IDS .and. I .ne. IDE-1 ) THEN !! not on boundary CYCLE test_i ENDIF ENDIF output_loop: DO k = KDS,KDE-1 desired_press = press_out(i,j,k) IF (desired_press .GT. press_in(i,j,LMIN)) THEN IF ((desired_press - press_in(i,j,LMIN)).LT. 10.) THEN ! 0.1 mb data_out(i,j,k) = data_in(i,j,LMIN) ELSE IF (extrapolate) THEN ! Extrapolate downward because desired P level is below ! the lowest level in our input data. Extrapolate using simple ! 1st derivative of value with respect to ln P for the bottom 2 ! input layers. ! Add a check to make sure we are not using the gradient of ! a very thin layer IF ( (press_in(i,j,LMIN)-press_in(i,j,LMIN+1)) .GT. 100.) THEN dlnp = log(press_in(i,j,LMIN))-log(press_in(i,j,LMIN+1)) dlnvaldlnp = ( log(data_in(i,j,LMIN)) - log(data_in(i,j,LMIN+1)) ) / dlnp ELSE dlnp = log(press_in(i,j,LMIN))-log(press_in(i,j,LMIN+2)) dlnvaldlnp = (log(data_in(i,j,LMIN)) - log(data_in(i,j,LMIN+2))) / dlnp ENDIF data_out(i,j,k) = exp(log(data_in(i,j,LMIN)) + dlnvaldlnp * & ( log(desired_press)-log(press_in(i,j,LMIN)))) ELSE data_out(i,j,k) = data_in(i,j,LMIN) ENDIF ENDIF ELSE IF (desired_press .LT. press_in(i,j,generic)) THEN IF ( (press_in(i,j,generic) - desired_press) .LT. 10.) THEN data_out(i,j,k) = data_in(i,j,generic) ELSE IF (extrapolate) THEN ! Extrapolate upward IF ((press_in(i,j,generic-1)-press_in(i,j,generic)).GT.50.) THEN dlnp =log(press_in(i,j,generic))-log(press_in(i,j,generic-1)) dlnvaldlnp=(log(data_in(i,j,generic))-log(data_in(i,j,generic-1)))/dlnp ELSE dlnp =log(press_in(i,j,generic))-log(press_in(i,j,generic-2)) dlnvaldlnp=(log(data_in(i,j,generic))-log(data_in(i,j,generic-2)))/dlnp ENDIF data_out(i,j,k) = exp(log(data_in(i,j,generic)) + & dlnvaldlnp * (log(desired_press)-log(press_in(i,j,generic)))) ELSE data_out(i,j,k) = data_in(i,j,generic) ENDIF ENDIF ELSE ! We can trap between two levels and linearly interpolate input_loop: DO kk = LMIN, generic-1 IF (desired_press .EQ. press_in(i,j,kk) )THEN data_out(i,j,k) = data_in(i,j,kk) EXIT input_loop ELSE IF ( (desired_press .LT. press_in(i,j,kk)) .AND. & (desired_press .GT. press_in(i,j,kk+1)) ) THEN ! do trapped in lnp dlnp = log(press_in(i,j,kk)) - log(press_in(i,j,kk+1)) dlnvaldlnp = (log(data_in(i,j,kk))-log(data_in(i,j,kk+1)))/dlnp data_out(i,j,k) = exp(log(data_in(i,j,kk+1))+ & dlnvaldlnp*(log(desired_press)-log(press_in(i,j,kk+1)))) EXIT input_loop ENDIF ENDDO input_loop ENDIF ENDDO output_loop ENDDO test_i ENDDO END SUBROUTINE interp_press2press_log !------------------------------------------------------------------- SUBROUTINE rh_to_mxrat (rh, t, p, q , wrt_liquid , & ids , ide , jds , jde , kds , kde , & ims , ime , jms , jme , kms , kme , & its , ite , jts , jte , kts , kte ) IMPLICIT NONE INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , & ims , ime , jms , jme , kms , kme , & its , ite , jts , jte , kts , kte LOGICAL , INTENT(IN) :: wrt_liquid ! REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(IN) :: p , t ! REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(INOUT) :: rh REAL , DIMENSION(ims:ime,jms:jme,kms:kme) , INTENT(IN) :: p , t REAL , DIMENSION(ims:ime,jms:jme,kms:kme) , INTENT(INOUT) :: rh REAL , DIMENSION(ims:ime,jms:jme,kms:kme) , INTENT(OUT) :: q ! Local vars INTEGER :: i , j , k REAL :: ew , q1 , t1 REAL, PARAMETER :: T_REF = 0.0 REAL, PARAMETER :: MW_AIR = 28.966 REAL, PARAMETER :: MW_VAP = 18.0152 REAL, PARAMETER :: A0 = 6.107799961 REAL, PARAMETER :: A1 = 4.436518521e-01 REAL, PARAMETER :: A2 = 1.428945805e-02 REAL, PARAMETER :: A3 = 2.650648471e-04 REAL, PARAMETER :: A4 = 3.031240396e-06 REAL, PARAMETER :: A5 = 2.034080948e-08 REAL, PARAMETER :: A6 = 6.136820929e-11 REAL, PARAMETER :: ES0 = 6.1121 REAL, PARAMETER :: C1 = 9.09718 REAL, PARAMETER :: C2 = 3.56654 REAL, PARAMETER :: C3 = 0.876793 REAL, PARAMETER :: EIS = 6.1071 REAL :: RHS REAL, PARAMETER :: TF = 273.16 REAL :: TK REAL :: ES REAL :: QS REAL, PARAMETER :: EPS = 0.622 REAL, PARAMETER :: SVP1 = 0.6112 REAL, PARAMETER :: SVP2 = 17.67 REAL, PARAMETER :: SVP3 = 29.65 REAL, PARAMETER :: SVPT0 = 273.15 ! This subroutine computes mixing ratio (q, kg/kg) from basic variables ! pressure (p, Pa), temperature (t, K) and relative humidity (rh, 1-100%). ! The reference temperature (t_ref, C) is used to describe the temperature ! at which the liquid and ice phase change occurs. DO k = kts , kte DO j = jts , MIN ( jde-1 , jte ) DO i = its , MIN (ide-1 , ite ) rh(i,j,k) = MIN ( MAX ( rh(i,j,k) , 1. ) , 100. ) END DO END DO END DO IF ( wrt_liquid ) THEN DO k = kts , kte DO j = jts , MIN ( jde-1 , jte ) DO i = its , MIN (ide-1 , ite ) es=svp1*10.*EXP(svp2*(t(i,j,k)-svpt0)/(t(i,j,k)-svp3)) qs=eps*es/(p(i,j,k)/100.-es) q(i,j,k)=MAX(.01*rh(i,j,k)*qs,0.0) END DO END DO END DO ELSE DO k = kts , kte DO j = jts , MIN ( jde-1 , jte ) DO i = its , MIN (ide-1 , ite ) t1 = t(i,j,k) - 273.16 ! Obviously dry. IF ( t1 .lt. -200. ) THEN q(i,j,k) = 0 ELSE ! First compute the ambient vapor pressure of water IF ( ( t1 .GE. t_ref ) .AND. ( t1 .GE. -47.) ) THEN ! liq phase ESLO ew = a0 + t1 * (a1 + t1 * (a2 + t1 * (a3 + t1 * (a4 + t1 * (a5 + t1 * a6))))) ELSE IF ( ( t1 .GE. t_ref ) .AND. ( t1 .LT. -47. ) ) then !liq phas poor ES ew = es0 * exp(17.67 * t1 / ( t1 + 243.5)) ELSE tk = t(i,j,k) rhs = -c1 * (tf / tk - 1.) - c2 * alog10(tf / tk) + & c3 * (1. - tk / tf) + alog10(eis) ew = 10. ** rhs END IF ! Now sat vap pres obtained compute local vapor pressure ew = MAX ( ew , 0. ) * rh(i,j,k) * 0.01 ! Now compute the specific humidity using the partial vapor ! pressures of water vapor (ew) and dry air (p-ew). The ! constants assume that the pressure is in hPa, so we divide ! the pressures by 100. q1 = mw_vap * ew q1 = q1 / (q1 + mw_air * (p(i,j,k)/100. - ew)) q(i,j,k) = q1 / (1. - q1 ) END IF END DO END DO END DO END IF END SUBROUTINE rh_to_mxrat !--=------------------------------------------------------------------ SUBROUTINE boundary_smooth(h, landmask, grid, nsmth , nrow & &, IDS,IDE,JDS,JDE,KDS,KDE & &, IMS,IME,JMS,JME,KMS,KME & &, ITS,ITE,JTS,JTE,KTS,KTE ) implicit none TYPE (domain) :: grid integer :: IDS,IDE,JDS,JDE,KDS,KDE integer :: IMS,IME,JMS,JME,KMS,KME integer :: ITS,ITE,JTS,JTE,KTS,KTE integer:: ihw(JDS:JDE-1),ihe(JDS:JDE-1),nsmth,nrow real:: h(IMS:IME,JMS:JME),landmask(IMS:IME,JMS:JME) real :: h_old(IMS:IME,JMS:JME) real :: hbms(IDS:IDE-1,JDS:JDE-1) real :: hse(IDS:IDE-1,JDS:JDE-1) real :: hne(IDS:IDE-1,JDS:JDE-1) integer :: IPS,IPE,JPS,JPE,KPS,KPE integer :: ihl, ihh, m2l, ibas,jmelin integer :: I,J,KS,IOFFSET,JSTART,JEND character (len=255) :: message ips=its ipe=ite jps=jts jpe=jte kps=kts kpe=kte do j= JTS,min(JTE,JDE-1) ihw(J)=-mod(J,2) ihe(j)=ihw(J)+1 end do do J=JTS,min(JTE,JDE-1) do I=ITS,min(ITE,IDE-1) hbms(I,J)=landmask(I,J) enddo enddo jmelin=(JDE-1)-nrow+1 ibas=nrow/2 m2l=mod(nrow,2) do j=jts,min(jte,jde-1) ihl=ibas+mod(j,2)+m2l*mod(J+1,2) ihh=(IDE-1)-ibas-m2l*mod(J+1,2) do i=its,min(ite,ide-1) if (I .ge. ihl .and. I .le. ihh .and. J .ge. nrow .and. J .le. jmelin) then hbms(I,J)=0. endif end do end do 634 format(30(f2.0,1x)) do KS=1,nsmth grid%ht_gc=h #ifdef DM_PARALLEL # include "HALO_NMM_MG.inc" #endif h=grid%ht_gc h_old=grid%ht_gc do J=JTS,min(JTE,JDE-1) do I=ITS, min(ITE,IDE-1) if (I .ge. (IDS+mod(J,2)) .and. J .gt. JDS .and. J .lt. JDE-1 .and. I .lt. IDE-1) then h(i,j)= ( h_old(i+ihe(j),j+1) + h_old(i+ihw(j),j-1) + h_old(i+ihe(j),j-1) + h_old(i+ihw(j),j+1) - & 4. *h_old(i,j) )*hbms(i,j)*0.125+h_old(i,j) endif enddo enddo ! special treatment for four corners if (hbms(1,1) .eq. 1 .and. ITS .le. 1 .and. JTS .le. 1) then h(1,1)=0.75*h(1,1)+0.125*h(1+ihe(1),2)+ & 0.0625*(h(2,1)+h(1,3)) endif if (hbms(IDE-1,1) .eq. 1 .and. ITE .ge. IDE-2 .and. JTS .le. 1) then h(IDE-1,1)=0.75*h(IDE-1,1)+0.125*h(IDE-1+ihw(1),2)+ & 0.0625*(h(IDE-1-1,1)+h(IDE-1,3)) endif if (hbms(1,JDE-1) .eq. 1 .and. ITS .le. 1 .and. JTE .ge. JDE-2) then h(1,JDE-1)=0.75*h(1,JDE-1)+0.125*h(1+ihe(JDE-1),JDE-1-1)+ & 0.0625*(h(2,JDE-1)+h(1,JDE-1-2)) endif if (hbms(IDE-1,JDE-1) .eq. 1 .and. ITE .ge. IDE-2 .and. JTE .ge. JDE-2) then h(IDE-1,JDE-1)=0.75*h(IDE-1,JDE-1)+0.125*h(IDE-1+ihw(JDE-1),JDE-1-1)+ & 0.0625*(h(IDE-1-1,JDE-1)+h(IDE-1,JDE-1-2)) endif do J=JMS,JME do I=IMS,IME grid%ht_gc(I,J)=h(I,J) enddo enddo #ifdef DM_PARALLEL # include "HALO_NMM_MG.inc" #endif do J=JMS,JME do I=IMS,IME h(I,J)=grid%ht_gc(I,J) enddo enddo ! S bound if (JTS .eq. JDS) then J=JTS do I=ITS,ITE if (I .ge. IDS+1 .and. I .le. IDE-2) then if (hbms(I,J) .eq. 1) then h(I,J)=0.75*h(I,J)+0.125*(h(I+ihw(J),J+1)+h(I+ihe(J),J+1)) endif endif enddo endif ! N bound if (JTE .eq. JDE) then J=JDE-1 write(message,*) 'DOING N BOUND SMOOTHING for J= ', J CALL wrf_debug(100,message) do I=ITS,min(ITE,IDE-1) if (hbms(I,J) .eq. 1 .and. I .ge. IDS+1 .and. I .le. IDE-2) then h(I,J)=0.75*h(I,J)+0.125*(h(I+ihw(J),J-1)+h(I+ihe(J),J-1)) endif enddo endif ! W bound if (ITS .eq. IDS) then I=ITS do J=JTS,min(JTE,JDE-1) if (hbms(I,J) .eq. 1 .and. J .ge. JDS+2 .and. J .le. JDE-3 .and. mod(J,2) .eq. 1) then h(I,J)=0.75*h(I,J)+0.125*(h(I+ihe(J),J+1)+h(I+ihe(J),J-1)) endif enddo endif ! E bound if (ITE .eq. IDE) then write(message,*) 'DOING E BOUND SMOOTHING for I= ', min(ITE,IDE-1) CALL wrf_debug(100,message) I=min(ITE,IDE-1) do J=JTS,min(JTE,JDE-1) if (hbms(I,J) .eq. 1 .and. J .ge. JDS+2 .and. J .le. JDE-3 .and. mod(J,2) .eq. 1) then h(I,J)=0.75*h(I,J)+0.125*(h(I+ihw(J),J+1)+h(I+ihw(J),J-1)) endif enddo endif enddo ! end ks loop do J=JMS,JME do I=IMS,IME grid%ht_gc(I,J)=h(I,J) enddo enddo #ifdef DM_PARALLEL #include "HALO_NMM_MG.inc" #endif do J=JMS,JME do I=IMS,IME h(I,J)=grid%ht_gc(I,J) enddo enddo ! extra smoothing along inner boundary if (JTS .eq. JDS) then if (ITE .eq. IDE) then IOFFSET=1 else IOFFSET=0 endif ! Southern Boundary do i=its,min(ITE,IDE-1)-IOFFSET h(i,2)=0.25*(h(i,1)+h(i+1,1)+ & h(i,3)+h(i+1,3)) enddo endif if (JTE .eq. JDE) then if (ITE .eq. IDE) then IOFFSET=1 else IOFFSET=0 endif do i=its,min(ITE,IDE-1)-IOFFSET h(i,(JDE-1)-1)=0.25*(h(i,(JDE-1)-2)+h(i+1,(JDE-1)-2)+ & h(i,JDE-1)+h(i+1,JDE-1)) enddo endif if (JTS .eq. 1) then JSTART=4 else JSTART=JTS+mod(JTS,2) ! needs to be even endif if (JTE .eq. JDE) then JEND=(JDE-1)-3 else JEND=JTE endif if (ITS .eq. IDS) then ! Western Boundary do j=JSTART,JEND,2 h(1,j)=0.25*(h(1,j-1)+h(2,j-1)+ & h(1,j+1)+h(2,j+1)) enddo endif if (ITE .eq. IDE) then ! Eastern Boundary do j=JSTART,JEND,2 h((IDE-1)-1,j)=0.25*(h((IDE-1)-1,j-1)+h((IDE-1),j-1)+ & h((IDE-1)-1,j+1)+h((IDE-1),j+1)) enddo endif END SUBROUTINE boundary_smooth !-------------------------------------------------------------------- SUBROUTINE monthly_interp_to_date ( field_in , date_str , field_out , & ids , ide , jds , jde , kds , kde , & ims , ime , jms , jme , kms , kme , & its , ite , jts , jte , kts , kte ) ! Linrarly in time interpolate data to a current valid time. The data is ! assumed to come in "monthly", valid at the 15th of every month. IMPLICIT NONE INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , & ims , ime , jms , jme , kms , kme , & its , ite , jts , jte , kts , kte CHARACTER (LEN=24) , INTENT(IN) :: date_str REAL , DIMENSION(ims:ime,jms:jme,12) , INTENT(IN) :: field_in REAL , DIMENSION(ims:ime, jms:jme) , INTENT(OUT) :: field_out ! Local vars INTEGER :: i , j , l INTEGER , DIMENSION(0:13) :: middle INTEGER :: target_julyr , target_julday , target_date INTEGER :: julyr , julday , int_month, next_month REAL :: gmt CHARACTER (LEN=4) :: yr CHARACTER (LEN=2) :: mon , day15 WRITE(day15,FMT='(I2.2)') 15 DO l = 1 , 12 WRITE(mon,FMT='(I2.2)') l CALL get_julgmt ( date_str(1:4)//'-'//mon//'-'//day15//'_'//'00:00:00.0000' , julyr , julday , gmt ) middle(l) = julyr*1000 + julday END DO l = 0 middle(l) = middle( 1) - 31 l = 13 middle(l) = middle(12) + 31 CALL get_julgmt ( date_str , target_julyr , target_julday , gmt ) target_date = target_julyr * 1000 + target_julday find_month : DO l = 0 , 12 IF ( ( middle(l) .LT. target_date ) .AND. ( middle(l+1) .GE. target_date ) ) THEN DO j = jts , MIN ( jde-1 , jte ) DO i = its , MIN (ide-1 , ite ) int_month = MOD ( l , 12 ) IF ( int_month .EQ. 0 ) int_month = 12 IF (int_month == 12) THEN next_month=1 ELSE next_month=int_month+1 ENDIF field_out(i,j) = ( field_in(i,j,next_month) * ( target_date - middle(l) ) + & field_in(i,j,int_month ) * ( middle(l+1) - target_date ) ) / & ( middle(l+1) - middle(l) ) END DO END DO EXIT find_month END IF END DO find_month END SUBROUTINE monthly_interp_to_date !--------------------------------------------------------------------- SUBROUTINE monthly_min_max ( field_in , field_min , field_max , & ids , ide , jds , jde , kds , kde , & ims , ime , jms , jme , kms , kme , & its , ite , jts , jte , kts , kte ) ! Plow through each month, find the max, min values for each i,j. IMPLICIT NONE INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , & ims , ime , jms , jme , kms , kme , & its , ite , jts , jte , kts , kte REAL , DIMENSION(ims:ime,jms:jme,12) , INTENT(IN) :: field_in REAL , DIMENSION(ims:ime, jms:jme) , INTENT(OUT) :: field_min , field_max ! Local vars INTEGER :: i , j , l REAL :: minner , maxxer DO j = jts , MIN(jde-1,jte) DO i = its , MIN(ide-1,ite) minner = field_in(i,j,1) maxxer = field_in(i,j,1) DO l = 2 , 12 IF ( field_in(i,j,l) .LT. minner ) THEN minner = field_in(i,j,l) END IF IF ( field_in(i,j,l) .GT. maxxer ) THEN maxxer = field_in(i,j,l) END IF END DO field_min(i,j) = minner field_max(i,j) = maxxer END DO END DO END SUBROUTINE monthly_min_max !----------------------------------------------------------------------- SUBROUTINE reverse_vert_coord ( field, start_z, end_z & &, IDS,IDE,JDS,JDE,KDS,KDE & &, IMS,IME,JMS,JME,KMS,KME & &, ITS,ITE,JTS,JTE,KTS,KTE ) IMPLICIT NONE INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , & ims , ime , jms , jme , kms , kme , & its , ite , jts , jte , kts , kte, & start_z, end_z REAL, INTENT(INOUT) :: field(IMS:IME,JMS:JME,end_z) ! local INTEGER :: I,J,L REAL, ALLOCATABLE :: dum3d(:,:,:) allocate(dum3d(IMS:IME,JMS:JME,end_z)) DO L=start_z,end_z DO J=jts,min(jte,jde-1) DO I=its,min(ite,ide-1) dum3d(I,J,L)=field(I,J,end_z-L+start_z) END DO END DO END DO DO L=start_z,end_z DO J=jts,min(jte,jde-1) DO I=its,min(ite,ide-1) field(I,J,L)=dum3d(I,J,L) END DO END DO END DO DEALLOCATE(dum3d) END SUBROUTINE reverse_vert_coord !-------------------------------------------------------------------- SUBROUTINE compute_nmm_levels(ninterface, ptop, eta_levels) USE module_model_constants IMPLICIT NONE character(len=132):: message integer :: ninterface,Lthick,L real, parameter:: gamma=.0065 real, parameter:: t_stand=288. real, parameter:: p_stand=101325. real :: maxdz_compute, ptop real :: plower,pupper,tlay, sum real :: eta_levels(ninterface) real, allocatable:: Z(:) real, allocatable:: deta_levels_spline(:) logical:: print_pbl_warn !---------------------------------------------------- allocate(Z(ninterface)) allocate(deta_levels_spline(ninterface-1)) CALL compute_eta_spline(ninterface-1,deta_levels_spline,ptop) sum=0. DO L=1,ninterface-1 sum=sum+deta_levels_spline(L) ENDDO eta_levels(1)=1.00 DO L=2,ninterface eta_levels(L)=eta_levels(L-1)-deta_levels_spline(L-1) ENDDO eta_levels(ninterface)=0.00 DO L=2,ninterface-1 eta_levels(L)=0.5*(eta_levels(L))+0.25*(eta_levels(L-1)+eta_levels(L+1)) ENDDO Z(1)=0. maxdz_compute=0. print_pbl_warn=.false. DO L=2,ninterface tlay=max( t_stand-gamma*Z(L-1), 216.5) plower=ptop+(p_stand-ptop)*eta_levels(L-1) pupper=ptop+(p_stand-ptop)*eta_levels(L) Z(L)=Z(L-1)+(tlay*r_d/g)*(log(plower)-log(pupper)) if (plower .gt. 85000. .and. pupper .lt. 85000. .and. L .lt. 10) then print_pbl_warn=.true. endif write(message,*) 'L, eta(l), pupper, Z(L): ', L, eta_levels(L),pupper,Z(L) CALL wrf_debug(100,message) if (Z(L)-Z(L-1) .gt. maxdz_compute) then Lthick=L endif maxdz_compute=max(maxdz_compute,Z(L)-Z(L-1)) ENDDO if (print_pbl_warn) then write(message,*) 'WARNING - PBL MAY BE POORLY RESOLVED WITH NUMBER OF VERTICAL LEVELS' CALL wrf_message(message) write(message,*) ' - CONSIDER INCREASING THE VERTICAL RESOLUTION' CALL wrf_message(message) endif write(message,*) 'thickest layer was: ', maxdz_compute , 'meters thick at level: ', Lthick CALL wrf_message(message) END SUBROUTINE compute_nmm_levels !--------------------------- SUBROUTINE compute_eta_spline(LM, dsg, ptop) IMPLICIT NONE real:: dsg(LM), ptop, sum, rsum real, allocatable:: xold(:),dold(:) real, allocatable:: xnew(:),sgm(:) real, allocatable:: pps(:),qqs(:),y2s(:) integer nlev,LM,L,KOLD IF (LM .ge. 46) THEN KOLD=9 allocate(xold(KOLD)) allocate(dold(KOLD)) xold(1)=.00 dold(1)=.006 xold(2)=.13 dold(2)=.009 xold(3)=.19 dold(3)=.012 xold(4)=.30 dold(4)=.036 xold(5)=.42 dold(5)=.041 xold(6)=.56 dold(6)=.040 xold(7)=.69 dold(7)=.018 if (ptop .ge. 2000.) then xold(8)=.90 dold(8)=.012 xold(9)=1.0 dold(9)=.006 else xold(8)=.90 dold(8)=.008 xold(9)=1.0 dold(9)=.003 endif ELSE KOLD=8 allocate(xold(KOLD)) allocate(dold(KOLD)) xold(1)=.00 dold(1)=.006 xold(2)=.18 dold(2)=.015 xold(3)=.32 dold(3)=.035 xold(4)=.50 dold(4)=.040 xold(5)=.68 dold(5)=.030 xold(6)=.75 dold(6)=.017 xold(7)=.85 dold(7)=.012 if (ptop .ge. 2000.) then xold(8)=1.0 dold(8)=.013 else xold(8)=1.0 dold(8)=.008 endif ENDIF allocate(xnew(lm)) allocate(sgm(lm+1)) allocate(pps(lm)) allocate(qqs(lm)) allocate(y2s(lm)) DO L=1,LM xnew(l)=float(l-1)/float(lm-1) ENDDO y2s=0. CALL spline(kold,xold,dold,y2s,lm,xnew,dsg,pps,qqs) sum=0. DO l=1,lm sum=sum+dsg(l) ENDDO rsum=1./sum sgm(1)=0. DO L=1,lm-1 dsg(l)=dsg(l)*rsum sgm(l+1)=sgm(l)+dsg(l) ENDDO sgm(lm+1)=1. dsg(lm)=sgm(lm+1)-sgm(lm) END SUBROUTINE compute_eta_spline ! ------------------------------------------------------------------- subroutine spline(NOLD,XOLD,YOLD,Y2,NNEW,XNEW,YNEW,P,q) ! ******************************************************************** ! * * ! * THIS IS A ONE-DIMENSIONAL CUBIC SPLINE FITTING ROUTINE * ! * PROGRAMED FOR A SMALL SCALAR MACHINE. * ! * * ! * PROGRAMER Z. JANJIC * ! * * ! * NOLD - NUMBER OF GIVEN VALUES OF THE FUNCTION. MUST BE GE 3. * ! * XOLD - LOCATIONS OF THE POINTS AT WHICH THE VALUES OF THE * ! * FUNCTION ARE GIVEN. MUST BE IN ASCENDING ORDER. * ! * YOLD - THE GIVEN VALUES OF THE FUNCTION AT THE POINTS XOLD. * ! * Y2 - THE SECOND DERIVATIVES AT THE POINTS XOLD. IF NATURAL * ! * SPLINE IS FITTED Y2(1)=0. AND Y2(NOLD)=0. MUST BE * ! * SPECIFIED. * ! * NNEW - NUMBER OF VALUES OF THE FUNCTION TO BE CALCULATED. * ! * XNEW - LOCATIONS OF THE POINTS AT WHICH THE VALUES OF THE * ! * FUNCTION ARE CALCULATED. XNEW(K) MUST BE GE XOLD(1) * ! * AND LE XOLD(NOLD). * ! * YNEW - THE VALUES OF THE FUNCTION TO BE CALCULATED. * ! * P, q - AUXILIARY VECTORS OF THE LENGTH NOLD-2. * ! * * ! ******************************************************************** ! ! LOG: ! ! JOVIC - July 2008 - fixed incorrectly dimensioned arrays, ! PYLE and do loop leading to out of bound array ! reference !------ ! ! PYLE - June 2007 - eliminated use of GO TO statements. ! !----------------------------------------------------------------------- IMPLICIT NONE !----------------------------------------------------------------------- INTEGER,INTENT(IN) :: NNEW,NOLD REAL,DIMENSION(NOLD),INTENT(IN) :: XOLD,YOLD REAL,DIMENSION(NNEW),INTENT(IN) :: XNEW REAL,DIMENSION(NNEW),INTENT(OUT) :: YNEW REAL,DIMENSION(NOLD+2),INTENT(INOUT) :: P,q,Y2 ! INTEGER :: K,K1,K2,KOLD,NOLDM1, K2_hold, K_hold REAL :: AK,BK,CK,DEN,DX,DXC,DXL,DXR,DYDXL,DYDXR & & ,RDX,RTDXC,X,XK,XSQ,Y2K,Y2KP1 !----------------------------------------------------------------------- NOLDM1=NOLD-1 DXL=XOLD(2)-XOLD(1) DXR=XOLD(3)-XOLD(2) DYDXL=(YOLD(2)-YOLD(1))/DXL DYDXR=(YOLD(3)-YOLD(2))/DXR RTDXC=0.5/(DXL+DXR) P(1)= RTDXC*(6.*(DYDXR-DYDXL)-DXL*Y2(1)) q(1)=-RTDXC*DXR K=3 first_loop: DO K=3,NOLD-1 DXL=DXR DYDXL=DYDXR DXR=XOLD(K+1)-XOLD(K) DYDXR=(YOLD(K+1)-YOLD(K))/DXR DXC=DXL+DXR DEN=1./(DXL*q(K-2)+DXC+DXC) P(K-1)= DEN*(6.*(DYDXR-DYDXL)-DXL*P(K-2)) q(K-1)=-DEN*DXR END DO first_loop DO K=NOLDM1,2,-1 Y2(K)=P(K-1)+q(K-1)*Y2(K+1) K_hold=K END DO K=K_hold !----------------------------------------------------------------------- second_loop: DO K1=1,NNEW XK=XNEW(K1) third_loop: DO K2=2,NOLD IF(XOLD(K2)>XK)THEN KOLD=K2-1 K2_hold=K2 exit third_loop ENDIF K2_hold=K2 END DO third_loop IF (XOLD(K2_hold) .le. XK) THEN YNEW(K1)=YOLD(NOLD) CYCLE second_loop ENDIF IF (K1 .eq. 1 .or. K .ne. KOLD) THEN K=KOLD Y2K=Y2(K) Y2KP1=Y2(K+1) DX=XOLD(K+1)-XOLD(K) RDX=1./DX AK=.1666667*RDX*(Y2KP1-Y2K) BK=0.5*Y2K CK=RDX*(YOLD(K+1)-YOLD(K))-.1666667*DX*(Y2KP1+Y2K+Y2K) ENDIF X=XK-XOLD(K) XSQ=X*X YNEW(K1)=AK*XSQ*X+BK*XSQ+CK*X+YOLD(K) END DO second_loop END SUBROUTINE SPLINE !-------------------------------------------------------------------- SUBROUTINE NMM_SH2O(IMS,IME,JMS,JME,ISTART,IM,JSTART,JM,& NSOIL,ISLTPK, & sm,sice,stc,smc,sh2o) !! INTEGER, PARAMETER:: NSOTYP=9 ! INTEGER, PARAMETER:: NSOTYP=16 INTEGER, PARAMETER:: NSOTYP=19 !!!!!!!!MAYBE??? REAL :: PSIS(NSOTYP),BETA(NSOTYP),SMCMAX(NSOTYP) REAL :: stc(IMS:IME,NSOIL,JMS:JME), & smc(IMS:IME,NSOIL,JMS:JME) REAL :: sh2o(IMS:IME,NSOIL,JMS:JME),sice(IMS:IME,JMS:JME),& sm(IMS:IME,JMS:JME) REAL :: HLICE,GRAV,T0,BLIM INTEGER :: ISLTPK(IMS:IME,JMS:JME) CHARACTER(LEN=255) :: message ! Constants used in cold start sh2o initialization DATA HLICE/3.335E5/,GRAV/9.81/,T0/273.15/ DATA BLIM/5.5/ ! DATA PSIS /0.04,0.62,0.47,0.14,0.10,0.26,0.14,0.36,0.04/ ! DATA BETA /4.26,8.72,11.55,4.74,10.73,8.17,6.77,5.25,4.26/ ! DATA SMCMAX /0.421,0.464,0.468,0.434,0.406, & ! 0.465,0.404,0.439,0.421/ !!! NOT SURE...PSIS=SATPSI, BETA=BB?? DATA PSIS /0.069, 0.036, 0.141, 0.759, 0.759, 0.355, & 0.135, 0.617, 0.263, 0.098, 0.324, 0.468, & 0.355, 0.000, 0.069, 0.036, 0.468, 0.069, 0.069 / DATA BETA/2.79, 4.26, 4.74, 5.33, 5.33, 5.25, & 6.66, 8.72, 8.17, 10.73, 10.39, 11.55, & 5.25, 0.00, 2.79, 4.26, 11.55, 2.79, 2.79 / DATA SMCMAX/0.339, 0.421, 0.434, 0.476, 0.476, 0.439, & 0.404, 0.464, 0.465, 0.406, 0.468, 0.468, & 0.439, 1.000, 0.200, 0.421, 0.468, 0.200, 0.339/ DO K=1,NSOIL DO J=JSTART,JM DO I=ISTART,IM !tst IF (smc(I,K,J) .gt. SMCMAX(ISLTPK(I,J))) then if (K .eq. 1) then write(message,*) 'I,J,reducing smc from ' ,I,J,smc(I,K,J), 'to ', SMCMAX(ISLTPK(I,J)) CALL wrf_debug(100,message) endif smc(I,K,J)=SMCMAX(ISLTPK(I,J)) ENDIF !tst IF ( (sm(I,J) .lt. 0.5) .and. (sice(I,J) .lt. 0.5) ) THEN IF (ISLTPK(I,J) .gt. 19) THEN WRITE(message,*) 'FORCING ISLTPK at : ', I,J CALL wrf_message(message) ISLTPK(I,J)=9 ELSEIF (ISLTPK(I,J) .le. 0) then WRITE(message,*) 'FORCING ISLTPK at : ', I,J CALL wrf_message(message) ISLTPK(I,J)=1 ENDIF ! cold start: determine liquid soil water content (sh2o) ! sh2o <= smc for t < 273.149K (-0.001C) IF (stc(I,K,J) .LT. 273.149) THEN ! first guess following explicit solution for Flerchinger Eqn from Koren ! et al, JGR, 1999, Eqn 17 (KCOUNT=0 in FUNCTION FRH2O). BX = BETA(ISLTPK(I,J)) IF ( BETA(ISLTPK(I,J)) .GT. BLIM ) BX = BLIM if ( GRAV*(-PSIS(ISLTPK(I,J))) .eq. 0 ) then write(message,*) 'TROUBLE' CALL wrf_message(message) write(message,*) 'I,J: ', i,J CALL wrf_message(message) write(message,*) 'grav, isltpk, psis(isltpk): ', grav,isltpk(I,J),& psis(isltpk(I,J)) CALL wrf_message(message) endif if (BX .eq. 0 .or. stc(I,K,J) .eq. 0) then write(message,*) 'TROUBLE -- I,J,BX, stc: ', I,J,BX,stc(I,K,J) CALL wrf_message(message) endif FK = (((HLICE/(GRAV*(-PSIS(ISLTPK(I,J)))))* & ((stc(I,K,J)-T0)/stc(I,K,J)))** & (-1/BX))*SMCMAX(ISLTPK(I,J)) IF (FK .LT. 0.02) FK = 0.02 sh2o(I,K,J) = MIN ( FK, smc(I,K,J) ) ! ---------------------------------------------------------------------- ! now use iterative solution for liquid soil water content using ! FUNCTION FRH2O (from the Eta "NOAH" land-surface model) with the ! initial guess for sh2o from above explicit first guess. sh2o(I,K,J)=FRH2O_init(stc(I,K,J),smc(I,K,J),sh2o(I,K,J), & SMCMAX(ISLTPK(I,J)),BETA(ISLTPK(I,J)), & PSIS(ISLTPK(I,J))) ELSE ! above freezing sh2o(I,K,J)=smc(I,K,J) ENDIF ELSE ! water point sh2o(I,K,J)=smc(I,K,J) ENDIF ! test on land/ice/sea if (sh2o(I,K,J) .gt. SMCMAX(ISLTPK(I,J))) then write(message,*) 'sh2o > THAN SMCMAX ', I,J,sh2o(I,K,J),SMCMAX(ISLTPK(I,J)),smc(I,K,J) CALL wrf_message(message) endif ENDDO ENDDO ENDDO END SUBROUTINE NMM_SH2O !------------------------------------------------------------------- FUNCTION FRH2O_init(TKELV,smc,sh2o,SMCMAX,B,PSIS) IMPLICIT NONE ! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC ! PURPOSE: CALCULATE AMOUNT OF SUPERCOOLED LIQUID SOIL WATER CONTENT ! IF TEMPERATURE IS BELOW 273.15K (T0). REQUIRES NEWTON-TYPE ITERATION ! TO SOLVE THE NONLINEAR IMPLICIT EQUATION GIVEN IN EQN 17 OF ! KOREN ET AL. (1999, JGR, VOL 104(D16), 19569-19585). ! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC ! ! New version (JUNE 2001): much faster and more accurate newton iteration ! achieved by first taking log of eqn cited above -- less than 4 ! (typically 1 or 2) iterations achieves convergence. Also, explicit ! 1-step solution option for special case of parameter Ck=0, which reduces ! the original implicit equation to a simpler explicit form, known as the ! ""Flerchinger Eqn". Improved handling of solution in the limit of ! freezing point temperature T0. ! ! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC ! ! INPUT: ! ! TKELV.........Temperature (Kelvin) ! smc...........Total soil moisture content (volumetric) ! sh2o..........Liquid soil moisture content (volumetric) ! SMCMAX........Saturation soil moisture content (from REDPRM) ! B.............Soil type "B" parameter (from REDPRM) ! PSIS..........Saturated soil matric potential (from REDPRM) ! ! OUTPUT: ! FRH2O.........supercooled liquid water content. ! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC REAL B REAL BLIM REAL BX REAL CK REAL DENOM REAL DF REAL DH2O REAL DICE REAL DSWL REAL ERROR REAL FK REAL FRH2O_init REAL GS REAL HLICE REAL PSIS REAL sh2o REAL smc REAL SMCMAX REAL SWL REAL SWLK REAL TKELV REAL T0 INTEGER NLOG INTEGER KCOUNT PARAMETER (CK=8.0) ! PARAMETER (CK=0.0) PARAMETER (BLIM=5.5) ! PARAMETER (BLIM=7.0) PARAMETER (ERROR=0.005) PARAMETER (HLICE=3.335E5) PARAMETER (GS = 9.81) PARAMETER (DICE=920.0) PARAMETER (DH2O=1000.0) PARAMETER (T0=273.15) ! ### LIMITS ON PARAMETER B: B < 5.5 (use parameter BLIM) #### ! ### SIMULATIONS SHOWED IF B > 5.5 UNFROZEN WATER CONTENT #### ! ### IS NON-REALISTICALLY HIGH AT VERY LOW TEMPERATURES #### ! ################################################################ ! BX = B IF ( B .GT. BLIM ) BX = BLIM ! ------------------------------------------------------------------ ! INITIALIZING ITERATIONS COUNTER AND ITERATIVE SOLUTION FLAG. NLOG=0 KCOUNT=0 ! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC ! C IF TEMPERATURE NOT SIGNIFICANTLY BELOW FREEZING (T0), sh2o = smc ! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC IF (TKELV .GT. (T0 - 1.E-3)) THEN FRH2O_init=smc ELSE ! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC IF (CK .NE. 0.0) THEN ! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC ! CCCCCCCCC OPTION 1: ITERATED SOLUTION FOR NONZERO CK CCCCCCCCCCC ! CCCCCCCCCCCC IN KOREN ET AL, JGR, 1999, EQN 17 CCCCCCCCCCCCCCCCC ! INITIAL GUESS FOR SWL (frozen content) SWL = smc-sh2o ! KEEP WITHIN BOUNDS. IF (SWL .GT. (smc-0.02)) SWL=smc-0.02 IF(SWL .LT. 0.) SWL=0. ! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC ! C START OF ITERATIONS ! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DO WHILE (NLOG .LT. 10 .AND. KCOUNT .EQ. 0) NLOG = NLOG+1 DF = ALOG(( PSIS*GS/HLICE ) * ( ( 1.+CK*SWL )**2. ) * & ( SMCMAX/(smc-SWL) )**BX) - ALOG(-(TKELV-T0)/TKELV) DENOM = 2. * CK / ( 1.+CK*SWL ) + BX / ( smc - SWL ) SWLK = SWL - DF/DENOM ! BOUNDS USEFUL FOR MATHEMATICAL SOLUTION. IF (SWLK .GT. (smc-0.02)) SWLK = smc - 0.02 IF(SWLK .LT. 0.) SWLK = 0. ! MATHEMATICAL SOLUTION BOUNDS APPLIED. DSWL=ABS(SWLK-SWL) SWL=SWLK ! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC ! CC IF MORE THAN 10 ITERATIONS, USE EXPLICIT METHOD (CK=0 APPROX.) ! CC WHEN DSWL LESS OR EQ. ERROR, NO MORE ITERATIONS REQUIRED. ! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC IF ( DSWL .LE. ERROR ) THEN KCOUNT=KCOUNT+1 END IF END DO ! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC ! C END OF ITERATIONS ! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC ! BOUNDS APPLIED WITHIN DO-BLOCK ARE VALID FOR PHYSICAL SOLUTION. FRH2O_init = smc - SWL ! CCCCCCCCCCCCCCCCCCCCCCCC END OPTION 1 CCCCCCCCCCCCCCCCCCCCCCCCCCC ENDIF IF (KCOUNT .EQ. 0) THEN ! Print*,'Flerchinger used in NEW version. Iterations=',NLOG ! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC ! CCCCC OPTION 2: EXPLICIT SOLUTION FOR FLERCHINGER EQ. i.e. CK=0 CCCCCCCC ! CCCCCCCCCCCCC IN KOREN ET AL., JGR, 1999, EQN 17 CCCCCCCCCCCCCCC FK=(((HLICE/(GS*(-PSIS)))*((TKELV-T0)/TKELV))**(-1/BX))*SMCMAX ! APPLY PHYSICAL BOUNDS TO FLERCHINGER SOLUTION IF (FK .LT. 0.02) FK = 0.02 FRH2O_init = MIN ( FK, smc ) ! CCCCCCCCCCCCCCCCCCCCCCCCC END OPTION 2 CCCCCCCCCCCCCCCCCCCCCCCCCC ENDIF ENDIF RETURN END FUNCTION FRH2O_init !-------------------------------------------------------------------- SUBROUTINE init_module_initialize END SUBROUTINE init_module_initialize !--------------------------------------------------------------------- !------------------------------------------------------------------------------- SUBROUTINE vortex ( ght_gc,rh_gc,t_gc,u_gc,v_gc,p_gc & &, ght_out,rh_out,t_out,u_out,v_out & &, ht_gc,tsk_gc,xice_gc & &, hlat_gc,hlon_gc,vlat_gc,vlon_gc & &, greenfrac_gc,albedo12m_gc,landusef_gc & &, soilctop_gc,soilcbot_gc & &, landusef_out,soilctop_out,soilcbot_out & &, num_veg_cat,num_soil_top_cat,num_soil_bot_cat & &, dx,internal_time_loop & &, start_z,end_z,sf_surface_physics & &, ids,ide,jds,jde,kds,kde & &, ims,ime,jms,jme,kms,kme & &, its,ite,jts,jte,kts,kte ) USE module_dm IMPLICIT NONE LOGICAL,EXTERNAL :: WRF_DM_ON_MONITOR INTEGER, INTENT(IN):: ids,ide,jds,jde,kds,kde, & ims,ime,jms,jme,kms,kme, & its,ite,jts,jte,kts,kte, & internal_time_loop,sf_surface_physics INTEGER, INTENT(IN):: num_veg_cat,num_soil_top_cat,num_soil_bot_cat,start_z, end_z REAL, DIMENSION(IMS:IME,JMS:JME,start_z:end_z),INTENT(INOUT) :: ght_gc,rh_gc,t_gc,u_gc,v_gc,p_gc REAL, DIMENSION(IMS:IME,JMS:JME,start_z:end_z),INTENT(INOUT) :: ght_out,rh_out,t_out,u_out,v_out REAL, DIMENSION(IMS:IME,JMS:JME),INTENT(INOUT) :: ht_gc,tsk_gc,xice_gc REAL, DIMENSION(IMS:IME,JMS:JME),INTENT(INOUT) :: hlat_gc,hlon_gc,vlat_gc,vlon_gc REAL, DIMENSION(ims:ime,jms:jme,12) , INTENT(INOUT) :: greenfrac_gc,albedo12m_gc REAL, DIMENSION(ims:ime,jms:jme,num_veg_cat), INTENT(INOUT) :: landusef_gc,landusef_out REAL, DIMENSION(ims:ime,jms:jme,num_soil_top_cat),INTENT(INOUT):: soilctop_gc,soilctop_out REAL, DIMENSION(ims:ime,jms:jme,num_soil_bot_cat),INTENT(INOUT):: soilcbot_gc,soilcbot_out REAL, INTENT(INOUT):: dx REAL, PARAMETER :: eradius=6371221.3 INTEGER :: i,j,l,iter,nmax,nres,nfilter INTEGER :: id0,jd0,id1,jd1 REAL :: factor,glon0,glat0,dist,cavlat,minval_gh REAL :: beta1,vmax,rmax,dx_km REAL, DIMENSION(start_z:end_z) :: vmax1,rmax1 REAL, DIMENSION(IMS:IME,JMS:JME) :: dist1,dir1,diffx,diffy REAL, DIMENSION(IMS:IME,JMS:JME,start_z:end_z) :: pert_ght,pert_rh,pert_t,pert_u,pert_v,pert_p,temp0 !New Vortex Related DATA beta1 /0.0/ ! 0.0 sets non-linear balance equation vortex. GFDL vortex for beta1 between 0.1 and 1.5 character*255 :: message REAL, ALLOCATABLE,DIMENSION(:,:) :: psc REAL, ALLOCATABLE,DIMENSION(:,:,:) :: u_temp,v_temp,t_temp,ght_temp,rh_temp !!BEGIN: LSM changes for LANDFALL: Subashini 7/27/2016 INTEGER :: imin,jmin,imax,jmax,proceed,VEG_ID,SOIL_ID LOGICAL :: mvland, logic_temp REAL :: s_temp NAMELIST / init_land /imin,jmin,imax,jmax,proceed,VEG_ID,SOIL_ID,mvland,logic_temp,s_temp !! END: LSM changes for LANDFALL : Subashini 7/27/2016 !---------------------------------------------------------------------------- ! PURPOSE: ! - HURRICANE VORTEX FILTERING ! - HURRICANE VORTEX INITIALIZATION ! ! This is gopal's doing !---------------------------------------------------------------------------- WRITE(0,*)'---------------- IDEAL CASE -------------------------' call wrf_debug(1,message) WRITE(message,*)'number of pressure levels',end_z call wrf_debug(1,message) WRITE(message,*)'number of vegcat levels',num_veg_cat call wrf_debug(1,message) WRITE(message,*)'number of soiltop levels',num_soil_top_cat call wrf_debug(1,message) WRITE(message,*)'number of soilbot levels',num_soil_bot_cat call wrf_debug(1,message) WRITE(message,*)'--------------- ght_gc --------------------------' ! call wrf_debug(1,message) ! WRITE(message,*)ght_gc(100,100,:) ! call wrf_debug(1,message) WRITE(message,*)'--------------- p_gc --------------------------' ! call wrf_debug(1,message) ! WRITE(message,*)p_gc(100,100,:) ! call wrf_debug(1,message) WRITE(message,*)'--------------- rh_gc ---------------------------' ! call wrf_debug(1,message) ! WRITE(message,*)rh_gc(100,100,:) ! call wrf_debug(1,message) WRITE(message,*)'--------------- t_gc --------------------------' ! call wrf_debug(1,message) ! WRITE(message,*)t_gc(100,100,:) ! call wrf_debug(1,message) WRITE(message,*)'------------------------------------------------' ! call wrf_debug(1,message) ! WRITE(message,*)'------------------------------------------------' ! call wrf_debug(1,message) ! ! SET UP IDEAL CONDITIONS ! !!BEGIN: LSM changes for LANDFALL: Subashini 7/27/2016 open (7,FILE='land.nml') read (UNIT=7, NML=init_land) close (UNIT=7) if (mvland .and. sf_surface_physics .ne. GFDLSLAB) then CALL wrf_error_fatal('LSM must be GFDLSLAB when mvland is true') endif do l = 1,num_veg_cat do j = jts, MIN(jte,jde-1) do i = its, MIN(ite,ide-1) ! create land patch that will move from W2E if(mvland .and. (i .ge. imin .and. i .le. imax .and. j .ge. jmin .and. j.le. jmax))then if(l==VEG_ID)THEN landusef_out(i,j,l)=1 ! barren land else landusef_out(i,j,l)=0 endif else ! original ocean world if(l==16)THEN landusef_out(i,j,l)=1 ! create ocean elsewhere else landusef_out(i,j,l)=0 endif endif ! enddo enddo enddo do l = 1,num_soil_top_cat do j = jts, MIN(jte,jde-1) do i = its, MIN(ite,ide-1) ! create land patch that will move from W2E if(mvland .and. (i .ge. imin .and. i .le. imax .and. j .ge. jmin .and. j.le. jmax))then if(l==SOIL_ID)THEN soilctop_out(i,j,l)=1 ! sandy soil else soilctop_out(i,j,l)=0 endif else ! original ocean world if(l==14)THEN soilctop_out(i,j,l)=1 ! create ocean everywhere else soilctop_out(i,j,l)=0 endif endif ! enddo enddo enddo do l = 1,num_soil_bot_cat do j = jts, MIN(jte,jde-1) do i = its, MIN(ite,ide-1) ! create land patch that will move from W2E if(mvland .and. (i .ge. imin .and. i .le. imax .and. j .ge. jmin .and. j.le. jmax))then if(l==SOIL_ID)THEN soilcbot_out(i,j,l)=1 ! sandy soil else soilcbot_out(i,j,l)=0 endif else ! original ocean world if(l==14)THEN soilcbot_out(i,j,l)=1 ! create ocean everywhere else soilcbot_out(i,j,l)=0 endif endif ! enddo enddo enddo landusef_gc=landusef_out !=landusef_gc soilcbot_gc=soilcbot_out !=soilcbot_gc soilctop_gc=soilctop_out !=soilctop_gc do j = jts, MIN(jte,jde-1) do i = its, MIN(ite,ide-1) xice_gc(i,j)=0. ht_gc(i,j)=0. ! uniform terrain ght_gc(i,j,1)=0. ! uniform gpm at level 1 ! create land patch that will move from W2E if(mvland .and. (i .ge. imin .and. i .le. imax .and. j .ge. jmin .and.j.le. jmax))then if(logic_temp .eq. .true.) then tsk_gc(i,j)=s_temp ! uniform land temperature or t_gc(i,j,1) for applying first level temperature else tsk_gc(i,j)= t_gc(i,j,1) endif else tsk_gc(i,j)=302.0 ! uniform SSTs endif enddo enddo !! END: LSM changes for LANDFALL : Subashini 7/27/2016 ! ! Make sure the GFDL analysis does not have temperature problems especially ! near the surface. I noticed some problems with the wps outputs (met_nmm file.) ! do j=jts,min(jte,jde-1) do i=its,min(ite,ide-1) if(t_gc(i,j,1) .le. 250.0)t_gc(i,j,1)=302.5 enddo enddo ! ! READ IN LAT-LONS FROM STORM MESSAGE FILE ! IF (WRF_DM_ON_MONITOR()) THEN !orig OPEN(21,file='../MESSAGES/storm.center',status='old') !repository messages OPEN(21,file='storm.center',status='old') read(21,*)glat0 read(21,*)glon0 close(21) ENDIF CALL wrf_dm_bcast_bytes (glat0,IWORDSIZE) CALL wrf_dm_bcast_bytes (glon0,IWORDSIZE) WRITE(message,*)'glat0 and glon0',glat0,glon0 call wrf_debug(1,message) ! ! ! FIND THE CENTER OF THE HURRICANE BASED ON THE STORM MESSAGE FILE ! id0 = -999 jd0 = -999 factor = 4.0*ATAN(1.)/180. dist = -1 do j=jts,min(jte,jde-1) do i=its,min(ite,ide-1) diffx(i,j) = (hlon_gc(i,j) - glon0)*factor diffy(i,j) = (hlat_gc(i,j) - glat0)*factor cavlat= cos((hlat_gc(i,j) + glat0)*0.5*factor) dist1(i,j) = eradius*sqrt(cavlat*cavlat*diffx(i,j)*diffx(i,j) + diffy(i,j)*diffy(i,j)) if(dist1(i,j).LE.dist .OR. dist.LT.0.)then dist = dist1(i,j) id0 = i jd0 = j endif enddo enddo ! ! SECONDARY SEARCH FOR THE CENTER ! id1=-999 jd1=-999 minval_gh=MINVAL(ght_gc(ids:ide-1,jds:jde-1,10)) write(message,*)'MIN VAL OF GH NEW=',minval_gh call wrf_debug(1,message) do j=jts,min(jte,jde-1) do i=its,min(ite,ide-1) if(minval_gh .EQ. ght_gc(i,j,10))then id1=i jd1=j endif enddo enddo ! jbao-put storm in center id1=ide/2 jd1=jde/2 ! end jbao-put storm in center IF(abs(ID0-ID1) .GE. 5 .OR. abs(JD0-JD1) .GE. 5)THEN ! call wrf_error_fatal("LAT LON INFO IN STROM MESSAGE FILE IN INCORRECT") write(message,*) 'WARNING - ACTUAL STORM LOCATION DIFFERENT FROM MODEL LOCATION' call wrf_debug(1,message) WRITE(message,*)'OLD CENTER',ID0,JD0 call wrf_debug(1,message) WRITE(message,*)'NEW CENTER',ID1,JD1 call wrf_debug(1,message) ELSE WRITE(message,*)'OLD CENTER',ID0,JD0 call wrf_debug(1,message) WRITE(message,*)'NEW CENTER',ID1,JD1 call wrf_debug(1,message) ENDIF ! ! RECOMPUTE RADIUS FROM THE CENTER BASED ON MODEL/ANALYSED STORM CENTER ! do j=jts,min(jte,jde-1) do i=its,min(ite,ide-1) diffx(i,j) = (hlon_gc(i,j) - hlon_gc(id1,jd1))*factor diffy(i,j) = (hlat_gc(i,j) - hlat_gc(id1,jd1))*factor cavlat= cos((hlat_gc(i,j) + hlat_gc(id1,jd1))*0.5*factor) dist1(i,j) = eradius*sqrt(cavlat*cavlat*diffx(i,j)*diffx(i,j) + diffy(i,j)*diffy(i,j)) enddo enddo ght_out=ght_gc t_out=t_gc rh_out=rh_gc u_out=u_gc v_out=v_gc ! ! Beta is the intensity parameter ! IF(BETA1 .NE. 0.)THEN ! Adjustment of the GFDL vortex ! ! DEFINE THE NUMBER OF RECRUSSIVE OPERATIONS (NMAX) REQUIRED TO ! REMOVE OR MODIFY THE INITIAL GRIB ANALYSIS. IF beta1 IS SET TO ! < 0.01, THE VORTEX IS COMPLETELY FILTERED OUT. ALSO, FILTERING ! RELATED TO COMPLETE VORTEX REMOVAL IS ONLY CALLED AT THE INITIAL ! TIME ASSUMING THAT THE VORTEX IS FAR AWAY FROM THE BOUNDARIES. ! nres=0.18/dx nfilter=1 ! if(nres==3.and. internal_time_loop .eq. 1)then nmax=9*nfilter else nmax=1*nfilter endif WRITE(message,*)'CALLING HBFILTER' call wrf_debug(1,message) CALL hbfilter ( ght_gc,ght_out & &, nmax,start_z,end_z & &, ids,ide,jds,jde & &, ims,ime,jms,jme & &, its,ite,jts,jte ) WRITE(message,*)'COMPLETED FILTERING ght_gc' call wrf_debug(1,message) CALL hbfilter ( t_gc,t_out & &, nmax,start_z,end_z & &, ids,ide,jds,jde & &, ims,ime,jms,jme & &, its,ite,jts,jte ) WRITE(message,*)'COMPLETED FILTERING t_gc' call wrf_debug(1,message) CALL hbfilter ( rh_gc,rh_out & &, nmax,start_z,end_z & &, ids,ide,jds,jde & &, ims,ime,jms,jme & &, its,ite,jts,jte ) WRITE(message,*)'COMPLETED FILTERING rh_gc' call wrf_debug(1,message) CALL hbfilter ( u_gc,u_out & &, nmax,start_z,end_z & &, ids,ide,jds,jde & &, ims,ime,jms,jme & &, its,ite,jts,jte ) WRITE(message,*)'COMPLETED FILTERING u_gc' call wrf_debug(1,message) CALL hbfilter ( v_gc,v_out & &, nmax,start_z,end_z & &, ids,ide,jds,jde & &, ims,ime,jms,jme & &, its,ite,jts,jte ) WRITE(message,*)'COMPLETED FILTERING v_gc' call wrf_debug(1,message) WRITE(message,*)'END OF HBFILTER' call wrf_debug(1,message) ! do l=start_z,end_z do j=jts,min(jte,jde-1) do i=its,min(ite,ide-1) pert_ght(i,j,l) = (ght_gc(i,j,l) - ght_out(i,j,l)) ! vortex only pert_t(i,j,l) = (t_gc(i,j,l) - t_out(i,j,l)) pert_rh(i,j,l) = (rh_gc(i,j,l) - rh_out(i,j,l)) pert_u(i,j,l) = (u_gc(i,j,l) - u_out(i,j,l)) pert_v(i,j,l) = (v_gc(i,j,l) - v_out(i,j,l)) enddo enddo enddo ! do l=start_z,end_z do j=jts,min(jte,jde-1) do i=its,min(ite,ide-1) ght_gc(i,j,l) = ght_out(i,j,l) + beta1*pert_ght(i,j,l) ! vortex only t_gc(i,j,l) = t_out(i,j,l) + beta1*pert_t(i,j,l) rh_gc(i,j,l) = rh_out(i,j,l) + beta1*pert_rh(i,j,l) u_gc(i,j,l) = u_out(i,j,l) + beta1*pert_u(i,j,l) v_gc(i,j,l) = v_out(i,j,l) + beta1*pert_v(i,j,l) enddo enddo enddo ! ! DIAGNOSTICS ght_out=0.;t_out=0.;rh_out=0.;u_out=0.;v_out=0. do l=start_z,end_z do j=jts,min(jte,jde-1) do i=its,min(ite,ide-1) ght_out(i,j,l)= beta1*pert_ght(i,j,l) t_out(i,j,l) = beta1*pert_t(i,j,l) rh_out(i,j,l) = beta1*pert_rh(i,j,l) u_out(i,j,l) = beta1*pert_u(i,j,l) v_out(i,j,l) = beta1*pert_v(i,j,l) enddo enddo enddo ELSE ! Balanced vortex option !jwb IF(internal_time_loop .eq. 1)THEN ! use balancing only at the initial time IF(internal_time_loop .le. 2)THEN ! use balancing only at the initial time ! ! INSERT NEW VORTEX USING THE BALANCING ALGORITHM DEVELOPED BY ZHIHUA ZENG & ! JIAN-WEN BAO. THIS ALGORITHM IS BASED ON INVERSE BALANCE EQUATION ! WRITE(message,*)'CALLING TCBOGUS' call wrf_debug(1,message) allocate(u_temp(ide,jde,end_z));allocate(v_temp(ide,jde,end_z)) allocate(t_temp(ide,jde,end_z));allocate(ght_temp(ide,jde,end_z)) allocate(rh_temp(ide,jde,end_z));allocate(psc(ide,jde)) dx_km=110000.0*dx ! dx in kilometers CALL tcbogus( u_temp,v_temp,t_temp,ght_temp,rh_temp,psc, & dx_km,id1,jd1, & start_z,end_z+1, & ids,ide,jds,jde, & ids,ide,jds,jde, & ids,ide,jds,jde ) ! DIAGNOSTICS AND REVERSAL OF LOOP ght_out=0.;t_out=0.;rh_out=0.;u_out=0.;v_out=0. do l=start_z,end_z do j=jts,min(jte,jde-1) do i=its,min(ite,ide-1) if(l .lt.end_z)ght_out(i,j,l+1)= ght_temp(i,j,end_z-l+1) t_out(i,j,l) = t_temp(i,j,end_z-l+1) rh_out(i,j,l) = rh_temp(i,j,end_z-l+1) u_out(i,j,l) = u_temp(i,j,end_z-l+1) v_out(i,j,l) = v_temp(i,j,end_z-l+1) enddo enddo enddo do l=start_z,end_z do j=jts,min(jte,jde-1) do i=its,min(ite,ide-1) ght_gc(i,j,l) = ght_out(i,j,l) t_gc(i,j,l) = t_out(i,j,l) rh_gc(i,j,l) = rh_out(i,j,l) u_gc(i,j,l) = u_out(i,j,l) v_gc(i,j,l) = v_out(i,j,l) if(l .eq. 1)p_gc(i,j,1)=psc(i,j) enddo enddo enddo WRITE(message,*)'--------------- new ght_gc --------------------------' ! call wrf_debug(1,message) ! WRITE(message,*)ght_gc(100,100,:) ! call wrf_debug(1,message) WRITE(message,*)'--------------- new p_gc --------------------------' ! call wrf_debug(1,message) ! WRITE(message,*)p_gc(100,100,:) ! call wrf_debug(1,message) WRITE(message,*)'--------------- new rh_gc ---------------------------' ! call wrf_debug(1,message) ! WRITE(message,*)rh_gc(100,100,:) ! call wrf_debug(1,message) WRITE(message,*)'--------------- new t_gc --------------------------' ! call wrf_debug(1,message) ! WRITE(message,*)t_gc(100,100,:) ! call wrf_debug(1,message) ! WRITE(message,*)'---------------------------------------------------' ! call wrf_debug(1,message) WRITE(0,*)'tcbogus completed' WRITE(message,*)'tcbogus completed' call wrf_debug(1,message) ENDIF ! for internal_time_loop .eq. 1 ENDIF ! endif for beta=0.- 1.5 END SUBROUTINE vortex !------------------------------------------------------------------------- ! SUBROUTINE hbfilter ! ------------------- ! input - field to be filtered ! output - filtered output ! nmax - recrussive application of filter ( suggested 3 for 9 km and 27 for 3 km) ! start_z,end_zstart - usually 1 to ktop ! ids,ide,jds,jde - domain index, namely,imin,imax,jmin,jmax ! ims,ime,jms,jme - memory index, for 1 processor set this to above values ! its,ite,jts,jte - tile index, for 1 processor set this to above values SUBROUTINE hbfilter ( input, output & &, nmax,start_z,end_z & &, ids,ide,jds,jde & &, ims,ime,jms,jme & &, its,ite,jts,jte ) IMPLICIT NONE INTEGER, INTENT(IN) :: ids,ide,jds,jde,ims,ime,jms,jme,its,ite,jts,jte INTEGER, INTENT(IN) :: start_z,end_z,nmax REAL :: dx REAL, DIMENSION(IMS:IME,JMS:JME,start_z:end_z),INTENT(IN) :: INPUT REAL, DIMENSION(IMS:IME,JMS:JME,start_z:end_z),INTENT(INOUT) :: OUTPUT ! local INTEGER :: I,J,K,N,JMIN,JMAX,IMIN,IMAX,NRES,REPEAT REAL :: PI REAL, DIMENSION(11) :: M DATA M /2,3,4,2,5,6,7,2,8,9,2/ INTEGER, PARAMETER :: NF=11 REAL, DIMENSION(NF) :: FK REAL, DIMENSION(IMS:IME,NF) :: XTU REAL, DIMENSION(JMS:JME,NF) :: YTU !------------------------------------------------------------------ ! ! Purpose: This routine filters and removes ! hurricane signals (Kurihara et al., 1993, MWR) ! ! Called from: vortex ! ! Original Code: Qingfu Liu (EMC) ! Modification history: This is gopal's doing !----------------------------------------------------------------- ! ! DOMAIN IMIN=1 !(1,ID0-40) IMAX=IDE-1 !min(IDE-1,ID0+40) JMIN=1 !max(1,JD0-67) JMAX=JDE-1 !min(JDE-1,JD0+67) ! INPUT: Local variable defined for testing PI=4.0* ATAN(1.0) DO K=start_z,end_z DO J=JMIN,JMAX DO I=IMIN,IMAX OUTPUT(I,J,K)= INPUT(I,J,K) !SIN(2.0*PI/(21-1)*(I-1)) ! INPUT(I,J,K) ENDDO ENDDO ENDDO ! DEFINE FILTER FUNCTION DO N=1,NF FK(N)=0.5/(1-COS(2.*PI/M(N))) ENDDO DO K=start_z,end_z !.. DO ZONAL FILTER DO REPEAT=1,nmax DO J=JMIN,JMAX-1 DO N=1,NF XTU(IMIN,N) = OUTPUT(IMIN,J,K) XTU(IMAX,N) = OUTPUT(IMAX,J,K) ENDDO DO I=IMIN+1,IMAX-1 XTU(I,1) = OUTPUT(I,J,K)+FK(1)*(OUTPUT(I-1,J,K)+OUTPUT(I+1,J,K)-2.*OUTPUT(I,J,K)) ENDDO DO N=2,NF DO I=IMIN+1,IMAX-1 XTU(I,N)=XTU(I,N-1)+FK(N)*(XTU(I-1,N-1)+XTU(I+1,N-1)-2.*XTU(I,N-1)) ENDDO ENDDO DO I=IMIN,IMAX-1 OUTPUT(I,J,K) = XTU(I,NF) ENDDO ENDDO ! J loop ENDDO ! End recrussive repeat !.. DO MERIDIONAL FILTER DO REPEAT=1,nmax DO I=IMIN,IMAX DO N=1,NF YTU(JMIN,N) = OUTPUT(I,JMIN,K) YTU(JMAX,N) = OUTPUT(I,JMAX,K) ENDDO DO J=JMIN+1,JMAX-1 YTU(J,1) = OUTPUT(I,J,K) + FK(1)*(OUTPUT(I,J-1,K) + OUTPUT(I,J+1,K) -2.*OUTPUT(I,J,K)) ENDDO DO N = 2,NF DO J = JMIN+1,JMAX-1 YTU(J,N) = YTU(J,N-1) + FK(N)*(YTU(J-1,N-1) + YTU(J+1,N-1) - 2.*YTU(J,N-1)) ENDDO ENDDO DO J = JMIN,JMAX-1 OUTPUT(I,J,K) = YTU(J,NF) ENDDO ENDDO ! I loop ENDDO ! End recrussive repeat ENDDO ! K loop RETURN END SUBROUTINE hbfilter !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++c ! This program performs the idealized vortex initialization for ARW_WRF c ! was written by Zhihua Zeng (Shanghai Typhoon Institute/CMA) and Jian-Wen c ! Bao (NOAA/PSD) in 2009. c ! c ! Contact Eamil: c ! zengzh@mail.typhoon.gov.cn or Jian-Wen.Bao@noaa.gov c ! c ! References: c ! Wang Y.,1995:On inverse balance equation in sigma-coordinates c ! for model initialization, MWR, 123, 482-488. c ! Zeng Z. et al,2009: Simulation of spiral rainbands and vertical c ! resolution, 10th annual WRF user's workshop,23-26 June, c ! Boulder,CO. c ! !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ! This program performs the vortex initialisation using a nonlinear ! balance equation in sigma coordinate following Wang (1995, MWR) ! but with some modifications. !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ subroutine tcbogus( uw,vw,tw,ghw,rhw,psc & &, dsc,id1,jd1 & &, start_z,end_z1 & &, ids,ide,jds,jde & &, ims,ime,jms,jme & &, its,ite,jts,jte ) implicit none include 'BALANCE_PARS.F' include 'BALANCE_COMS.F' INTEGER, INTENT(IN):: start_z,end_z1,id1,jd1 INTEGER, INTENT(IN):: ids,ide,jds,jde,ims,ime,jms,jme,its,ite,jts,jte ! REAL, INTENT(IN) :: dsc REAL, DIMENSION(IMS:IME,JMS:JME), INTENT(INOUT) :: psc REAL, DIMENSION(IMS:IME,JMS:JME,start_z:end_z1),INTENT(INOUT) :: uw,vw,tw,ghw,rhw ! local variables INTEGER, PARAMETER :: nv=30 INTEGER :: i,j,k,km,kd,indx,jtc,itc INTEGER :: ji,ii,itime,jjt,nwest REAL :: f0,beta0,beta,pia,betk,txc,tyc,sigma REAL :: tk,t6,rh,rm,vm0,rm0,ts REAL, DIMENSION(JMS:JME) :: fc,fe REAL, DIMENSION(start_z:end_z1) :: up1,pp1,pp1o REAL, DIMENSION(start_z:end_z1-1) :: up,tvb REAL, DIMENSION(IMS:IME,JMS:JME) :: betkc REAL, DIMENSION(IMS:IME,JMS:JME,start_z:end_z1-1) :: uc,vc,tc,rqc,ppc REAL, DIMENSION(IMS:IME,JMS:JME,start_z:end_z1-1) :: tvbc,ghc,rhc character*255 message ! common /comp1/tk(kmx),t6(kmx1),rh(kmx),rm,vm0,ts km=end_z1-1 ! !******************************************************************* open(22,file='sigma.d',form='formatted') do k=1,end_z1 read(22,*) kd, sig(k) enddo 5 format(2x,i2,4x,f10.6) close(22) do k=1,km sig1(k)=0.5*(sig(k)+sig(k+1)) enddo ! open(20,file='input.d',form='formatted',status='unknown') read(20,1) indx read(20,2) nenv read(20,3) vm0 read(20,3) rm0 ! read(20,4) tsf0 close(20) 1 format(i1) 2 format(i2) 3 format(f5.1) 4 format(f8.2) ! rm=rm0*1000.0 ! ! write(6,*) vm0,rm0 ! write(6,*) nenv ! do k=1,km up(k)=0.0 enddo do k=1,end_z1 up1(k)=0.0 enddo jjt=0.0 ! nwest=-(ide-1)/2+25 jjt=max(jjt,nwest) write(6,*) jjt ! txc=id1 !(float(ide)+1.0)/3.0 !2.0 tyc=jd1 !(float(jde)+1.0)/3.0 !2.0 write(6,*)'txc= ', txc,'tyc= ',tyc ! call tem(tk,t6,rh,km,end_z1,pp1) ! ts=t6(end_z1) write(6,*) ts !------------------------------------- !--to calculate the Coriolis parameter !------------------------------------- pia=4.0*atan(1.0)/180.0 f0=14.584e-5*sin(pia*fi0) beta=14.584e-5*cos(pia*fi0)/6.3712e6 if(indx.eq.0) beta=0.0 do i=1,jde fc(i)=f0+beta*(float(i)-tyc)*dsc enddo !--------------------------------------------------- call env(tc,psc,fc,ide,jde,dsc,tyc,up1,km,end_z1) !--------------------------------------------------- do k=1,km do i=1,jde do j=1,ide tc(j,i,k)=tk(k)+tc(j,i,k) enddo enddo enddo ! call result(uc,vc,tc,rqc,psc,fc,ide,jde,dsc,tyc,txc, & up1,up,km,end_z1,beta) write(message,*)'out of result' call wrf_debug(1,message) ! itime=0.0 ! do k=1,km qb(k)=rqc(1,1,k) enddo ! do i=1,jde do k=1,km do j=1,ide ppc(j,i,k)=sig1(k)*psc(j,i) enddo enddo enddo ! !-------------------------------------- !--to calculate the height of the u,v,t !-------------------------------------- do k=1,km tvb(k)=tk(k)*(1.0+0.608*qb(k)) enddo betk=(rgas/grv)*log(1.0/sig1(km)) poz(km)=betk*tvb(km) ! write(10,*) poz(km) do k=km-1,1,-1 betk=0.5*(rgas/grv)*log(sig1(k+1)/sig1(k)) poz(k)=poz(k+1)+betk*(tvb(k)+tvb(k+1)) ! write(10,*) poz(k) enddo ! write(message,*)'before interpolation' call wrf_debug(1,message) call interpz(uc,vc,tc,ppc,rqc,ide,jde,km) ! if(vm0.gt.1.0) then ! !-------------------------------------- do k=1,km do i=1,jde do j=1,ide rhc(j,i,k)=rh(k) fe(i)=0.0 enddo enddo enddo !-------------------------------------- do k=1,km do i=1,jde do j=1,ide tvbc(j,i,k)=tc(j,i,k)*(1.0+0.608*rqc(j,i,k)) enddo enddo enddo do i=1,jde do j=1,ide betkc(j,i)=(rgas/grv)*log(0.01*psc(j,i)/pp1(km)) !! psc is hPa ghc(j,i,km)=betkc(j,i)*tvbc(j,i,km) enddo enddo do k=km-1,1,-1 do i=1,jde do j=1,ide ! write(63,*)k,i,j,pp1(k+1),pp1(k) if(pp1(k) .gt. 0.)then betkc(j,i)=0.5*(rgas/grv)*log(pp1(k+1)/pp1(k)) else betkc(j,i)=0.0 endif ghc(j,i,k)=ghc(j,i,k+1)+betkc(j,i)*(tvbc(j,i,k)+tvbc(j,i,k+1)) enddo enddo enddo ! write(message,*)'---------- test write out -----------' call wrf_debug(1,message) do k=1,km do i=1,jde do j=1,ide uw(j,i,k)=uc(j,i,k) vw(j,i,k)=vc(j,i,k) tw(j,i,k)=tc(j,i,k) ghw(j,i,k)=ghc(j,i,k) rhw(j,i,k)=rhc(j,i,k) if(i .eq. 10 .and. j .eq. 10) then write(message,*)k,uw(j,i,k),vw(j,i,k),ghw(j,i,k),tw(j,i,k) call wrf_debug(1,message) endif if(i .eq. 100 .and. j .eq. 98) then write(message,*)k,uw(j,i,k),vw(j,i,k),ghw(j,i,k),tw(j,i,k) call wrf_debug(1,message) endif enddo enddo enddo endif ! return end subroutine tcbogus ! !----------------------------------------------------------------- ! subroutine result(u,v,t,rq,ps,f,lq,lp,ds,ty0,tx0,up1, & up,km,end_z1,beta) ! !----------------------------------------------------------------- implicit none include 'BALANCE_PARS.F' include 'BALANCE_COMS.F' integer i,j,k,lq,lp,km,end_z1,nnt,jti,iti real ds,ds2,err,uee,tx0,ty0,beta,ajac,vor,forc1,forc2,forc3 real ees,rqs,sgg,reh,ta,pp,alpi real u(lq,lp,km),v(lq,lp,km),t(lq,lp,km),rq(lq,lp,km),ps(lq,lp) real alp(lq,lp),tt(lq,lp,end_z1),tt2(2:km),xx(lq,lp),forc(lq,lp) real ux(lq,lp),vy(lq,lp),tt1(lq,lp,km),gg1(lq,lp,end_z1),f(lp) real up1(end_z1),up(km) real tk,t6,rh,rm,vm0,ts character*255 message common /comp1/tk(kmx),t6(kmx1),rh(kmx),rm,vm0,ts ! ds2=ds*ds err=3.0e-6 do i=1,lp do j=1,lq forc(j,i)=0.0 alp(j,i)=log(ps(j,i)) enddo enddo ! uee=up1(end_z1) WRITE(message,*)'CALLING VORTEX 1' call wrf_debug(1,message) call vortex1(ux,vy,1.0,ds,lq,lp,ty0,tx0,uee) do 1 i=2,lp-1 do 1 j=2,lq-1 ajac=((ux(j+1,i)-ux(j-1,i))*(vy(j,i+1)-vy(j,i-1))-(ux(j,i+1) & -ux(j,i-1))*(vy(j+1,i)-vy(j-1,i)))/(2.0*ds2) vor=f(i)*(vy(j+1,i)-vy(j-1,i)-ux(j,i+1)+ux(j,i-1))/(2.0*ds) forc(j,i)=ds2*(ajac+vor-beta*ux(j,i))/(rgas*ts) 1 continue call posn(alp,forc,5.0e-10,1,lq,lp) do i=1,lp do j=1,lq ps(j,i)=exp(alp(j,i)) enddo enddo do 2 i=1,lp do 2 j=1,lq do 2 k=1,km t(j,i,k)=t(j,i,k)-tk(k) 2 continue ! do 5 k=1,end_z1 ! sgg=sig(k) uee=up1(k) WRITE(message,*)'CALLING VORTEX 2',k call wrf_debug(1,message) call vortex1(ux,vy,sgg,ds,lq,lp,ty0,tx0,uee) do 3 i=2,lp-1 do 3 j=2,lq-1 ajac=((ux(j+1,i)-ux(j-1,i))*(vy(j,i+1)-vy(j,i-1))-(ux(j,i+1) & -ux(j,i-1))*(vy(j+1,i)-vy(j-1,i)))/(2.0*ds2) vor=f(i)*(vy(j+1,i)-vy(j-1,i)-ux(j,i+1)+ux(j,i-1))/(2.0*ds) gg1(j,i,k)=ajac+vor-beta*ux(j,i) 3 continue ! 5 continue ! nnt=0 6 nnt=nnt+1 ! do 8 j=1,lq do 8 i=1,lp do k=2,km tt2(k)=t(j,i,k-1)+(t(j,i,k)-t(j,i,k-1))*(sig(k)-sig1(k-1)) & /(sig1(k)-sig1(k-1)) enddo do k=1,km if(k.eq.1) then tt(j,i,k)=(tt2(k+1)-t(j,i,k))/log(sig(k+1)/sig1(k)) else if(k.eq.km) then tt(j,i,k)=(t(j,i,k)-tt2(k))/log(sig1(k)/sig(k)) else tt(j,i,k)=(tt2(k+1)-tt2(k))/log(sig(k+1)/sig(k)) endif enddo 8 continue ! do 15 k=1,km ! do 9 j=2,lq-1 do 9 i=2,lp-1 alpi=alp(j-1,i)+alp(j,i-1)-4.0*alp(j,i)+alp(j+1,i)+alp(j,i+1) if(k.eq.1) then forc1=(t6(k+1)-tk(k))*alpi/log(sig(k+1)/sig1(k)) forc3=(gg1(j,i,k+1)-gg1(j,i,k+2))*ds2/ & (rgas*log(sig(k+2)/sig(k+1))) else forc1=(t6(k+1)-t6(k))*alpi/log(sig(k+1)/sig(k)) forc3=(gg1(j,i,k)-gg1(j,i,k+1))*ds2/(rgas*log(sig(k+1)/sig(k))) endif forc2=(tt(j+1,i,k)+tt(j,i,k))*(alp(j+1,i)-alp(j,i))-(tt(j,i,k)+ & tt(j-1,i,k))*(alp(j,i)-alp(j-1,i))+(tt(j,i+1,k)+tt(j,i,k))*(alp & (j,i+1)-alp(j,i))-(tt(j,i,k)+tt(j,i-1,k))*(alp(j,i)-alp(j,i-1)) forc(j,i)=forc1+0.5*forc2+forc3 9 continue do 10 j=1,lq do 10 i=1,lp xx(j,i)=t(j,i,k) 10 continue ! call posn(xx,forc,err,nnt,lq,lp) ! do 12 j=1,lq do 12 i=1,lp 12 t(j,i,k)=xx(j,i) ! 15 continue ! if(nnt.ne.4) goto 6 ! jti=int(tx0+0.5) iti=int(ty0+0.5) ! do k=1,km write(6,*) t(jti,iti,k) enddo ! do 18 k=1,km do 18 i=1,lp do 18 j=1,lq t(j,i,k)=t(j,i,k)+tk(k) reh=amin1(1.0,rh(k)*0.01) if(nphy.eq.0) reh=0.0 pp=sig1(k)*ps(j,i) ta=t(j,i,k) ees=611.2*exp(cs2*(ta-cs3)/(ta-cs4)) ees=amin1(0.5*pp,ees) rqs=0.622*ees/(pp-ees) rq(j,i,k)=reh*rqs 18 continue ! do 20 k=1,km sgg=sig1(k) uee=up(k) WRITE(message,*)'CALLING VORTEX 3',k call wrf_debug(1,message) call vortex1(ux,vy,sgg,ds,lq,lp,ty0,tx0,uee) do 20 i=1,lp do 20 j=1,lq u(j,i,k)=ux(j,i) v(j,i,k)=vy(j,i) 20 continue ! return end subroutine result ! !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ! subroutine vortex1(ux,vy,sgg,ds,lq,lp,ty0,tx0,uee) ! !---------------------------------------------------------------- implicit none include 'BALANCE_PARS.F' character*255 message integer i,j,lq,lp real ahem,pi,sgt,sgb,sgg,ds,ty0,tx0,uee,vm,b,r0,rr,rd,ty,tx real ux(lq,lp),vy(lq,lp) real tk,t6,rh,rm,vm0,ts,vt common /comp1/tk(kmx),t6(kmx1),rh(kmx),rm,vm0,ts ! WRITE(message,*)'INSIDE VORTEX' call wrf_debug(1,message) ahem=1.0 if(fi0.lt.0.0) ahem=-1.0 pi=2.0*atan(1.0) sgt=0.15 sgb=0.95 if(sgg.ge.sgb) then vm=vm0 else if(sgg.le.sgt) then vm=0.0 else vm=vm0*sin(pi*(sgg-sgt)/(sgb-sgt)) endif ! b=1.0 r0=10.0*rm do i=1,lp ty=(float(i)-ty0)*ds do j=1,lq tx=(float(j)-tx0)*2.0*ds ! factor 2.0: gopal's doing for E grid rr=sqrt(ty*ty+tx*tx) rd=rr/rm vt=vm/rm*(exp((1.0-rd**b)/b)-abs(rr-rm)/(r0-rm) & *exp((1.0-(r0/rm)**b)/b)) if(rr.ge.r0) vt=0.0 ux(j,i)=-vt*ty*ahem+uee vy(j,i)=vt*tx*ahem enddo enddo ! return end subroutine vortex1 ! !++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ! subroutine posn(xx,forc,err,nnt,lq,lp) ! !-------------------------------------------------------- implicit none include 'BALANCE_PARS.F' include 'BALANCE_COMS.F' ! integer i,j,nnt,lq,lp real xx(lq,lp),forc(lq,lp),err real*8 rji,ane,arm,anorm real*8 tem(lq,lp),wk(lq,lp),fc(lq,lp) ! do i=1,lp do j=1,lq tem(j,i)=xx(j,i) fc(j,i)=forc(j,i) enddo enddo ! if(nnt.eq.1.and.nenv.ne.0) then do 1 i=2,lp-1 do 1 j=2,lq-1 tem(j,i)=tem(j,i)-0.25*fc(j,i) 1 continue endif ! nt=0 2 nt=nt+1 ! anorm=0.0d0 ! do 10 i=2,lp-1 do 10 j=2,lq-1 rji=fc(j,i)-(tem(j-1,i)+tem(j,i-1)+ & tem(j+1,i)+tem(j,i+1)-4.0d0*tem(j,i)) ane=-0.25d0*rji wk(j,i)=tem(j,i)+ane arm=abs(ane) if(arm.gt.anorm) anorm=arm 10 continue ! do 20 i=2,lp-1 do 20 j=2,lq-1 20 tem(j,i)=wk(j,i) ! !Gopal's doing on June 17th, 2015. The line below was erroneously removed from the code. if(anorm.gt.err) goto 2 write(6,*) nt,anorm,err ! do i=1,lp do j=1,lq xx(j,i)=tem(j,i) enddo enddo ! return end subroutine posn ! !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ! subroutine tem(ttk,ttk1,rh,km,end_z1,pp1) ! !---------------------------------------------------------------- implicit none include 'BALANCE_PARS.F' include 'BALANCE_COMS.F' integer km,end_z1,nv,k,n1 parameter(nv=30) real ppk(km),ttk(km),ttk1(end_z1),tt(nv),pp(nv),rh(km) real pp1(end_z1) real rho(nv),t1,t2,p1,p2,pk,r1,r2,b,rk ! ! Wilis Island sounding (January) ! open(30,file='sound.d',form='formatted',status='old') do k=1,nv read(30,*) pp(k),tt(k),rho(k) enddo ! 2 format(1x,f7.1,7x,f7.2,8x,f4.1) close(30) ! do k=1,nv tt(k)=tt(k)+273.15 enddo ! do k=1,km ppk(k)=sig1(k)*pse/100.0 pp1(k)=sig(k)*pse/100.0 enddo pp1(end_z1)=pse/100.0 ! do k=1,km pk=ppk(k) ! do n1=1,nv-1 if(pp(n1).le.pk.and.pp(n1+1).gt.pk) then t1=tt(n1) t2=tt(n1+1) p1=pp(n1) p2=pp(n1+1) endif enddo write(6,*) ' p1, p2=',p1,p2 ! b=(t2-t1)/(p2-p1) ttk(k)=t1+b*(pk-p1) enddo ! do k=2,end_z1 pk=pp1(k) ! do n1=1,nv-1 if(pp(n1).le.pk.and.pp(n1+1).gt.pk) then t1=tt(n1) t2=tt(n1+1) p1=pp(n1) p2=pp(n1+1) endif enddo ! b=(t2-t1)/(p2-p1) ttk1(k)=t1+b*(pk-p1) enddo ttk1(1)=ttk(1) ! do k=1,km pk=ppk(k) ! do n1=1,nv-1 if(pp(n1).le.pk.and.pp(n1+1).gt.pk) then r1=rho(n1) r2=rho(n1+1) p1=pp(n1) p2=pp(n1+1) endif enddo ! b=(r2-r1)/(p2-p1) rk=r1+b*(pk-p1) rh(k)=rk enddo ! open(23,file='pp.d',form='formatted') do k=1,end_z1 write(23,3) k,sig(k),pp1(k),ttk1(k)-273.15 enddo ! do k=1,km write(23,4) k,sig1(k),ppk(k),ttk(k)-273.15,rh(k) enddo 3 format(2x,'k=',i2,', sig=',f8.4,', p=',f8.3,', t=',f7.2) 4 format(2x,'k=',i2,', sig=',f8.4,', p=',f8.3,', t=',f7.2,', rh=', & f6.3) close(23) ! return end subroutine tem subroutine env(t,ps,f,lq,lp,ds,ty0,up1,km,end_z1) ! !------------------------------------------------------------ ! This subroutine performs initialization for a zonal flow ! using a nonlinear balance equation in sigma coordinate. !------------------------------------------------------------ implicit none include 'BALANCE_PARS.F' include 'BALANCE_COMS.F' integer i,j,k,lq,lp,km,end_z1,nnt real ds,ty0,forc1,forc2,forc3 real alp(lp),tt(lp,end_z1),xx(lp),forc(lp) real alpy(lp),f(lp),t(lq,lp,km),ps(lq,lp) real tt2(2:km),up1(end_z1) real tk,t6,rh,rm,vm0,ts common /comp1/tk(kmx),t6(kmx1),rh(kmx),rm,vm0,ts ! it0=int(ty0+0.5) ! do 1 k=1,km do 1 i=1,lp do 1 j=1,lq 1 t(j,i,k)=0.0 ! do i=1,lp alpy(i)=-f(i)*up1(end_z1)/(rgas*ts) enddo alp(it0)=log(pse) do i=it0-1,1,-1 alp(i)=alp(i+1)-0.5*(alpy(i)+alpy(i+1))*ds enddo do i=it0+1,lp alp(i)=alp(i-1)+0.5*(alpy(i)+alpy(i-1))*ds enddo ! do 2 j=1,lq do 2 i=1,lp ps(j,i)=exp(alp(i)) 2 continue ! nnt=0 6 nnt=nnt+1 ! do 8 i=1,lp do k=2,km tt2(k)=t(1,i,k-1)+(t(1,i,k)-t(1,i,k-1))*(sig(k)-sig1(k-1)) & /(sig1(k)-sig1(k-1)) enddo do k=1,km if(k.eq.1) then tt(i,k)=(tt2(k+1)-t(1,i,k))/log(sig(k+1)/sig1(k)) else if(k.eq.km) then tt(i,k)=(t(1,i,k)-tt2(k))/log(sig1(k)/sig(k)) else tt(i,k)=(tt2(k+1)-tt2(k))/log(sig(k+1)/sig(k)) endif enddo 8 continue ! do 15 k=1,km ! Bug fix for idealized simulation by Lin Zhu. forc2. The bug is associated with ! thermal wind relationship in sheared environment. do 9 i=1,lp forc1=tt(i,k)*alpy(i) if(k.eq.1) then forc2= - f(i)*(up1(3)-up1(2))/(rgas*log(sig(3)/sig(2))) forc3=(t6(2)-tk(1))*alpy(i)/log(sig(2)/sig1(1)) else forc2= - f(i)*(up1(3)-up1(2))/(rgas*log(sig(3)/sig(2))) forc3=(t6(k+1)-t6(k))*alpy(i)/log(sig(k+1)/sig(k)) endif forc(i)=forc1+forc2+forc3 9 continue ! xx(it0)=0.0 do i=it0-1,1,-1 xx(i)=xx(i+1)-0.5*(forc(i)+forc(i+1))*ds enddo do i=it0+1,lp xx(i)=xx(i-1)+0.5*(forc(i)+forc(i-1))*ds enddo ! do i=1,lp do j=1,lq t(j,i,k)=xx(i) enddo enddo ! ! 15 continue ! if(nnt.ne.6) goto 6 ! return end subroutine env ! !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++ real function splin(x,a,b,c,d) implicit none real a,b,c,d real x(4) splin=(-a*x(1)+b*x(2)+c*x(3)-d*x(4))/81.0 return end function splin ! !------------------------------------------------------ ! subroutine interpz(u,v,t,pp,qv,lq,lp,km) ! !****************************************************** implicit none include 'BALANCE_PARS.F' include 'BALANCE_COMS.F' integer i,j,k,lq,lp,km real ees,pa,ta,rqs,betk real u(lq,lp,km),v(lq,lp,km),qv(lq,lp,km),t(lq,lp,km) real tv(lq,lp,km),pp(lq,lp,km),pzz(lq,lp,km) !*************************************************** ! do i=1,lp do k=1,km do j=1,lq tv(j,i,k)=t(j,i,k)*(1.0+0.608*qv(j,i,k)) enddo enddo enddo ! do i=1,lp do j=1,lq betk=(rgas/grv)*log(1.0/sig1(km)) pzz(j,i,km)=betk*tv(j,i,km) enddo enddo do i=1,lp do k=km-1,1,-1 do j=1,lq betk=0.5*(rgas/grv)*log(sig1(k+1)/sig1(k)) pzz(j,i,k)=pzz(j,i,k+1)+betk*(tv(j,i,k)+tv(j,i,k+1)) enddo enddo enddo ! do i=1,lp do k=1,km do j=1,lq pa=pp(j,i,k) ta=t(j,i,k) ees=611.2*exp(cs2*(ta-cs3)/(ta-cs4)) rqs=0.622*ees/(pa-ees) qv(j,i,k)=qv(j,i,k)/rqs enddo enddo enddo ! call interps(u,pzz,poz,lq,lp,km) call interps(v,pzz,poz,lq,lp,km) call interps(t,pzz,poz,lq,lp,km) call interps(pp,pzz,poz,lq,lp,km) call interps(qv,pzz,poz,lq,lp,km) ! do i=1,lp do k=1,km do j=1,lq pa=pp(j,i,k) ta=t(j,i,k) ees=611.2*exp(cs2*(ta-cs3)/(ta-cs4)) rqs=0.622*ees/(pa-ees) qv(j,i,k)=qv(j,i,k)*rqs enddo enddo enddo ! return end subroutine interpz ! !------------------------------------------------- ! this subroutine performs vertical interpotation !------------------------------------------------- ! subroutine interps(as,zz,poz,lq,lp,km) ! !------------------------------------------------- implicit none integer i,j,k,lq,lp,km,n1,nn real a1,a2,sg,sg1,sg2,b real as(lq,lp,km),poz(km),zz(lq,lp,km),wk(lq,lp,km) ! do 100 i=1,lp do 100 j=1,lq ! do 100 k=1,km sg=poz(k) if(sg.le.zz(j,i,km)) then n1=km-1 else if(sg.gt.zz(j,i,1)) then n1=1 else do nn=1,km-1 if(sg.gt.zz(j,i,nn+1).and.sg.le.zz(j,i,nn)) then n1=nn goto 20 endif enddo 20 continue endif a1=as(j,i,n1+1) a2=as(j,i,n1) sg1=zz(j,i,n1+1) sg2=zz(j,i,n1) sg2=zz(j,i,n1) b=(a2-a1)/(sg2-sg1) wk(j,i,k)=a1+b*(sg-sg1) ! 100 continue ! do 200 i=1,lp do 200 k=1,km do 200 j=1,lq as(j,i,k)=wk(j,i,k) 200 continue ! return end subroutine interps #if ( HWRF == 1 ) ! compute earth lat-lons for before interpolations. This is gopal's doing for ocean coupling !============================================================================================ SUBROUTINE EARTH_LATLON_hwrf ( HLAT,HLON,VLAT,VLON, & !Earth lat,lon at H and V points DLMD1,DPHD1,WBD1,SBD1, & !input res,west & south boundaries, CENTRAL_LAT,CENTRAL_LON, & ! central lat,lon, all in degrees IDS,IDE,JDS,JDE,KDS,KDE, & IMS,IME,JMS,JME,KMS,KME, & ITS,ITE,JTS,JTE,KTS,KTE ) !============================================================================ ! IMPLICIT NONE INTEGER, INTENT(IN ) :: IDS,IDE,JDS,JDE,KDS,KDE INTEGER, INTENT(IN ) :: IMS,IME,JMS,JME,KMS,KME INTEGER, INTENT(IN ) :: ITS,ITE,JTS,JTE,KTS,KTE REAL, INTENT(IN ) :: DLMD1,DPHD1,WBD1,SBD1 REAL, INTENT(IN ) :: CENTRAL_LAT,CENTRAL_LON REAL, DIMENSION(IMS:IME,JMS:JME), INTENT(OUT) :: HLAT,HLON,VLAT,VLON ! local INTEGER,PARAMETER :: KNUM=SELECTED_REAL_KIND(13) INTEGER :: I,J REAL(KIND=KNUM) :: WB,SB,DLM,DPH,TPH0,STPH0,CTPH0 REAL(KIND=KNUM) :: TDLM,TDPH,TLMH,TLMV,TLMH0,TLMV0,TPHH,TPHV,DTR REAL(KIND=KNUM) :: STPH,CTPH,STPV,CTPV,PI_2 REAL(KIND=KNUM) :: SPHH,CLMH,FACTH,SPHV,CLMV,FACTV REAL(KIND=KNUM), DIMENSION(IMS:IME,JMS:JME) :: GLATH,GLONH,GLATV,GLONV !------------------------------------------------------------------------- ! PI_2 = ACOS(0.) DTR = PI_2/90. WB = WBD1 * DTR ! WB: western boundary in radians SB = SBD1 * DTR ! SB: southern boundary in radians DLM = DLMD1 * DTR ! DLM: dlamda in radians DPH = DPHD1 * DTR ! DPH: dphi in radians TDLM = DLM + DLM ! TDLM: 2.0*dlamda TDPH = DPH + DPH ! TDPH: 2.0*DPH ! For earth lat lon only TPH0 = CENTRAL_LAT*DTR ! TPH0: central lat in radians STPH0 = SIN(TPH0) CTPH0 = COS(TPH0) DO J = JTS,MIN(JTE,JDE) !-1) ! H./ This loop takes care of zig-zag ! ! \.H starting points along j TLMH0 = WB - TDLM + MOD(J+1,2) * DLM ! ./ TLMH (rotated lats at H points) TLMV0 = WB - TDLM + MOD(J,2) * DLM ! H (//ly for V points) TPHH = SB + (J-1)*DPH ! TPHH (rotated lons at H points) are simple trans. TPHV = TPHH ! TPHV (rotated lons at V points) are simple trans. STPH = SIN(TPHH) CTPH = COS(TPHH) STPV = SIN(TPHV) CTPV = COS(TPHV) ! .H DO I = ITS,MIN(ITE,IDE) !-1) ! / TLMH = TLMH0 + I*TDLM ! \.H .U .H ! !H./ ----><---- SPHH = CTPH0 * STPH + STPH0 * CTPH * COS(TLMH) ! DLM + DLM GLATH(I,J)=ASIN(SPHH) ! GLATH: Earth Lat in radians CLMH = CTPH*COS(TLMH)/(COS(GLATH(I,J))*CTPH0) & - TAN(GLATH(I,J))*TAN(TPH0) IF(CLMH .GT. 1.) CLMH = 1.0 IF(CLMH .LT. -1.) CLMH = -1.0 FACTH = 1. IF(TLMH .GT. 0.) FACTH = -1. GLONH(I,J) = -CENTRAL_LON*DTR + FACTH*ACOS(CLMH) ENDDO DO I = ITS,MIN(ITE,IDE) !-1) TLMV = TLMV0 + I*TDLM SPHV = CTPH0 * STPV + STPH0 * CTPV * COS(TLMV) GLATV(I,J) = ASIN(SPHV) CLMV = CTPV*COS(TLMV)/(COS(GLATV(I,J))*CTPH0) & - TAN(GLATV(I,J))*TAN(TPH0) IF(CLMV .GT. 1.) CLMV = 1. IF(CLMV .LT. -1.) CLMV = -1. FACTV = 1. IF(TLMV .GT. 0.) FACTV = -1. GLONV(I,J) = -CENTRAL_LON*DTR + FACTV*ACOS(CLMV) ENDDO ENDDO ! Conversion to degrees (may not be required, eventually) DO J = JTS, MIN(JTE,JDE) !-1) DO I = ITS, MIN(ITE,IDE) !-1) HLAT(I,J) = GLATH(I,J) / DTR HLON(I,J)= -GLONH(I,J)/DTR IF(HLON(I,J) .GT. 180.) HLON(I,J) = HLON(I,J) - 360. IF(HLON(I,J) .LT. -180.) HLON(I,J) = HLON(I,J) + 360. ! VLAT(I,J) = GLATV(I,J) / DTR VLON(I,J) = -GLONV(I,J) / DTR IF(VLON(I,J) .GT. 180.) VLON(I,J) = VLON(I,J) - 360. IF(VLON(I,J) .LT. -180.) VLON(I,J) = VLON(I,J) + 360. ENDDO ENDDO END SUBROUTINE EARTH_LATLON_hwrf SUBROUTINE G2T2H_hwrf( SM,HRES_SM, & ! output grid index and weights HLAT,HLON, & ! target (nest) input lat lon in degrees DLMD1,DPHD1,WBD1,SBD1, & ! parent res, west and south boundaries CENTRAL_LAT,CENTRAL_LON, & ! parent central lat,lon, all in degrees P_IDE,P_JDE,P_IMS,P_IME,P_JMS,P_JME, & ! parent imax and jmax IDS,IDE,JDS,JDE,KDS,KDE, & ! target (nest) dIMEnsions IMS,IME,JMS,JME,KMS,KME, & ITS,ITE,JTS,JTE,KTS,KTE ) ! !*** Tom Black - Initial Version !*** Gopal - Revised Version for WRF (includes coincident grid points) !*** !*** GIVEN PARENT CENTRAL LAT-LONS, RESOLUTION AND WESTERN AND SOUTHERN BOUNDARY, !*** AND THE NESTED GRID LAT-LONS AT H POINTS, THIS ROUTINE FIRST LOCATES THE !*** INDICES,IIH,JJH, OF THE PARENT DOMAIN'S H POINTS THAT LIES CLOSEST TO THE !*** h POINTS OF THE NESTED DOMAIN ! !============================================================================ ! IMPLICIT NONE INTEGER, INTENT(IN ) :: IDS,IDE,JDS,JDE,KDS,KDE INTEGER, INTENT(IN ) :: IMS,IME,JMS,JME,KMS,KME INTEGER, INTENT(IN ) :: ITS,ITE,JTS,JTE,KTS,KTE INTEGER, INTENT(IN ) :: P_IDE,P_JDE,P_IMS,P_IME,P_JMS,P_JME REAL, INTENT(IN ) :: DLMD1,DPHD1,WBD1,SBD1 REAL, INTENT(IN ) :: CENTRAL_LAT,CENTRAL_LON REAL, DIMENSION(P_IMS:P_IME,P_JMS:P_JME), INTENT(IN) :: SM REAL, DIMENSION(IMS:IME,JMS:JME), INTENT(IN) :: HLAT,HLON REAL, DIMENSION(IMS:IME,JMS:JME), INTENT(OUT) :: HRES_SM ! local INTEGER,PARAMETER :: KNUM=SELECTED_REAL_KIND(13) INTEGER :: IMT,JMT,N2R,MK,K,I,J,DSLP0,DSLOPE,N INTEGER :: NROW,NCOL,KROWS REAL(KIND=KNUM) :: X,Y,Z,TLAT,TLON REAL(KIND=KNUM) :: PI_2,D2R,R2D,GLAT,GLON,DPH,DLM,TPH0,TLM0,WB,SB REAL(KIND=KNUM) :: ROW,COL,SLP0,TLATHC,TLONHC,DENOM,SLOPE REAL(KIND=KNUM) :: TLAT1,TLAT2,TLON1,TLON2,DLM1,DLM2,DLM3,DLM4,D1,D2 REAL(KIND=KNUM) :: DLA1,DLA2,DLA3,DLA4,S1,R1,DS1,AN1,AN2,AN3 ! Q REAL(KIND=KNUM) :: DL1,DL2,DL3,DL4,DL1I,DL2I,DL3I,DL4I,SUMDL,TLONO,TLATO REAL(KIND=KNUM) :: DTEMP REAL , DIMENSION(IMS:IME,JMS:JME) :: TLATHX,TLONHX INTEGER, DIMENSION(IMS:IME,JMS:JME) :: KOUTB REAL SUM,AMAXVAL REAL, DIMENSION (4, ims:ime, jms:jme ) :: NBWGT LOGICAL FLIP REAL, DIMENSION(IMS:IME,JMS:JME) :: HBWGT1,HBWGT2,HBWGT3,HBWGT4 INTEGER, DIMENSION(IMS:IME,JMS:JME) :: IIH,JJH character*255 message !------------------------------------------------------------------------------- IMT=2*P_IDE-2 ! parent i dIMEnsions JMT=P_JDE/2 ! parent j dIMEnsions PI_2=ACOS(0.) D2R=PI_2/90. R2D=1./D2R DPH=DPHD1*D2R DLM=DLMD1*D2R TPH0= CENTRAL_LAT*D2R TLM0=-CENTRAL_LON*D2R ! NOTE THE MINUS HERE WB=WBD1*D2R ! CONVERT NESTED GRID H POINTS FROM GEODETIC SB=SBD1*D2R SLP0=DPHD1/DLMD1 DSLP0=NINT(R2D*ATAN(SLP0)) DS1=SQRT(DPH*DPH+DLM*DLM) ! Q AN1=ASIN(DLM/DS1) AN2=ASIN(DPH/DS1) DO J = JTS,MIN(JTE,JDE) !-1) DO I = ITS,MIN(ITE,IDE) !-1) !*** !*** LOCATE TARGET h POINTS (HLAT AND HLON) ON THE PARENT DOMAIN AND !*** DETERMINE THE INDICES IN TERMS OF THE PARENT DOMAIN. FIRST !*** CONVERT NESTED GRID h POINTS FROM GEODETIC TO TRANSFORMED !*** COORDINATE ON THE PARENT GRID ! GLAT=HLAT(I,J)*D2R GLON= (360. - HLON(I,J))*D2R X=COS(TPH0)*COS(GLAT)*COS(GLON-TLM0)+SIN(TPH0)*SIN(GLAT) Y=-COS(GLAT)*SIN(GLON-TLM0) Z=COS(TPH0)*SIN(GLAT)-SIN(TPH0)*COS(GLAT)*COS(GLON-TLM0) TLAT=R2D*ATAN(Z/SQRT(X*X+Y*Y)) TLON=R2D*ATAN(Y/X) ! ROW=TLAT/DPHD1+JMT ! JMT IS THE CENTRAL ROW OF THE PARENT DOMAIN COL=TLON/DLMD1+P_IDE-1 ! (P_IDE-1) IS THE CENTRAL COLUMN OF THE PARENT DOMAIN NROW=INT(ROW + 0.001) ! ROUND-OFF IS AVOIDED WITHOUT USING NINT ON PURPOSE NCOL=INT(COL + 0.001) TLAT=TLAT*D2R TLON=TLON*D2R ! WRITE(60,*)'============================================================' ! WRITE(60,*)' ','i=',i,'j=',j !*** !*** !*** FIRST CONSIDER THE SITUATION WHERE THE POINT h IS AT !*** !*** V H !*** !*** !*** h !*** H V !*** !*** THEN LOCATE THE NEAREST H POINT ON THE PARENT GRID !*** IF(MOD(NROW,2).EQ.1.AND.MOD(NCOL,2).EQ.1.OR. & MOD(NROW,2).EQ.0.AND.MOD(NCOL,2).EQ.0)THEN TLAT1=(NROW-JMT)*DPH TLAT2=TLAT1+DPH TLON1=(NCOL-(P_IDE-1))*DLM TLON2=TLON1+DLM DLM1=TLON-TLON1 DLM2=TLON-TLON2 ! D1=ACOS(COS(TLAT)*COS(TLAT1)*COS(DLM1)+SIN(TLAT)*SIN(TLAT1)) ! D2=ACOS(COS(TLAT)*COS(TLAT2)*COS(DLM2)+SIN(TLAT)*SIN(TLAT2)) DTEMP=MIN(1.0_KNUM,COS(TLAT)*COS(TLAT1)*COS(DLM1)+SIN(TLAT)*SIN(TLAT1)) D1=ACOS(DTEMP) DTEMP=MIN(1.0_KNUM,COS(TLAT)*COS(TLAT2)*COS(DLM2)+SIN(TLAT)*SIN(TLAT2)) D2=ACOS(DTEMP) IF(D1.GT.D2)THEN NROW=NROW+1 ! FIND THE NEAREST H ROW NCOL=NCOL+1 ! FIND THE NEAREST H COLUMN ENDIF ! WRITE(60,*)'NEAREST PARENT IS:','col=',COL,'row=',ROW,'ncol=',NCOL,'nrow=',NROW ELSE !*** !*** NOW CONSIDER THE SITUATION WHERE THE POINT h IS AT !*** !*** H V !*** !*** !*** h !*** V H !*** !*** THEN LOCATE THE NEAREST H POINT ON THE PARENT GRID !*** !*** TLAT1=(NROW+1-JMT)*DPH TLAT2=TLAT1-DPH TLON1=(NCOL-(P_IDE-1))*DLM TLON2=TLON1+DLM DLM1=TLON-TLON1 DLM2=TLON-TLON2 ! D1=ACOS(COS(TLAT)*COS(TLAT1)*COS(DLM1)+SIN(TLAT)*SIN(TLAT1)) ! D2=ACOS(COS(TLAT)*COS(TLAT2)*COS(DLM2)+SIN(TLAT)*SIN(TLAT2)) DTEMP=MIN(1.0_KNUM,COS(TLAT)*COS(TLAT1)*COS(DLM1)+SIN(TLAT)*SIN(TLAT1)) D1=ACOS(DTEMP) DTEMP=MIN(1.0_KNUM,COS(TLAT)*COS(TLAT2)*COS(DLM2)+SIN(TLAT)*SIN(TLAT2)) D2=ACOS(DTEMP) IF(D1.LT.D2)THEN NROW=NROW+1 ! FIND THE NEAREST H ROW ELSE NCOL=NCOL+1 ! FIND THE NEAREST H COLUMN ENDIF ! WRITE(60,*)'NEAREST PARENT IS:','col=',COL,'row=',ROW,'ncol=',NCOL,'nrow=',NROW ENDIF KROWS=((NROW-1)/2)*IMT IF(MOD(NROW,2).EQ.1)THEN K=KROWS+(NCOL+1)/2 ELSE K=KROWS+P_IDE-1+NCOL/2 ENDIF !*** !*** WE NOW KNOW THAT THE INNER GRID POINT IN QUESTION IS !*** NEAREST TO THE CENTER K AS SEEN BELOW. WE MUST FIND !*** WHICH OF THE FOUR H-BOXES (OF WHICH THIS H POINT IS !*** A VERTEX) SURROUNDS THE INNER GRID h POINT IN QUESTION. !*** !** !*** H !*** !*** !*** !*** H V H !*** !*** !*** h !*** H V H V H !*** !*** !*** !*** H V H !*** !*** !*** !*** H !*** !*** !*** FIND THE SLOPE OF THE LINE CONNECTING h AND THE CENTER H. !*** N2R=K/IMT MK=MOD(K,IMT) ! IF(MK.EQ.0)THEN TLATHC=SB+(2*N2R-1)*DPH ELSE TLATHC=SB+(2*N2R+(MK-1)/(P_IDE-1))*DPH ENDIF ! IF(MK.LE.(P_IDE-1))THEN TLONHC=WB+2*(MK-1)*DLM ELSE TLONHC=WB+(2*(MK-(P_IDE-1))-1)*DLM ENDIF ! !*** EXECUTE CAUTION IF YOU NEED TO CHANGE THESE CONDITIONS. SINCE WE ARE !*** DEALING WITH SLOPES TO GENERATE DIAMOND SHAPE H BOXES, WE NEED TO BE !*** CAREFUL HERE ! IF(ABS(TLON-TLONHC) .LE. 1.E-4)TLONHC=TLON IF(ABS(TLAT-TLATHC) .LE. 1.E-4)TLATHC=TLAT DENOM=(TLON-TLONHC) ! !*** !***STORE THE LOCATION OF THE WESTERNMOST VERTEX OF THE H-BOX ON !***THE OUTER GRID THAT SURROUNDS THE h POINT ON THE INNER GRID. !*** !*** COINCIDENT CONDITIONS IF(DENOM.EQ.0.0)THEN IF(TLATHC.EQ.TLAT)THEN KOUTB(I,J)=K IIH(I,J) = NCOL JJH(I,J) = NROW TLATHX(I,J)=TLATHC TLONHX(I,J)=TLONHC HBWGT1(I,J)=1.0 HBWGT2(I,J)=0.0 HBWGT3(I,J)=0.0 HBWGT4(I,J)=0.0 ! WRITE(60,*)'TRIVIAL SOLUTION' ELSE ! SAME LONGITUDE BUT DIFFERENT LATS ! IF(TLATHC .GT. TLAT)THEN ! NESTED POINT SOUTH OF PARENT KOUTB(I,J)=K-(P_IDE-1) IIH(I,J) = NCOL-1 JJH(I,J) = NROW-1 TLATHX(I,J)=TLATHC-DPH TLONHX(I,J)=TLONHC-DLM ! WRITE(60,*)'VANISHING SLOPE, -ve: TLATHC-DPH, TLONHC-DLM' ELSE ! NESTED POINT NORTH OF PARENT KOUTB(I,J)=K+(P_IDE-1)-1 IIH(I,J) = NCOL-1 JJH(I,J) = NROW+1 TLATHX(I,J)=TLATHC+DPH TLONHX(I,J)=TLONHC-DLM ! WRITE(60,*)'VANISHING SLOPE, +ve: TLATHC+DPH, TLONHC-DLM' ENDIF !*** !*** !*** 4 !*** !*** h !*** 1 2 !*** !*** 3 !*** DL 1-4 ARE THE ANGULAR DISTANCES FROM h TO EACH VERTEX TLATO=TLATHX(I,J) TLONO=TLONHX(I,J) DLM1=TLON-TLONO DLA1=TLAT-TLATO ! Q ! DL1=ACOS(COS(TLAT)*COS(TLATO)*COS(DLM1)+SIN(TLAT)*SIN(TLATO)) ! Q DL1=SQRT(DLM1*DLM1+DLA1*DLA1) ! Q ! TLATO=TLATHX(I,J) TLONO=TLONHX(I,J)+2.*DLM DLM2=TLON-TLONO DLA2=TLAT-TLATO ! Q ! DL2=ACOS(COS(TLAT)*COS(TLATO)*COS(DLM2)+SIN(TLAT)*SIN(TLATO)) ! Q DL2=SQRT(DLM2*DLM2+DLA2*DLA2) ! Q ! TLATO=TLATHX(I,J)-DPH TLONO=TLONHX(I,J)+DLM DLM3=TLON-TLONO DLA3=TLAT-TLATO ! Q ! DL3=ACOS(COS(TLAT)*COS(TLATO)*COS(DLM3)+SIN(TLAT)*SIN(TLATO)) ! Q DL3=SQRT(DLM3*DLM3+DLA3*DLA3) ! Q TLATO=TLATHX(I,J)+DPH TLONO=TLONHX(I,J)+DLM DLM4=TLON-TLONO DLA4=TLAT-TLATO ! Q ! DL4=ACOS(COS(TLAT)*COS(TLATO)*COS(DLM4)+SIN(TLAT)*SIN(TLATO)) ! Q DL4=SQRT(DLM4*DLM4+DLA4*DLA4) ! Q ! THE BILINEAR WEIGHTS !*** !*** AN3=ATAN2(DLA1,DLM1) ! Q R1=DL1*SIN(AN2-AN3)/SIN(2.*AN1) S1=DL1*SIN(2.*PI_2-2*AN1-AN2+AN3)/SIN(2.*AN1) R1=R1/DS1 S1=S1/DS1 DL1I=(1.-R1)*(1.-S1) DL2I=R1*S1 DL3I=R1*(1.-S1) DL4I=(1.-R1)*S1 ! HBWGT1(I,J)=DL1I HBWGT2(I,J)=DL2I HBWGT3(I,J)=DL3I HBWGT4(I,J)=DL4I ! ENDIF ELSE ! !*** NON-COINCIDENT POINTS ! SLOPE=(TLAT-TLATHC)/DENOM DSLOPE=NINT(R2D*ATAN(SLOPE)) IF(DSLOPE.LE.DSLP0.AND.DSLOPE.GE.-DSLP0)THEN IF(TLON.GT.TLONHC)THEN ! IF(TLONHC.GE.-WB-DLM)CALL wrf_error_fatal("1H:NESTED DOMAIN TOO CLOSE TO THE BOUNDARY OF PARENT") KOUTB(I,J)=K IIH(I,J) = NCOL JJH(I,J) = NROW TLATHX(I,J)=TLATHC TLONHX(I,J)=TLONHC ! WRITE(60,*)'HERE WE GO1: TLATHC, TLONHC' ELSE ! IF(TLONHC.LE.WB+DLM)CALL wrf_error_fatal("2H:NESTED DOMAIN TOO CLOSE TO THE BOUNDARY OF PARENT") KOUTB(I,J)=K-1 IIH(I,J) = NCOL-2 JJH(I,J) = NROW TLATHX(I,J)=TLATHC TLONHX(I,J)=TLONHC -2.*DLM ! WRITE(60,*)'HERE WE GO2: TLATHC, TLONHC -2.*DLM' ENDIF ! ELSEIF(DSLOPE.GT.DSLP0)THEN IF(TLON.GT.TLONHC)THEN ! IF(TLATHC.GE.-SB-DPH)CALL wrf_error_fatal("3H:NESTED DOMAIN TOO CLOSE TO THE BOUNDARY OF PARENT") KOUTB(I,J)=K+(P_IDE-1)-1 IIH(I,J) = NCOL-1 JJH(I,J) = NROW+1 TLATHX(I,J)=TLATHC+DPH TLONHX(I,J)=TLONHC-DLM ! WRITE(60,*)'HERE WE GO3: TLATHC+DPH, TLONHC-DLM' ELSE ! IF(TLATHC.LE.SB+DPH)CALL wrf_error_fatal("4H:NESTED DOMAIN TOO CLOSE TO THE BOUNDARY OF PARENT") KOUTB(I,J)=K-(P_IDE-1) IIH(I,J) = NCOL-1 JJH(I,J) = NROW-1 TLATHX(I,J)=TLATHC-DPH TLONHX(I,J)=TLONHC-DLM ! WRITE(60,*)'HERE WE GO4: TLATHC-DPH, TLONHC-DLM' ENDIF ! ELSEIF(DSLOPE.LT.-DSLP0)THEN IF(TLON.GT.TLONHC)THEN ! IF(TLATHC.LE.SB+DPH)CALL wrf_error_fatal("5H:NESTED DOMAIN TOO CLOSE TO THE BOUNDARY OF PARENT") KOUTB(I,J)=K-(P_IDE-1) IIH(I,J) = NCOL-1 JJH(I,J) = NROW-1 TLATHX(I,J)=TLATHC-DPH TLONHX(I,J)=TLONHC-DLM ! WRITE(60,*)'HERE WE GO5: TLATHC-DPH, TLONHC-DLM' ELSE ! IF(TLATHC.GE.-SB-DPH)CALL wrf_error_fatal("6H:NESTED DOMAIN TOO CLOSE TO THE BOUNDARY OF PARENT") KOUTB(I,J)=K+(P_IDE-1)-1 IIH(I,J) = NCOL-1 JJH(I,J) = NROW+1 TLATHX(I,J)=TLATHC+DPH TLONHX(I,J)=TLONHC-DLM ! WRITE(60,*)'HERE WE GO6: TLATHC+DPH, TLONHC-DLM' ENDIF ENDIF ! !*** NOW WE WILL MOVE AS FOLLOWS: !*** !*** !*** 4 !*** !*** !*** !*** h !*** 1 2 !*** !*** !*** !*** !*** 3 !*** !*** !*** !*** DL 1-4 ARE THE ANGULAR DISTANCES FROM h TO EACH VERTEX TLATO=TLATHX(I,J) TLONO=TLONHX(I,J) DLM1=TLON-TLONO DLA1=TLAT-TLATO ! Q ! DL1=ACOS(COS(TLAT)*COS(TLATO)*COS(DLM1)+SIN(TLAT)*SIN(TLATO)) ! Q DL1=SQRT(DLM1*DLM1+DLA1*DLA1) ! Q ! TLATO=TLATHX(I,J) ! redundant computations TLONO=TLONHX(I,J)+2.*DLM DLM2=TLON-TLONO DLA2=TLAT-TLATO ! Q ! DL2=ACOS(COS(TLAT)*COS(TLATO)*COS(DLM2)+SIN(TLAT)*SIN(TLATO)) ! Q DL2=SQRT(DLM2*DLM2+DLA2*DLA2) ! Q ! TLATO=TLATHX(I,J)-DPH TLONO=TLONHX(I,J)+DLM DLM3=TLON-TLONO DLA3=TLAT-TLATO ! Q ! DL3=ACOS(COS(TLAT)*COS(TLATO)*COS(DLM3)+SIN(TLAT)*SIN(TLATO)) ! Q DL3=SQRT(DLM3*DLM3+DLA3*DLA3) ! Q ! TLATO=TLATHX(I,J)+DPH TLONO=TLONHX(I,J)+DLM DLM4=TLON-TLONO DLA4=TLAT-TLATO ! Q ! DL4=ACOS(COS(TLAT)*COS(TLATO)*COS(DLM4)+SIN(TLAT)*SIN(TLATO)) ! Q DL4=SQRT(DLM4*DLM4+DLA4*DLA4) ! Q ! THE BILINEAR WEIGHTS !*** AN3=ATAN2(DLA1,DLM1) ! Q R1=DL1*SIN(AN2-AN3)/SIN(2.*AN1) S1=DL1*SIN(2.*PI_2-2*AN1-AN2+AN3)/SIN(2.*AN1) R1=R1/DS1 S1=S1/DS1 DL1I=(1.-R1)*(1.-S1) DL2I=R1*S1 DL3I=R1*(1.-S1) DL4I=(1.-R1)*S1 ! HBWGT1(I,J)=DL1I HBWGT2(I,J)=DL2I HBWGT3(I,J)=DL3I HBWGT4(I,J)=DL4I ! ENDIF ! !*** FINALLY STORE IIH IN TERMS OF E-GRID INDEX ! IIH(I,J)=NINT(0.5*IIH(I,J)) ENDDO ENDDO ! !*** EXTENSION TO NEAREST NEIGHBOR ! DO J = JTS,MIN(JTE,JDE) !-1) DO I = ITS,MIN(ITE,IDE) !-1) NBWGT(1,I,J)=HBWGT1(I,J) NBWGT(2,I,J)=HBWGT2(I,J) NBWGT(3,I,J)=HBWGT3(I,J) NBWGT(4,I,J)=HBWGT4(I,J) ENDDO ENDDO DO J = JTS,MIN(JTE,JDE) !-1) DO I = ITS,MIN(ITE,IDE) !-1) AMAXVAL=0. DO N=1,4 AMAXVAL=amax1(NBWGT(N,I,J),AMAXVAL) ENDDO ! FLIP=.TRUE. SUM=0.0 DO N=1,4 IF(AMAXVAL .EQ. NBWGT(N,I,J) .AND. FLIP)THEN NBWGT(N,I,J)=1.0 FLIP=.FALSE. ELSE NBWGT(N,I,J)=0.0 ENDIF SUM=SUM+NBWGT(N,I,J) IF(SUM .GT. 1.0)CALL wrf_error_fatal ( "horizontal interp error - interp_hnear_nmm" ) ENDDO IF((NBWGT(1,I,J)+NBWGT(2,I,J)+NBWGT(3,I,J)+NBWGT(4,I,J)) .NE. 1)THEN WRITE(message,*)'------------------------------------------------------------------------' call wrf_debug(1,message) WRITE(message,*)'FATAL: SOMETHING IS WRONG WITH THE WEIGHTS IN module_initialize_real.F' call wrf_debug(1,message) WRITE(message,*)'------------------------------------------------------------------------' call wrf_debug(1,message) STOP ENDIF ! WRITE(66,*)I,J,NBWGT(1,I,J),NBWGT(2,I,J),NBWGT(3,I,J),NBWGT(4,I,J) ENDDO ENDDO DO J=MAX(3,JTS),MIN(JTE,JDE) !-1) DO I=MAX(3,ITS),MIN(ITE,IDE) !-1) IF(MOD(JJH(I,J),2) .NE. 0)THEN ! 1,3,5,7 HRES_SM(I,J) = NBWGT(1,I,J)*SM(IIH(I,J),JJH(I,J) ) & + NBWGT(2,I,J)*SM(IIH(I,J)+1, JJH(I,J) ) & + NBWGT(3,I,J)*SM(IIH(I,J), JJH(I,J)-1) & + NBWGT(4,I,J)*SM(IIH(I,J), JJH(I,J)+1) ! WRITE(68,*)I,J,SM(IIH(I,J),JJH(I,J)),SM(IIH(I,J)+1, JJH(I,J)), & ! SM(IIH(I,J), JJH(I,J)-1),SM(IIH(I,J), JJH(I,J)+1),HRES_SM(I,J) ELSE HRES_SM(I,J) = NBWGT(1,I,J)*SM(IIH(I,J), JJH(I,J) ) & + NBWGT(2,I,J)*SM(IIH(I,J)+1, JJH(I,J) ) & + NBWGT(3,I,J)*SM(IIH(I,J)+1, JJH(I,J)-1) & + NBWGT(4,I,J)*SM(IIH(I,J)+1, JJH(I,J)+1) ! WRITE(68,*)I,J,SM(IIH(I,J),JJH(I,J)),SM(IIH(I,J)+1, JJH(I,J)), & ! SM(IIH(I,J)+1, JJH(I,J)-1),SM(IIH(I,J)+1, JJH(I,J)+1),HRES_SM(I,J) ENDIF ENDDO ENDDO ! Boundary treatment in J direction DO J=MAX(3,JTS),MIN(JTE,JDE) HRES_SM(2,J)=HRES_SM(3,J) HRES_SM(1,J)=HRES_SM(2,J) END DO ! Boundary treatment in J direction and 4 corners DO I=ITS,MIN(ITE,IDE) HRES_SM(I,2)=HRES_SM(I,3) HRES_SM(I,1)=HRES_SM(I,2) END DO RETURN END SUBROUTINE G2T2H_hwrf !======================================================================================== ! end gopal's doing for ocean coupling !============================================================================================ #endif END MODULE module_initialize_ideal