!REAL:MODEL_LAYER:INITIALIZATION #ifndef VERT_UNIT ! This MODULE holds the routines which are used to perform various initializations ! for the individual domains, specifically for the Eulerian, mass-based coordinate. !----------------------------------------------------------------------- !****MARS: modified May 2007 MODULE module_initialize USE module_bc USE module_configure USE module_domain USE module_io_domain USE module_model_constants USE module_state_description USE module_timing USE module_soil_pre USE module_date_time #ifdef DM_PARALLEL USE module_dm #endif REAL , SAVE :: p_top_save INTEGER :: internal_time_loop 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 :: dyn_opt INTEGER :: idum1, idum2 CALL nl_get_dyn_opt ( 1, dyn_opt ) CALL set_scalar_indices_from_config ( head_grid%id , idum1, idum2 ) IF ( dyn_opt .eq. 1 & .or. dyn_opt .eq. 2 & .or. dyn_opt .eq. 3 & ) THEN CALL init_domain_rk( grid & ! #include "em_actual_new_args.inc" ! ) ELSE WRITE(0,*)' init_domain: unknown or unimplemented dyn_opt = ',dyn_opt CALL wrf_error_fatal ( 'ERROR-dyn_opt-wrong-in-namelist' ) ENDIF END SUBROUTINE init_domain !------------------------------------------------------------------- SUBROUTINE init_domain_rk ( grid & ! #include "em_dummy_new_args.inc" ! ) USE module_optional_si_input IMPLICIT NONE ! Input space and data. No gridded meteorological data has been stored, though. ! TYPE (domain), POINTER :: grid TYPE (domain) :: grid #include "em_dummy_new_decl.inc" 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 :: loop , num_seaice_changes 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 INTEGER :: ns ! Local data INTEGER :: error REAL :: p_surf, p_level REAL :: cof1, cof2 REAL :: qvf , qvf1 , qvf2 , pd_surf REAL :: p00 , t00 , a REAL :: hold_znw LOGICAL :: were_bad LOGICAL :: stretch_grid, dry_sounding, debug INTEGER IICOUNT REAL :: p_top_requested , temp INTEGER :: num_metgrid_levels REAL , DIMENSION(max_eta) :: eta_levels REAL :: max_dz ! INTEGER , PARAMETER :: nl_max = 1000 ! REAL , DIMENSION(nl_max) :: grid%em_dn integer::oops1,oops2 REAL :: zap_close_levels INTEGER :: force_sfc_in_vinterp INTEGER :: interp_type , lagrange_order LOGICAL :: lowest_lev_from_sfc LOGICAL :: we_have_tavgsfc INTEGER :: lev500 , loop_count REAL :: zl , zu , pl , pu , z500 , dz500 , tvsfc , dpmu !-- Carsel and Parrish [1988] REAL , DIMENSION(100) :: lqmi !****MARS INTEGER :: sizegcm, kold, knew,inew,jnew REAL :: pa, indic, p1, p2, pn REAL, ALLOCATABLE, DIMENSION (:,:,:) :: sig, ap, bp, box REAL :: randnum REAL :: tiso REAl :: yeah, yeahc, yeah2 !****MARS INTEGER :: hypsometric_opt = 1 ! classic LOGICAL :: interp_theta = .true. ! classic !INTEGER :: hypsometric_opt = 2 ! Wee et al. 2012 correction !LOGICAL :: interp_theta = .false. ! Wee et al. 2012 correction REAL :: pfu, pfd, phm REAL :: tpot #ifdef DM_PARALLEL # include "em_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 ; ! note that tile is entire patch its = grid%sp32 ; ite = grid%ep32 ; ! note that tile is entire patch jts = grid%sp33 ; jte = grid%ep33 ; ! note that 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 ; ! note that tile is entire patch jts = grid%sp32 ; jte = grid%ep32 ; ! note that tile is entire patch kts = grid%sp33 ; kte = grid%ep33 ; ! note that 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 ; ! note that tile is entire patch kts = grid%sp32 ; kte = grid%ep32 ; ! note that tile is entire patch jts = grid%sp33 ; jte = grid%ep33 ; ! note that tile is entire patch END SELECT CALL model_to_grid_config_rec ( grid%id , model_config_rec , config_flags ) ! 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%step_number = 0 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 ! To define the base state, we call a USER MODIFIED routine to set the three ! necessary constants: p00 (sea level pressure, Pa), t00 (sea level temperature, K), ! and A (temperature difference, from 1000 mb to 300 mb, K). CALL const_module_initialize ( p00 , t00 , a , tiso ) #if 0 !KLUDGE, this is for testing only if ( flag_metgrid .eq. 1 ) then read (20+grid%id) grid%em_ht_gc read (20+grid%id) grid%em_xlat_gc read (20+grid%id) grid%em_xlong_gc read (20+grid%id) msft read (20+grid%id) msfu read (20+grid%id) msfv read (20+grid%id) f read (20+grid%id) e read (20+grid%id) sina read (20+grid%id) cosa read (20+grid%id) grid%landmask read (20+grid%id) grid%landusef read (20+grid%id) grid%soilctop read (20+grid%id) grid%soilcbot read (20+grid%id) grid%vegcat read (20+grid%id) grid%soilcat else write (20+grid%id) grid%em_ht write (20+grid%id) grid%em_xlat write (20+grid%id) grid%em_xlong write (20+grid%id) msft write (20+grid%id) msfu write (20+grid%id) msfv write (20+grid%id) f write (20+grid%id) e write (20+grid%id) sina write (20+grid%id) cosa write (20+grid%id) grid%landmask write (20+grid%id) grid%landusef write (20+grid%id) grid%soilctop write (20+grid%id) grid%soilcbot write (20+grid%id) grid%vegcat write (20+grid%id) grid%soilcat endif #endif ! Is there any vertical interpolation to do? The "old" data comes in on the correct ! vertical locations already. IF ( flag_metgrid .EQ. 1 ) THEN ! <----- START OF VERTICAL INTERPOLATION PART ----> ! Variables that are named differently between SI and WPS. DO j = jts, MIN(jte,jde-1) DO i = its, MIN(ite,ide-1) grid%tmn(i,j) = grid%em_tmn_gc(i,j) grid%xlat(i,j) = grid%em_xlat_gc(i,j) grid%xlong(i,j) = grid%em_xlong_gc(i,j) grid%ht(i,j) = grid%em_ht_gc(i,j) !!****MARS grid%m_tsurf(i,j) = grid%em_tsk_gc(i,j) grid%m_albedo(i,j) = grid%em_albedo_gcm_gc(i,j) grid%m_ti(i,j) = grid%em_therm_inert_gc(i,j) grid%slpx(i,j) = grid%em_slpx_gc(i,j) grid%slpy(i,j) = grid%em_slpy_gc(i,j) grid%m_emiss(i,j) = grid%st000010(i,j) grid%m_co2ice(i,j) = grid%st010040(i,j) grid%m_h2oice(i,j) = grid%sm100200(i,j) grid%m_q2(i,:,j) = 0. !! one more security ... co2ice cannot be negative IF (grid%m_co2ice(i,j) .lt. 0.) grid%m_co2ice(i,j)=0. IF (grid%m_h2oice(i,j) .lt. 0.) grid%m_h2oice(i,j)=0. DO k = 1, config_flags%num_soil_layers grid%m_tsoil(i,k,j)=grid%em_tsoil_gc(i,k+1,j) !!ici k+1, because em_tsoil_gc dim is num_metgrid_levels !! ENDDO !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! #ifdef NEWPHYS grid%m_wstar(i,j)=0. grid%m_fluxrad(i,j)=0. grid%m_z0(i,j) = 0. grid%m_z0(i,j) = grid%em_z0_gc(i,j)*0.01 !! in cm in surface.nc but in m in physiq.F !IF (config_flags%init_Z0 .ne. 0.) THEN ! grid%z0 = grid%z0*0. + config_flags%init_Z0 !ENDIF ! here, that bit is necessary for new soil model ! IF (config_flags%init_TI .ne. 0.) THEN grid%m_ti = grid%m_ti*0. + config_flags%init_TI print *, 'constant thermal inertia ', config_flags%init_TI ENDIF DO k = 1, config_flags%num_soil_layers grid%m_isoil(i,k,j)=grid%em_isoil_gc(i,k+1,j) grid%m_dsoil(i,k,j)=grid%em_dsoil_gc(i,k+1,j) ENDDO DO k = 1, config_flags%num_soil_layers !!!!!!!!!!!!!!!!! DONE in soil_setting.F IF (grid%m_dsoil(i,k,j) == -999.) THEN !! old soil depths (or) no info in files grid%m_dsoil(i,k,j) = sqrt(887.75/3.14)*((2.**(k-0.5))-1.) * grid%m_ti(i,j) / wvolcapa !!! ATTENTION il faut interpoler si le nombre de niveaux change !!! voir soil_setting.F (olddepthdef=.true. ; interpol=.true.) !!! mais: en meso-echelle on a juste a prendre le mm nombre de niveaux que le GCM ENDIF IF (grid%m_isoil(i,k,j) == -999.) THEN !! old soil model (or) no 3D thermal inertia grid%m_isoil(i,k,j) = grid%m_ti(i,j) ELSE IF (grid%m_dsoil(i,k,j) .le. sqrt(88775./3.14) * grid%m_ti(i,j) / wvolcapa) THEN grid%m_isoil(i,k,j) = grid%m_ti(i,j) !! if depth < skin depth, we use hi-res TI ELSE !! if depth > skin depth, we use low-res (GCM) TI !! except for a transition layer !! EM: and, well, it would be wrong to sum up TI values !! EM: (cf. last page of soil model technical document) IF (grid%m_dsoil(i,k-1,j) .le. sqrt(88775./3.14) * grid%m_ti(i,j) / wvolcapa) THEN grid%m_isoil(i,k,j) = & sqrt( & ( grid%m_dsoil(i,k+1,j) - grid%m_dsoil(i,k-1,j) ) & / & ( ( (grid%m_dsoil(i,k,j) - grid%m_dsoil(i,k-1,j)) & / (grid%m_isoil(i,k-1,j)*grid%m_isoil(i,k-1,j)) ) & + & ( (grid%m_dsoil(i,k+1,j) - grid%m_dsoil(i,k,j)) & / (grid%m_isoil(i,k+1,j)*grid%m_isoil(i,k+1,j)) ) & ) & ) ENDIF !! grid%m_isoil(i,k-1,j) was changed at previous step to value grid%m_ti(i,j) !! grid%m_isoil(i,k+1,j) is defined to large-scale value grid%em_isoil_gc ENDIF ENDIF IF (grid%m_tsoil(i,k,j) .lt. 20.) THEN !!! une securite pour les anciens diagfi qui n'ont que 10 niveaux IF (k .ne. 1) grid%m_tsoil(i,k,j) = grid%m_tsoil(i,k-1,j) ENDIF !!!!!!!!!!!!!!!!! DONE in soil_setting.F ENDDO #endif !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! grid%m_gw(i,1,j)=grid%st040100(i,j) !!ZMEA grid%m_gw(i,2,j)=grid%st100200(i,j) !!ZSTD grid%m_gw(i,3,j)=grid%sm000010(i,j) !!ZSIG grid%m_gw(i,4,j)=grid%sm010040(i,j) !!ZGAM grid%m_gw(i,5,j)=grid%sm040100(i,j) !!ZTHE END DO END DO !!****MARS !!****MARS !! User-defined constants initialisations !! defined by the namelist entries !! init_TI : fixed value for thermal inertia !! init_AL : fixed value for albedo !! init_U : fixed value for zonal wind !! init_V : fixed value for meridional wind !! init_WX & init_WY : fixed wind profile taken at these coordinates !! init_MU : multiply zonal wind by a constant !! init_MV : multiply meridional wind by a constant !! init_LES : LES mode (LOGICAL) IF (config_flags%init_TI .ne. 0.) THEN !DO j = jts, MIN(jte,jde-1) !DO i = its, MIN(ite,ide-1) ! grid%m_ti(i,j) = config_flags%init_TI !ENDDO !ENDDO grid%m_ti = grid%m_ti*0. + config_flags%init_TI print *, 'constant thermal inertia ', config_flags%init_TI !yeah2=0. !yeahc=0. !DO j = jts, MIN(jte,jde-1) !DO i = its, MIN(ite,ide-1) !yeah2 = grid%tsk(i,j) + yeah2 !yeahc = yeahc + 1. !ENDDO !ENDDO !print *, 'constant skin temperature ', yeah2 / yeahc !DO j = jts, MIN(jte,jde-1) !DO i = its, MIN(ite,ide-1) !grid%tsk(i,j) = yeah2 / yeahc !ENDDO !ENDDO ! ! DO k = 1, config_flags%num_soil_layers !yeah=0. !yeahc=0. !DO j = jts, MIN(jte,jde-1) !DO i = its, MIN(ite,ide-1) !yeah = grid%m_tsoil(i,k,j) + yeah !yeahc = yeahc + 1. !ENDDO !ENDDO !print *, 'constant soil temperature ', k, yeah / yeahc !DO j = jts, MIN(jte,jde-1) !DO i = its, MIN(ite,ide-1) !grid%m_tsoil(i,k,j) = yeah / yeahc !ENDDO !ENDDO ! ENDDO ENDIF IF (config_flags%init_AL .ne. 0.) THEN grid%m_albedo = grid%m_albedo*0. + config_flags%init_AL print *, 'constant albedo ', config_flags%init_AL ENDIF IF ( (config_flags%init_WX .ne. 0) .and. (config_flags%init_WY .ne. 0) ) THEN DO j = jts, MIN(jte,jde-1) DO i = its, MIN(ite,ide-1) grid%em_u_gc(i,:,j)=grid%em_u_gc(config_flags%init_WX,:,config_flags%init_WY) ! zonal wind grid%em_v_gc(i,:,j)=grid%em_v_gc(config_flags%init_WX,:,config_flags%init_WY) ! meridional wind ENDDO ENDDO !! FIX for the STAGGERED SPECIFICITY grid%em_u_gc(MIN(ite,ide-1)+1,:,:)=grid%em_u_gc(MIN(ite,ide-1),:,:) grid%em_v_gc(:,:,MIN(jte,jde-1)+1)=grid%em_v_gc(:,:,MIN(jte,jde-1)) !! CHECK print *, 'wind profile' print *, 'took at ...', config_flags%init_WX, config_flags%init_WY print *, '--zonal' print *, grid%em_u_gc(config_flags%init_WX,:,config_flags%init_WY) print *, '--meridional' print *, grid%em_v_gc(config_flags%init_WX,:,config_flags%init_WY) ENDIF IF (config_flags%init_MU .ne. 0.) THEN grid%em_u_gc = grid%em_u_gc*config_flags%init_MU print *, 'multiply zonal wind ', config_flags%init_MU ENDIF IF (config_flags%init_MV .ne. 0.) THEN grid%em_v_gc = grid%em_v_gc*config_flags%init_MV print *, 'multiply meridional wind ', config_flags%init_MV ENDIF IF (config_flags%init_LES) THEN print *, '*** LES MODE ***' print *, 'setting uniform values and profiles' print *, 'u', grid%em_u_gc(its+1,:,jts+1) print *, 'v', grid%em_v_gc(its+1,:,jts+1) print *, 't', grid%em_t_gc(its+1,:,jts+1) print *, 'p', grid%em_rh_gc(its+1,:,jts+1) print *, 'geop', grid%em_ght_gc(its+1,:,jts+1) print *, 'albedo', grid%m_albedo(its+1,jts+1) print *, 'thermal inertia', grid%m_ti(its+1,jts+1) print *, 'topography', grid%ht(its+1,jts+1) print *, 'toposoil', grid%toposoil(its+1,jts+1) print *, 'surface temperature', grid%m_tsurf(its+1,jts+1) print *, 'surface pressure', grid%psfc(its+1,jts+1), grid%em_psfc_gc(its+1,jts+1) DO j = jts, MIN(jte,jde-1) DO i = its, MIN(ite,ide-1) grid%em_u_gc(i,:,j)=grid%em_u_gc(its+1,:,jts+1) grid%em_v_gc(i,:,j)=grid%em_v_gc(its+1,:,jts+1) grid%em_t_gc(i,:,j)=grid%em_t_gc(its+1,:,jts+1) grid%em_rh_gc(i,:,j)=grid%em_rh_gc(its+1,:,jts+1) grid%em_ght_gc(i,:,j) = grid%em_ght_gc(its+1,:,jts+1) grid%m_albedo(i,j) = grid%m_albedo(its+1,jts+1) grid%m_ti(i,j) = grid%m_ti(its+1,jts+1) grid%ht(i,j) = grid%ht(its+1,jts+1) grid%toposoil(i,j) = grid%toposoil(its+1,jts+1) grid%m_tsurf(i,j) = grid%m_tsurf(its+1,jts+1) grid%psfc(i,j) = grid%psfc(its+1,jts+1) grid%em_psfc_gc(i,j) = grid%em_psfc_gc(its+1,jts+1) grid%slpx(i,j) = 0. grid%slpy(i,j) = 0. grid%m_emiss(i,j) = 0.95 grid%m_co2ice(i,j) = 0. grid%m_h2oice(i,j) = 0. grid%m_tsoil(i,:,j)=grid%m_tsoil(its+1,:,jts+1) !!! grid%m_isoil(i,:,j)=grid%m_isoil(its+1,:,jts+1) grid%m_dsoil(i,:,j)=grid%m_dsoil(its+1,:,jts+1) !! T.Michaels trick to break symmetry CALL RANDOM_NUMBER(randnum) grid%em_t_gc(i,1,j)=grid%em_t_gc(its+1,1,jts+1) + 0.1*2.*(0.5-randnum) CALL RANDOM_NUMBER(randnum) grid%em_t_gc(i,2,j)=grid%em_t_gc(its+1,2,jts+1) + 0.1*2.*(0.5-randnum) !CALL RANDOM_NUMBER(randnum) !grid%em_t_gc(i,3,j)=grid%em_t_gc(its+1,3,jts+1) + 0.1*2.*(0.5-randnum) !CALL RANDOM_NUMBER(randnum) !grid%em_t_gc(i,4,j)=grid%em_t_gc(its+1,4,jts+1) + 0.1*2.*(0.5-randnum) !CALL RANDOM_NUMBER(randnum) !grid%em_t_gc(i,5,j)=grid%em_t_gc(its+1,5,jts+1) + 0.1*2.*(0.5-randnum) ENDDO ENDDO ENDIF IF (config_flags%init_U .ne. 0.) THEN grid%em_u_gc = grid%em_u_gc*0. + config_flags%init_U print *, 'constant zonal wind ', config_flags%init_U ENDIF IF (config_flags%init_V .ne. 0.) THEN grid%em_v_gc = grid%em_v_gc*0. + config_flags%init_V print *, 'constant meridional wind ', config_flags%init_V ENDIF !!!!!!!!!!!!!!!!!!! !!! READ PROFILE !! !!!!!!!!!!!!!!!!!!! ! !open(unit=10,file='input_sounding',form='formatted',status='old') !rewind(10) !read(10,*) grid%em_u_gc(1,:,1) !!****MARS !! !! case with idealized topography !! !!CALL ideal_topo ( grid%ht , 2000., 6., & !CALL ideal_topo ( grid%ht , 2000., 3., & ! ids , ide , jds , jde , kds , kde , & ! ims , ime , jms , jme , kms , kme , & ! its , ite , jts , jte , kts , kte ) !!****MARS ! If we have any input low-res surface pressure, we store it. !!****MARS !!fix pour être certain d'être avec les bons flag print *,flag_psfc print *,flag_soilhgt print *,flag_metgrid flag_psfc=1 flag_soilhgt=1 flag_metgrid=1 !!**** TODO: trouver quand même pourquoi donne 0 :) pa=999999. !!****MARS IF ( flag_psfc .EQ. 1 ) THEN DO j = jts, MIN(jte,jde-1) DO i = its, MIN(ite,ide-1) grid%em_psfc_gc(i,j) = grid%psfc(i,j) !!!****MARS: em_p_gc is only a way to count vertical levels in WPS :) !!!****MARS: is filled here with real pressure levels grid%em_p_gc(i,:,j) = grid%em_rh_gc(i,:,j) !!!****MARS grid%em_p_gc(i,1,j) = grid%psfc(i,j) !!!****MARS IF (pa .gt. grid%em_p_gc(i,1,j)) pa=grid%em_p_gc(i,1,j) !!!****MARS END DO END DO END IF print *, 'found minimum pressure (Pa) :',pa !!!!****MARS !!!!****MARS !!!! define new hybrid coordinate levels !!!! with transition level between sigma and pressure !!!! lower than input data ! ! ! !--get vertical size of the GCM input array ! sizegcm=SIZE(grid%em_rh_gc(1,:,1)) ! ALLOCATE(sig(MIN(ite,ide-1)-its+1,sizegcm, MIN(jte,jde-1)-jts+1)) ! ALLOCATE(ap(MIN(ite,ide-1)-its+1,sizegcm, MIN(jte,jde-1)-jts+1)) ! ALLOCATE(bp(MIN(ite,ide-1)-its+1,sizegcm, MIN(jte,jde-1)-jts+1)) ! !ALLOCATE(box(MIN(ite,ide-1)-its+1,sizegcm, MIN(jte,jde-1)-jts+1)) ! ! ! ! !--define sigma levels, ! !--then derive new sigma levels, and new pressure levels ! DO j = jts, MIN(jte,jde-1) ! DO i = its, MIN(ite,ide-1) ! ! ! old sigma levels ! sig(i,:,j)=grid%em_p_gc(i,:,j)/grid%em_psfc_gc(i,j) !! sig(i,:,j)=grid%em_p_gc(20,:,20)/grid%em_psfc_gc(20,20) ! ! new pressure levels ! ! - pressure_new = ap_new + bp_new * ps_gcm ! ! - bp_new is converging more rapidly than bp ! ! ... while conserving the same structure near the surface ! ! ! ! NB: grid%zap_close_levels ne sert pas dans vert_interp_old :) ! ! NB: peut donc servir pour préciser une constante reelle ! ! NB: qui permet de rehausser la zone de transition ! ! ! bp(i,:,j)=sqrt(sqrt(exp(1.-1./(sig(i,:,j)**4)))) ! ap(i,:,j)=pa*exp(-grid%zap_close_levels/10.)*(sig(i,:,j)-bp(i,:,j)) ! grid%em_rh_gc(i,:,j)=ap(i,:,j)+bp(i,:,j)*grid%em_psfc_gc(i,j) ! ! ! avoid extrapolation at the top ! ! -- the last level is thus unsignificant ! grid%em_p_gc(i,sizegcm,j)=grid%em_p_gc(i,sizegcm,j)/100. !! grid%em_p_gc(i,sizegcm,j)=grid%em_p_gc(i,sizegcm,j)/10000. ! ! ENDDO ! ENDDO ! ! ! ! ! !-- check that the biggest differences are higher ! print *, 'sigma levels' ! print *, sig(its+1,:,jts+1) ! print *, 'old pressure levels' ! print *, grid%em_p_gc(its+1,:,jts+1) ! print *, 'new pressure levels' ! print *, grid%em_rh_gc(its+1,:,jts+1) ! !!print *, 't_gc', SIZE(grid%em_t_gc(1,:,1)) !!print *, 'p_gc', SIZE(grid%em_p_gc(1,:,1)) !!print *, 't_2', SIZE(grid%em_t_2(1,:,1)) !!print *, 'rh_gc', SIZE(grid%em_rh_gc(1,:,1)) ! ! !!-------- !!-- interpolate on the new pressure levels !!-------- ! ! DO j = jts, MIN(jte,jde-1) ! DO i = its, MIN(ite,ide-1) ! !DO knew = 1,sizegcm ! loop on each level of the new grid ! ! DO kold = 1,sizegcm-1 ! find the two enclosing levels ! ! ! indic becomes negative when the two levels are found ! indic=(grid%em_p_gc(i,kold,j)-grid%em_rh_gc(i,knew,j))& ! *(grid%em_p_gc(i,kold+1,j)-grid%em_rh_gc(i,knew,j)) ! ! ! 1. the two levels are found - define p1,p2,pn and exit the loop ! IF (indic < 0.) THEN ! !IF ((i == its) .AND. (j == jts)) THEN !just a check ! ! print *, 'new level', grid%em_rh_gc(i,knew,j) ! ! print *, 'interp levels', grid%em_p_gc(i,kold,j), & ! ! grid%em_p_gc(i,kold+1,j) ! !ENDIF ! p1 = ALOG(grid%em_p_gc(i,kold,j)) ! p2 = ALOG(grid%em_p_gc(i,kold+1,j)) ! pn = ALOG(grid%em_rh_gc(i,knew,j)) ! EXIT ! ! ! 2. must handle the case (usually close to the surface) ! ! of similar new/old levels - then exit with the right kold value ! ELSE IF (1-abs(grid%em_rh_gc(i,knew,j)/grid%em_p_gc(i,kold,j)) .lt. 1e-8) THEN ! !print *,grid%em_p_gc(i,kold,j),grid%em_rh_gc(i,knew,j) ! EXIT ! ELSE IF (1-abs(grid%em_rh_gc(i,knew,j)/grid%em_p_gc(i,kold+1,j)) .lt. 1e-8) THEN ! !print *,grid%em_p_gc(i,kold+1,j),grid%em_rh_gc(i,knew,j) ! kold=kold+1 ! EXIT ! ! ! 3. continue looping if the two levels are not found .... ! ENDIF ! ENDDO ! ! ! this is an initialization, useful for case 2 (and erased just below if case 1) ! grid%em_t_2(i,knew,j)= grid%em_t_gc(i,kold,j) ! grid%em_u_2(i,knew,j)= grid%em_u_gc(i,kold,j) ! grid%em_v_2(i,knew,j)= grid%em_v_gc(i,kold,j) ! ! ! case 1: OK, in the previous loop, the two levels were found, and stored in p1 and p2 ! ! ... thus interpolation can be performed ! IF (indic < 0.) THEN ! grid%em_t_2(i,knew,j)= ( grid%em_t_gc(i,kold,j) * ( p2 - pn ) + & ! grid%em_t_gc(i,kold+1,j) * ( pn - p1 ) ) / & ! ( p2 - p1 ) ! grid%em_u_2(i,knew,j)= ( grid%em_u_gc(i,kold,j) * ( p2 - pn ) + & ! grid%em_u_gc(i,kold+1,j) * ( pn - p1 ) ) / & ! ( p2 - p1 ) ! grid%em_v_2(i,knew,j)= ( grid%em_v_gc(i,kold,j) * ( p2 - pn ) + & ! grid%em_v_gc(i,kold+1,j) * ( pn - p1 ) ) / & ! ( p2 - p1 ) ! ENDIF ! ! !ENDDO ! ! ENDDO ! ENDDO !grid%em_u_2(MIN(ite,ide-1)+1,:,:)=grid%em_u_2(MIN(ite,ide-1),:,:) !grid%em_v_2(:,:,MIN(jte,jde-1)+1)=grid%em_v_2(:,:,MIN(jte,jde-1)) !!-------- !!-- end - interpolate on the new pressure levels !!-------- ! !-- interpolate on the new pressure levels ! CALL vert_interp_old ( grid%em_t_gc , & ! --- interpolate this field ! grid%em_p_gc, & ! --- with coordinates ! grid%em_t_2, & ! --- to obtain the new field ! grid%em_rh_gc, & ! --- on coordinates ! sizegcm, & ! 'T', & ! --- no staggering (will be done later) ! 2, & ! --- log p interpolation ! 1, & ! --- (0) lagrange_order ! .false., & ! --- (0) lowest_lev_from_sfc ! 0., & ! --- (0) zap_close_levels ! 0, & ! --- (0) force_sfc_in_vinterp ! ids , ide , jds , jde , kds , kde , & ! ims , ime , jms , jme , kms , kme , & ! its , ite , jts , jte , kts , kte ) ! CALL vert_interp_old ( grid%em_u_gc , & ! --- interpolate this field ! grid%em_p_gc, & ! --- with coordinates ! grid%em_u_2 , & ! --- to obtain the new field ! grid%em_rh_gc, & ! --- on coordinates ! sizegcm, & ! 'U', & ! --- no staggering (will be done later) ! 2, & ! --- log p interpolation ! 1, & ! --- (0) lagrange_order ! .false., & ! --- (0) lowest_lev_from_sfc ! 0., & ! --- (0) zap_close_levels ! 0, & ! --- (0) force_sfc_in_vinterp ! ids , ide , jds , jde , kds , kde , & ! ims , ime , jms , jme , kms , kme , & ! its , ite , jts , jte , kts , kte ) ! CALL vert_interp_old ( grid%em_v_gc , & ! --- interpolate this field ! grid%em_p_gc, & ! --- with coordinates ! grid%em_v_2 , & ! --- to obtain the new field ! grid%em_rh_gc, & ! --- on coordinates ! sizegcm, & ! 'V', & ! --- no staggering (will be done later) ! 2, & ! --- log p interpolation ! 1, & ! --- (0) lagrange_order ! .false., & ! --- (0) lowest_lev_from_sfc ! 0., & ! --- (0) zap_close_levels ! 0, & ! --- (0) force_sfc_in_vinterp ! !ids , ide , jds , jde , kds , sizegcm , & ! !ims , ime , jms , jme , kms , sizegcm , & ! !its , ite , jts , jte , kts , sizegcm ) ! ids , ide , jds , jde , kds , kde , & ! ims , ime , jms , jme , kms , kme , & ! its , ite , jts , jte , kts , kte ) ! ! ! !-- save the new field and the new pressure coordinates ! !-- these will be regarded now as the inputs from the GCM ! grid%em_t_gc=grid%em_t_2 ! grid%em_t_2(:,:,:)=0. ! grid%em_u_gc=grid%em_u_2 ! grid%em_u_2(:,:,:)=0. ! grid%em_v_gc=grid%em_v_2 ! grid%em_v_2(:,:,:)=0. ! grid%em_p_gc=grid%em_rh_gc ! grid%em_rh_gc(:,:,:)=0. !!!!****MARS !!!****MARS ! If we have the low-resolution surface elevation, stick that in the ! "input" locations of the 3d height. We still have the "hi-res" topo ! stuck in the grid%em_ht array. The grid%landmask if test is required as some sources ! have ZERO elevation over water (thank you very much). IF ( flag_soilhgt .EQ. 1) THEN DO j = jts, MIN(jte,jde-1) DO i = its, MIN(ite,ide-1) IF ( grid%landmask(i,j) .GT. 0.5 ) THEN grid%em_ght_gc(i,1,j) = grid%toposoil(i,j) grid%em_ht_gc(i,j)= grid%toposoil(i,j) END IF END DO END DO END IF ! Assign surface fields with original input values. If this is hybrid data, ! the values are not exactly representative. However - this is only for ! plotting purposes and such at the 0h of the forecast, so we are not all that ! worried. !****MARS ! DO j = jts, min(jde-1,jte) ! DO i = its, min(ide,ite) ! grid%u10(i,j)=grid%em_u_gc(i,1,j) ! END DO ! END DO ! ! DO j = jts, min(jde,jte) ! DO i = its, min(ide-1,ite) ! grid%v10(i,j)=grid%em_v_gc(i,1,j) ! END DO ! END DO !****MARS ! DO j = jts, min(jde-1,jte) ! DO i = its, min(ide-1,ite) ! grid%t2(i,j)=grid%em_t_gc(i,1,j) ! END DO ! END DO ! The number of vertical levels in the input data. There is no staggering for ! different variables. num_metgrid_levels = grid%num_metgrid_levels ! The requested ptop for real data cases. p_top_requested = grid%p_top_requested ! Compute the top pressure, grid%p_top. For isobaric data, this is just the ! top level. For the generalized vertical coordinate data, we find the ! max pressure on the top level. We have to be careful of two things: ! 1) the value has to be communicated, 2) the value can not increase ! at subsequent times from the initial value. IF ( internal_time_loop .EQ. 1 ) THEN CALL find_p_top ( grid%em_p_gc , grid%p_top , & ids , ide , jds , jde , 1 , num_metgrid_levels , & ims , ime , jms , jme , 1 , num_metgrid_levels , & its , ite , jts , jte , 1 , num_metgrid_levels ) !! ^---- equivalent to: !!grid%ptop=MINVAL(grid%em_p_gc(:,:,:)) !!!!obsolete !print *,'ptop GCM',grid%em_rh_gc(2,1,2) !IF (grid%em_rh_gc(2,1,2) == 0) THEN ! print *,'ptop cannot be 0' ! stop !ENDIF !grid%p_top=grid%em_rh_gc(2,1,2) !!!!obsolete #ifdef DM_PARALLEL grid%p_top = wrf_dm_max_real ( grid%p_top ) #endif ! Compare the requested grid%p_top with the value available from the input data. print *,'p_top_requested = ',p_top_requested print *,'allowable grid%p_top in data = ',grid%p_top IF ( p_top_requested .LT. grid%p_top ) THEN CALL wrf_error_fatal ( 'p_top_requested < grid%p_top possible from data' ) END IF ! The grid%p_top valus is the max of what is available from the data and the ! requested value. We have already compared <, so grid%p_top is directly set to ! the value in the namelist. grid%p_top = p_top_requested ! For subsequent times, we have to remember what the grid%p_top for the first ! time was. Why? If we have a generalized vert coordinate, the grid%p_top value ! could fluctuate. p_top_save = grid%p_top ELSE CALL find_p_top ( grid%em_p_gc , grid%p_top , & ids , ide , jds , jde , 1 , num_metgrid_levels , & ims , ime , jms , jme , 1 , num_metgrid_levels , & its , ite , jts , jte , 1 , num_metgrid_levels ) #ifdef DM_PARALLEL grid%p_top = wrf_dm_max_real ( grid%p_top ) #endif IF ( grid%p_top .GT. p_top_save ) THEN print *,'grid%p_top from last time period = ',p_top_save print *,'grid%p_top from this time period = ',grid%p_top CALL wrf_error_fatal ( 'grid%p_top > previous value' ) END IF grid%p_top = p_top_save ENDIF !****MARS !****MARS print *,'ptop GCM', grid%p_top print *,'sample: pressure at its jts' print *,grid%em_p_gc(its,:,jts) !****MARS !****MARS !****MARS: useless !****MARS: ! ! Get the monthly values interpolated to the current date for the traditional monthly ! ! fields of green-ness fraction and background albedo. ! ! CALL monthly_interp_to_date ( grid%em_greenfrac , 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%em_albedo12m , 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%em_greenfrac , 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. ! 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) ! 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 ! END DO ! END DO !!****MARS: !!****MARS: useless !!****MARS: ! ! Compute the mixing ratio from the input relative humidity. ! ! IF ( flag_qv .NE. 1 ) THEN ! CALL rh_to_mxrat (grid%em_rh_gc, grid%em_t_gc, grid%em_p_gc, grid%em_qv_gc , .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 ) ! END IF !!****MARS: !!grid%em_rh_gc are GCM equivalent eta_levels !!****MARS ! Two ways to get the surface pressure. 1) If we have the low-res input surface ! pressure and the low-res topography, then we can do a simple hydrostatic ! relation. 2) Otherwise we compute the surface pressure from the sea-level ! pressure. ! Note that on output, grid%em_psfc is now hi-res. The low-res surface pressure and ! elevation are grid%em_psfc_gc and grid%em_ht_gc (same as grid%em_ght_gc(k=1)). !!****MARS: switch off this option !!****MARS: --> cf sfcprs2 and geopotential function at 500mb ! IF ( config_flags%adjust_heights ) THEN ! we_have_tavgsfc = ( flag_tavgsfc == 1 ) ! ELSE ! we_have_tavgsfc = .FALSE. ! END IF !****MARS: we_have_tavgsfc = .FALSE. !****MARS: hi-res psfc is done if the flag 'sfcp_to_sfcp' is active (recommended) IF ( ( flag_psfc .EQ. 1 ) .AND. ( flag_soilhgt .EQ. 1 ) .AND. & ( config_flags%sfcp_to_sfcp ) ) THEN print *,'compute psfc from hi-res topography' CALL sfcprs2(grid%em_t_gc, grid%em_qv_gc, grid%em_ght_gc, grid%em_psfc_gc, grid%ht, & grid%em_tavgsfc, grid%em_p_gc, grid%psfc, we_have_tavgsfc, & ids , ide , jds , jde , 1 , num_metgrid_levels , & ims , ime , jms , jme , 1 , num_metgrid_levels , & its , ite , jts , jte , 1 , num_metgrid_levels ) !****MARS: here, in reality, grid%em_p_gc is not used !****MARS: no sea-level pressure inputs possible ! ELSE ! CALL sfcprs (grid%em_t_gc, grid%em_qv_gc, grid%em_ght_gc, grid%em_pslv_gc, grid%ht, & ! grid%em_tavgsfc, grid%em_p_gc, grid%psfc, we_have_tavgsfc, & ! ids , ide , jds , jde , 1 , num_metgrid_levels , & ! ims , ime , jms , jme , 1 , num_metgrid_levels , & ! its , ite , jts , jte , 1 , num_metgrid_levels ) !****MARS: no sea-level pressure inputs possible ! If we have no input surface pressure, we'd better stick something in there. IF ( flag_psfc .NE. 1 ) THEN DO j = jts, MIN(jte,jde-1) DO i = its, MIN(ite,ide-1) grid%em_psfc_gc(i,j) = grid%psfc(i,j) grid%em_p_gc(i,1,j) = grid%psfc(i,j) END DO END DO END IF END IF !!!****MARS: !!!****MARS: old stuff !!! grid%em_p_gc is needed ... so it is computed from eta_gcm ! !print *,'computing pressure levels for input data...' ! ! !! pressure is computed from eta_gcm and hi-res topography ! DO j = jts, MIN(jte,jde-1) ! DO i = its, MIN(ite,ide-1) !!!psfc ou em_psfc_gc ? em_psfc_gc, sinon c'est faux et déclenche instabilités !grid%em_p_gc(i,:,j)=grid%em_rh_gc(i,:,j)*(grid%em_psfc_gc(i,j)-grid%em_rh_gc(2,1,2))+grid%em_rh_gc(2,1,2) !grid%em_p_gc(i,1,j)=grid%em_psfc_gc(i,j) ! ! ! END DO ! END DO !! !!****MARS: !! Integrate the mixing ratio to get the vapor pressure. ! !CALL integ_moist ( grid%em_qv_gc , grid%em_p_gc , grid%em_pd_gc , grid%em_t_gc , grid%em_ght_gc , grid%em_intq_gc , & ! ids , ide , jds , jde , 1 , num_metgrid_levels , & ! ims , ime , jms , jme , 1 , num_metgrid_levels , & ! its , ite , jts , jte , 1 , num_metgrid_levels ) !!****MARS !!****MARS !! and now, convert the GCM sigma levels into WRF sigma levels using hi-res surface pressure !!DO j = jts , MIN ( jde-1 , jte ) !!DO i = its , MIN (ide-1 , ite ) !! !! grid%em_pd_gc(i,:,j)=ap(i,:,j)+bp(i,:,j)*grid%psfc(i,j) !! !!END DO !!END DO IF ( planet == "mars" ) then !--get vertical size of the GCM input array and allocate new stuff sizegcm=SIZE(grid%em_rh_gc(1,:,1)) ALLOCATE(sig(MIN(ite,ide-1)-its+1,sizegcm, MIN(jte,jde-1)-jts+1)) !ALLOCATE(ap(MIN(ite,ide-1)-its+1,sizegcm, MIN(jte,jde-1)-jts+1)) ALLOCATE(bp(MIN(ite,ide-1)-its+1,sizegcm, MIN(jte,jde-1)-jts+1)) DO j = jts , MIN ( jde-1 , jte ) DO i = its , MIN (ide-1 , ite ) !!! Define old sigma levels for each column sig(i,:,j)=grid%em_p_gc(i,:,j)/grid%em_psfc_gc(i,j) !!! Compute new sigma levels from old sigma levels with GCM (low-res) and WRF (hi-res) surface pressure !!! (dimlevs,sigma_gcm, ps_gcm, ps_hr, sigma_hr) CALL build_sigma_hr(sizegcm,sig(i,:,j),grid%em_psfc_gc(i,j),grid%psfc(i,j),bp(i,:,j)) !!! Calculate new pressure levels grid%em_pd_gc(i,:,j)=bp(i,:,j)*grid%psfc(i,j) END DO END DO DEALLOCATE(sig) DEALLOCATE(bp) !!****MARS who knows... grid%em_rh_gc(:,:,:)=0. !!****MARS !grid%em_pd_gc=grid%em_p_gc !!****MARS ELSE ! VENUS !! Compute the difference between the dry, total surface pressure (input) and the !! dry top pressure (constant). ! CALL p_dts ( grid%em_mu0 , grid%em_intq_gc , grid%psfc , grid%p_top , & ids , ide , jds , jde , 1 , num_metgrid_levels , & ims , ime , jms , jme , 1 , num_metgrid_levels , & its , ite , jts , jte , 1 , num_metgrid_levels ) ENDIF IF ( planet == "mars" ) then !!****MARS DO j = jts , MIN ( jde-1 , jte ) DO i = its , MIN (ide-1 , ite ) grid%em_mu0(i,j) = grid%psfc(i,j) - grid%p_top END DO END DO !!****MARS ELSE ! VENUS !! Compute the dry, hydrostatic surface pressure. ! CALL p_dhs ( grid%em_pdhs , grid%ht , p00 , t00 , a , & ids , ide , jds , jde , kds , kde , & ims , ime , jms , jme , kms , kme , & its , ite , jts , jte , kts , kte ) ENDIF !!****MARS: voir remarques dans la routine !!****MARS: dry hydrostatic pressure comes from the GCM ... ! DO j = jts , MIN ( jde-1 , jte ) ! DO i = its , MIN (ide-1 , ite ) ! grid%em_pdhs(i,j) = grid%psfc(i,j) ! END DO ! END DO !!****MARS: em_pdhs ne sert qu'ici ! ! Compute the eta levels if not defined already. !!TODO: pb when ptop<1Pa IF ( grid%em_znw(1) .NE. 1.0 ) THEN eta_levels(1:kde) = model_config_rec%eta_levels(1:kde) max_dz = model_config_rec%max_dz !!****MARS IF (grid%force_sfc_in_vinterp == 0) grid%force_sfc_in_vinterp = 8 !!default choice !!****MARS CALL compute_eta ( grid%em_znw , & eta_levels , max_eta , max_dz , & grid%force_sfc_in_vinterp, & !!ne sert pas par ailleurs grid%p_top , g , p00 , cvpm , a , r_d , cp , t00 , p1000mb , t0 , & tiso, & ids , ide , jds , jde , kds , kde , & ims , ime , jms , jme , kms , kme , & its , ite , jts , jte , kts , kte ) END IF IF ( interp_theta ) THEN ! The input field is temperature, we want potential temp. !****MARS: here em_p_gc is really needed ! CALL t_to_theta ( grid%em_t_gc , grid%em_p_gc , p00 , & ids , ide , jds , jde , 1 , num_metgrid_levels , & ims , ime , jms , jme , 1 , num_metgrid_levels , & its , ite , jts , jte , 1 , num_metgrid_levels ) ENDIF ! On the eta surfaces, compute the dry pressure = mu eta, stored in ! grid%em_pb, since it is a pressure, and we don't need another kms:kme 3d ! array floating around. The grid%em_pb array is re-computed as the base pressure ! later after the vertical interpolations are complete. CALL p_dry ( grid%em_mu0 , grid%em_znw , grid%p_top , grid%em_pb , & ids , ide , jds , jde , kds , kde , & ims , ime , jms , jme , kms , kme , & its , ite , jts , jte , kts , kte ) print *, 'test sample' print *, grid%em_pb(its+10,:,jts+10) print *, 'test sample 2' print *, grid%em_pb(its,:,jts) !****MARS !****MARS: old stuff !****MARS !!! and now eta levels from the GCM are computed with the WRF ptop and GCM psfc !!! and em_pb is filled with WRF eta levels to prepare interpolation !print *,'computing eta levels for input data...' ! DO j = jts, MIN(jte,jde-1) ! DO i = its, MIN(ite,ide-1) ! !!grid%em_psfc_gc: pb en haut!!!! !!!!valeurs plus grandes que 1 et extrapolation !!grid%em_p_gc(i,:,j)=(grid%em_p_gc(i,:,j)-grid%p_top)/(grid%psfc(i,j)-grid%p_top) !!!!utile si l'on est proche de la surface, mais pb plus haut ! !grid%em_p_gc(i,:,j)=(grid%em_p_gc(i,:,j)-grid%p_top)/(grid%em_psfc_gc(i,j)-grid%p_top) !grid%em_pb(i,:,j)=grid%em_znw(:) ! !! !!!!manage negative values !!DO k=1,num_metgrid_levels !! grid%em_p_gc(i,k,j)=MAX(0.,grid%em_p_gc(i,k,j)) !!END DO !! ! ! END DO ! END DO !! !!print *,'sample: eta GCM at its jts' !!print *,grid%em_p_gc(its,:,jts) !!print *,'sample: eta WRF at its jts' !!print *,grid%em_pb(its,:,jts) !! !!print *,grid%em_p_gc(:,2,:) !!print *, 'yeah yeah' !!grid%em_pd_gc(:,:,:)=grid%em_p_gc(:,:,:) !! !****MARS !****MARS: old stuff !****MARS ! All of the vertical interpolations are done in dry-pressure space. The ! input data has had the moisture removed (grid%em_pd_gc). The target levels (grid%em_pb) ! had the vapor pressure removed from the surface pressure, then they were ! scaled by the eta levels. interp_type = grid%interp_type lagrange_order = grid%lagrange_order lowest_lev_from_sfc = grid%lowest_lev_from_sfc zap_close_levels = grid%zap_close_levels force_sfc_in_vinterp = grid%force_sfc_in_vinterp !!****MARS: normalement c'est vert_interp !!****MARS: mais les résultats sont trop discontinus > retour à une !!****MARS: interpolation plus classique CALL vert_interp_old ( grid%em_qv_gc , grid%em_pd_gc , moist(:,:,:,P_QV) , grid%em_pb , & num_metgrid_levels , 'Q' , & interp_type , lagrange_order , lowest_lev_from_sfc , & zap_close_levels , force_sfc_in_vinterp , & ids , ide , jds , jde , kds , kde , & ims , ime , jms , jme , kms , kme , & its , ite , jts , jte , kts , kte ) ! Depending on the setting of interp_theta = T/F, t_gc is is either theta Xor ! temperature, and that means that the t_2 field is also the associated field. ! It is better to interpolate temperature and potential temperature in LOG(p), ! regardless of requested default. !!****MARS: normalement c'est vert_interp CALL vert_interp_old ( grid%em_t_gc , grid%em_pd_gc , grid%em_t_2 , grid%em_pb , & num_metgrid_levels , 'T' , & interp_type , lagrange_order , lowest_lev_from_sfc , & zap_close_levels , force_sfc_in_vinterp , & ids , ide , jds , jde , kds , kde , & ims , ime , jms , jme , kms , kme , & its , ite , jts , jte , kts , kte ) IF ( .NOT. interp_theta ) THEN !! correction Wee et al. 2012 !! first interpolate temperature (see above) !! then interpolate pressure !! and in the end compute potential temperature !! scalar just an intermediate thing for interpolated pressure !! -- it is reinitialized afterwards !! It is better to interpolate pressure in p regardless default options CALL vert_interp_old ( grid%em_p_gc , grid%em_pd_gc , scalar(:,:,:,1) , grid%em_pb , & num_metgrid_levels , 'T' , & 1, lagrange_order , lowest_lev_from_sfc , & zap_close_levels , force_sfc_in_vinterp , & ids , ide , jds , jde , kds , kde , & ims , ime , jms , jme , kms , kme , & its , ite , jts , jte , kts , kte ) CALL t_to_theta ( grid%em_t_2 , scalar(:,:,:,1) , p00 , & ids , ide , jds , jde , kds , kde , & ims , ime , jms , jme , kms , kme , & its , ite , jts , jte , kts , kte ) scalar(:,:,:,1) = 0. END IF !!!!!!****MARS****!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!! ****MARS MARS MARS for order in rank for each option, check in REGISTRY.EM !!!!!! [a little bit too hardcoded here unfortunately... but e.g. P_QH2O must be known] !!!!!! [there is a more flexible option with flags in metgrid.tbl, TBD?] !!!!!! NB: real_em.F must also be modified !!!!!! NB2: qvapor is not used to avoid collision with earth-related calculations if ( ( config_flags%mars == 2 ) .OR. ( config_flags%mars == 3 ) ) then print *, '**** INTERPOLATE DUSTQ **** RANK 2 in SCALAR' CALL vert_interp_old ( grid%em_dustq_gc , grid%em_pd_gc , scalar(:,:,:,2) , grid%em_pb , & num_metgrid_levels , 'Q' , & interp_type , lagrange_order , lowest_lev_from_sfc , & zap_close_levels , force_sfc_in_vinterp , & ids , ide , jds , jde , kds , kde , & ims , ime , jms , jme , kms , kme , & its , ite , jts , jte , kts , kte ) endif if ( ( config_flags%mars == 3 ) ) then print *, '**** INTERPOLATE DUSTN **** RANK 3 in SCALAR' CALL vert_interp_old ( grid%em_dustn_gc , grid%em_pd_gc , scalar(:,:,:,3) , grid%em_pb , & num_metgrid_levels , 'Q' , & interp_type , lagrange_order , lowest_lev_from_sfc , & zap_close_levels , force_sfc_in_vinterp , & ids , ide , jds , jde , kds , kde , & ims , ime , jms , jme , kms , kme , & its , ite , jts , jte , kts , kte ) endif if ( ( config_flags%mars == 1 ) .OR. ( config_flags%mars == 11 ) .OR. ( config_flags%mars == 12 ) .OR. ( config_flags%mars == 32 ) ) then print *, '**** INTERPOLATE HV **** RANK 2 in SCALAR' !print *, size(scalar(0,0,0,:)), P_QH2O, P_QH2O_ICE CALL vert_interp_old ( grid%em_hv_gc , grid%em_pd_gc , scalar(:,:,:,2) , grid%em_pb , & num_metgrid_levels , 'Q' , & interp_type , lagrange_order , lowest_lev_from_sfc , & zap_close_levels , force_sfc_in_vinterp , & ids , ide , jds , jde , kds , kde , & ims , ime , jms , jme , kms , kme , & its , ite , jts , jte , kts , kte ) print *, '**** INTERPOLATE HI **** RANK 3 in SCALAR' CALL vert_interp_old ( grid%em_hi_gc , grid%em_pd_gc , scalar(:,:,:,3) , grid%em_pb , & num_metgrid_levels , 'Q' , & interp_type , lagrange_order , lowest_lev_from_sfc , & zap_close_levels , force_sfc_in_vinterp , & ids , ide , jds , jde , kds , kde , & ims , ime , jms , jme , kms , kme , & its , ite , jts , jte , kts , kte ) endif #ifdef NEWPHYS if ( (config_flags%mars == 10) ) then print *, '**** INTERPOLATE CO2 **** RANK 2 in SCALAR' CALL vert_interp_old ( grid%em_co2_gc , grid%em_pd_gc , scalar(:,:,:,2) , grid%em_pb , & num_metgrid_levels , 'Q' , & interp_type , lagrange_order , lowest_lev_from_sfc , & zap_close_levels , force_sfc_in_vinterp , & ids , ide , jds , jde , kds , kde , & ims , ime , jms , jme , kms , kme , & its , ite , jts , jte , kts , kte ) endif if ( (config_flags%mars == 11) .OR. (config_flags%mars == 12) .OR. (config_flags%mars == 32) ) then print *, '**** INTERPOLATE DUSTQ **** RANK 4 in SCALAR' CALL vert_interp_old ( grid%em_dustq_gc , grid%em_pd_gc , scalar(:,:,:,4) , grid%em_pb , & num_metgrid_levels , 'Q' , & interp_type , lagrange_order , lowest_lev_from_sfc , & zap_close_levels , force_sfc_in_vinterp , & ids , ide , jds , jde , kds , kde , & ims , ime , jms , jme , kms , kme , & its , ite , jts , jte , kts , kte ) print *, '**** INTERPOLATE DUSTN **** RANK 5 in SCALAR' CALL vert_interp_old ( grid%em_dustn_gc , grid%em_pd_gc , scalar(:,:,:,5) , grid%em_pb , & num_metgrid_levels , 'Q' , & interp_type , lagrange_order , lowest_lev_from_sfc , & zap_close_levels , force_sfc_in_vinterp , & ids , ide , jds , jde , kds , kde , & ims , ime , jms , jme , kms , kme , & its , ite , jts , jte , kts , kte ) endif if ( (config_flags%mars == 12) .OR. (config_flags%mars == 32) ) then print *, '**** INTERPOLATE CCNQ **** RANK 6 in SCALAR' CALL vert_interp_old ( grid%em_ccnq_gc , grid%em_pd_gc , scalar(:,:,:,6) , grid%em_pb , & num_metgrid_levels , 'Q' , & interp_type , lagrange_order , lowest_lev_from_sfc , & zap_close_levels , force_sfc_in_vinterp , & ids , ide , jds , jde , kds , kde , & ims , ime , jms , jme , kms , kme , & its , ite , jts , jte , kts , kte ) print *, '**** INTERPOLATE CCNN **** RANK 7 in SCALAR' CALL vert_interp_old ( grid%em_ccnn_gc , grid%em_pd_gc , scalar(:,:,:,7) , grid%em_pb , & num_metgrid_levels , 'Q' , & interp_type , lagrange_order , lowest_lev_from_sfc , & zap_close_levels , force_sfc_in_vinterp , & ids , ide , jds , jde , kds , kde , & ims , ime , jms , jme , kms , kme , & its , ite , jts , jte , kts , kte ) endif if ( (config_flags%mars == 32) ) then print *, '**** INTERPOLATE CO2 **** RANK 8 in SCALAR' CALL vert_interp_old ( grid%em_co2_gc , grid%em_pd_gc , scalar(:,:,:,8) , grid%em_pb , & num_metgrid_levels , 'Q' , & interp_type , lagrange_order , lowest_lev_from_sfc , & zap_close_levels , force_sfc_in_vinterp , & ids , ide , jds , jde , kds , kde , & ims , ime , jms , jme , kms , kme , & its , ite , jts , jte , kts , kte ) print *, '**** set other CO2 tracers to 0' scalar(:,:,:,9) = 0. scalar(:,:,:,10) = 0. scalar(:,:,:,11) = 0. endif #endif !#ifdef NEWPHYS !!VENUS photochemistry !if ( config_flags%mars == 34 ) then ! print*,'grid%em_qco2_gc',grid%em_qco2_gc(0,:,0) ! CALL vert_interp_old ( grid%em_qco2_gc , grid%em_pd_gc , scalar(:,:,:,2) , grid%em_pb , & ! num_metgrid_levels , 'Q' , & ! interp_type , lagrange_order , lowest_lev_from_sfc , & ! zap_close_levels , force_sfc_in_vinterp , & ! ids , ide , jds , jde , kds , kde , & ! ims , ime , jms , jme , kms , kme , & ! its , ite , jts , jte , kts , kte ) ! ! CALL vert_interp_old ( grid%em_qco_gc , grid%em_pd_gc , scalar(:,:,:,3), grid%em_pb , & ! num_metgrid_levels , 'Q' , & ! interp_type , lagrange_order ,lowest_lev_from_sfc , & ! zap_close_levels , force_sfc_in_vinterp , & ! ids , ide , jds , jde , kds , kde , & ! ims , ime , jms , jme , kms , kme , & ! its , ite , jts , jte , kts , kte ) ! ! CALL vert_interp_old ( grid%em_qh2_gc , grid%em_pd_gc , scalar(:,:,:,4),grid%em_pb , & ! num_metgrid_levels , 'Q' , & ! interp_type , lagrange_order,lowest_lev_from_sfc , & ! zap_close_levels , force_sfc_in_vinterp , & ! ids , ide , jds , jde , kds , kde , & ! ims , ime , jms , jme , kms , kme , & ! its , ite , jts , jte , kts , kte ) ! ! CALL vert_interp_old ( grid%em_qh2o_gc , grid%em_pd_gc , scalar(:,:,:,5),grid%em_pb , & ! num_metgrid_levels , 'Q' , & ! interp_type , lagrange_order,lowest_lev_from_sfc , & ! zap_close_levels , force_sfc_in_vinterp , & ! ids , ide , jds , jde , kds , kde , & ! ims , ime , jms , jme , kms , kme , & ! its , ite , jts , jte , kts , kte ) ! ! CALL vert_interp_old ( grid%em_qo1d_gc , grid%em_pd_gc , scalar(:,:,:,6),grid%em_pb , & ! num_metgrid_levels , 'Q' , & ! interp_type , lagrange_order,lowest_lev_from_sfc , & ! zap_close_levels , force_sfc_in_vinterp , & ! ids , ide , jds , jde , kds , kde , & ! ims , ime , jms , jme , kms , kme , & ! its , ite , jts , jte , kts , kte ) ! ! CALL vert_interp_old ( grid%em_qo_gc , grid%em_pd_gc , scalar(:,:,:,7),grid%em_pb , & ! num_metgrid_levels , 'Q' , & ! interp_type , lagrange_order,lowest_lev_from_sfc , & ! zap_close_levels , force_sfc_in_vinterp , & ! ids , ide , jds , jde , kds , kde , & ! ims , ime , jms , jme , kms , kme , & ! its , ite , jts , jte , kts , kte ) ! ! CALL vert_interp_old ( grid%em_qo2_gc , grid%em_pd_gc , scalar(:,:,:,8),grid%em_pb , & ! num_metgrid_levels , 'Q' , & ! interp_type , lagrange_order,lowest_lev_from_sfc , & ! zap_close_levels , force_sfc_in_vinterp , & ! ids , ide , jds , jde , kds , kde , & ! ims , ime , jms , jme , kms , kme , & ! its , ite , jts , jte , kts , kte ) ! ! CALL vert_interp_old ( grid%em_qo2dg_gc , grid%em_pd_gc , scalar(:,:,:,9),grid%em_pb , & ! num_metgrid_levels , 'Q' , & ! interp_type , lagrange_order,lowest_lev_from_sfc , & ! zap_close_levels , force_sfc_in_vinterp , & ! ids , ide , jds , jde , kds , kde , & ! ims , ime , jms , jme , kms , kme , & ! its , ite , jts , jte , kts , kte ) ! ! CALL vert_interp_old ( grid%em_qo3_gc , grid%em_pd_gc , scalar(:,:,:,10),grid%em_pb , & ! num_metgrid_levels , 'Q' , & ! interp_type , lagrange_order,lowest_lev_from_sfc , & ! zap_close_levels , force_sfc_in_vinterp , & ! ids , ide , jds , jde , kds , kde , & ! ims , ime , jms , jme , kms , kme , & ! its , ite , jts , jte , kts , kte ) ! ! CALL vert_interp_old ( grid%em_qh_gc , grid%em_pd_gc , scalar(:,:,:,11),grid%em_pb , & ! num_metgrid_levels , 'Q' , & ! interp_type , lagrange_order,lowest_lev_from_sfc , & ! zap_close_levels , force_sfc_in_vinterp , & ! ids , ide , jds , jde , kds , kde , & ! ims , ime , jms , jme , kms , kme , & ! its , ite , jts , jte , kts , kte ) ! ! CALL vert_interp_old ( grid%em_qoh_gc , grid%em_pd_gc , scalar(:,:,:,12), grid%em_pb , & ! num_metgrid_levels , 'Q' , & ! interp_type , lagrange_order,lowest_lev_from_sfc , & ! zap_close_levels , force_sfc_in_vinterp , & ! ids , ide , jds , jde , kds , kde , & ! ims , ime , jms , jme , kms , kme , & ! its , ite , jts , jte , kts , kte ) ! ! CALL vert_interp_old ( grid%em_qho2_gc , grid%em_pd_gc , scalar(:,:,:,13),grid%em_pb , & ! num_metgrid_levels , 'Q' , & ! interp_type , lagrange_order,lowest_lev_from_sfc , & ! zap_close_levels , force_sfc_in_vinterp , & ! ids , ide , jds , jde , kds , kde , & ! ims , ime , jms , jme , kms , kme , & ! its , ite , jts , jte , kts , kte ) ! ! CALL vert_interp_old ( grid%em_qh2o2_gc , grid%em_pd_gc , scalar(:,:,:,14),grid%em_pb , & ! num_metgrid_levels , 'Q' , & ! interp_type , lagrange_order,lowest_lev_from_sfc , & ! zap_close_levels , force_sfc_in_vinterp , & ! ids , ide , jds , jde , kds , kde , & ! ims , ime , jms , jme , kms , kme , & ! its , ite , jts , jte , kts , kte ) ! ! CALL vert_interp_old ( grid%em_qcl_gc , grid%em_pd_gc , scalar(:,:,:,15),grid%em_pb , & ! num_metgrid_levels , 'Q' , & ! interp_type , lagrange_order ,lowest_lev_from_sfc , & ! zap_close_levels , force_sfc_in_vinterp , & ! ids , ide , jds , jde , kds , kde , & ! ims , ime , jms , jme , kms , kme , & ! its , ite , jts , jte , kts , kte ) ! ! CALL vert_interp_old ( grid%em_qclo_gc , grid%em_pd_gc , scalar(:,:,:,16),grid%em_pb , & ! num_metgrid_levels , 'Q' , & ! interp_type , lagrange_order,lowest_lev_from_sfc , & ! zap_close_levels , force_sfc_in_vinterp , & ! ids , ide , jds , jde , kds , kde , & ! ims , ime , jms , jme , kms , kme , & ! its , ite , jts , jte , kts , kte ) ! ! CALL vert_interp_old ( grid%em_qcl2_gc , grid%em_pd_gc , scalar(:,:,:,17),grid%em_pb , & ! num_metgrid_levels , 'Q' , & ! interp_type , lagrange_order,lowest_lev_from_sfc , & ! zap_close_levels , force_sfc_in_vinterp , & ! ids , ide , jds , jde , kds , kde , & ! ims , ime , jms , jme , kms , kme , & ! its , ite , jts , jte , kts , kte ) ! ! CALL vert_interp_old ( grid%em_qhcl_gc , grid%em_pd_gc , scalar(:,:,:,18),grid%em_pb , & ! num_metgrid_levels , 'Q' , & ! interp_type , lagrange_order,lowest_lev_from_sfc , & ! zap_close_levels , force_sfc_in_vinterp , & ! ids , ide , jds , jde , kds , kde , & ! ims , ime , jms , jme , kms , kme , & ! its , ite , jts , jte , kts , kte ) ! ! CALL vert_interp_old ( grid%em_qhocl_gc , grid%em_pd_gc , scalar(:,:,:,19),grid%em_pb , & ! num_metgrid_levels , 'Q' , & ! interp_type , lagrange_order,lowest_lev_from_sfc , & ! zap_close_levels , force_sfc_in_vinterp , & ! ids , ide , jds , jde , kds , kde , & ! ims , ime , jms , jme , kms , kme , & ! its , ite , jts , jte , kts , kte ) ! ! CALL vert_interp_old ( grid%em_qclco_gc , grid%em_pd_gc , scalar(:,:,:,20),grid%em_pb , & ! num_metgrid_levels , 'Q' , & ! interp_type , lagrange_order,lowest_lev_from_sfc , & ! zap_close_levels , force_sfc_in_vinterp , & ! ids , ide , jds , jde , kds , kde , & ! ims , ime , jms , jme , kms , kme , & ! its , ite , jts , jte , kts , kte ) ! ! CALL vert_interp_old ( grid%em_qclco3_gc , grid%em_pd_gc , scalar(:,:,:,21),grid%em_pb , & ! num_metgrid_levels , 'Q' , & ! interp_type , lagrange_order,lowest_lev_from_sfc , & ! zap_close_levels , force_sfc_in_vinterp , & ! ids , ide , jds , jde , kds , kde , & ! ims , ime , jms , jme , kms , kme , & ! its , ite , jts , jte , kts , kte ) ! ! CALL vert_interp_old ( grid%em_qcocl2_gc , grid%em_pd_gc , scalar(:,:,:,22),grid%em_pb , & ! num_metgrid_levels , 'Q' , & ! interp_type , lagrange_order,lowest_lev_from_sfc , & ! zap_close_levels , force_sfc_in_vinterp , & ! ids , ide , jds , jde , kds , kde , & ! ims , ime , jms , jme , kms , kme , & ! its , ite , jts , jte , kts , kte ) ! ! CALL vert_interp_old ( grid%em_qs_gc , grid%em_pd_gc , scalar(:,:,:,23), grid%em_pb , & ! num_metgrid_levels , 'Q' , & ! interp_type , lagrange_order ,lowest_lev_from_sfc , & ! zap_close_levels , force_sfc_in_vinterp , & ! ids , ide , jds , jde , kds , kde , & ! ims , ime , jms , jme , kms , kme , & ! its , ite , jts , jte , kts , kte ) ! ! CALL vert_interp_old ( grid%em_qso_gc , grid%em_pd_gc , scalar(:,:,:,24),grid%em_pb , & ! num_metgrid_levels , 'Q' , & ! interp_type , lagrange_order,lowest_lev_from_sfc , & ! zap_close_levels , force_sfc_in_vinterp , & ! ids , ide , jds , jde , kds , kde , & ! ims , ime , jms , jme , kms , kme , & ! its , ite , jts , jte , kts , kte ) ! ! CALL vert_interp_old ( grid%em_qso2_gc , grid%em_pd_gc , scalar(:,:,:,25), grid%em_pb , & ! num_metgrid_levels , 'Q' , & ! interp_type , lagrange_order ,lowest_lev_from_sfc , & ! zap_close_levels , force_sfc_in_vinterp , & ! ids , ide , jds , jde , kds , kde , & ! ims , ime , jms , jme , kms , kme , & ! its , ite , jts , jte , kts , kte ) ! ! CALL vert_interp_old ( grid%em_qo3_gc , grid%em_pd_gc , scalar(:,:,:,26),grid%em_pb , & ! num_metgrid_levels , 'Q' , & ! interp_type , lagrange_order ,lowest_lev_from_sfc , & ! zap_close_levels , force_sfc_in_vinterp , & ! ids , ide , jds , jde , kds , kde , & ! ims , ime , jms , jme , kms , kme , & ! its , ite , jts , jte , kts , kte ) ! ! CALL vert_interp_old ( grid%em_qs2o2_gc , grid%em_pd_gc , scalar(:,:,:,27),grid%em_pb , & ! num_metgrid_levels , 'Q' , & ! interp_type , lagrange_order,lowest_lev_from_sfc , & ! zap_close_levels , force_sfc_in_vinterp , & ! ids , ide , jds , jde , kds , kde , & ! ims , ime , jms , jme , kms , kme , & ! its , ite , jts , jte , kts , kte ) ! ! CALL vert_interp_old ( grid%em_qocs_gc , grid%em_pd_gc , scalar(:,:,:,28),grid%em_pb , & ! num_metgrid_levels , 'Q' , & ! interp_type , lagrange_order,lowest_lev_from_sfc , & ! zap_close_levels , force_sfc_in_vinterp , & ! ids , ide , jds , jde , kds , kde , & ! ims , ime , jms , jme , kms , kme , & ! its , ite , jts , jte , kts , kte ) ! ! CALL vert_interp_old ( grid%em_qhso3_gc , grid%em_pd_gc , scalar(:,:,:,29),grid%em_pb , & ! num_metgrid_levels , 'Q' , & ! interp_type , lagrange_order,lowest_lev_from_sfc , & ! zap_close_levels , force_sfc_in_vinterp , & ! ids , ide , jds , jde , kds , kde , & ! ims , ime , jms , jme , kms , kme , & ! its , ite , jts , jte , kts , kte ) ! ! CALL vert_interp_old ( grid%em_qh2so4_gc , grid%em_pd_gc , scalar(:,:,:,30),grid%em_pb , & ! num_metgrid_levels , 'Q' , & ! interp_type , lagrange_order,lowest_lev_from_sfc , & ! zap_close_levels , force_sfc_in_vinterp , & ! ids , ide , jds , jde , kds , kde , & ! ims , ime , jms , jme , kms , kme , & ! its , ite , jts , jte , kts , kte ) ! ! CALL vert_interp_old ( grid%em_qs2_gc , grid%em_pd_gc , scalar(:,:,:,31),grid%em_pb , & ! num_metgrid_levels , 'Q' , & ! interp_type , lagrange_order,lowest_lev_from_sfc , & ! zap_close_levels , force_sfc_in_vinterp , & ! ids , ide , jds , jde , kds , kde , & ! ims , ime , jms , jme , kms , kme , & ! its , ite , jts , jte , kts , kte ) ! ! CALL vert_interp_old ( grid%em_qclso2_gc , grid%em_pd_gc , scalar(:,:,:,32),grid%em_pb , & ! num_metgrid_levels , 'Q' , & ! interp_type , lagrange_order,lowest_lev_from_sfc , & ! zap_close_levels , force_sfc_in_vinterp , & ! ids , ide , jds , jde , kds , kde , & ! ims , ime , jms , jme , kms , kme , & ! its , ite , jts , jte , kts , kte ) ! ! CALL vert_interp_old ( grid%em_qoscl_gc , grid%em_pd_gc , scalar(:,:,:,33),grid%em_pb , & ! num_metgrid_levels , 'Q' , & ! interp_type , lagrange_order,lowest_lev_from_sfc , & ! zap_close_levels , force_sfc_in_vinterp , & ! ids , ide , jds , jde , kds , kde , & ! ims , ime , jms , jme , kms , kme , & ! its , ite , jts , jte , kts , kte ) ! ! CALL vert_interp_old ( grid%em_qh2oliq_gc , grid%em_pd_gc , scalar(:,:,:,34),grid%em_pb , & ! num_metgrid_levels , 'Q' , & ! interp_type , lagrange_order,lowest_lev_from_sfc , & ! zap_close_levels , force_sfc_in_vinterp , & ! ids , ide , jds , jde , kds , kde , & ! ims , ime , jms , jme , kms , kme , & ! its , ite , jts , jte , kts , kte ) ! ! CALL vert_interp_old ( grid%em_qh2so4liq_gc , grid%em_pd_gc , scalar(:,:,:,35),grid%em_pb , & ! num_metgrid_levels , 'Q' , & ! interp_type , lagrange_order,lowest_lev_from_sfc , & ! zap_close_levels , force_sfc_in_vinterp , & ! ids , ide , jds , jde , kds , kde , & ! ims , ime , jms , jme , kms , kme , & ! its , ite , jts , jte , kts , kte ) ! !endif !#endif !!! we want any scalar (i.e. tracer) to be positive !!! and because of interpolation it is possible that negative values occur... WHERE( scalar < 0. ) scalar = 0. !!!!!!****MARS****!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! #if 0 ! Uncomment the Registry entries to activate these. This adds ! noticeably to the allocated space for the model. IF ( flag_qr .EQ. 1 ) THEN DO im = PARAM_FIRST_SCALAR, num_3d_m IF ( im .EQ. P_QR ) THEN CALL vert_interp_old ( qr_gc , grid%em_pd_gc , moist(:,:,:,P_QR) , grid%em_pb , & num_metgrid_levels , 'Q' , & interp_type , lagrange_order , lowest_lev_from_sfc , & zap_close_levels , force_sfc_in_vinterp , & ids , ide , jds , jde , kds , kde , & ims , ime , jms , jme , kms , kme , & its , ite , jts , jte , kts , kte ) END IF END DO END IF IF ( flag_qc .EQ. 1 ) THEN DO im = PARAM_FIRST_SCALAR, num_3d_m IF ( im .EQ. P_QC ) THEN CALL vert_interp_old ( qc_gc , grid%em_pd_gc , moist(:,:,:,P_QC) , grid%em_pb , & num_metgrid_levels , 'Q' , & interp_type , lagrange_order , lowest_lev_from_sfc , & zap_close_levels , force_sfc_in_vinterp , & ids , ide , jds , jde , kds , kde , & ims , ime , jms , jme , kms , kme , & its , ite , jts , jte , kts , kte ) END IF END DO END IF IF ( flag_qi .EQ. 1 ) THEN DO im = PARAM_FIRST_SCALAR, num_3d_m IF ( im .EQ. P_QI ) THEN CALL vert_interp_old ( qi_gc , grid%em_pd_gc , moist(:,:,:,P_QI) , grid%em_pb , & num_metgrid_levels , 'Q' , & interp_type , lagrange_order , lowest_lev_from_sfc , & zap_close_levels , force_sfc_in_vinterp , & ids , ide , jds , jde , kds , kde , & ims , ime , jms , jme , kms , kme , & its , ite , jts , jte , kts , kte ) END IF END DO END IF IF ( flag_qs .EQ. 1 ) THEN DO im = PARAM_FIRST_SCALAR, num_3d_m IF ( im .EQ. P_QS ) THEN CALL vert_interp_old ( qs_gc , grid%em_pd_gc , moist(:,:,:,P_QS) , grid%em_pb , & num_metgrid_levels , 'Q' , & interp_type , lagrange_order , lowest_lev_from_sfc , & zap_close_levels , force_sfc_in_vinterp , & ids , ide , jds , jde , kds , kde , & ims , ime , jms , jme , kms , kme , & its , ite , jts , jte , kts , kte ) END IF END DO END IF IF ( flag_qg .EQ. 1 ) THEN DO im = PARAM_FIRST_SCALAR, num_3d_m IF ( im .EQ. P_QG ) THEN CALL vert_interp_old ( qg_gc , grid%em_pd_gc , moist(:,:,:,P_QG) , grid%em_pb , & num_metgrid_levels , 'Q' , & interp_type , lagrange_order , lowest_lev_from_sfc , & zap_close_levels , force_sfc_in_vinterp , & ids , ide , jds , jde , kds , kde , & ims , ime , jms , jme , kms , kme , & its , ite , jts , jte , kts , kte ) END IF END DO END IF #endif #ifdef DM_PARALLEL ips = its ; ipe = ite ; jps = jts ; jpe = jte ; kps = kts ; kpe = kte ! For the U and V vertical interpolation, we need the pressure defined ! at both the locations for the horizontal momentum, which we get by ! averaging two pressure values (i and i-1 for U, j and j-1 for V). The ! pressure field on input (grid%em_pd_gc) and the pressure of the new coordinate ! (grid%em_pb) are both communicated with an 8 stencil. # include "HALO_EM_VINTERP_UV_1.inc" #endif !!****MARS: normalement c'est vert_interp CALL vert_interp_old ( grid%em_u_gc , grid%em_pd_gc , grid%em_u_2, grid%em_pb , & num_metgrid_levels , 'U' , & interp_type , lagrange_order , lowest_lev_from_sfc , & zap_close_levels , force_sfc_in_vinterp , & ids , ide , jds , jde , kds , kde , & ims , ime , jms , jme , kms , kme , & its , ite , jts , jte , kts , kte ) !!****MARS: normalement c'est vert_interp CALL vert_interp_old ( grid%em_v_gc , grid%em_pd_gc , grid%em_v_2, grid%em_pb , & num_metgrid_levels , 'V' , & interp_type , lagrange_order , lowest_lev_from_sfc , & zap_close_levels , force_sfc_in_vinterp , & ids , ide , jds , jde , kds , kde , & ims , ime , jms , jme , kms , kme , & its , ite , jts , jte , kts , kte ) !!****MARS !!****MARS !! !! old obsolete method !! ------------------- !! !!! and now eta levels from the GCM are computed with the WRF ptop and GCM psfc !!! and em_pb is filled with WRF eta levels to prepare interpolation !print *,'computing eta levels for input data...' ! ! DO j = jts, MIN(jte,jde-1) ! DO i = its, MIN(ite,ide-1) ! !! grid%em_psfc_gc: pb en haut!!!! !!!!valeurs plus grandes que 1 et extrapolation !! grid%em_p_gc(i,:,j)=(grid%em_p_gc(i,:,j)-grid%p_top)/(grid%psfc(i,j)-grid%p_top) !!!!utile si l'on est proche de la surface, mais pb plus haut ! !grid%em_pd_gc(i,:,j)=(grid%em_p_gc(i,:,j)-grid%p_top)/(grid%em_psfc_gc(i,j)-grid%p_top) !grid%em_pb(i,:,j)=grid%em_znw(:) ! !! !!!!manage negative values !!DO k=1,num_metgrid_levels !! grid%em_p_gc(i,k,j)=MAX(0.,grid%em_p_gc(i,k,j)) !!END DO !! ! ! END DO ! END DO ! !print *,'sample: eta GCM at its jts' !print *,grid%em_pd_gc(its,:,jts) !print *,'sample: eta WRF at its jts' !print *,grid%em_pb(its,:,jts) !! !!!****MARS ! ! ! !!!!****MARS !!!! !!!!grid%force_sfc_in_vinterp ne sert pas dans vert_interp_old :) !!!!peut donc servir pour préciser le nombre de niveaux !!!!pris à partir de l'interpolation eta ! !IF (grid%force_sfc_in_vinterp .NE. 0) THEN ! ! !!!save in an array that is now unused ! !!!the previously performed pressure interpolation ! grid%em_qv_gc(:,:,:)=grid%em_t_2(:,:,:) ! ! ! !!!perform interpolation on eta levels ! print *, 'interpolate on eta levels for near-surface fields' ! CALL vert_interp_old ( grid%em_t_gc , grid%em_pd_gc , grid%em_t_2, grid%em_pb , & ! num_metgrid_levels , 'T' , & ! interp_type , lagrange_order , lowest_lev_from_sfc ,& ! zap_close_levels , force_sfc_in_vinterp , & ! ids , ide , jds , jde , kds , kde , & ! ims , ime , jms , jme , kms , kme , & ! its , ite , jts , jte , kts , kte ) ! ! !!!take the first layers from the eta interpolation ! print *, 'the first ', & ! grid%force_sfc_in_vinterp, & ! 'layers will be taken from eta interpolation' ! grid%em_qv_gc(:,1:grid%force_sfc_in_vinterp,:)=grid%em_t_2(:,1:grid%force_sfc_in_vinterp,:) ! ! !!!fix the possible little discontinuity at the limit ! !!!between the two interpolation methods ! grid%em_qv_gc(:,grid%force_sfc_in_vinterp+1,:)= & ! 0.5*(grid%em_t_2(:,grid%force_sfc_in_vinterp,:) + & !!eta interpolation below ! grid%em_qv_gc(:,grid%force_sfc_in_vinterp+2,:)) !!pressure interpolation above ! ! ! !!!assign the final result in t_2 ! grid%em_t_2(:,:,:)=grid%em_qv_gc(:,:,:) ! grid%em_qv_gc(:,:,:)=0. ! ! !ELSE ! ! !ENDIF !!****MARS !!****MARS END IF ! <----- END OF VERTICAL INTERPOLATION PART ----> !****MARS: no need ! ! Protect against bad grid%em_tsk values over water by supplying grid%sst (if it is ! ! available, and if the grid%sst is reasonable). ! ! 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%sst(i,j) .GT. 200. ) .AND. ( grid%sst(i,j) .LT. 350. ) ) THEN ! grid%tsk(i,j) = grid%sst(i,j) ! ENDIF ! END DO ! END DO ! ! ! Save the grid%em_tsk field for later use in the sea ice surface temperature ! ! for the Noah LSM scheme. ! ! DO j = jts, MIN(jte,jde-1) ! DO i = its, MIN(ite,ide-1) ! grid%tsk_save(i,j) = grid%tsk(i,j) ! END DO ! END DO ! !!****MARS: no need ! ! Take the data from the input file and store it in the variables that ! ! use the WRF naming and ordering conventions. ! ! DO j = jts, MIN(jte,jde-1) ! DO i = its, MIN(ite,ide-1) ! IF ( grid%snow(i,j) .GE. 10. ) then ! grid%snowc(i,j) = 1. ! ELSE ! grid%snowc(i,j) = 0.0 ! END IF ! END DO ! END DO ! ! ! Set flag integers for presence of snowh and soilw fields ! ! grid%ifndsnowh = flag_snowh ! IF (num_sw_levels_input .GE. 1) THEN ! grid%ifndsoilw = 1 ! ELSE ! grid%ifndsoilw = 0 ! END IF ! !****MARS: no need ! ! We require input data for the various LSM schemes. ! ! enough_data : SELECT CASE ( model_config_rec%sf_surface_physics(grid%id) ) ! ! CASE (LSMSCHEME) ! IF ( num_st_levels_input .LT. 2 ) THEN ! CALL wrf_error_fatal ( 'Not enough soil temperature data for Noah LSM scheme.') ! END IF ! ! CASE (RUCLSMSCHEME) ! IF ( num_st_levels_input .LT. 2 ) THEN ! CALL wrf_error_fatal ( 'Not enough soil temperature data for RUC LSM scheme.') ! END IF ! ! END SELECT enough_data ! ! ! For sf_surface_physics = 1, we want to use close to a 30 cm value ! ! for the bottom level of the soil temps. ! ! fix_bottom_level_for_temp : SELECT CASE ( model_config_rec%sf_surface_physics(grid%id) ) ! ! CASE (SLABSCHEME) ! IF ( flag_tavgsfc .EQ. 1 ) THEN ! DO j = jts , MIN(jde-1,jte) ! DO i = its , MIN(ide-1,ite) ! grid%tmn(i,j) = grid%em_tavgsfc(i,j) ! END DO ! END DO ! ELSE IF ( flag_st010040 .EQ. 1 ) THEN ! DO j = jts , MIN(jde-1,jte) ! DO i = its , MIN(ide-1,ite) ! grid%tmn(i,j) = grid%st010040(i,j) ! END DO ! END DO ! ELSE IF ( flag_st000010 .EQ. 1 ) THEN ! DO j = jts , MIN(jde-1,jte) ! DO i = its , MIN(ide-1,ite) ! grid%tmn(i,j) = grid%st000010(i,j) ! END DO ! END DO ! ELSE IF ( flag_soilt020 .EQ. 1 ) THEN ! DO j = jts , MIN(jde-1,jte) ! DO i = its , MIN(ide-1,ite) ! grid%tmn(i,j) = grid%soilt020(i,j) ! END DO ! END DO ! ELSE IF ( flag_st007028 .EQ. 1 ) THEN ! DO j = jts , MIN(jde-1,jte) ! DO i = its , MIN(ide-1,ite) ! grid%tmn(i,j) = grid%st007028(i,j) ! END DO ! END DO ! ELSE ! CALL wrf_debug ( 0 , 'No 10-40 cm, 0-10 cm, 7-28, or 20 cm soil temperature data for grid%em_tmn') ! CALL wrf_debug ( 0 , 'Using 1 degree static annual mean temps' ) ! END IF ! ! CASE (LSMSCHEME) ! ! CASE (RUCLSMSCHEME) ! ! END SELECT fix_bottom_level_for_temp ! ! ! Adjustments for the seaice field PRIOR to the grid%tslb computations. This is ! ! is for the 5-layer scheme. ! ! num_veg_cat = SIZE ( grid%landusef , DIM=2 ) ! num_soil_top_cat = SIZE ( grid%soilctop , DIM=2 ) ! num_soil_bot_cat = SIZE ( grid%soilcbot , DIM=2 ) ! CALL nl_get_seaice_threshold ( grid%id , grid%seaice_threshold ) ! CALL nl_get_isice ( grid%id , grid%isice ) ! CALL nl_get_iswater ( grid%id , grid%iswater ) ! CALL adjust_for_seaice_pre ( grid%xice , grid%landmask , grid%tsk , grid%ivgtyp , grid%vegcat , grid%lu_index , & ! grid%xland , grid%landusef , grid%isltyp , grid%soilcat , grid%soilctop , & ! grid%soilcbot , grid%tmn , & ! grid%seaice_threshold , & ! num_veg_cat , num_soil_top_cat , num_soil_bot_cat , & ! grid%iswater , grid%isice , & ! model_config_rec%sf_surface_physics(grid%id) , & ! ids , ide , jds , jde , kds , kde , & ! ims , ime , jms , jme , kms , kme , & ! its , ite , jts , jte , kts , kte ) ! ! ! 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 , DIM=2 ) ! num_soil_top_cat = SIZE ( grid%soilctop , DIM=2 ) ! num_soil_bot_cat = SIZE ( grid%soilcbot , DIM=2 ) ! ! 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) ) ! ! ! Make all the veg/soil parms the same so as not to confuse the developer. ! ! 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 ! ! END IF ! ! ! Land use assignment. ! ! DO j = jts, MIN(jde-1,jte) ! DO i = its, MIN(ide-1,ite) ! grid%lu_index(i,j) = grid%ivgtyp(i,j) ! IF ( grid%lu_index(i,j) .NE. model_config_rec%iswater(grid%id) ) THEN ! grid%landmask(i,j) = 1 ! grid%xland(i,j) = 1 ! ELSE ! grid%landmask(i,j) = 0 ! grid%xland(i,j) = 2 ! END IF ! END DO ! END DO ! ! ! 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_soilhgt .EQ. 1 ) THEN ! adjust_soil : SELECT CASE ( model_config_rec%sf_surface_physics(grid%id) ) ! ! CASE ( SLABSCHEME , LSMSCHEME , RUCLSMSCHEME ) ! CALL adjust_soil_temp_new ( grid%tmn , model_config_rec%sf_surface_physics(grid%id) , & ! grid%tsk , grid%ht , grid%toposoil , grid%landmask , flag_soilhgt , & ! grid%st000010 , grid%st010040 , grid%st040100 , grid%st100200 , grid%st010200 , & ! flag_st000010 , flag_st010040 , flag_st040100 , flag_st100200 , flag_st010200 , & ! grid%st000007 , grid%st007028 , grid%st028100 , grid%st100255 , & ! flag_st000007 , flag_st007028 , flag_st028100 , flag_st100255 , & ! grid%soilt000 , grid%soilt005 , grid%soilt020 , grid%soilt040 , grid%soilt160 , & ! grid%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 ) ! ! END SELECT adjust_soil ! END IF ! ! ! Fix grid%em_tmn and grid%em_tsk. ! ! fix_tsk_tmn : SELECT CASE ( model_config_rec%sf_surface_physics(grid%id) ) ! ! CASE ( SLABSCHEME , LSMSCHEME , RUCLSMSCHEME ) ! 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%sst(i,j) .GT. 240. ) .AND. ( grid%sst(i,j) .LT. 350. ) ) THEN ! grid%tmn(i,j) = grid%sst(i,j) ! grid%tsk(i,j) = grid%sst(i,j) ! ELSE IF ( grid%landmask(i,j) .LT. 0.5 ) THEN ! grid%tmn(i,j) = grid%tsk(i,j) ! END IF ! END DO ! END DO ! END SELECT fix_tsk_tmn ! ! ! Is the grid%em_tsk reasonable? ! !!**** MARS DO j = jts, MIN(jde-1,jte) DO i = its, MIN(ide-1,ite) !!grid%tsk(i,j)=200 grid%tmn(i,j)=0 grid%sst(i,j)=0 !!no use on Mars!! grid%tslb(i,:,j)=0 !!tslb is 3D field END DO END DO !!**** MARS ! IF ( internal_time_loop .NE. 1 ) THEN ! DO j = jts, MIN(jde-1,jte) ! DO i = its, MIN(ide-1,ite) ! IF ( grid%tsk(i,j) .LT. 170 .or. grid%tsk(i,j) .GT. 400. ) THEN ! grid%tsk(i,j) = grid%em_t_2(i,1,j) ! END IF ! END DO ! END DO ! ELSE ! DO j = jts, MIN(jde-1,jte) ! DO i = its, MIN(ide-1,ite) ! IF ( grid%tsk(i,j) .LT. 170 .or. grid%tsk(i,j) .GT. 400. ) THEN ! print *,'error in the grid%em_tsk' ! print *,'i,j=',i,j ! print *,'grid%landmask=',grid%landmask(i,j) ! print *,'grid%tsk, grid%sst, grid%tmn=',grid%tsk(i,j),grid%sst(i,j),grid%tmn(i,j) ! if(grid%tmn(i,j).gt.170. .and. grid%tmn(i,j).lt.400.)then ! grid%tsk(i,j)=grid%tmn(i,j) ! else if(grid%sst(i,j).gt.170. .and. grid%sst(i,j).lt.400.)then ! grid%tsk(i,j)=grid%sst(i,j) ! else ! CALL wrf_error_fatal ( 'grid%em_tsk unreasonable' ) ! end if ! END IF ! END DO ! END DO ! END IF ! ! ! Is the grid%em_tmn reasonable? ! ! DO j = jts, MIN(jde-1,jte) ! DO i = its, MIN(ide-1,ite) ! IF ( ( ( grid%tmn(i,j) .LT. 170. ) .OR. ( grid%tmn(i,j) .GT. 400. ) ) & ! .AND. ( grid%landmask(i,j) .GT. 0.5 ) ) THEN ! IF ( model_config_rec%sf_surface_physics(grid%id) .NE. LSMSCHEME ) THEN ! print *,'error in the grid%em_tmn' ! print *,'i,j=',i,j ! print *,'grid%landmask=',grid%landmask(i,j) ! print *,'grid%tsk, grid%sst, grid%tmn=',grid%tsk(i,j),grid%sst(i,j),grid%tmn(i,j) ! END IF ! ! if(grid%tsk(i,j).gt.170. .and. grid%tsk(i,j).lt.400.)then ! grid%tmn(i,j)=grid%tsk(i,j) ! else if(grid%sst(i,j).gt.170. .and. grid%sst(i,j).lt.400.)then ! grid%tmn(i,j)=grid%sst(i,j) ! else ! CALL wrf_error_fatal ( 'grid%em_tmn unreasonable' ) ! endif ! END IF ! END DO ! END DO ! ! interpolate_soil_tmw : SELECT CASE ( model_config_rec%sf_surface_physics(grid%id) ) ! ! CASE ( SLABSCHEME , LSMSCHEME , RUCLSMSCHEME ) ! CALL process_soil_real ( grid%tsk , grid%tmn , & ! grid%landmask , grid%sst , & ! st_input , sm_input , sw_input , st_levels_input , sm_levels_input , sw_levels_input , & ! grid%zs , grid%dzs , grid%tslb , grid%smois , 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 ) ! ! END SELECT interpolate_soil_tmw ! ! ! 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 /) !! 0.004, 0.065, 0.020, 0.004, 0.008 /) ! has extra levels for playa, lava, and white sand ! ! ! 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. ! ! IF ( domain_ClockIsStartTime(grid) ) THEN ! 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%tslb(i,1,j) .gt. 200 ) .and. & ! ( grid%tslb(i,1,j) .lt. 400 ) .and. ( grid%smois(i,1,j) .lt. 0.005 ) ) then ! print *,'Noah -> Noah: bad soil moisture at i,j = ',i,j,grid%smois(i,:,j) ! iicount = iicount + 1 ! grid%smois(i,:,j) = 0.005 ! END IF ! END DO ! END DO ! IF ( iicount .GT. 0 ) THEN ! print *,'Noah -> Noah: total number of small soil moisture locations = ',iicount ! END IF ! ELSE IF ( FLAG_SOILM000 .EQ. 1 ) THEN ! DO j = jts, MIN(jde-1,jte) ! DO i = its, MIN(ide-1,ite) ! grid%smois(i,:,j) = grid%smois(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%tslb(i,1,j) .gt. 200 ) .and. & ! ( grid%tslb(i,1,j) .lt. 400 ) .and. ( grid%smois(i,1,j) .lt. 0.005 ) ) then ! print *,'RUC -> Noah: bad soil moisture at i,j = ',i,j,grid%smois(i,:,j) ! iicount = iicount + 1 ! grid%smois(i,:,j) = 0.005 ! END IF ! END DO ! END DO ! IF ( iicount .GT. 0 ) THEN ! print *,'RUC -> Noah: total number of small soil moisture locations = ',iicount ! 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%smois(i,:,j) = MAX ( grid%smois(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 ! END IF ! ! ! Is the grid%tslb reasonable? ! ! IF ( internal_time_loop .NE. 1 ) THEN ! DO j = jts, MIN(jde-1,jte) ! DO ns = 1 , model_config_rec%num_soil_layers ! DO i = its, MIN(ide-1,ite) ! IF ( grid%tslb(i,ns,j) .LT. 170 .or. grid%tslb(i,ns,j) .GT. 400. ) THEN ! grid%tslb(i,ns,j) = grid%em_t_2(i,1,j) ! grid%smois(i,ns,j) = 0.3 ! END IF ! END DO ! END DO ! END DO ! ELSE ! DO j = jts, MIN(jde-1,jte) ! DO i = its, MIN(ide-1,ite) ! IF ( ( ( grid%tslb(i,1,j) .LT. 170. ) .OR. ( grid%tslb(i,1,j) .GT. 400. ) ) .AND. & ! ( grid%landmask(i,j) .GT. 0.5 ) ) THEN ! IF ( ( model_config_rec%sf_surface_physics(grid%id) .NE. LSMSCHEME ) .AND. & ! ( model_config_rec%sf_surface_physics(grid%id) .NE. RUCLSMSCHEME ) ) THEN ! print *,'error in the grid%tslb' ! print *,'i,j=',i,j ! print *,'grid%landmask=',grid%landmask(i,j) ! print *,'grid%tsk, grid%sst, grid%tmn=',grid%tsk(i,j),grid%sst(i,j),grid%tmn(i,j) ! print *,'grid%tslb = ',grid%tslb(i,:,j) ! print *,'old grid%smois = ',grid%smois(i,:,j) ! grid%smois(i,1,j) = 0.3 ! grid%smois(i,2,j) = 0.3 ! grid%smois(i,3,j) = 0.3 ! grid%smois(i,4,j) = 0.3 ! END IF ! ! IF ( (grid%tsk(i,j).GT.170. .AND. grid%tsk(i,j).LT.400.) .AND. & ! (grid%tmn(i,j).GT.170. .AND. grid%tmn(i,j).LT.400.) ) THEN ! fake_soil_temp : SELECT CASE ( model_config_rec%sf_surface_physics(grid%id) ) ! CASE ( SLABSCHEME ) ! DO ns = 1 , model_config_rec%num_soil_layers ! grid%tslb(i,ns,j) = ( grid%tsk(i,j)*(3.0 - grid%zs(ns)) + & ! grid%tmn(i,j)*(0.0 - grid%zs(ns)) ) /(3.0 - 0.0) ! END DO ! CASE ( LSMSCHEME , RUCLSMSCHEME ) ! CALL wrf_error_fatal ( 'Assigning constant soil moisture, bad idea') ! DO ns = 1 , model_config_rec%num_soil_layers ! grid%tslb(i,ns,j) = ( grid%tsk(i,j)*(3.0 - grid%zs(ns)) + & ! grid%tmn(i,j)*(0.0 - grid%zs(ns)) ) /(3.0 - 0.0) ! END DO ! END SELECT fake_soil_temp ! else if(grid%tsk(i,j).gt.170. .and. grid%tsk(i,j).lt.400.)then ! CALL wrf_error_fatal ( 'grid%tslb unreasonable 1' ) ! DO ns = 1 , model_config_rec%num_soil_layers ! grid%tslb(i,ns,j)=grid%tsk(i,j) ! END DO ! else if(grid%sst(i,j).gt.170. .and. grid%sst(i,j).lt.400.)then ! CALL wrf_error_fatal ( 'grid%tslb unreasonable 2' ) ! DO ns = 1 , model_config_rec%num_soil_layers ! grid%tslb(i,ns,j)=grid%sst(i,j) ! END DO ! else if(grid%tmn(i,j).gt.170. .and. grid%tmn(i,j).lt.400.)then ! CALL wrf_error_fatal ( 'grid%tslb unreasonable 3' ) ! DO ns = 1 , model_config_rec%num_soil_layers ! grid%tslb(i,ns,j)=grid%tmn(i,j) ! END DO ! else ! CALL wrf_error_fatal ( 'grid%tslb unreasonable 4' ) ! endif ! END IF ! END DO ! END DO ! END IF ! ! ! Adjustments for the seaice field AFTER the grid%tslb computations. This is ! ! is for the Noah LSM scheme. ! ! num_veg_cat = SIZE ( grid%landusef , DIM=2 ) ! num_soil_top_cat = SIZE ( grid%soilctop , DIM=2 ) ! num_soil_bot_cat = SIZE ( grid%soilcbot , DIM=2 ) ! CALL nl_get_seaice_threshold ( grid%id , grid%seaice_threshold ) ! CALL nl_get_isice ( grid%id , grid%isice ) ! CALL nl_get_iswater ( grid%id , grid%iswater ) ! CALL adjust_for_seaice_post ( grid%xice , grid%landmask , grid%tsk , grid%tsk_save , & ! grid%ivgtyp , grid%vegcat , grid%lu_index , & ! grid%xland , grid%landusef , grid%isltyp , grid%soilcat , & ! grid%soilctop , & ! grid%soilcbot , grid%tmn , grid%vegfra , & ! grid%tslb , grid%smois , grid%sh2o , & ! grid%seaice_threshold , & ! num_veg_cat , num_soil_top_cat , num_soil_bot_cat , & ! model_config_rec%num_soil_layers , & ! grid%iswater , grid%isice , & ! model_config_rec%sf_surface_physics(grid%id) , & ! ids , ide , jds , jde , kds , kde , & ! ims , ime , jms , jme , kms , kme , & ! its , ite , jts , jte , kts , kte ) ! ! ! Let us make sure (again) that the grid%landmask and the veg/soil categories match. ! !oops1=0 !oops2=0 ! DO j = jts, MIN(jde-1,jte) ! DO i = its, MIN(ide-1,ite) ! IF ( ( ( grid%landmask(i,j) .LT. 0.5 ) .AND. & ! ( grid%ivgtyp(i,j) .NE. config_flags%iswater .OR. grid%isltyp(i,j) .NE. 14 ) ) .OR. & ! ( ( grid%landmask(i,j) .GT. 0.5 ) .AND. & ! ( grid%ivgtyp(i,j) .EQ. config_flags%iswater .OR. grid%isltyp(i,j) .EQ. 14 ) ) ) THEN ! IF ( grid%tslb(i,1,j) .GT. 1. ) THEN !oops1=oops1+1 ! grid%ivgtyp(i,j) = 5 ! grid%isltyp(i,j) = 8 ! grid%landmask(i,j) = 1 ! grid%xland(i,j) = 1 ! ELSE IF ( grid%sst(i,j) .GT. 1. ) THEN !oops2=oops2+1 ! grid%ivgtyp(i,j) = config_flags%iswater ! grid%isltyp(i,j) = 14 ! grid%landmask(i,j) = 0 ! grid%xland(i,j) = 2 ! ELSE ! print *,'the grid%landmask and soil/veg cats do not match' ! print *,'i,j=',i,j ! print *,'grid%landmask=',grid%landmask(i,j) ! print *,'grid%ivgtyp=',grid%ivgtyp(i,j) ! print *,'grid%isltyp=',grid%isltyp(i,j) ! print *,'iswater=', config_flags%iswater ! print *,'grid%tslb=',grid%tslb(i,:,j) ! print *,'grid%sst=',grid%sst(i,j) ! CALL wrf_error_fatal ( 'mismatch_landmask_ivgtyp' ) ! END IF ! END IF ! END DO ! END DO !if (oops1.gt.0) then !print *,'points artificially set to land : ',oops1 !endif !if(oops2.gt.0) then !print *,'points artificially set to water: ',oops2 !endif !! fill grid%sst array with grid%em_tsk if missing in real input (needed for time-varying grid%sst in wrf) ! DO j = jts, MIN(jde-1,jte) ! DO i = its, MIN(ide-1,ite) ! IF ( flag_sst .NE. 1 ) THEN ! grid%sst(i,j) = grid%tsk(i,j) ! ENDIF ! END DO ! END DO ! From the full level data, we can get the half levels, reciprocals, and layer ! thicknesses. These are all defined at half level locations, so one less level. ! We allow the vertical coordinate to *accidently* come in upside down. We want ! the first full level to be the ground surface. ! Check whether grid%em_znw (full level) data are truly full levels. If not, we need to adjust them ! to be full levels. ! in this test, we check if grid%em_znw(1) is neither 0 nor 1 (within a tolerance of 10**-5) were_bad = .false. IF ( ( (grid%em_znw(1).LT.(1-1.E-5) ) .OR. ( grid%em_znw(1).GT.(1+1.E-5) ) ).AND. & ( (grid%em_znw(1).LT.(0-1.E-5) ) .OR. ( grid%em_znw(1).GT.(0+1.E-5) ) ) ) THEN were_bad = .true. print *,'Your grid%em_znw input values are probably half-levels. ' print *,grid%em_znw print *,'WRF expects grid%em_znw values to be full levels. ' print *,'Adjusting now to full levels...' ! We want to ignore the first value if it's negative IF (grid%em_znw(1).LT.0) THEN grid%em_znw(1)=0 END IF DO k=2,kde grid%em_znw(k)=2*grid%em_znw(k)-grid%em_znw(k-1) END DO END IF ! Let's check our changes IF ( ( ( grid%em_znw(1) .LT. (1-1.E-5) ) .OR. ( grid%em_znw(1) .GT. (1+1.E-5) ) ).AND. & ( ( grid%em_znw(1) .LT. (0-1.E-5) ) .OR. ( grid%em_znw(1) .GT. (0+1.E-5) ) ) ) THEN print *,'The input grid%em_znw height values were half-levels or erroneous. ' print *,'Attempts to treat the values as half-levels and change them ' print *,'to valid full levels failed.' CALL wrf_error_fatal("bad grid%em_znw values from input files") ELSE IF ( were_bad ) THEN print *,'...adjusted. grid%em_znw array now contains full eta level values. ' ENDIF IF ( grid%em_znw(1) .LT. grid%em_znw(kde) ) THEN DO k=1, kde/2 hold_znw = grid%em_znw(k) grid%em_znw(k)=grid%em_znw(kde+1-k) grid%em_znw(kde+1-k)=hold_znw END DO END IF DO k=1, kde-1 grid%em_dnw(k) = grid%em_znw(k+1) - grid%em_znw(k) grid%em_rdnw(k) = 1./grid%em_dnw(k) grid%em_znu(k) = 0.5*(grid%em_znw(k+1)+grid%em_znw(k)) END DO ! Now the same sort of computations with the half eta levels, even ANOTHER ! level less than the one above. DO k=2, kde-1 grid%em_dn(k) = 0.5*(grid%em_dnw(k)+grid%em_dnw(k-1)) grid%em_rdn(k) = 1./grid%em_dn(k) grid%em_fnp(k) = .5* grid%em_dnw(k )/grid%em_dn(k) grid%em_fnm(k) = .5* grid%em_dnw(k-1)/grid%em_dn(k) END DO ! Scads of vertical coefficients. cof1 = (2.*grid%em_dn(2)+grid%em_dn(3))/(grid%em_dn(2)+grid%em_dn(3))*grid%em_dnw(1)/grid%em_dn(2) cof2 = grid%em_dn(2) /(grid%em_dn(2)+grid%em_dn(3))*grid%em_dnw(1)/grid%em_dn(3) grid%cf1 = grid%em_fnp(2) + cof1 grid%cf2 = grid%em_fnm(2) - cof1 - cof2 grid%cf3 = cof2 grid%cfn = (.5*grid%em_dnw(kde-1)+grid%em_dn(kde-1))/grid%em_dn(kde-1) grid%cfn1 = -.5*grid%em_dnw(kde-1)/grid%em_dn(kde-1) ! Inverse grid distances. grid%rdx = 1./config_flags%dx grid%rdy = 1./config_flags%dy ! Some of the many weird geopotential initializations that we'll see today: grid%em_ph0 is total, ! and grid%em_ph_2 is a perturbation from the base state geopotential. We set the base geopotential ! at the lowest level to terrain elevation * gravity. DO j=jts,jte DO i=its,ite grid%em_ph0(i,1,j) = grid%ht(i,j) * g grid%em_ph_2(i,1,j) = 0. END DO END DO ! Base state potential temperature and inverse density (alpha = 1/rho) from ! the half eta levels and the base-profile surface pressure. Compute 1/rho ! from equation of state. The potential temperature is a perturbation from t0. DO j = jts, MIN(jte,jde-1) DO i = its, MIN(ite,ide-1) !****MARS !TODO: etudier si une meilleure formule n'existe pas pour Mars !TODO: mais il s'agit juste d'un etat de base ... !****MARS ! Base state pressure is a function of eta level and terrain, only, plus ! the hand full of constants: p00 (sea level pressure, Pa), t00 (sea level ! temperature, K), and A (temperature difference, from 1000 mb to 300 mb, K). !!****MARS !!ici il s'agit de definir un etat de base, de reference !!- on ne peut prendre le profil de temperature du modele !! qui conduit a des instabilites !! grid%em_t_init(i,k,j)=grid%em_t_2(i,k,j) - t0 est a eviter donc. !!- pour la pression de surface, aucune information !! sur un profil de temperature variable et non equilibre !! ne doit transparaitre !! p_surf = grid%psfc(i,j) pourquoi pas ... mais t y est utilisee ... !! !!>> l'etat de base ne doit dependre "geographiquement" que de la topographie !! !!****MARS p_surf = p00 * EXP ( -t00/a + ( (t00/a)**2 - 2.*g*grid%ht(i,j)/a/r_d ) **0.5 ) DO k = 1, kte-1 grid%em_php(i,k,j) = grid%em_znw(k)*(p_surf - grid%p_top) + grid%p_top ! temporary, full lev base pressure grid%em_pb(i,k,j) = grid%em_znu(k)*(p_surf - grid%p_top) + grid%p_top ! temp = MAX ( 200., t00 + A*LOG(grid%em_pb(i,k,j)/p00) ) ! temp = t00 + A*LOG(grid%em_pb(i,k,j)/p00) !!! MODIF WRFV3.1 - parameter tiso !!! have to change as well in start_em temp = MAX ( tiso, t00 + A*LOG(grid%em_pb(i,k,j)/p00) ) IF (( i .EQ. its ) .AND. ( j .EQ. jts )) print *, temp, k !!! MODIF WRFV3.1 - parameter tiso IF (planet .eq. "mars" ) THEN grid%em_t_init(i,k,j) = temp*(p00/grid%em_pb(i,k,j))**(r_d/cp) - t0 ELSE grid%em_t_init(i,k,j) = (temp**nu + nu*(TT00**nu)*log((p00/grid%em_pb(i,k,j))**rcp))**(1/nu) -t0 ENDIF grid%em_alb(i,k,j) = (r_d/p1000mb)*(grid%em_t_init(i,k,j)+t0)*(grid%em_pb(i,k,j)/p1000mb)**cvpm END DO ! Base state mu is defined as base state surface pressure minus grid%p_top grid%em_mub(i,j) = p_surf - grid%p_top ! Dry surface pressure is defined as the following (this mu is from the input file ! computed from the dry pressure). Here the dry pressure is just reconstituted. pd_surf = grid%em_mu0(i,j) + grid%p_top ! Integrate base geopotential, starting at terrain elevation. This assures that ! the base state is in exact hydrostatic balance with respect to the model equations. ! This field is on full levels. grid%em_phb(i,1,j) = grid%ht(i,j) * g IF (hypsometric_opt == 1) THEN DO k = 2,kte grid%em_phb(i,k,j) = grid%em_phb(i,k-1,j) - grid%em_dnw(k-1)*grid%em_mub(i,j)*grid%em_alb(i,k-1,j) END DO ELSE IF (hypsometric_opt == 2) THEN DO k = 2,kte pfu = grid%em_mub(i,j)*grid%em_znw(k) + grid%p_top pfd = grid%em_mub(i,j)*grid%em_znw(k-1) + grid%p_top phm = grid%em_mub(i,j)*grid%em_znu(k-1) + grid%p_top grid%em_phb(i,k,j) = grid%em_phb(i,k-1,j) + grid%em_alb(i,k-1,j)*phm*LOG(pfd/pfu) END DO END IF END DO END DO ! Fill in the outer rows and columns to allow us to be sloppy. IF ( ite .EQ. ide ) THEN i = ide DO j = jts, MIN(jde-1,jte) grid%em_mub(i,j) = grid%em_mub(i-1,j) grid%em_mu_2(i,j) = grid%em_mu_2(i-1,j) DO k = 1, kte-1 grid%em_pb(i,k,j) = grid%em_pb(i-1,k,j) grid%em_t_init(i,k,j) = grid%em_t_init(i-1,k,j) grid%em_alb(i,k,j) = grid%em_alb(i-1,k,j) END DO DO k = 1, kte grid%em_phb(i,k,j) = grid%em_phb(i-1,k,j) END DO END DO END IF IF ( jte .EQ. jde ) THEN j = jde DO i = its, ite grid%em_mub(i,j) = grid%em_mub(i,j-1) grid%em_mu_2(i,j) = grid%em_mu_2(i,j-1) DO k = 1, kte-1 grid%em_pb(i,k,j) = grid%em_pb(i,k,j-1) grid%em_t_init(i,k,j) = grid%em_t_init(i,k,j-1) grid%em_alb(i,k,j) = grid%em_alb(i,k,j-1) END DO DO k = 1, kte grid%em_phb(i,k,j) = grid%em_phb(i,k,j-1) END DO END DO END IF ! Compute the perturbation dry pressure (grid%em_mub + grid%em_mu_2 + ptop = dry grid%em_psfc). DO j = jts, min(jde-1,jte) DO i = its, min(ide-1,ite) grid%em_mu_2(i,j) = grid%em_mu0(i,j) - grid%em_mub(i,j) END DO END DO ! Fill in the outer rows and columns to allow us to be sloppy. IF ( ite .EQ. ide ) THEN i = ide DO j = jts, MIN(jde-1,jte) grid%em_mu_2(i,j) = grid%em_mu_2(i-1,j) END DO END IF IF ( jte .EQ. jde ) THEN j = jde DO i = its, ite grid%em_mu_2(i,j) = grid%em_mu_2(i,j-1) END DO END IF lev500 = 0 DO j = jts, min(jde-1,jte) DO i = its, min(ide-1,ite) ! Assign the potential temperature (perturbation from t0) and qv on all the mass ! point locations. DO k = 1 , kde-1 grid%em_t_2(i,k,j) = grid%em_t_2(i,k,j) - t0 END DO !!--------------------------------------------------------------- !!****MARS: no 500mb adjustment needed !!****MARS: must keep however the hydrostatic equation integration performed in this loop ! !!****MARS: the DO WHILE loop is deactivated, since we will always be in the case !!****MARS: ... of "ELSE dpmu = 0." !!--------------------------------------------------------------- ! dpmu = 10001. ! loop_count = 0 ! ! DO WHILE ( ( ABS(dpmu) .GT. 10. ) .AND. & ! ( loop_count .LT. 5 ) ) ! ! loop_count = loop_count + 1 ! Integrate the hydrostatic equation (from the RHS of the bigstep vertical momentum ! equation) down from the top to get the pressure perturbation. First get the pressure ! perturbation, moisture, and inverse density (total and perturbation) at the top-most level. k = kte-1 qvf1 = 0.5*(moist(i,k,j,P_QV)+moist(i,k,j,P_QV)) qvf2 = 1./(1.+qvf1) qvf1 = qvf1*qvf2 grid%em_p(i,k,j) = - 0.5*(grid%em_mu_2(i,j)+qvf1*grid%em_mub(i,j))/grid%em_rdnw(k)/qvf2 qvf = 1. + rvovrd*moist(i,k,j,P_QV) grid%em_alt(i,k,j) = (r_d/p1000mb)*(grid%em_t_2(i,k,j)+t0)*qvf& *(((grid%em_p(i,k,j)+grid%em_pb(i,k,j))/p1000mb)**cvpm) grid%em_al(i,k,j) = grid%em_alt(i,k,j) - grid%em_alb(i,k,j) ! Now, integrate down the column to compute the pressure perturbation, and diagnose the two ! inverse density fields (total and perturbation). DO k=kte-2,1,-1 qvf1 = 0.5*(moist(i,k,j,P_QV)+moist(i,k+1,j,P_QV)) qvf2 = 1./(1.+qvf1) qvf1 = qvf1*qvf2 grid%em_p(i,k,j) = grid%em_p(i,k+1,j) - (grid%em_mu_2(i,j) + qvf1*grid%em_mub(i,j))/qvf2/grid%em_rdn(k+1) qvf = 1. + rvovrd*moist(i,k,j,P_QV) grid%em_alt(i,k,j) = (r_d/p1000mb)*(grid%em_t_2(i,k,j)+t0)*qvf* & (((grid%em_p(i,k,j)+grid%em_pb(i,k,j))/p1000mb)**cvpm) grid%em_al(i,k,j) = grid%em_alt(i,k,j) - grid%em_alb(i,k,j) END DO ! This is the hydrostatic equation used in the model after the small timesteps. In ! the model, grid%em_al (inverse density) is computed from the geopotential. IF (hypsometric_opt == 1) THEN DO k = 2,kte grid%em_ph_2(i,k,j) = grid%em_ph_2(i,k-1,j) - & grid%em_dnw(k-1) * ( (grid%em_mub(i,j)+grid%em_mu_2(i,j))*grid%em_al(i,k-1,j) & + grid%em_mu_2(i,j)*grid%em_alb(i,k-1,j) ) grid%em_ph0(i,k,j) = grid%em_ph_2(i,k,j) + grid%em_phb(i,k,j) END DO ELSE IF (hypsometric_opt == 2) THEN ! Alternative hydrostatic eq.: dZ = -al*p*dLOG(p), where p is dry pressure. ! Note that al*p approximates Rd*T and dLOG(p) does z. ! Here T varies mostly linear with z, the first-order integration produces better result. PRINT*,"WEE ET AL. 2012 CORRECTION." grid%em_ph_2(i,1,j) = grid%em_phb(i,1,j) DO k = 2,kte pfu = grid%em_mu0(i,j)*grid%em_znw(k) + grid%p_top pfd = grid%em_mu0(i,j)*grid%em_znw(k-1) + grid%p_top phm = grid%em_mu0(i,j)*grid%em_znu(k-1) + grid%p_top grid%em_ph_2(i,k,j) = grid%em_ph_2(i,k-1,j) + grid%em_alt(i,k-1,j)*phm*LOG(pfd/pfu) END DO DO k = 1,kte grid%em_ph_2(i,k,j) = grid%em_ph_2(i,k,j) - grid%em_phb(i,k,j) END DO END IF ! ! Adjust the column pressure so that the computed 500 mb height is close to the ! ! input value (of course, not when we are doing hybrid input). ! ! IF ( ( flag_metgrid .EQ. 1 ) .AND. ( i .EQ. its ) .AND. ( j .EQ. jts ) ) THEN ! DO k = 1 , num_metgrid_levels ! IF ( ABS ( grid%em_p_gc(i,k,j) - 50000. ) .LT. 1. ) THEN ! lev500 = k ! EXIT ! END IF ! END DO ! END IF ! ! ! We only do the adjustment of height if we have the input data on pressure ! ! surfaces, and folks have asked to do this option. ! ! IF ( ( flag_metgrid .EQ. 1 ) .AND. & ! ( config_flags%adjust_heights ) .AND. & ! ( lev500 .NE. 0 ) ) THEN ! ! DO k = 2 , kte-1 ! ! ! Get the pressures on the full eta levels (grid%em_php is defined above as ! ! the full-lev base pressure, an easy array to use for 3d space). ! ! pl = grid%em_php(i,k ,j) + & ! ( grid%em_p(i,k-1 ,j) * ( grid%em_znw(k ) - grid%em_znu(k ) ) + & ! grid%em_p(i,k ,j) * ( grid%em_znu(k-1 ) - grid%em_znw(k ) ) ) / & ! ( grid%em_znu(k-1 ) - grid%em_znu(k ) ) ! pu = grid%em_php(i,k+1,j) + & ! ( grid%em_p(i,k-1+1,j) * ( grid%em_znw(k +1) - grid%em_znu(k+1) ) + & ! grid%em_p(i,k +1,j) * ( grid%em_znu(k-1+1) - grid%em_znw(k+1) ) ) / & ! ( grid%em_znu(k-1+1) - grid%em_znu(k+1) ) ! ! ! If these pressure levels trap 500 mb, use them to interpolate ! ! to the 500 mb level of the computed height. !!**** PB on MARS .... ? ! IF ( ( pl .GE. 50000. ) .AND. ( pu .LT. 50000. ) ) THEN ! zl = ( grid%em_ph_2(i,k ,j) + grid%em_phb(i,k ,j) ) / g ! zu = ( grid%em_ph_2(i,k+1,j) + grid%em_phb(i,k+1,j) ) / g ! ! z500 = ( zl * ( LOG(50000.) - LOG(pu ) ) + & ! zu * ( LOG(pl ) - LOG(50000.) ) ) / & ! ( LOG(pl) - LOG(pu) ) !! z500 = ( zl * ( (50000.) - (pu ) ) + & !! zu * ( (pl ) - (50000.) ) ) / & !! ( (pl) - (pu) ) ! ! ! Compute the difference of the 500 mb heights (computed minus input), and ! ! then the change in grid%em_mu_2. The grid%em_php is still full-levels, base pressure. ! ! dz500 = z500 - grid%em_ght_gc(i,lev500,j) ! tvsfc = ((grid%em_t_2(i,1,j)+t0)*((grid%em_p(i,1,j)+grid%em_php(i,1,j))/p1000mb)**(r_d/cp)) * & ! (1.+0.6*moist(i,1,j,P_QV)) ! dpmu = ( grid%em_php(i,1,j) + grid%em_p(i,1,j) ) * EXP ( g * dz500 / ( r_d * tvsfc ) ) ! dpmu = dpmu - ( grid%em_php(i,1,j) + grid%em_p(i,1,j) ) ! grid%em_mu_2(i,j) = grid%em_mu_2(i,j) - dpmu ! EXIT ! END IF ! ! END DO ! ELSE ! dpmu = 0. ! END IF ! ! END DO END DO END DO !!****MARS: we use WPS ! ! ! If this is data from the SI, then we probably do not have the original ! ! surface data laying around. Note that these are all the lowest levels ! ! of the respective 3d arrays. For surface pressure, we assume that the ! ! vertical gradient of grid%em_p prime is zilch. This is not all that important. ! ! These are filled in so that the various plotting routines have something ! ! to play with at the initial time for the model. ! ! IF ( flag_metgrid .NE. 1 ) THEN ! DO j = jts, min(jde-1,jte) ! DO i = its, min(ide,ite) ! grid%u10(i,j)=grid%em_u_2(i,1,j) ! END DO ! END DO ! ! DO j = jts, min(jde,jte) ! DO i = its, min(ide-1,ite) ! grid%v10(i,j)=grid%em_v_2(i,1,j) ! END DO ! END DO ! ! DO j = jts, min(jde-1,jte) ! DO i = its, min(ide-1,ite) ! p_surf = p00 * EXP ( -t00/a + ( (t00/a)**2 - 2.*g*grid%ht(i,j)/a/r_d ) **0.5 ) ! grid%psfc(i,j)=p_surf + grid%em_p(i,1,j) ! grid%q2(i,j)=moist(i,1,j,P_QV) ! grid%th2(i,j)=grid%em_t_2(i,1,j)+300. ! grid%t2(i,j)=grid%th2(i,j)*(((grid%em_p(i,1,j)+grid%em_pb(i,1,j))/p00)**(r_d/cp)) ! END DO ! END DO ! ! ! If this data is from WPS, then we have previously assigned the surface ! ! data for u, v, and t. If we have an input qv, welp, we assigned that one, ! ! too. Now we pick up the left overs, and if RH came in - we assign the ! ! mixing ratio. ! ! ELSE IF ( flag_metgrid .EQ. 1 ) THEN ! !!****MARS: we use WPS DO j = jts, min(jde-1,jte) DO i = its, min(ide-1,ite) p_surf = p00 * EXP ( -t00/a + ( (t00/a)**2 - 2.*g*grid%ht(i,j)/a/r_d ) **0.5 ) ! recompute the value of surface pressure as calculated by sfcprs2 grid%psfc(i,j)=p_surf + grid%em_p(i,1,j) !!grid%th2 is used for other purpose !grid%th2(i,j)=grid%t2(i,j)*(p00/(grid%em_p(i,1,j)+grid%em_pb(i,1,j)))**(r_d/cp) grid%th2(i,j)=0. !!TODO TODO TODO - waiting for an input END DO END DO !!NB: q2 is used for other purpose ... !IF ( flag_qv .NE. 1 ) THEN ! DO j = jts, min(jde-1,jte) ! DO i = its, min(ide-1,ite) ! grid%q2(i,j)=moist(i,1,j,P_QV) ! END DO ! END DO !END IF !!NB: q2 is used for other purpose ... ! END IF !!!!MARS !!!! !!!! useful for history files @ first step !!!! grid%em_phtot = grid%em_ph0 grid%em_ptot = grid%em_p + grid%em_pb print *, 'OK OK OK OK' !!!! !!!!MARS ips = its ; ipe = ite ; jps = jts ; jpe = jte ; kps = kts ; kpe = kte #ifdef DM_PARALLEL # include "HALO_EM_INIT_1.inc" # include "HALO_EM_INIT_2.inc" # include "HALO_EM_INIT_3.inc" # include "HALO_EM_INIT_4.inc" # include "HALO_EM_INIT_5.inc" #endif RETURN END SUBROUTINE init_domain_rk !--------------------------------------------------------------------- SUBROUTINE const_module_initialize ( p00 , t00 , a, tiso ) USE module_configure IMPLICIT NONE ! For the real-data-cases only. REAL , INTENT(OUT) :: p00 , t00 , a, tiso CALL nl_get_base_pres ( 1 , p00 ) CALL nl_get_base_temp ( 1 , t00 ) CALL nl_get_base_lapse ( 1 , a ) CALL nl_get_tiso ( 1 , tiso ) END SUBROUTINE const_module_initialize !------------------------------------------------------------------- SUBROUTINE rebalance_driver ( grid ) IMPLICIT NONE TYPE (domain) :: grid CALL rebalance( grid & ! #include "em_actual_new_args.inc" ! ) END SUBROUTINE rebalance_driver !--------------------------------------------------------------------- SUBROUTINE rebalance ( grid & ! #include "em_dummy_new_args.inc" ! ) IMPLICIT NONE TYPE (domain) :: grid #include "em_dummy_new_decl.inc" TYPE (grid_config_rec_type) :: config_flags REAL :: p_surf , pd_surf, p_surf_int , pb_int , ht_hold REAL :: qvf , qvf1 , qvf2 REAL :: p00 , t00 , a, tiso, temp1, temp2 REAL , DIMENSION(:,:,:) , ALLOCATABLE :: t_init_int ! Local domain indices and counters. INTEGER :: num_veg_cat , num_soil_top_cat , num_soil_bot_cat 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 REAL :: pfu, pfd, phm INTEGER :: hypsometric_opt = 1 ! classic !INTEGER :: hypsometric_opt = 2 ! Wee et al. 2012 correction #ifdef DM_PARALLEL # include "em_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 ; ! note that tile is entire patch its = grid%sp32 ; ite = grid%ep32 ; ! note that tile is entire patch jts = grid%sp33 ; jte = grid%ep33 ; ! note that 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 ; ! note that tile is entire patch jts = grid%sp32 ; jte = grid%ep32 ; ! note that tile is entire patch kts = grid%sp33 ; kte = grid%ep33 ; ! note that 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 ; ! note that tile is entire patch kts = grid%sp32 ; kte = grid%ep32 ; ! note that tile is entire patch jts = grid%sp33 ; jte = grid%ep33 ; ! note that tile is entire patch END SELECT ALLOCATE ( t_init_int(ims:ime,kms:kme,jms:jme) ) ! Some of the many weird geopotential initializations that we'll see today: grid%em_ph0 is total, ! and grid%em_ph_2 is a perturbation from the base state geopotential. We set the base geopotential ! at the lowest level to terrain elevation * gravity. DO j=jts,jte DO i=its,ite grid%em_ph0(i,1,j) = grid%ht_fine(i,j) * g grid%em_ph_2(i,1,j) = 0. END DO END DO ! To define the base state, we call a USER MODIFIED routine to set the three ! necessary constants: p00 (sea level pressure, Pa), t00 (sea level temperature, K), ! and A (temperature difference, from 1000 mb to 300 mb, K). CALL const_module_initialize ( p00 , t00 , a , tiso ) ! Base state potential temperature and inverse density (alpha = 1/rho) from ! the half eta levels and the base-profile surface pressure. Compute 1/rho ! from equation of state. The potential temperature is a perturbation from t0. DO j = jts, MIN(jte,jde-1) DO i = its, MIN(ite,ide-1) ! Base state pressure is a function of eta level and terrain, only, plus ! the hand full of constants: p00 (sea level pressure, Pa), t00 (sea level ! temperature, K), and A (temperature difference, from 1000 mb to 300 mb, K). ! The fine grid terrain is ht_fine, the interpolated is grid%em_ht. p_surf = p00 * EXP ( -t00/a + ( (t00/a)**2 - 2.*g*grid%ht_fine(i,j)/a/r_d ) **0.5 ) p_surf_int = p00 * EXP ( -t00/a + ( (t00/a)**2 - 2.*g*grid%ht(i,j) /a/r_d ) **0.5 ) DO k = 1, kte-1 grid%em_pb(i,k,j) = grid%em_znu(k)*(p_surf - grid%p_top) + grid%p_top pb_int = grid%em_znu(k)*(p_surf_int - grid%p_top) + grid%p_top ! grid%em_t_init(i,k,j) = (t00 + A*LOG(grid%em_pb(i,k,j)/p00))*(p00/grid%em_pb(i,k,j))**(r_d/cp) - t0 ! t_init_int(i,k,j)= (t00 + A*LOG(pb_int /p00))*(p00/pb_int )**(r_d/cp) - t0 temp1 = MAX(tiso,t00+A*LOG(grid%em_pb(i,k,j)/p00)) temp2 = MAX(tiso,t00+A*LOG( pb_int/p00)) IF (planet .eq. "mars" ) THEN grid%em_t_init(i,k,j) = temp1*(p00/grid%em_pb(i,k,j))**(r_d/cp) - t0 t_init_int(i,k,j) = temp2*(p00/pb_int )**(r_d/cp) - t0 ELSE grid%em_t_init(i,k,j) = (temp1**nu + nu*(TT00**nu)*log((p00/grid%em_pb(i,k,j))**(rcp)))**(1/nu) - t0 t_init_int(i,k,j) = (temp2**nu + nu*(TT00**nu)*log((p00/pb_int)**(rcp)))**(1/nu) - t0 ENDIF grid%em_alb(i,k,j) = (r_d/p1000mb)*(grid%em_t_init(i,k,j)+t0)*(grid%em_pb(i,k,j)/p1000mb)**cvpm END DO ! Base state mu is defined as base state surface pressure minus grid%p_top grid%em_mub(i,j) = p_surf - grid%p_top ! Dry surface pressure is defined as the following (this mu is from the input file ! computed from the dry pressure). Here the dry pressure is just reconstituted. pd_surf = ( grid%em_mub(i,j) + grid%em_mu_2(i,j) ) + grid%p_top ! Integrate base geopotential, starting at terrain elevation. This assures that ! the base state is in exact hydrostatic balance with respect to the model equations. ! This field is on full levels. grid%em_phb(i,1,j) = grid%ht_fine(i,j) * g IF (hypsometric_opt == 1) THEN DO k = 2,kte grid%em_phb(i,k,j) = grid%em_phb(i,k-1,j) - grid%em_dnw(k-1)*grid%em_mub(i,j)*grid%em_alb(i,k-1,j) END DO ELSE IF (hypsometric_opt == 2) THEN DO k = 2,kte pfu = grid%em_mub(i,j)*grid%em_znw(k) + grid%p_top pfd = grid%em_mub(i,j)*grid%em_znw(k-1) + grid%p_top phm = grid%em_mub(i,j)*grid%em_znu(k-1) + grid%p_top grid%em_phb(i,k,j) = grid%em_phb(i,k-1,j) + grid%em_alb(i,k-1,j)*phm*LOG(pfd/pfu) END DO END IF END DO END DO ! Replace interpolated terrain with fine grid values. DO j = jts, MIN(jte,jde-1) DO i = its, MIN(ite,ide-1) grid%ht(i,j) = grid%ht_fine(i,j) END DO END DO ! Perturbation fields. DO j = jts, min(jde-1,jte) DO i = its, min(ide-1,ite) ! The potential temperature is THETAnest = THETAinterp + ( TBARnest - TBARinterp) DO k = 1 , kde-1 grid%em_t_2(i,k,j) = grid%em_t_2(i,k,j) + ( grid%em_t_init(i,k,j) - t_init_int(i,k,j) ) END DO ! Integrate the hydrostatic equation (from the RHS of the bigstep vertical momentum ! equation) down from the top to get the pressure perturbation. First get the pressure ! perturbation, moisture, and inverse density (total and perturbation) at the top-most level. k = kte-1 qvf1 = 0.5*(moist(i,k,j,P_QV)+moist(i,k,j,P_QV)) qvf2 = 1./(1.+qvf1) qvf1 = qvf1*qvf2 grid%em_p(i,k,j) = - 0.5*(grid%em_mu_2(i,j)+qvf1*grid%em_mub(i,j))/grid%em_rdnw(k)/qvf2 qvf = 1. + rvovrd*moist(i,k,j,P_QV) grid%em_alt(i,k,j) = (r_d/p1000mb)*(grid%em_t_2(i,k,j)+t0)*qvf* & (((grid%em_p(i,k,j)+grid%em_pb(i,k,j))/p1000mb)**cvpm) grid%em_al(i,k,j) = grid%em_alt(i,k,j) - grid%em_alb(i,k,j) ! Now, integrate down the column to compute the pressure perturbation, and diagnose the two ! inverse density fields (total and perturbation). DO k=kte-2,1,-1 qvf1 = 0.5*(moist(i,k,j,P_QV)+moist(i,k+1,j,P_QV)) qvf2 = 1./(1.+qvf1) qvf1 = qvf1*qvf2 grid%em_p(i,k,j) = grid%em_p(i,k+1,j) - (grid%em_mu_2(i,j) + qvf1*grid%em_mub(i,j))/qvf2/grid%em_rdn(k+1) qvf = 1. + rvovrd*moist(i,k,j,P_QV) grid%em_alt(i,k,j) = (r_d/p1000mb)*(grid%em_t_2(i,k,j)+t0)*qvf* & (((grid%em_p(i,k,j)+grid%em_pb(i,k,j))/p1000mb)**cvpm) grid%em_al(i,k,j) = grid%em_alt(i,k,j) - grid%em_alb(i,k,j) END DO ! This is the hydrostatic equation used in the model after the small timesteps. In ! the model, grid%al (inverse density) is computed from the geopotential. IF (hypsometric_opt == 1) THEN DO k = 2,kte grid%em_ph_2(i,k,j) = grid%em_ph_2(i,k-1,j) - & grid%em_dnw(k-1) * ( (grid%em_mub(i,j)+grid%em_mu_2(i,j))*grid%em_al(i,k-1,j) & + grid%em_mu_2(i,j)*grid%em_alb(i,k-1,j) ) grid%em_ph0(i,k,j) = grid%em_ph_2(i,k,j) + grid%em_phb(i,k,j) END DO ELSE IF (hypsometric_opt == 2) THEN ! Alternative hydrostatic eq.: dZ = -al*p*dLOG(p), where p is dry pressure. ! Note that al*p approximates Rd*T and dLOG(p) does z. ! Here T varies mostly linear with z, the first-order integration produces better result. grid%em_ph_2(i,1,j) = grid%em_phb(i,1,j) DO k = 2,kte pfu = grid%em_mu0(i,j)*grid%em_znw(k) + grid%p_top pfd = grid%em_mu0(i,j)*grid%em_znw(k-1) + grid%p_top phm = grid%em_mu0(i,j)*grid%em_znu(k-1) + grid%p_top grid%em_ph_2(i,k,j) = grid%em_ph_2(i,k-1,j) + grid%em_alt(i,k-1,j)*phm*LOG(pfd/pfu) END DO DO k = 1,kte grid%em_ph_2(i,k,j) = grid%em_ph_2(i,k,j) - grid%em_phb(i,k,j) END DO END IF END DO END DO DEALLOCATE ( t_init_int ) ips = its ; ipe = ite ; jps = jts ; jpe = jte ; kps = kts ; kpe = kte #ifdef DM_PARALLEL # include "HALO_EM_INIT_1.inc" # include "HALO_EM_INIT_2.inc" # include "HALO_EM_INIT_3.inc" # include "HALO_EM_INIT_4.inc" # include "HALO_EM_INIT_5.inc" #endif END SUBROUTINE rebalance !--------------------------------------------------------------------- RECURSIVE SUBROUTINE find_my_parent ( grid_ptr_in , grid_ptr_out , id_i_am , id_wanted , found_the_id ) USE module_domain TYPE(domain) , POINTER :: grid_ptr_in , grid_ptr_out TYPE(domain) , POINTER :: grid_ptr_sibling INTEGER :: id_wanted , id_i_am LOGICAL :: found_the_id found_the_id = .FALSE. grid_ptr_sibling => grid_ptr_in DO WHILE ( ASSOCIATED ( grid_ptr_sibling ) ) IF ( grid_ptr_sibling%grid_id .EQ. id_wanted ) THEN found_the_id = .TRUE. grid_ptr_out => grid_ptr_sibling RETURN ELSE IF ( grid_ptr_sibling%num_nests .GT. 0 ) THEN grid_ptr_sibling => grid_ptr_sibling%nests(1)%ptr CALL find_my_parent ( grid_ptr_sibling , grid_ptr_out , id_i_am , id_wanted , found_the_id ) ELSE grid_ptr_sibling => grid_ptr_sibling%sibling END IF END DO END SUBROUTINE find_my_parent #endif !--------------------------------------------------------------------- #ifdef VERT_UNIT !This is a main program for a small unit test for the vertical interpolation. program vint implicit none integer , parameter :: ij = 3 integer , parameter :: keta = 30 integer , parameter :: kgen =20 integer :: ids , ide , jds , jde , kds , kde , & ims , ime , jms , jme , kms , kme , & its , ite , jts , jte , kts , kte integer :: generic real , dimension(1:ij,kgen,1:ij) :: fo , po real , dimension(1:ij,1:keta,1:ij) :: fn_calc , fn_interp , pn integer, parameter :: interp_type = 1 ! 2 ! integer, parameter :: lagrange_order = 2 ! 1 integer :: lagrange_order logical, parameter :: lowest_lev_from_sfc = .FALSE. ! .TRUE. real , parameter :: zap_close_levels = 500. ! 100. integer, parameter :: force_sfc_in_vinterp = 0 ! 6 integer :: k ids = 1 ; ide = ij ; jds = 1 ; jde = ij ; kds = 1 ; kde = keta ims = 1 ; ime = ij ; jms = 1 ; jme = ij ; kms = 1 ; kme = keta its = 1 ; ite = ij ; jts = 1 ; jte = ij ; kts = 1 ; kte = keta generic = kgen print *,' ' print *,'------------------------------------' print *,'UNIT TEST FOR VERTICAL INTERPOLATION' print *,'------------------------------------' print *,' ' do lagrange_order = 1 , 2 print *,' ' print *,'------------------------------------' print *,'Lagrange Order = ',lagrange_order print *,'------------------------------------' print *,' ' call fillitup ( fo , po , fn_calc , pn , & ids , ide , jds , jde , kds , kde , & ims , ime , jms , jme , kms , kme , & its , ite , jts , jte , kts , kte , & generic , lagrange_order ) print *,' ' print *,'Level Pressure Field' print *,' (Pa) (generic)' print *,'------------------------------------' print *,' ' do k = 1 , generic write (*,fmt='(i2,2x,f12.3,1x,g15.8)' ) & k,po(2,k,2),fo(2,k,2) end do print *,' ' call vert_interp ( fo , po , fn_interp , pn , & generic , 'T' , & interp_type , lagrange_order , lowest_lev_from_sfc , & zap_close_levels , force_sfc_in_vinterp , & ids , ide , jds , jde , kds , kde , & ims , ime , jms , jme , kms , kme , & its , ite , jts , jte , kts , kte ) print *,'Multi-Order Interpolator' print *,'------------------------------------' print *,' ' print *,'Level Pressure Field Field Field' print *,' (Pa) Calc Interp Diff' print *,'------------------------------------' print *,' ' do k = kts , kte-1 write (*,fmt='(i2,2x,f12.3,1x,3(g15.7))' ) & k,pn(2,k,2),fn_calc(2,k,2),fn_interp(2,k,2),fn_calc(2,k,2)-fn_interp(2,k,2) end do call vert_interp_old ( fo , po , fn_interp , pn , & generic , 'T' , & interp_type , lagrange_order , lowest_lev_from_sfc , & zap_close_levels , force_sfc_in_vinterp , & ids , ide , jds , jde , kds , kde , & ims , ime , jms , jme , kms , kme , & its , ite , jts , jte , kts , kte ) print *,'Linear Interpolator' print *,'------------------------------------' print *,' ' print *,'Level Pressure Field Field Field' print *,' (Pa) Calc Interp Diff' print *,'------------------------------------' print *,' ' do k = kts , kte-1 write (*,fmt='(i2,2x,f12.3,1x,3(g15.7))' ) & k,pn(2,k,2),fn_calc(2,k,2),fn_interp(2,k,2),fn_calc(2,k,2)-fn_interp(2,k,2) end do end do end program vint subroutine wrf_error_fatal (string) character (len=*) :: string print *,string stop end subroutine wrf_error_fatal subroutine fillitup ( fo , po , fn , pn , & ids , ide , jds , jde , kds , kde , & ims , ime , jms , jme , kms , kme , & its , ite , jts , jte , kts , kte , & generic , lagrange_order ) implicit none integer , intent(in) :: ids , ide , jds , jde , kds , kde , & ims , ime , jms , jme , kms , kme , & its , ite , jts , jte , kts , kte integer , intent(in) :: generic , lagrange_order real , dimension(ims:ime,generic,jms:jme) , intent(out) :: fo , po real , dimension(ims:ime,kms:kme,jms:jme) , intent(out) :: fn , pn integer :: i , j , k real , parameter :: piov2 = 3.14159265358 / 2. k = 1 do j = jts , jte do i = its , ite po(i,k,j) = 102000. end do end do do k = 2 , generic do j = jts , jte do i = its , ite po(i,k,j) = ( 5000. * ( 1 - (k-1) ) + 100000. * ( (k-1) - (generic-1) ) ) / (1. - real(generic-1) ) end do end do end do if ( lagrange_order .eq. 1 ) then do k = 1 , generic do j = jts , jte do i = its , ite fo(i,k,j) = po(i,k,j) ! fo(i,k,j) = sin(po(i,k,j) * piov2 / 102000. ) end do end do end do else if ( lagrange_order .eq. 2 ) then do k = 1 , generic do j = jts , jte do i = its , ite fo(i,k,j) = (((po(i,k,j)-5000.)/102000.)*((102000.-po(i,k,j))/102000.))*102000. ! fo(i,k,j) = sin(po(i,k,j) * piov2 / 102000. ) end do end do end do end if !!!!!!!!!!!! do k = kts , kte do j = jts , jte do i = its , ite pn(i,k,j) = ( 5000. * ( 0 - (k-1) ) + 102000. * ( (k-1) - (kte-1) ) ) / (-1. * real(kte-1) ) end do end do end do do k = kts , kte-1 do j = jts , jte do i = its , ite pn(i,k,j) = ( pn(i,k,j) + pn(i,k+1,j) ) /2. end do end do end do if ( lagrange_order .eq. 1 ) then do k = kts , kte-1 do j = jts , jte do i = its , ite fn(i,k,j) = pn(i,k,j) ! fn(i,k,j) = sin(pn(i,k,j) * piov2 / 102000. ) end do end do end do else if ( lagrange_order .eq. 2 ) then do k = kts , kte-1 do j = jts , jte do i = its , ite fn(i,k,j) = (((pn(i,k,j)-5000.)/102000.)*((102000.-pn(i,k,j))/102000.))*102000. ! fn(i,k,j) = sin(pn(i,k,j) * piov2 / 102000. ) end do end do end do end if end subroutine fillitup #endif !--------------------------------------------------------------------- SUBROUTINE vert_interp ( fo , po , fnew , pnu , & generic , var_type , & interp_type , lagrange_order , lowest_lev_from_sfc , & zap_close_levels , force_sfc_in_vinterp , & ids , ide , jds , jde , kds , kde , & ims , ime , jms , jme , kms , kme , & its , ite , jts , jte , kts , kte ) ! Vertically interpolate the new field. The original field on the original ! pressure levels is provided, and the new pressure surfaces to interpolate to. IMPLICIT NONE INTEGER , INTENT(IN) :: interp_type , lagrange_order LOGICAL , INTENT(IN) :: lowest_lev_from_sfc REAL , INTENT(IN) :: zap_close_levels INTEGER , INTENT(IN) :: force_sfc_in_vinterp INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , & ims , ime , jms , jme , kms , kme , & its , ite , jts , jte , kts , kte INTEGER , INTENT(IN) :: generic CHARACTER (LEN=1) :: var_type REAL , DIMENSION(ims:ime,generic,jms:jme) , INTENT(IN) :: fo , po REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(IN) :: pnu REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(OUT) :: fnew REAL , DIMENSION(ims:ime,generic,jms:jme) :: forig , porig REAL , DIMENSION(ims:ime,kms:kme,jms:jme) :: pnew ! Local vars INTEGER :: i , j , k , ko , kn , k1 , k2 , ko_1 , ko_2 , knext INTEGER :: istart , iend , jstart , jend , kstart , kend INTEGER , DIMENSION(ims:ime,kms:kme ) :: k_above , k_below INTEGER , DIMENSION(ims:ime ) :: ks INTEGER , DIMENSION(ims:ime ) :: ko_above_sfc INTEGER :: count , zap , kst LOGICAL :: any_below_ground REAL :: p1 , p2 , pn, hold REAL , DIMENSION(1:generic) :: ordered_porig , ordered_forig REAL , DIMENSION(kts:kte) :: ordered_pnew , ordered_fnew !****MARS !big problems ... discontinuity in the interpolated fields ... print *, '25/05/2007: decided to use simple linear interpolations' print *, 'use that one at your own risk' !stop !****MARS ! Horiontal loop bounds for different variable types. IF ( var_type .EQ. 'U' ) THEN istart = its iend = ite jstart = jts jend = MIN(jde-1,jte) kstart = kts kend = kte-1 DO j = jstart,jend DO k = 1,generic DO i = MAX(ids+1,its) , MIN(ide-1,ite) porig(i,k,j) = ( po(i,k,j) + po(i-1,k,j) ) * 0.5 END DO END DO IF ( ids .EQ. its ) THEN DO k = 1,generic porig(its,k,j) = po(its,k,j) END DO END IF IF ( ide .EQ. ite ) THEN DO k = 1,generic porig(ite,k,j) = po(ite-1,k,j) END DO END IF DO k = kstart,kend DO i = MAX(ids+1,its) , MIN(ide-1,ite) pnew(i,k,j) = ( pnu(i,k,j) + pnu(i-1,k,j) ) * 0.5 END DO END DO IF ( ids .EQ. its ) THEN DO k = kstart,kend pnew(its,k,j) = pnu(its,k,j) END DO END IF IF ( ide .EQ. ite ) THEN DO k = kstart,kend pnew(ite,k,j) = pnu(ite-1,k,j) END DO END IF END DO ELSE IF ( var_type .EQ. 'V' ) THEN istart = its iend = MIN(ide-1,ite) jstart = jts jend = jte kstart = kts kend = kte-1 DO i = istart,iend DO k = 1,generic DO j = MAX(jds+1,jts) , MIN(jde-1,jte) porig(i,k,j) = ( po(i,k,j) + po(i,k,j-1) ) * 0.5 END DO END DO IF ( jds .EQ. jts ) THEN DO k = 1,generic porig(i,k,jts) = po(i,k,jts) END DO END IF IF ( jde .EQ. jte ) THEN DO k = 1,generic porig(i,k,jte) = po(i,k,jte-1) END DO END IF DO k = kstart,kend DO j = MAX(jds+1,jts) , MIN(jde-1,jte) pnew(i,k,j) = ( pnu(i,k,j) + pnu(i,k,j-1) ) * 0.5 END DO END DO IF ( jds .EQ. jts ) THEN DO k = kstart,kend pnew(i,k,jts) = pnu(i,k,jts) END DO END IF IF ( jde .EQ. jte ) THEN DO k = kstart,kend pnew(i,k,jte) = pnu(i,k,jte-1) END DO END IF END DO ELSE IF ( ( var_type .EQ. 'W' ) .OR. ( var_type .EQ. 'Z' ) ) THEN istart = its iend = MIN(ide-1,ite) jstart = jts jend = MIN(jde-1,jte) kstart = kts kend = kte DO j = jstart,jend DO k = 1,generic DO i = istart,iend porig(i,k,j) = po(i,k,j) END DO END DO DO k = kstart,kend DO i = istart,iend pnew(i,k,j) = pnu(i,k,j) END DO END DO END DO ELSE IF ( ( var_type .EQ. 'T' ) .OR. ( var_type .EQ. 'Q' ) ) THEN istart = its iend = MIN(ide-1,ite) jstart = jts jend = MIN(jde-1,jte) kstart = kts kend = kte-1 DO j = jstart,jend DO k = 1,generic DO i = istart,iend porig(i,k,j) = po(i,k,j) END DO END DO DO k = kstart,kend DO i = istart,iend pnew(i,k,j) = pnu(i,k,j) END DO END DO END DO ELSE istart = its iend = MIN(ide-1,ite) jstart = jts jend = MIN(jde-1,jte) kstart = kts kend = kte-1 DO j = jstart,jend DO k = 1,generic DO i = istart,iend porig(i,k,j) = po(i,k,j) END DO END DO DO k = kstart,kend DO i = istart,iend pnew(i,k,j) = pnu(i,k,j) END DO END DO END DO END IF DO j = jstart , jend ! The lowest level is the surface. Levels 2 through "generic" are supposed to ! be "bottom-up". Flip if they are not. This is based on the input pressure ! array. IF ( porig(its,2,j) .LT. porig(its,generic,j) ) THEN DO kn = 2 , ( generic + 1 ) / 2 DO i = istart , iend hold = porig(i,kn,j) porig(i,kn,j) = porig(i,generic+2-kn,j) porig(i,generic+2-kn,j) = hold forig(i,kn,j) = fo (i,generic+2-kn,j) forig(i,generic+2-kn,j) = fo (i,kn,j) END DO DO i = istart , iend forig(i,1,j) = fo (i,1,j) END DO END DO ELSE DO kn = 1 , generic DO i = istart , iend forig(i,kn,j) = fo (i,kn,j) END DO END DO END IF ! Skip all of the levels below ground in the original data based upon the surface pressure. ! The ko_above_sfc is the index in the pressure array that is above the surface. If there ! are no levels underground, this is index = 2. The remaining levels are eligible for use ! in the vertical interpolation. DO i = istart , iend ko_above_sfc(i) = -1 END DO DO ko = kstart+1 , kend DO i = istart , iend IF ( ko_above_sfc(i) .EQ. -1 ) THEN IF ( porig(i,1,j) .GT. porig(i,ko,j) ) THEN ko_above_sfc(i) = ko END IF END IF END DO END DO ! Piece together columns of the original input data. Pass the vertical columns to ! the iterpolator. DO i = istart , iend ! If the surface value is in the middle of the array, three steps: 1) do the ! values below the ground (this is just to catch the occasional value that is ! inconsistently below the surface based on input data), 2) do the surface level, then ! 3) add in the levels that are above the surface. For the levels next to the surface, ! we check to remove any levels that are "too close". When building the column of input ! pressures, we also attend to the request for forcing the surface analysis to be used ! in a few lower eta-levels. ! How many levels have we skipped in the input column. zap = 0 ! Fill in the column from up to the level just below the surface with the input ! presssure and the input field (orig or old, which ever). For an isobaric input ! file, this data is isobaric. IF ( ko_above_sfc(i) .GT. 2 ) THEN count = 1 DO ko = 2 , ko_above_sfc(i)-1 ordered_porig(count) = porig(i,ko,j) ordered_forig(count) = forig(i,ko,j) count = count + 1 END DO ! Make sure the pressure just below the surface is not "too close", this ! will cause havoc with the higher order interpolators. In case of a "too close" ! instance, we toss out the offending level (NOT the surface one) by simply ! decrementing the accumulating loop counter. IF ( ordered_porig(count-1) - porig(i,1,j) .LT. zap_close_levels ) THEN count = count -1 zap = 1 END IF ! Add in the surface values. ordered_porig(count) = porig(i,1,j) ordered_forig(count) = forig(i,1,j) count = count + 1 ! A usual way to do the vertical interpolation is to pay more attention to the ! surface data. Why? Well it has about 20x the density as the upper air, so we ! hope the analysis is better there. We more strongly use this data by artificially ! tossing out levels above the surface that are beneath a certain number of prescribed ! eta levels at this (i,j). The "zap" value is how many levels of input we are ! removing, which is used to tell the interpolator how many valid values are in ! the column. The "count" value is the increment to the index of levels, and is ! only used for assignments. IF ( force_sfc_in_vinterp .GT. 0 ) THEN ! Get the pressure at the eta level. We want to remove all input pressure levels ! between the level above the surface to the pressure at this eta surface. That ! forces the surface value to be used through the selected eta level. Keep track ! of two things: the level to use above the eta levels, and how many levels we are ! skipping. knext = ko_above_sfc(i) find_level : DO ko = ko_above_sfc(i) , generic IF ( porig(i,ko,j) .LE. pnew(i,force_sfc_in_vinterp,j) ) THEN knext = ko exit find_level ELSE zap = zap + 1 END IF END DO find_level ! No request for special interpolation, so we just assign the next level to use ! above the surface as, ta da, the first level above the surface. I know, wow. ELSE knext = ko_above_sfc(i) END IF ! One more time, make sure the pressure just above the surface is not "too close", this ! will cause havoc with the higher order interpolators. In case of a "too close" ! instance, we toss out the offending level above the surface (NOT the surface one) by simply ! incrementing the loop counter. Here, count-1 is the surface level and knext is either ! the next level up OR it is the level above the prescribed number of eta surfaces. IF ( ordered_porig(count-1) - porig(i,knext,j) .LT. zap_close_levels ) THEN kst = knext+1 zap = zap + 1 ELSE kst = knext END IF DO ko = kst , generic ordered_porig(count) = porig(i,ko,j) ordered_forig(count) = forig(i,ko,j) count = count + 1 END DO ! This is easy, the surface is the lowest level, just stick them in, in this order. OK, ! there are a couple of subtleties. We have to check for that special interpolation that ! skips some input levels so that the surface is used for the lowest few eta levels. Also, ! we must macke sure that we still do not have levels that are "too close" together. ELSE ! Initialize no input levels have yet been removed from consideration. zap = 0 ! The surface is the lowest level, so it gets set right away to location 1. ordered_porig(1) = porig(i,1,j) ordered_forig(1) = forig(i,1,j) ! We start filling in the array at loc 2, as in just above the level we just stored. count = 2 ! Are we forcing the interpolator to skip valid input levels so that the ! surface data is used through more levels? Essentially as above. IF ( force_sfc_in_vinterp .GT. 0 ) THEN knext = 2 find_level2: DO ko = 2 , generic IF ( porig(i,ko,j) .LE. pnew(i,force_sfc_in_vinterp,j) ) THEN knext = ko exit find_level2 ELSE zap = zap + 1 END IF END DO find_level2 ELSE knext = 2 END IF ! Fill in the data above the surface. The "knext" index is either the one ! just above the surface OR it is the index associated with the level that ! is just above the pressure at this (i,j) of the top eta level that is to ! be directly impacted with the surface level in interpolation. DO ko = knext , generic IF ( ordered_porig(count-1) - porig(i,ko,j) .LT. zap_close_levels ) THEN zap = zap + 1 CYCLE END IF ordered_porig(count) = porig(i,ko,j) ordered_forig(count) = forig(i,ko,j) count = count + 1 END DO END IF ! Now get the column of the "new" pressure data. So, this one is easy. DO kn = kstart , kend ordered_pnew(kn) = pnew(i,kn,j) END DO ! The polynomials are either in pressure or LOG(pressure). IF ( interp_type .EQ. 1 ) THEN CALL lagrange_setup ( var_type , & ordered_porig , ordered_forig , generic-zap , lagrange_order , & ordered_pnew , ordered_fnew , kend-kstart+1 ,i,j) ELSE CALL lagrange_setup ( var_type , & LOG(ordered_porig(1:generic-zap)) , ordered_forig , generic-zap , lagrange_order , & LOG(ordered_pnew(kstart:kend)) , ordered_fnew , kend-kstart+1 ,i,j) END IF ! Save the computed data. DO kn = kstart , kend fnew(i,kn,j) = ordered_fnew(kn) END DO ! There may have been a request to have the surface data from the input field ! to be assigned as to the lowest eta level. This assumes thin layers (usually ! the isobaric original field has the surface from 2-m T and RH, and 10-m U and V). IF ( lowest_lev_from_sfc ) THEN fnew(i,1,j) = forig(i,ko_above_sfc(i)-1,j) END IF END DO END DO END SUBROUTINE vert_interp !--------------------------------------------------------------------- SUBROUTINE vert_interp_old ( forig , po , fnew , pnu , & generic , var_type , & interp_type , lagrange_order , lowest_lev_from_sfc , & zap_close_levels , force_sfc_in_vinterp , & ids , ide , jds , jde , kds , kde , & ims , ime , jms , jme , kms , kme , & its , ite , jts , jte , kts , kte ) ! Vertically interpolate the new field. The original field on the original ! pressure levels is provided, and the new pressure surfaces to interpolate to. IMPLICIT NONE INTEGER , INTENT(IN) :: interp_type , lagrange_order LOGICAL , INTENT(IN) :: lowest_lev_from_sfc REAL , INTENT(IN) :: zap_close_levels INTEGER , INTENT(IN) :: force_sfc_in_vinterp INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , & ims , ime , jms , jme , kms , kme , & its , ite , jts , jte , kts , kte INTEGER , INTENT(IN) :: generic CHARACTER (LEN=1) :: var_type ! REAL , DIMENSION(ims:ime,generic,jms:jme) , INTENT(IN) :: forig , po !****MARS !error with g95 and warning with pgf90 REAL , DIMENSION(ims:ime,generic,jms:jme) , INTENT(IN) :: po REAL , DIMENSION(ims:ime,generic,jms:jme) , INTENT(INOUT) :: forig REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(IN) :: pnu REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(OUT) :: fnew REAL , DIMENSION(ims:ime,generic,jms:jme) :: porig REAL , DIMENSION(ims:ime,kms:kme,jms:jme) :: pnew ! Local vars INTEGER :: i , j , k , ko , kn , k1 , k2 , ko_1 , ko_2 INTEGER :: istart , iend , jstart , jend , kstart , kend INTEGER , DIMENSION(ims:ime,kms:kme ) :: k_above , k_below INTEGER , DIMENSION(ims:ime ) :: ks INTEGER , DIMENSION(ims:ime ) :: ko_above_sfc LOGICAL :: any_below_ground REAL :: p1 , p2 , pn !****MARS integer vert_extrap integer kn_save vert_extrap = 0 kn_save = 0 !****MARS ! Horizontal loop bounds for different variable types. IF ( var_type .EQ. 'U' ) THEN istart = its iend = ite jstart = jts jend = MIN(jde-1,jte) kstart = kts kend = kte-1 DO j = jstart,jend DO k = 1,generic DO i = MAX(ids+1,its) , MIN(ide-1,ite) porig(i,k,j) = ( po(i,k,j) + po(i-1,k,j) ) * 0.5 END DO END DO IF ( ids .EQ. its ) THEN DO k = 1,generic porig(its,k,j) = po(its,k,j) END DO END IF IF ( ide .EQ. ite ) THEN DO k = 1,generic porig(ite,k,j) = po(ite-1,k,j) END DO END IF DO k = kstart,kend DO i = MAX(ids+1,its) , MIN(ide-1,ite) pnew(i,k,j) = ( pnu(i,k,j) + pnu(i-1,k,j) ) * 0.5 END DO END DO IF ( ids .EQ. its ) THEN DO k = kstart,kend pnew(its,k,j) = pnu(its,k,j) END DO END IF IF ( ide .EQ. ite ) THEN DO k = kstart,kend pnew(ite,k,j) = pnu(ite-1,k,j) END DO END IF END DO ELSE IF ( var_type .EQ. 'V' ) THEN istart = its iend = MIN(ide-1,ite) jstart = jts jend = jte kstart = kts kend = kte-1 DO i = istart,iend DO k = 1,generic DO j = MAX(jds+1,jts) , MIN(jde-1,jte) porig(i,k,j) = ( po(i,k,j) + po(i,k,j-1) ) * 0.5 END DO END DO IF ( jds .EQ. jts ) THEN DO k = 1,generic porig(i,k,jts) = po(i,k,jts) END DO END IF IF ( jde .EQ. jte ) THEN DO k = 1,generic porig(i,k,jte) = po(i,k,jte-1) END DO END IF DO k = kstart,kend DO j = MAX(jds+1,jts) , MIN(jde-1,jte) pnew(i,k,j) = ( pnu(i,k,j) + pnu(i,k,j-1) ) * 0.5 END DO END DO IF ( jds .EQ. jts ) THEN DO k = kstart,kend pnew(i,k,jts) = pnu(i,k,jts) END DO END IF IF ( jde .EQ. jte ) THEN DO k = kstart,kend pnew(i,k,jte) = pnu(i,k,jte-1) END DO END IF END DO ELSE IF ( ( var_type .EQ. 'W' ) .OR. ( var_type .EQ. 'Z' ) ) THEN istart = its iend = MIN(ide-1,ite) jstart = jts jend = MIN(jde-1,jte) kstart = kts kend = kte DO j = jstart,jend DO k = 1,generic DO i = istart,iend porig(i,k,j) = po(i,k,j) END DO END DO DO k = kstart,kend DO i = istart,iend pnew(i,k,j) = pnu(i,k,j) END DO END DO END DO ELSE IF ( ( var_type .EQ. 'T' ) .OR. ( var_type .EQ. 'Q' ) ) THEN istart = its iend = MIN(ide-1,ite) jstart = jts jend = MIN(jde-1,jte) kstart = kts kend = kte-1 DO j = jstart,jend DO k = 1,generic DO i = istart,iend porig(i,k,j) = po(i,k,j) END DO END DO DO k = kstart,kend DO i = istart,iend pnew(i,k,j) = pnu(i,k,j) END DO END DO END DO ELSE istart = its iend = MIN(ide-1,ite) jstart = jts jend = MIN(jde-1,jte) kstart = kts kend = kte-1 DO j = jstart,jend DO k = 1,generic DO i = istart,iend porig(i,k,j) = po(i,k,j) END DO END DO DO k = kstart,kend DO i = istart,iend pnew(i,k,j) = pnu(i,k,j) END DO END DO END DO END IF DO j = jstart , jend ! Skip all of the levels below ground in the original data based upon the surface pressure. ! The ko_above_sfc is the index in the pressure array that is above the surface. If there ! are no levels underground, this is index = 2. The remaining levels are eligible for use ! in the vertical interpolation. DO i = istart , iend ko_above_sfc(i) = -1 END DO DO ko = kstart+1 , kend DO i = istart , iend IF ( ko_above_sfc(i) .EQ. -1 ) THEN IF ( porig(i,1,j) .GT. porig(i,ko,j) ) THEN ko_above_sfc(i) = ko !!****MARS !!old stuff !! !! Pressure level may be OK, however data from the diagfi is possibly missing !IF (forig(i,ko,j) .EQ. -1.0e+30) THEN ! ko_above_sfc(i) = -1 !END IF ! !! Once the right start level is found, check that it is OK ! !! >> first column should be 1e30 or so, second column should be a realistic value ! !IF ( ko_above_sfc(i) .NE. -1 ) THEN ! ! print *, 'verif', forig(i,ko-1,j), forig(i,ko,j), forig(i,ko+1,j), ko ! !END IF !! !!****MARS END IF END IF END DO END DO ! Initialize interpolation location. These are the levels in the original pressure ! data that are physically below and above the targeted new pressure level. DO kn = kts , kte DO i = its , ite k_above(i,kn) = -1 k_below(i,kn) = -2 END DO END DO ! Starting location is no lower than previous found location. This is for O(n logn) ! and not O(n^2), where n is the number of vertical levels to search. DO i = its , ite ks(i) = 1 END DO ! Find trapping layer for interpolation. The kn index runs through all of the "new" ! levels of data. DO kn = kstart , kend DO i = istart , iend ! For each "new" level (kn), we search to find the trapping levels in the "orig" ! data. Most of the time, the "new" levels are the eta surfaces, and the "orig" ! levels are the input pressure levels. found_trap_above : DO ko = ks(i) , generic-1 ! Because we can have levels in the interpolation that are not valid, ! let's toss out any candidate orig pressure values that are below ground ! based on the surface pressure. If the level =1, then this IS the surface ! level, so we HAVE to keep that one, but maybe not the ones above. If the ! level (ks) is NOT=1, then we have to just CYCLE our loop to find a legit ! below-pressure value. If we are not below ground, then we choose two ! neighboring levels to test whether they surround the new pressure level. ! The input trapping levels that we are trying is the surface and the first valid ! level above the surface. IF ( ( ko .LT. ko_above_sfc(i) ) .AND. ( ko .EQ. 1 ) ) THEN ko_1 = ko ko_2 = ko_above_sfc(i) !!****MARS !!old remark: the possible issue is fixed later in the code ... !!****MARS ! The "below" level is underground, cycle until we get to a valid pressure ! above ground. ELSE IF ( ( ko .LT. ko_above_sfc(i) ) .AND. ( ko .NE. 1 ) ) THEN CYCLE found_trap_above ! The "below" level is above the surface, so we are in the clear to test these ! two levels out. ELSE ko_1 = ko ko_2 = ko+1 END IF ! The test of the candidate levels: "below" has to have a larger pressure, and ! "above" has to have a smaller pressure. ! OK, we found the correct two surrounding levels. The locations are saved for use in the ! interpolation. IF ( ( porig(i,ko_1,j) .GE. pnew(i,kn,j) ) .AND. & ( porig(i,ko_2,j) .LT. pnew(i,kn,j) ) ) THEN k_above(i,kn) = ko_2 k_below(i,kn) = ko_1 ks(i) = ko_1 EXIT found_trap_above ! What do we do is we need to extrapolate the data underground? This happens when the ! lowest pressure that we have is physically "above" the new target pressure. Our ! actions depend on the type of variable we are interpolating. ELSE IF ( porig(i,1,j) .LT. pnew(i,kn,j) ) THEN !!****MARS !!old stuff !!check: values are usually quite close !print *,porig(i,1,j),pnew(i,kn,j) !!****MARS ! For horizontal winds and moisture, we keep a constant value under ground. IF ( ( var_type .EQ. 'U' ) .OR. & ( var_type .EQ. 'V' ) .OR. & ( var_type .EQ. 'Q' ) ) THEN k_above(i,kn) = 1 ks(i) = 1 ! For temperature and height, we extrapolate the data. Hopefully, we are not ! extrapolating too far. For pressure level input, the eta levels are always ! contained within the surface to p_top levels, so no extrapolation is ever ! required. ELSE IF ( ( var_type .EQ. 'Z' ) .OR. & ( var_type .EQ. 'T' ) ) THEN k_above(i,kn) = ko_above_sfc(i) k_below(i,kn) = 1 ks(i) = 1 !!!****MARS !!old stuff !k_above(i,kn) = 1 !ks(i) = 1 !!!"Hopefully, we are not extrapolating too far" !!!>> true on Mars ?? !!!****MARS ! Just a catch all right now. ELSE k_above(i,kn) = 1 ks(i) = 1 END IF EXIT found_trap_above ! The other extrapolation that might be required is when we are going above the ! top level of the input data. Usually this means we chose a P_PTOP value that ! was inappropriate, and we should stop and let someone fix this mess. ELSE IF ( porig(i,generic,j) .GT. pnew(i,kn,j) ) THEN print *,'data is too high, try a lower p_top' print *,'pnew=',pnew(i,kn,j),'i',i,'j',j,'kn',kn print *,'pnew=',pnew(i,:,j) print *,'porig=',porig(i,:,j) CALL wrf_error_fatal ('requested p_top is higher than input data, lower p_top') END IF END DO found_trap_above END DO END DO ! Linear vertical interpolation. DO kn = kstart , kend DO i = istart , iend IF ( k_above(i,kn) .EQ. 1 ) THEN !!!****MARS !!old stuff !!!ne doit pas arriver avec la temperature si l'on definit bien le champ au sol !IF (forig(i,1,j) .EQ. -1.0e+30) THEN ! print *,'no data here - surface - var is ...',var_type,i,j,1 ! print *,'setting to first level with data...',ko_above_sfc(i),porig(i,ko_above_sfc(i),j) ! forig(i,1,j) = forig(i,ko_above_sfc(i),j) ! !IF ( ( var_type .EQ. 'U' ) .OR. & ! ! ( var_type .EQ. 'V' ) .OR. & ! ! ( var_type .EQ. 'Q' ) ) THEN ! ! print *,'zero wind at the ground' ! ! forig(i,1,j) = 0 ! !ENDIF ! IF (forig(i,1,j) .EQ. -1.0e+30) THEN ! print *,'well ... are you sure ?' ! stop ! ENDIF !END IF !!!****MARS fnew(i,kn,j) = forig(i,1,j) ELSE k2 = MAX ( k_above(i,kn) , 2) k1 = MAX ( k_below(i,kn) , 1) IF ( k1 .EQ. k2 ) THEN CALL wrf_error_fatal ( 'identical values in the interp, bad for divisions' ) END IF !!!****MARS !!old stuff !IF (forig(i,k2,j) .EQ. -1.0e+30) THEN ! print *,'no data here - level above - you_d better stop',i,j,k2 ! stop !END IF !IF (forig(i,k1,j) .EQ. -1.0e+30) THEN ! print *,'no data here - level below - var is ...',var_type,i,j,k1 ! print *,'setting to first level with data...',ko_above_sfc(i),porig(i,ko_above_sfc(i),j) ! forig(i,k1,j) = forig(i,ko_above_sfc(i),j) ! !!!VERIFIER QUE LA TEMPERATURE AU SOL N'EST PAS CONCERNEE ! !!!(montagnes=sources locales de chaleur) ! !!!normalement, pas de souci, et lors de l'exécution rien ne s'affiche !END IF !!!****MARS IF ( interp_type .EQ. 1 ) THEN p1 = porig(i,k1,j) p2 = porig(i,k2,j) pn = pnew(i,kn,j) ELSE IF ( interp_type .EQ. 2 ) THEN p1 = ALOG(porig(i,k1,j)) p2 = ALOG(porig(i,k2,j)) pn = ALOG(pnew(i,kn,j)) END IF IF ( ( p1-pn) * (p2-pn) > 0. ) THEN ! CALL wrf_error_fatal ( 'both trapping pressures are on the same side of the new pressure' ) ! CALL wrf_debug ( 0 , 'both trapping pressures are on the same side of the new pressure' ) !!!****MARS vert_extrap = vert_extrap + 1 !print *, 'extrapolate', pnew(i,kn,j)-porig(i,k1,j), 'for WRF level', kn IF (kn_save < kn) kn_save=kn !!!****MARS END IF fnew(i,kn,j) = ( forig(i,k1,j) * ( p2 - pn ) + & forig(i,k2,j) * ( pn - p1 ) ) / & ( p2 - p1 ) END IF END DO END DO search_below_ground : DO kn = kstart , kend any_below_ground = .FALSE. DO i = istart , iend IF ( k_above(i,kn) .EQ. 1 ) THEN fnew(i,kn,j) = forig(i,1,j) any_below_ground = .TRUE. END IF END DO IF ( .NOT. any_below_ground ) THEN EXIT search_below_ground END IF END DO search_below_ground ! There may have been a request to have the surface data from the input field ! to be assigned as to the lowest eta level. This assumes thin layers (usually ! the isobaric original field has the surface from 2-m T and RH, and 10-m U and V). DO i = istart , iend IF ( lowest_lev_from_sfc ) THEN fnew(i,1,j) = forig(i,ko_above_sfc(i),j) END IF END DO END DO print *,'VERT EXTRAP = ', vert_extrap print *,'finished with ... ', var_type print *,'max WRF eta level where extrap. occurs: ',kn_save END SUBROUTINE vert_interp_old !--------------------------------------------------------------------- SUBROUTINE lagrange_setup ( var_type , all_x , all_y , all_dim , n , target_x , target_y , target_dim ,i,j) ! We call a Lagrange polynomial interpolator. The parallel concerns are put off as this ! is initially set up for vertical use. The purpose is an input column of pressure (all_x), ! and the associated pressure level data (all_y). These are assumed to be sorted (ascending ! or descending, no matter). The locations to be interpolated to are the pressures in ! target_x, probably the new vertical coordinate values. The field that is output is the ! target_y, which is defined at the target_x location. Mostly we expect to be 2nd order ! overlapping polynomials, with only a single 2nd order method near the top and bottom. ! When n=1, this is linear; when n=2, this is a second order interpolator. IMPLICIT NONE CHARACTER (LEN=1) :: var_type INTEGER , INTENT(IN) :: all_dim , n , target_dim REAL, DIMENSION(all_dim) , INTENT(IN) :: all_x , all_y REAL , DIMENSION(target_dim) , INTENT(IN) :: target_x REAL , DIMENSION(target_dim) , INTENT(OUT) :: target_y ! Brought in for debug purposes, all of the computations are in a single column. INTEGER , INTENT(IN) :: i,j ! Local vars REAL , DIMENSION(n+1) :: x , y REAL :: target_y_1 , target_y_2 LOGICAL :: found_loc INTEGER :: loop , loc_center_left , loc_center_right , ist , iend , target_loop IF ( all_dim .LT. n+1 ) THEN print *,'all_dim = ',all_dim print *,'order = ',n print *,'i,j = ',i,j print *,'p array = ',all_x print *,'f array = ',all_y print *,'p target= ',target_x CALL wrf_error_fatal ( 'troubles, the interpolating order is too large for this few input values' ) END IF IF ( n .LT. 1 ) THEN CALL wrf_error_fatal ( 'pal, linear is about as low as we go' ) END IF ! Loop over the list of target x and y values. DO target_loop = 1 , target_dim ! Find the two trapping x values, and keep the indices. found_loc = .FALSE. find_trap : DO loop = 1 , all_dim -1 IF ( ( target_x(target_loop) - all_x(loop) ) * ( target_x(target_loop) - all_x(loop+1) ) .LE. 0.0 ) THEN loc_center_left = loop loc_center_right = loop+1 found_loc = .TRUE. !****MARS: check if no errors here !print *,'interpolating ... ',var_type ! print *,'i,j = ',i,j ! print *,'target pressure and value = ',target_x(target_loop),target_y(target_loop) ! DO loop = 1 , all_dim ! print *,'column of pressure and value = ',all_x(loop),all_y(loop) ! END DO !END IF !****MARS EXIT find_trap END IF END DO find_trap IF ( ( .NOT. found_loc ) .AND. ( target_x(target_loop) .GT. all_x(1) ) ) THEN IF ( var_type .EQ. 'T' ) THEN write(6,fmt='(A,2i5,2f11.3)') & ' --> extrapolating TEMPERATURE near sfc: i,j,psfc, p target = ',& i,j,all_x(1),target_x(target_loop) target_y(target_loop) = ( all_y(1) * ( target_x(target_loop) - all_x(2) ) + & all_y(2) * ( all_x(1) - target_x(target_loop) ) ) / & ( all_x(1) - all_x(2) ) ELSE !write(6,fmt='(A,2i5,2f11.3)') & !' --> extrapolating zero gradient near sfc: i,j,psfc, p target = ',& !i,j,all_x(1),target_x(target_loop) target_y(target_loop) = all_y(1) END IF CYCLE ELSE IF ( .NOT. found_loc ) THEN !****MARS: normally, no errors here (otherwise, keep this part commented ?) print *, var_type print *,'i,j = ',i,j print *,'target pressure and value = ',target_x(target_loop),target_y(target_loop) DO loop = 1 , all_dim print *,'column of pressure and value = ',all_x(loop),all_y(loop) END DO CALL wrf_error_fatal ( 'troubles, could not find trapping x locations' ) !****MARS: end of 'keep this part commented' END IF ! Even or odd order? We can put the value in the middle if this is ! an odd order interpolator. For the even guys, we'll do it twice ! and shift the range one index, then get an average. IF ( MOD(n,2) .NE. 0 ) THEN IF ( ( loc_center_left -(((n+1)/2)-1) .GE. 1 ) .AND. & ( loc_center_right+(((n+1)/2)-1) .LE. all_dim ) ) THEN ist = loc_center_left -(((n+1)/2)-1) iend = iend + n CALL lagrange_interp ( all_x(ist:iend) , all_y(ist:iend) , n , target_x(target_loop) , target_y(target_loop) ) ELSE IF ( .NOT. found_loc ) THEN CALL wrf_error_fatal ( 'I doubt this will happen, I will only do 2nd order for now' ) END IF END IF ELSE IF ( MOD(n,2) .EQ. 0 ) THEN IF ( ( loc_center_left -(((n )/2)-1) .GE. 1 ) .AND. & ( loc_center_right+(((n )/2) ) .LE. all_dim ) .AND. & ( loc_center_left -(((n )/2) ) .GE. 1 ) .AND. & ( loc_center_right+(((n )/2)-1) .LE. all_dim ) ) THEN ist = loc_center_left -(((n )/2)-1) iend = ist + n CALL lagrange_interp ( all_x(ist:iend) , all_y(ist:iend) , n , target_x(target_loop) , target_y_1 ) ist = loc_center_left -(((n )/2) ) iend = ist + n CALL lagrange_interp ( all_x(ist:iend) , all_y(ist:iend) , n , target_x(target_loop) , target_y_2 ) target_y(target_loop) = ( target_y_1 + target_y_2 ) * 0.5 ELSE IF ( ( loc_center_left -(((n )/2)-1) .GE. 1 ) .AND. & ( loc_center_right+(((n )/2) ) .LE. all_dim ) ) THEN ist = loc_center_left -(((n )/2)-1) iend = ist + n CALL lagrange_interp ( all_x(ist:iend) , all_y(ist:iend) , n , target_x(target_loop) , target_y(target_loop) ) ELSE IF ( ( loc_center_left -(((n )/2) ) .GE. 1 ) .AND. & ( loc_center_right+(((n )/2)-1) .LE. all_dim ) ) THEN ist = loc_center_left -(((n )/2) ) iend = ist + n CALL lagrange_interp ( all_x(ist:iend) , all_y(ist:iend) , n , target_x(target_loop) , target_y(target_loop) ) ELSE CALL wrf_error_fatal ( 'unauthorized area, you should not be here' ) END IF END IF END DO END SUBROUTINE lagrange_setup !--------------------------------------------------------------------- SUBROUTINE lagrange_interp ( x , y , n , target_x , target_y ) ! Interpolation using Lagrange polynomials. ! P(x) = f(x0)Ln0(x) + ... + f(xn)Lnn(x) ! where Lnk(x) = (x -x0)(x -x1)...(x -xk-1)(x -xk+1)...(x -xn) ! --------------------------------------------- ! (xk-x0)(xk-x1)...(xk-xk-1)(xk-xk+1)...(xk-xn) IMPLICIT NONE INTEGER , INTENT(IN) :: n REAL , DIMENSION(0:n) , INTENT(IN) :: x , y REAL , INTENT(IN) :: target_x REAL , INTENT(OUT) :: target_y ! Local vars INTEGER :: i , k REAL :: numer , denom , Px REAL , DIMENSION(0:n) :: Ln Px = 0. DO i = 0 , n numer = 1. denom = 1. DO k = 0 , n IF ( k .EQ. i ) CYCLE numer = numer * ( target_x - x(k) ) denom = denom * ( x(i) - x(k) ) END DO Ln(i) = y(i) * numer / denom Px = Px + Ln(i) END DO target_y = Px END SUBROUTINE lagrange_interp #ifndef VERT_UNIT !--------------------------------------------------------------------- SUBROUTINE p_dry ( mu0 , eta , pdht , pdry , & ids , ide , jds , jde , kds , kde , & ims , ime , jms , jme , kms , kme , & its , ite , jts , jte , kts , kte ) ! Compute reference pressure and the reference mu. 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) , INTENT(IN) :: mu0 REAL , DIMENSION( kms:kme ) , INTENT(IN) :: eta REAL :: pdht REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(OUT) :: pdry ! Local vars INTEGER :: i , j , k REAL , DIMENSION( kms:kme ) :: eta_h DO k = kts , kte-1 eta_h(k) = ( eta(k) + eta(k+1) ) * 0.5 END DO DO j = jts , MIN ( jde-1 , jte ) DO k = kts , kte-1 DO i = its , MIN (ide-1 , ite ) pdry(i,k,j) = eta_h(k) * mu0(i,j) + pdht END DO END DO END DO END SUBROUTINE p_dry !--------------------------------------------------------------------- SUBROUTINE p_dts ( pdts , intq , psfc , p_top , & ids , ide , jds , jde , kds , kde , & ims , ime , jms , jme , kms , kme , & its , ite , jts , jte , kts , kte ) ! Compute difference between the dry, total surface pressure and the top pressure. IMPLICIT NONE INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , & ims , ime , jms , jme , kms , kme , & its , ite , jts , jte , kts , kte REAL , INTENT(IN) :: p_top REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: psfc REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: intq REAL , DIMENSION(ims:ime,jms:jme) , INTENT(OUT) :: pdts ! Local vars INTEGER :: i , j , k DO j = jts , MIN ( jde-1 , jte ) DO i = its , MIN (ide-1 , ite ) pdts(i,j) = psfc(i,j) - intq(i,j) - p_top END DO END DO END SUBROUTINE p_dts !--------------------------------------------------------------------- SUBROUTINE p_dhs ( pdhs , ht , p0 , t0 , a , & ids , ide , jds , jde , kds , kde , & ims , ime , jms , jme , kms , kme , & its , ite , jts , jte , kts , kte ) ! Compute dry, hydrostatic surface pressure. 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) , INTENT(IN) :: ht REAL , DIMENSION(ims:ime, jms:jme) , INTENT(OUT) :: pdhs REAL , INTENT(IN) :: p0 , t0 , a ! Local vars INTEGER :: i , j , k !****MARS .... REAL , PARAMETER :: Rd = 192. REAL , PARAMETER :: g = 3.72 print *,'compute dry, hydrostatic surface pressure' !****MARS .... DO j = jts , MIN ( jde-1 , jte ) DO i = its , MIN (ide-1 , ite ) pdhs(i,j) = p0 * EXP ( -t0/a + SQRT ( (t0/a)**2 - 2. * g * ht(i,j)/(a * Rd) ) ) END DO END DO !****MARS !****MARS cette formule est-elle juste sur Mars ? !****MARS >> a premiere vue, ne donne pas de resultats absurdes !****TODO: il y a peut etre meilleur ! !****MARS !print *,pdhs !stop END SUBROUTINE p_dhs !--------------------------------------------------------------------- SUBROUTINE find_p_top ( p , p_top , & ids , ide , jds , jde , kds , kde , & ims , ime , jms , jme , kms , kme , & its , ite , jts , jte , kts , kte ) ! Find the largest pressure in the top level. This is our p_top. We are ! assuming that the top level is the location where the pressure is a minimum ! for each column. In cases where the top surface is not isobaric, a ! communicated value must be shared in the calling routine. Also in cases ! where the top surface is not isobaric, care must be taken that the new ! maximum pressure is not greater than the previous value. This test is ! also handled in the calling routine. IMPLICIT NONE INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , & ims , ime , jms , jme , kms , kme , & its , ite , jts , jte , kts , kte REAL :: p_top REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(IN) :: p ! Local vars INTEGER :: i , j , k, min_lev i = its j = jts p_top = p(i,2,j) min_lev = 2 DO k = 2 , kte IF ( p_top .GT. p(i,k,j) ) THEN p_top = p(i,k,j) min_lev = k END IF END DO k = min_lev p_top = p(its,k,jts) DO j = jts , MIN ( jde-1 , jte ) DO i = its , MIN (ide-1 , ite ) p_top = MAX ( p_top , p(i,k,j) ) END DO END DO END SUBROUTINE find_p_top !--------------------------------------------------------------------- SUBROUTINE t_to_theta ( t , p , p00 , & ids , ide , jds , jde , kds , kde , & ims , ime , jms , jme , kms , kme , & its , ite , jts , jte , kts , kte ) ! Compute dry, hydrostatic surface pressure. IMPLICIT NONE INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , & ims , ime , jms , jme , kms , kme , & its , ite , jts , jte , kts , kte REAL , INTENT(IN) :: p00 REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(IN) :: p REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(INOUT) :: t ! Local vars INTEGER :: i , j , k !****MARS warning warning hardcoded !!!! ! REAL , PARAMETER :: Rd = 192. ! REAL , PARAMETER :: Cp = 844.6 REAL , PARAMETER :: Rd = 191. REAL , PARAMETER :: Cp = 744.5 !****MARS DO j = jts , MIN ( jde-1 , jte ) DO k = kts , kte DO i = its , MIN (ide-1 , ite ) t(i,k,j) = t(i,k,j) * ( p00 / p(i,k,j) ) ** (Rd / Cp) END DO END DO END DO END SUBROUTINE t_to_theta !--------------------------------------------------------------------- SUBROUTINE integ_moist ( q_in , p_in , pd_out , t_in , ght_in , intq , & ids , ide , jds , jde , kds , kde , & ims , ime , jms , jme , kms , kme , & its , ite , jts , jte , kts , kte ) ! Integrate the moisture field vertically. Mostly used to get the total ! vapor pressure, which can be subtracted from the total pressure to get ! the dry pressure. 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,kms:kme,jms:jme) , INTENT(IN) :: q_in , p_in , t_in , ght_in REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(OUT) :: pd_out REAL , DIMENSION(ims:ime, jms:jme) , INTENT(OUT) :: intq ! Local vars INTEGER :: i , j , k INTEGER , DIMENSION(ims:ime) :: level_above_sfc REAL , DIMENSION(ims:ime,jms:jme) :: psfc , tsfc , qsfc, zsfc REAL , DIMENSION(ims:ime,kms:kme) :: q , p , t , ght, pd REAL :: rhobar , qbar , dz REAL :: p1 , p2 , t1 , t2 , q1 , q2 , z1, z2 LOGICAL :: upside_down !****MARS REAL , PARAMETER :: Rd = 192. REAL , PARAMETER :: g = 3.72 !****MARS ! Get a surface value, always the first level of a 3d field. DO j = jts , MIN ( jde-1 , jte ) DO i = its , MIN (ide-1 , ite ) psfc(i,j) = p_in(i,kts,j) tsfc(i,j) = t_in(i,kts,j) qsfc(i,j) = q_in(i,kts,j) zsfc(i,j) = ght_in(i,kts,j) END DO END DO IF ( p_in(its,kts+1,jts) .LT. p_in(its,kte,jts) ) THEN upside_down = .TRUE. ELSE upside_down = .FALSE. END IF DO j = jts , MIN ( jde-1 , jte ) ! Initialize the integrated quantity of moisture to zero. DO i = its , MIN (ide-1 , ite ) intq(i,j) = 0. END DO IF ( upside_down ) THEN DO i = its , MIN (ide-1 , ite ) p(i,kts) = p_in(i,kts,j) t(i,kts) = t_in(i,kts,j) q(i,kts) = q_in(i,kts,j) ght(i,kts) = ght_in(i,kts,j) DO k = kts+1,kte p(i,k) = p_in(i,kte+2-k,j) t(i,k) = t_in(i,kte+2-k,j) q(i,k) = q_in(i,kte+2-k,j) ght(i,k) = ght_in(i,kte+2-k,j) END DO END DO ELSE DO i = its , MIN (ide-1 , ite ) DO k = kts,kte p(i,k) = p_in(i,k ,j) t(i,k) = t_in(i,k ,j) q(i,k) = q_in(i,k ,j) ght(i,k) = ght_in(i,k ,j) END DO END DO END IF ! Find the first level above the ground. If all of the levels are above ground, such as ! a terrain following lower coordinate, then the first level above ground is index #2. DO i = its , MIN (ide-1 , ite ) level_above_sfc(i) = -1 IF ( p(i,kts+1) .LT. psfc(i,j) ) THEN level_above_sfc(i) = kts+1 ELSE find_k : DO k = kts+1,kte-1 IF ( ( p(i,k )-psfc(i,j) .GE. 0. ) .AND. & ( p(i,k+1)-psfc(i,j) .LT. 0. ) ) THEN level_above_sfc(i) = k+1 EXIT find_k END IF END DO find_k IF ( level_above_sfc(i) .EQ. -1 ) THEN print *,'i,j = ',i,j print *,'p = ',p(i,:) print *,'p sfc = ',psfc(i,j) CALL wrf_error_fatal ( 'Could not find level above ground') END IF END IF END DO DO i = its , MIN (ide-1 , ite ) ! Account for the moisture above the ground. pd(i,kte) = p(i,kte) DO k = kte-1,level_above_sfc(i),-1 rhobar = ( p(i,k ) / ( Rd * t(i,k ) ) + & p(i,k+1) / ( Rd * t(i,k+1) ) ) * 0.5 qbar = ( q(i,k ) + q(i,k+1) ) * 0.5 dz = ght(i,k+1) - ght(i,k) intq(i,j) = intq(i,j) + g * qbar * rhobar / (1. + qbar) * dz pd(i,k) = p(i,k) - intq(i,j) END DO ! Account for the moisture between the surface and the first level up. IF ( ( p(i,level_above_sfc(i)-1)-psfc(i,j) .GE. 0. ) .AND. & ( p(i,level_above_sfc(i) )-psfc(i,j) .LT. 0. ) .AND. & ( level_above_sfc(i) .GT. kts ) ) THEN p1 = psfc(i,j) p2 = p(i,level_above_sfc(i)) t1 = tsfc(i,j) t2 = t(i,level_above_sfc(i)) q1 = qsfc(i,j) q2 = q(i,level_above_sfc(i)) z1 = zsfc(i,j) z2 = ght(i,level_above_sfc(i)) rhobar = ( p1 / ( Rd * t1 ) + & p2 / ( Rd * t2 ) ) * 0.5 qbar = ( q1 + q2 ) * 0.5 dz = z2 - z1 IF ( dz .GT. 0.1 ) THEN intq(i,j) = intq(i,j) + g * qbar * rhobar / (1. + qbar) * dz END IF ! Fix the underground values. DO k = level_above_sfc(i)-1,kts+1,-1 pd(i,k) = p(i,k) - intq(i,j) END DO END IF pd(i,kts) = psfc(i,j) - intq(i,j) END DO IF ( upside_down ) THEN DO i = its , MIN (ide-1 , ite ) pd_out(i,kts,j) = pd(i,kts) DO k = kts+1,kte pd_out(i,kte+2-k,j) = pd(i,k) END DO END DO ELSE DO i = its , MIN (ide-1 , ite ) DO k = kts,kte pd_out(i,k,j) = pd(i,k) END DO END DO END IF END DO !!!****MARS: no water vapor pressure !! DO k = level_above_sfc(i)-1,kts+1,-1 !! pd(i,k) = p(i,k) !! END DO !! pd(i,kts) = psfc(i,j) !!!****MARS END SUBROUTINE integ_moist !--------------------------------------------------------------------- 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,kms:kme,jms:jme) , INTENT(OUT) :: q ! Local vars INTEGER :: i , j , k REAL :: ew , q1 , t1 !****MARS .... regler si besoin .... !****MARS 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 !****MARS !****MARS ! 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 j = jts , MIN ( jde-1 , jte ) DO k = kts , kte DO i = its , MIN (ide-1 , ite ) rh(i,k,j) = MIN ( MAX ( rh(i,k,j) , 1. ) , 100. ) END DO END DO END DO IF ( wrt_liquid ) THEN DO j = jts , MIN ( jde-1 , jte ) DO k = kts , kte DO i = its , MIN (ide-1 , ite ) es=svp1*10.*EXP(svp2*(t(i,k,j)-svpt0)/(t(i,k,j)-svp3)) qs=eps*es/(p(i,k,j)/100.-es) q(i,k,j)=MAX(.01*rh(i,k,j)*qs,0.0) END DO END DO END DO ELSE DO j = jts , MIN ( jde-1 , jte ) DO k = kts , kte DO i = its , MIN (ide-1 , ite ) t1 = t(i,k,j) - 273.16 ! Obviously dry. IF ( t1 .lt. -200. ) THEN q(i,k,j) = 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,k,j) 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,k,j) * 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,k,j)/100. - ew)) q(i,k,j) = q1 / (1. - q1 ) END IF END DO END DO END DO END IF !!****MARS !!TODO: change once tracers are activated ? !q=0. !!****MARS END SUBROUTINE rh_to_mxrat !--------------------------------------------------------------------- SUBROUTINE compute_eta ( znw , & eta_levels , max_eta , max_dz , & fixedpbl, & p_top , g , p00 , cvpm , a , r_d , cp , t00 , p1000mb , t0 , & tiso, & ids , ide , jds , jde , kds , kde , & ims , ime , jms , jme , kms , kme , & its , ite , jts , jte , kts , kte ) ! Compute eta levels, either using given values from the namelist (hardly ! a computation, yep, I know), or assuming a constant dz above the PBL, ! knowing p_top and the number of eta levels. IMPLICIT NONE INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , & ims , ime , jms , jme , kms , kme , & its , ite , jts , jte , kts , kte REAL , INTENT(IN) :: max_dz REAL , INTENT(IN) :: p_top , g , p00 , cvpm , a , r_d , cp , t00 , p1000mb , t0 INTEGER , INTENT(IN) :: max_eta REAL , DIMENSION (max_eta) , INTENT(IN) :: eta_levels REAL , DIMENSION (kts:kte) , INTENT(OUT) :: znw ! Local vars INTEGER :: k REAL :: mub , t_init , p_surf , pb, ztop, ztop_pbl , dz , temp REAL , DIMENSION(kts:kte) :: dnw INTEGER , PARAMETER :: prac_levels = 17 INTEGER :: loop , loop1 REAL , DIMENSION(prac_levels) :: znw_prac , znu_prac , dnw_prac REAL , DIMENSION(kts:kte) :: alb , phb REAL :: z_scale REAL, INTENT(IN) :: tiso !****MARS !****MARS INTEGER :: fixedpbl ! usually, 8 first layers are fixed ! change this parameter if the top is very ! low print *, 'check Mars: p_top , g , p00 , cvpm , a , r_d , cp , t00 , p1000mb , t0' print *, p_top , g , p00 , cvpm , a , r_d , cp , t00 , p1000mb , t0 !-----solution alternative: definir dans la namelist les niveaux verticaux !****MARS !****MARS ! Gee, do the eta levels come in from the namelist? IF ( ABS(eta_levels(1)+1.) .GT. 0.0000001 ) THEN IF ( ( ABS(eta_levels(1 )-1.) .LT. 0.0000001 ) .AND. & ( ABS(eta_levels(kde)-0.) .LT. 0.0000001 ) ) THEN DO k = kds+1 , kde-1 znw(k) = eta_levels(k) END DO znw( 1) = 1. znw(kde) = 0. ELSE !CALL wrf_error_fatal ( 'First eta level should be 1.0 and the last 0.0 in namelist' ) print *, 'ok that s bad so I read the file, got it' !!MARS !!MARS open(unit=12,file='levels',form='formatted',status='old') rewind(12) DO k = kds, kde-1 read(12,*) znw(k) write(6,*) 'read level ', k, znw(k) ENDDO close(12) znw( 1) = 1. znw(kde) = 0. ! z_scale = .40 ! DO k=1, kde ! znw(k) = (exp(-(k-1)/float(kde-1)/z_scale) - exp(-1./z_scale))/ & ! (1.-exp(-1./z_scale)) ! ENDDO ! znw(1) = 1.0000 ! znw(2) = 0.9995 ! znw(3) = 0.9980 ! znw(4) = 0.9950 ! znw(5) = 0.9850 ! znw(6) = 0.9700 ! znw(7) = 0.9400 ! znw(8) = 0.9000 !!MARS !!MARS END IF !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! p_surf=p00 print *, 'prescribed levels' DO k = 1, kde pb = znw(k) * (p_surf - p_top) + p_top print *, 'level', k, & ', pressure (Pa)', pb, & ', logp height (m)', -10000.*log(pb/p00) END DO !mub = p_surf - p_top !DO k = 1, kde-1 ! pb = (znw(k)+znw(k+1))*0.5 * (p_surf - p_top) + p_top ! !temp = MAX ( 200., t00 + A*LOG(pb/p00) ) ! temp = t00 + A*LOG(pb/p00) ! t_init = temp*(p00/pb)**(r_d/cp) - t0 ! alb(k) = (r_d/p1000mb)*(t_init+t0)*(pb/p1000mb)**cvpm !END DO !phb(1) = 0. !DO k = 2,kde ! phb(k) = phb(k-1) - (znw(k)-znw(k-1)) * mub*alb(k-1) !END DO !ztop = phb(kde)/g !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Compute eta levels assuming a constant delta z above the PBL. ELSE ! Compute top of the atmosphere with some silly levels. We just want to ! integrate to get a reasonable value for ztop. We use the planned PBL-esque ! levels, and then just coarse resolution above that. We know p_top, and we ! have the base state vars. p_surf = p00 ! znw_prac = (/ 1.000 , 0.993 , 0.983 , 0.970 , 0.954 , 0.934 , 0.909 , & ! 0.88 , 0.8 , 0.7 , 0.6 , 0.5 , 0.4 , 0.3 , 0.2 , 0.1 , 0.0 /) !****MARS !****MARS ! on Mars, this is important to correctly resolve the surface ! -- levels were changed to get closer to the surface ! -- values were chosen as done typically in LMD GCM simulations !TODO: better repartition ? ! znw_prac = (/ 1.000 , & ! 0.9999 , & !1m ! 0.9995 , & !5m ! 0.9980 , & !20m ! 0.9950 , & !55m ! 0.9850 , & !166m ! 0.9550 , & !504m 0.9700 , & !334m 0.9400 , & !676m ! 0.9000 , & ! 0.8 , 0.7 , 0.6 , 0.5 , 0.4 , 0.3 , 0.2 , 0.1 , 0.0 /) znw_prac = (/ 1.000 , & 0.9995 , & !5m 0.9980 , & !20m 0.9950 , & !55m 0.9850 , & !166m 0.9700 , & !334m 0.9400 , & !676m 0.9000 , & 0.8 , 0.7 , 0.6 , 0.5 , 0.4 , 0.3 , 0.2 , 0.1 , 0.0 /) !****MARS !****MARS DO k = 1 , prac_levels - 1 znu_prac(k) = ( znw_prac(k) + znw_prac(k+1) ) * 0.5 dnw_prac(k) = znw_prac(k+1) - znw_prac(k) END DO DO k = 1, prac_levels-1 pb = znu_prac(k)*(p_surf - p_top) + p_top !! temp = MAX ( 200., t00 + A*LOG(pb/p00) ) ! temp = t00 + A*LOG(pb/p00) temp = MAX ( tiso, t00 + A*LOG(pb/p00) ) IF (planet .eq. "mars" ) THEN t_init = temp*(p00/pb)**(r_d/cp) - t0 ELSE t_init = (temp**nu + nu*(TT00**nu)*log((p00/pb)**(rcp)))**(1/nu) - t0 ENDIF alb(k) = (r_d/p1000mb)*(t_init+t0)*(pb/p1000mb)**cvpm END DO ! Base state mu is defined as base state surface pressure minus p_top mub = p_surf - p_top ! Integrate base geopotential, starting at terrain elevation. phb(1) = 0. DO k = 2,prac_levels phb(k) = phb(k-1) - dnw_prac(k-1)*mub*alb(k-1) END DO ! So, now we know the model top in meters. Get the average depth above the PBL ! of each of the remaining levels. We are going for a constant delta z thickness. ztop = phb(prac_levels) / g ztop_pbl = phb(fixedpbl) / g dz = ( ztop - ztop_pbl ) / REAL ( kde - fixedpbl ) ! Standard levels near the surface so no one gets in trouble. DO k = 1 , fixedpbl znw(k) = znw_prac(k) END DO ! Using d phb(k)/ d eta(k) = -mub * alb(k), eqn 2.9 ! Skamarock et al, NCAR TN 468. Use full levels, so ! use twice the thickness. DO k = fixedpbl, kte-1 pb = znw(k) * (p_surf - p_top) + p_top !! temp = MAX ( 200., t00 + A*LOG(pb/p00) ) ! temp = t00 + A*LOG(pb/p00) temp = MAX ( tiso, t00 + A*LOG(pb/p00) ) IF (planet .eq. "mars" ) THEN t_init = temp*(p00/pb)**(r_d/cp) - t0 ELSE t_init = (temp**nu + nu*(TT00**nu)*log((p00/pb)**(rcp)))**(1/nu) -t0 ENDIF alb(k) = (r_d/p1000mb)*(t_init+t0)*(pb/p1000mb)**cvpm znw(k+1) = znw(k) - dz*g / ( mub*alb(k) ) END DO znw(kte) = 0.000 ! There is some iteration. We want the top level, ztop, to be ! consistent with the delta z, and we want the half level values ! to be consistent with the eta levels. The inner loop to 10 gets ! the eta levels very accurately, but has a residual at the top, due ! to dz changing. We reset dz five times, and then things seem OK. DO loop1 = 1 , 5 DO loop = 1 , 10 DO k = fixedpbl, kte-1 pb = (znw(k)+znw(k+1))*0.5 * (p_surf - p_top) + p_top !! temp = MAX ( 200., t00 + A*LOG(pb/p00) ) ! temp = t00 + A*LOG(pb/p00) temp = MAX ( tiso, t00 + A*LOG(pb/p00) ) IF (planet .eq. "mars" ) THEN t_init = temp*(p00/pb)**(r_d/cp) - t0 ELSE t_init = (temp**nu + nu*(TT00**nu)*log((p00/pb)**(rcp)))**(1/nu) -t0 ENDIF alb(k) = (r_d/p1000mb)*(t_init+t0)*(pb/p1000mb)**cvpm znw(k+1) = znw(k) - dz*g / ( mub*alb(k) ) !!****MARS !!attention 'base_lapse' ne doit pas etre trop grand !!sinon ... des NaN car temperatures negatives en haut !IF ( ( loop1 .EQ. 5 ) .AND. ( loop .EQ. 10 ) ) THEN ! IF (k .EQ. 8) THEN ! print *, 'p,t,z,k' ! END IF ! print *, pb,temp,znw(k+1),k !END IF !****MARS END DO IF ( ( loop1 .EQ. 5 ) .AND. ( loop .EQ. 10 ) ) THEN print *,'Converged znw(kte) should be 0.0 = ',znw(kte) END IF znw(kte) = 0.000 END DO ! Here is where we check the eta levels values we just computed. DO k = 1, kde-1 pb = (znw(k)+znw(k+1))*0.5 * (p_surf - p_top) + p_top !! temp = MAX ( 200., t00 + A*LOG(pb/p00) ) ! temp = t00 + A*LOG(pb/p00) temp = MAX ( tiso, t00 + A*LOG(pb/p00) ) IF (planet .eq. "mars" ) THEN t_init = temp*(p00/pb)**(r_d/cp) - t0 ELSE t_init = (temp**nu + nu*(TT00**nu)*log((p00/pb)**(rcp)))**(1/nu) -t0 ENDIF alb(k) = (r_d/p1000mb)*(t_init+t0)*(pb/p1000mb)**cvpm END DO phb(1) = 0. DO k = 2,kde phb(k) = phb(k-1) - (znw(k)-znw(k-1)) * mub*alb(k-1) END DO ! Reset the model top and the dz, and iterate. ztop = phb(kde)/g ztop_pbl = phb(fixedpbl)/g dz = ( ztop - ztop_pbl ) / REAL ( kde - fixedpbl ) END DO ! ****MARS print *, 'eta_levels= ', znw ! Display the computed levels print *,'WRF levels are:' print *,'z (m) = ',phb(1)/g do k = 2 ,kte print *,'z (m) and dz (m) = ',phb(k)/g,(phb(k)-phb(k-1))/g !! little check of the repartition if (k>2) then if ((phb(k)-2.*phb(k-1)+phb(k-2))/g < -1.e-2) then print *, 'problem on the repartition' print *, '>> try to decrease force_sfc_in_vinterp (<8)' print *, '>> or increase model top (i.e. lower ptop)' print *, (phb(k)-2.*phb(k-1)+phb(k-2))/g stop endif endif end do ! ****MARS IF ( dz .GT. max_dz ) THEN print *,'z (m) = ',phb(1)/g do k = 2 ,kte print *,'z (m) and dz (m) = ',phb(k)/g,(phb(k)-phb(k-1))/g end do print *,'dz (m) above fixed eta levels = ',dz print *,'namelist max_dz (m) = ',max_dz print *,'namelist p_top (Pa) = ',p_top CALL wrf_debug ( 0, 'You need one of three things:' ) CALL wrf_debug ( 0, '1) More eta levels to reduce the dz: e_vert' ) CALL wrf_debug ( 0, '2) A lower p_top so your total height is reduced: p_top_requested') CALL wrf_debug ( 0, '3) Increase the maximum allowable eta thickness: max_dz') CALL wrf_debug ( 0, 'All are namelist options') CALL wrf_error_fatal ( 'dz above fixed eta levels is too large') END IF END IF END SUBROUTINE compute_eta !--------------------------------------------------------------------- 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,12,jms:jme) , 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,1,j) maxxer = field_in(i,1,j) DO l = 2 , 12 IF ( field_in(i,l,j) .LT. minner ) THEN minner = field_in(i,l,j) END IF IF ( field_in(i,l,j) .GT. maxxer ) THEN maxxer = field_in(i,l,j) END IF END DO field_min(i,j) = minner field_max(i,j) = maxxer END DO END DO END SUBROUTINE monthly_min_max !--------------------------------------------------------------------- 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,12,jms:jme) , 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 , month1 , month2 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 = l IF ( ( int_month .EQ. 0 ) .OR. ( int_month .EQ. 12 ) ) THEN month1 = 12 month2 = 1 ELSE month1 = int_month month2 = month1 + 1 END IF field_out(i,j) = ( field_in(i,month2,j) * ( target_date - middle(l) ) + & field_in(i,month1,j) * ( 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 sfcprs (t, q, height, pslv, ter, avgsfct, p, & psfc, ez_method, & ids , ide , jds , jde , kds , kde , & ims , ime , jms , jme , kms , kme , & its , ite , jts , jte , kts , kte ) ! Computes the surface pressure using the input height, ! temperature and q (already computed from relative ! humidity) on p surfaces. Sea level pressure is used ! to extrapolate a first guess. IMPLICIT NONE !****MARS: ok not used REAL , PARAMETER :: Rd = 192. REAL , PARAMETER :: Cp = 844.6 REAL, PARAMETER :: g = 3.72 REAL, PARAMETER :: pconst = 610. !****MARS !****MARS .... to be changed if used REAL, PARAMETER :: gamma = 6.5E-3 REAL, PARAMETER :: TC = 273.15 + 17.5 REAL, PARAMETER :: gammarg = gamma * Rd / g REAL, PARAMETER :: rov2 = Rd / 2. !****MARS .... to be changed if used INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , & ims , ime , jms , jme , kms , kme , & its , ite , jts , jte , kts , kte LOGICAL , INTENT ( IN ) :: ez_method REAL , DIMENSION (ims:ime,kms:kme,jms:jme) , INTENT(IN ):: t, q, height, p REAL , DIMENSION (ims:ime, jms:jme) , INTENT(IN ):: pslv , ter, avgsfct REAL , DIMENSION (ims:ime, jms:jme) , INTENT(OUT):: psfc INTEGER :: i INTEGER :: j INTEGER :: k INTEGER , DIMENSION (its:ite,jts:jte) :: k500 , k700 , k850 LOGICAL :: l1 LOGICAL :: l2 LOGICAL :: l3 LOGICAL :: OK REAL :: gamma78 ( its:ite,jts:jte ) REAL :: gamma57 ( its:ite,jts:jte ) REAL :: ht ( its:ite,jts:jte ) REAL :: p1 ( its:ite,jts:jte ) REAL :: t1 ( its:ite,jts:jte ) REAL :: t500 ( its:ite,jts:jte ) REAL :: t700 ( its:ite,jts:jte ) REAL :: t850 ( its:ite,jts:jte ) REAL :: tfixed ( its:ite,jts:jte ) REAL :: tsfc ( its:ite,jts:jte ) REAL :: tslv ( its:ite,jts:jte ) ! We either compute the surface pressure from a time averaged surface temperature ! (what we will call the "easy way"), or we try to remove the diurnal impact on the ! surface temperature (what we will call the "other way"). Both are essentially ! corrections to a sea level pressure with a high-resolution topography field. !****MARS .... !****MARS .... the mean sea level method is abandoned print *, 'no sea level pressure on Mars, please' stop !****MARS .... IF ( ez_method ) THEN DO j = jts , MIN(jde-1,jte) DO i = its , MIN(ide-1,ite) psfc(i,j) = pslv(i,j) * ( 1.0 + gamma * ter(i,j) / avgsfct(i,j) ) ** ( - g / ( Rd * gamma ) ) END DO END DO ELSE ! Find the locations of the 850, 700 and 500 mb levels. k850 = 0 ! find k at: P=850 k700 = 0 ! P=700 k500 = 0 ! P=500 i = its j = jts DO k = kts+1 , kte IF (NINT(p(i,k,j)) .EQ. 85000) THEN k850(i,j) = k ELSE IF (NINT(p(i,k,j)) .EQ. 70000) THEN k700(i,j) = k ELSE IF (NINT(p(i,k,j)) .EQ. 50000) THEN k500(i,j) = k END IF END DO IF ( ( k850(i,j) .EQ. 0 ) .OR. ( k700(i,j) .EQ. 0 ) .OR. ( k500(i,j) .EQ. 0 ) ) THEN DO j = jts , MIN(jde-1,jte) DO i = its , MIN(ide-1,ite) psfc(i,j) = pslv(i,j) * ( 1.0 + gamma * ter(i,j) / t(i,1,j) ) ** ( - g / ( Rd * gamma ) ) END DO END DO RETURN #if 0 ! Possibly it is just that we have a generalized vertical coord, so we do not ! have the values exactly. Do a simple assignment to a close vertical level. DO j = jts , MIN(jde-1,jte) DO i = its , MIN(ide-1,ite) DO k = kts+1 , kte-1 IF ( ( p(i,k,j) - 85000. ) * ( p(i,k+1,j) - 85000. ) .LE. 0.0 ) THEN k850(i,j) = k END IF IF ( ( p(i,k,j) - 70000. ) * ( p(i,k+1,j) - 70000. ) .LE. 0.0 ) THEN k700(i,j) = k END IF IF ( ( p(i,k,j) - 50000. ) * ( p(i,k+1,j) - 50000. ) .LE. 0.0 ) THEN k500(i,j) = k END IF END DO END DO END DO ! If we *still* do not have the k levels, punt. I mean, we did try. OK = .TRUE. DO j = jts , MIN(jde-1,jte) DO i = its , MIN(ide-1,ite) IF ( ( k850(i,j) .EQ. 0 ) .OR. ( k700(i,j) .EQ. 0 ) .OR. ( k500(i,j) .EQ. 0 ) ) THEN OK = .FALSE. PRINT '(A)','(i,j) = ',i,j,' Error in finding p level for 850, 700 or 500 hPa.' DO K = kts+1 , kte PRINT '(A,I3,A,F10.2,A)','K = ',k,' PRESSURE = ',p(i,k,j),' Pa' END DO PRINT '(A)','Expected 850, 700, and 500 mb values, at least.' END IF END DO END DO IF ( .NOT. OK ) THEN CALL wrf_error_fatal ( 'wrong pressure levels' ) END IF #endif ! We are here if the data is isobaric and we found the levels for 850, 700, ! and 500 mb right off the bat. ELSE DO j = jts , MIN(jde-1,jte) DO i = its , MIN(ide-1,ite) k850(i,j) = k850(its,jts) k700(i,j) = k700(its,jts) k500(i,j) = k500(its,jts) END DO END DO END IF ! The 850 hPa level of geopotential height is called something special. DO j = jts , MIN(jde-1,jte) DO i = its , MIN(ide-1,ite) ht(i,j) = height(i,k850(i,j),j) END DO END DO ! The variable ht is now -ter/ht(850 hPa). The plot thickens. DO j = jts , MIN(jde-1,jte) DO i = its , MIN(ide-1,ite) ht(i,j) = -ter(i,j) / ht(i,j) END DO END DO ! Make an isothermal assumption to get a first guess at the surface ! pressure. This is to tell us which levels to use for the lapse ! rates in a bit. DO j = jts , MIN(jde-1,jte) DO i = its , MIN(ide-1,ite) psfc(i,j) = pslv(i,j) * (pslv(i,j) / p(i,k850(i,j),j)) ** ht(i,j) END DO END DO ! Get a pressure more than pconst Pa above the surface - p1. The ! p1 is the top of the level that we will use for our lapse rate ! computations. DO j = jts , MIN(jde-1,jte) DO i = its , MIN(ide-1,ite) IF ( ( psfc(i,j) - 95000. ) .GE. 0. ) THEN p1(i,j) = 85000. ELSE IF ( ( psfc(i,j) - 70000. ) .GE. 0. ) THEN p1(i,j) = psfc(i,j) - pconst ELSE p1(i,j) = 50000. END IF END DO END DO ! Compute virtual temperatures for k850, k700, and k500 layers. Now ! you see why we wanted Q on pressure levels, it all is beginning ! to make sense. DO j = jts , MIN(jde-1,jte) DO i = its , MIN(ide-1,ite) t850(i,j) = t(i,k850(i,j),j) * (1. + 0.608 * q(i,k850(i,j),j)) t700(i,j) = t(i,k700(i,j),j) * (1. + 0.608 * q(i,k700(i,j),j)) t500(i,j) = t(i,k500(i,j),j) * (1. + 0.608 * q(i,k500(i,j),j)) END DO END DO ! Compute lapse rates between these three levels. These are ! environmental values for each (i,j). DO j = jts , MIN(jde-1,jte) DO i = its , MIN(ide-1,ite) gamma78(i,j) = ALOG(t850(i,j) / t700(i,j)) / ALOG (p(i,k850(i,j),j) / p(i,k700(i,j),j) ) gamma57(i,j) = ALOG(t700(i,j) / t500(i,j)) / ALOG (p(i,k700(i,j),j) / p(i,k500(i,j),j) ) END DO END DO DO j = jts , MIN(jde-1,jte) DO i = its , MIN(ide-1,ite) IF ( ( psfc(i,j) - 95000. ) .GE. 0. ) THEN t1(i,j) = t850(i,j) ELSE IF ( ( psfc(i,j) - 85000. ) .GE. 0. ) THEN t1(i,j) = t700(i,j) * (p1(i,j) / (p(i,k700(i,j),j))) ** gamma78(i,j) ELSE IF ( ( psfc(i,j) - 70000. ) .GE. 0.) THEN t1(i,j) = t500(i,j) * (p1(i,j) / (p(i,k500(i,j),j))) ** gamma57(i,j) ELSE t1(i,j) = t500(i,j) ENDIF END DO END DO ! From our temperature way up in the air, we extrapolate down to ! the sea level to get a guess at the sea level temperature. DO j = jts , MIN(jde-1,jte) DO i = its , MIN(ide-1,ite) tslv(i,j) = t1(i,j) * (pslv(i,j) / p1(i,j)) ** gammarg END DO END DO ! The new surface temperature is computed from the with new sea level ! temperature, just using the elevation and a lapse rate. This lapse ! rate is -6.5 K/km. DO j = jts , MIN(jde-1,jte) DO i = its , MIN(ide-1,ite) tsfc(i,j) = tslv(i,j) - gamma * ter(i,j) END DO END DO ! A small correction to the sea-level temperature, in case it is too warm. DO j = jts , MIN(jde-1,jte) DO i = its , MIN(ide-1,ite) tfixed(i,j) = tc - 0.005 * (tsfc(i,j) - tc) ** 2 END DO END DO DO j = jts , MIN(jde-1,jte) DO i = its , MIN(ide-1,ite) l1 = tslv(i,j) .LT. tc l2 = tsfc(i,j) .LE. tc l3 = .NOT. l1 IF ( l2 .AND. l3 ) THEN tslv(i,j) = tc ELSE IF ( ( .NOT. l2 ) .AND. l3 ) THEN tslv(i,j) = tfixed(i,j) END IF END DO END DO ! Finally, we can get to the surface pressure. DO j = jts , MIN(jde-1,jte) DO i = its , MIN(ide-1,ite) p1(i,j) = - ter(i,j) * g / ( rov2 * ( tsfc(i,j) + tslv(i,j) ) ) psfc(i,j) = pslv(i,j) * EXP ( p1(i,j) ) END DO END DO END IF ! Surface pressure and sea-level pressure are the same at sea level. ! DO j = jts , MIN(jde-1,jte) ! DO i = its , MIN(ide-1,ite) ! IF ( ABS ( ter(i,j) ) .LT. 0.1 ) THEN ! psfc(i,j) = pslv(i,j) ! END IF ! END DO ! END DO END SUBROUTINE sfcprs !--------------------------------------------------------------------- SUBROUTINE sfcprs2(t, q, height, psfc_in, ter, avgsfct, p, & psfc, ez_method, & ids , ide , jds , jde , kds , kde , & ims , ime , jms , jme , kms , kme , & its , ite , jts , jte , kts , kte ) ! Computes the surface pressure using the input height, ! temperature and q (already computed from relative ! humidity) on p surfaces. Sea level pressure is used ! to extrapolate a first guess. IMPLICIT NONE !****MARS: beware, hardcoded !!! ! REAL , PARAMETER :: Rd = 192. REAL, PARAMETER :: Rd = 191. REAL, PARAMETER :: g = 3.72 !****MARS INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , & ims , ime , jms , jme , kms , kme , & its , ite , jts , jte , kts , kte LOGICAL , INTENT ( IN ) :: ez_method REAL , DIMENSION (ims:ime,kms:kme,jms:jme) , INTENT(IN ):: t, q, height, p REAL , DIMENSION (ims:ime, jms:jme) , INTENT(IN ):: psfc_in , ter, avgsfct REAL , DIMENSION (ims:ime, jms:jme) , INTENT(OUT):: psfc INTEGER :: i INTEGER :: j INTEGER :: k REAL :: tv_sfc_avg , tv_sfc , del_z ! Compute the new surface pressure from the old surface pressure, and a ! known change in elevation at the surface. !****MARS: as is done in MCD/pres0 with the MOLA topography :) !!--------- !! del_z = diff in surface topo, lo-res vs hi-res !grid%em_ght_gc - grid%ht !!--------- !!* em_ght_gc: surface geopotential height from the GCM !!* ht: hi-res altimetry ! psfc = psfc_in * exp ( g del_z / (Rd Tv_sfc ) ) !!--------- IF ( ez_method ) THEN !! !!****MARS: 'ez_method' is 'we_have_tavgsfc', hard-coded as false !! DO j = jts , MIN(jde-1,jte) DO i = its , MIN(ide-1,ite) tv_sfc_avg = avgsfct(i,j) * (1. + 0.608 * q(i,1,j)) del_z = height(i,1,j) - ter(i,j) psfc(i,j) = psfc_in(i,j) * EXP ( g * del_z / ( Rd * tv_sfc_avg ) ) END DO END DO ELSE !! !!****MARS .... here is what is done for Mars !! DO j = jts , MIN(jde-1,jte) DO i = its , MIN(ide-1,ite) ! tv_sfc = t(i,1,j) * (1. + 0.608 * q(i,1,j)) !!****MARS: 0.608 >> nonsense on Mars tv_sfc = t(i,1,j) !!****MARS .... changer pour t_1km - 7e couche GCM !!****MARS .... spiga et al. (2007) tv_sfc = t(i,8,j) del_z = height(i,1,j) - ter(i,j) psfc(i,j) = psfc_in(i,j) * EXP ( g * del_z / ( Rd * tv_sfc ) ) !****MARS !****MARS .... which temperature is used in the Laplace formula ? !!****MARS: hardcoded as 220K (t0) !!****MARS: pas une enorme influence !psfc(i,j) = psfc_in(i,j) * EXP ( g * del_z / ( Rd * 220 ) ) ! !****MARS .... check of the altimetry differences ! print *,del_z, tv_sfc END DO END DO print *, '1 km temperatures - max' print *, MAXVAL(t(:,8,:)) END IF END SUBROUTINE sfcprs2 !--------------------------------------------------------------------- SUBROUTINE init_module_initialize END SUBROUTINE init_module_initialize !--------------------------------------------------------------------- SUBROUTINE constante3(field, field_custom, & ids , ide , jds , jde , kds , kde , & ims , ime , jms , jme , kms , kme , & its , ite , jts , jte , kts , kte ) IMPLICIT NONE REAL :: field_custom REAL, DIMENSION (ims:ime,kms:kme,jms:jme), INTENT(INOUT):: field INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , & ims , ime , jms , jme , kms , kme , & its , ite , jts , jte , kts , kte !!****MARS: set the 3D field to a constant value field(:,:,:)=field_custom END SUBROUTINE constante3 !--------------------------------------------------------------------- SUBROUTINE constante2(field, field_custom, & ids , ide , jds , jde , kds , kde , & ims , ime , jms , jme , kms , kme , & its , ite , jts , jte , kts , kte ) IMPLICIT NONE REAL :: field_custom REAL, DIMENSION (ims:ime,jms:jme), INTENT(INOUT):: field INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , & ims , ime , jms , jme , kms , kme , & its , ite , jts , jte , kts , kte !!****MARS: set the 3D field to a constant value field(:,:)=field_custom END SUBROUTINE constante2 !--------------------------------------------------------------------- subroutine build_sigma_hr(dimlevs,sigma_gcm,ps_gcm,ps_hr,sigma_hr) !,p_pgcm) implicit none ! include "constants_mcd.inc" !--------------------------------------- ! written by E. Millour and F. Forget ! Mars Climate Database v4.2 ! see DDD page 27 and following !--------------------------------------- INTEGER , INTENT(IN) :: dimlevs ! inputs real sigma_gcm(dimlevs) ! GCM sigma levels real ps_gcm ! GCM surface pressure real ps_hr ! High res surface pressure ! outputs real sigma_hr(dimlevs) ! High res sigma levels ! real p_pgcm(dimlevs) ! high res to GCM pressure ratios ! local variables integer l real x ! lower layer compression (-0.9