!``````````````````````````````````````````````````````````````````` subroutine toptinit_user(n2, n3, igd, topt, xt, yt) !+ Where user may customize topography values !~ 29.08.2012 (last modified; TIM) use univ_kind_defs_mod, only: f4 use main_memory_mod, only: gd, initialization_type, run_toptinit_user implicit none !=== Argument declarations: integer, intent(IN) :: n2, n3, igd real(f4), dimension(n2), intent(IN) :: xt real(f4), dimension(n3), intent(IN) :: yt real(f4), dimension(n2,n3), intent(INOUT) :: topt !== Local declarations: ! This subroutine is the intended location for a user to customize TOPT, ! the surface topography array. It is called after all other types of ! initialization of this field, so this subroutine has the last word. ! By default the subroutine makes no change to the field. The commented ! lines below serve as a template for user-designed changes; the example ! shown is the Witch of Agnesi mountain, a common test case. Note that ! this routine is called for each grid separately, so attention to the ! current value of NG in this routine may be required by the user: real(f4) :: half_width, height, halfw2, r integer :: i, j real(f4) :: diameter, depth, height_ag, cent_xt, cent_yt, c1, c2, c3 integer :: cent_i, cent_j !=== Executable statements: if (.not.run_toptinit_user) return ! Not desired, currently ! Parabolic crater: diameter = 3500. ! (measured from the top of one rim to the other) [m] depth = 0.12*diameter**0.96 ! Total depth (top of rim to bottom of floor); [m] ! functional form from Matias_etal (LPSC 2000) height_ag = depth/4. ! Height (above surrounding plain) of top of rim [m] cent_i = n2/2 ! Center grid point coords cent_j = max(n3/2, 1) ! cent_xt = xt(cent_i) cent_yt = yt(cent_j) c1 = 4.*depth/(diameter*diameter) c2 = height_ag-depth c3 = diameter*0.5-sqrt(-c2/c1) do j=1,n3 do i=1,n2 r = sqrt((xt(i)-cent_xt)**2+(yt(j)-cent_yt)**2) if (abs(r) <= diameter*0.5) then topt(i,j) = c1*r*r+c2 else if (abs(r) <= diameter*0.5+c3) then r = diameter-r topt(i,j) = c1*r*r+c2 else topt(i,j) = 0. end if end do end do !if(igd == 1) then !if (initial == 1) then ! half_width = 15000. ! height = 1000. ! halfw2 = half_width*half_width ! ! do j=1,n3 ! do i=1,n2 ! r = sqrt(xt(i)**2+yt(j)**2) ! topt(i,j) = height*halfw2/(halfw2+r**2) !xt(i)**2) ! ! enddo ! enddo !elseif(igd == 2) then ! ! Whatever you want here ! !end if !*** UNCOMMENT the below statement to zero out topography on all grids: !topt(1:n2,1:n3) = 0. end subroutine toptinit_user !``````````````````````````````````````````````````````````````````` subroutine sfcinit_user(n1, n2, n3, n4, n5, n6, igd, vnam, glat, glon, & TI_static, albd_static) !+ Where user may customize SBPS characteristics' values !~ 29.08.2012 (last modified; TIM) use univ_kind_defs_mod, only: f4 use main_memory_mod, only: a, initialization_type, TI_val, run_sfcinit_user implicit none !=== Argument declarations: integer, intent(IN) :: n1, n2, n3, n4, n5, n6, igd character(len=*), intent(IN) :: vnam real(f4), dimension(n2,n3), intent(IN) :: glat, glon real(f4), dimension(n2,n3,n5), intent(INOUT) :: TI_static, albd_static !== Local declarations: ! This subroutine is the intended location for a user to customize any ! sub/super surface (SBPS) characteristic array. It is called after all ! other types of initialization of this field, so this subroutine has the ! last word. By default the subroutine makes no change to the field. The ! commented lines below serve as a template for user-designed changes; ! Note that this routine is called for each grid separately, so attention to ! the current value of NG in this routine may be required by the user: ! !------------------------------------------------------------------------ ! *** IMPORTANT ***: Do not alter SBPSTMCD here, as it will not be used ! properly; instead, alter TISTAT appropriately. !------------------------------------------------------------------------ integer :: i, j, k, ip !=== Executable statements: if (.not.run_sfcinit_user) return ! Not desired, currently if (vnam == 'FILEDATA') then if (initialization_type == 2) then ! do ip=1,n5 ! Loop over all patches ! do j=1,n3 ! do i=1,n2 ! if(glat(i,j) < -79.9499) then ! South polar TI kludge ! TI_static(i,j,ip) = 100. ! elseif(glat(i,j) > 79.9499) then ! North polar TI kludge ! TI_static(i,j,ip) = TI_val ! endif ! enddo ! enddo ! enddo end if else if (vnam == 'NOFILE') then if (igd == 1) then !***** WHATEVER ****** end if end if end subroutine sfcinit_user !``````````````````````````````````````````````````````````````````` subroutine LES_iperturb(m1,m2,m3,thp,lev) !~ 19.11.2008 (last modified; TIM) !+ Slightly randomly perturb thermal LES init field (so not so homogeneous) use univ_kind_defs_mod, only: f4 implicit none !=== Argument declarations: integer, intent(IN) :: m1,m2,m3,lev real(f4), dimension(m1,m2,m3), intent(INOUT) :: thp !== Local declarations: real(f4), dimension(m2,m3) :: harvest !=== Executable statements: call random_number(harvest) ! Get a horizontal grid of numbers between 0 & 1 thp(lev,:,:) = thp(lev,:,:)+(harvest(:,:)/5._f4-0.1_f4) ! Add random perturbation ! to desired theta level (+/- 0.1 K) end subroutine LES_iperturb !``````````````````````````````````````````````````````````````````` subroutine tracer_init(n1,n2,n3,ng,tID,tracer,topo,rtgt,zt,thp,this_is_a_restart) !~ 19.11.2008 (last modified; TIM) ! 11.May.2016 (last modified; JPLA; including thp potential temperature variable) !+ Initialize some surface quantities (has last word) use univ_kind_defs_mod, only: f4 use univ_math_constants_mod, only: ZERO_f4 implicit none !=== Argument declarations: integer, intent(IN) :: n1,n2,n3,ng,tID real(f4), dimension(n2,n3), intent(IN) :: topo,rtgt real(f4), dimension(n1), intent(IN) :: zt logical, intent(IN) :: this_is_a_restart real(f4), dimension(n1,n2,n3), intent(INOUT) :: tracer,thp ! Is actually: a(isclpn(X,ng)) !== Local declarations: integer :: i,j,ie,ib,je,jb,k real(f4) :: az,z_agl !=== Executable statements: if(this_is_a_restart) then do j=1,n3 do i=1,n2 do k=1,n1 tracer(k,i,j) = ZERO_f4 if (ng <= 3 .and. tID == 3) tracer(k,i,j) = 1. enddo enddo enddo !tracer = 0. ! Make sure everything else is zeroed out if(ng == 4) then ! Only want to do this on grid 4 jb = 40 ! frames Gale crater je = 60 ! ib = 37 ! ie = 55 ! select case(tID) case (1) jb = 40 ! je = 60 ! ib = 37 ! ie = 55 ! do j=jb,je do i=ib,ie ! tracer if az <= 12km and z_AGL > 200m do k=1,n1 z_agl = zt(k)*rtgt(i,j) if(z_agl <= 500.) then tracer(k,i,j) = 1. endif enddo enddo enddo case (2) jb = 40 ! je = 60 ! ib = 37 ! ie = 55 ! do j=jb,je do i=ib,ie ! tracer if z_AGL > 500m and <= 2000m do k=1,n1 z_agl = zt(k)*rtgt(i,j) if(z_agl > 500. .and. z_agl <= 2000.) then tracer(k,i,j) = 1. endif enddo enddo enddo case (3) jb = 40 ! je = 60 ! ib = 37 ! ie = 55 ! do j=jb,je do i=ib,ie ! tracer if az <= 12km and z_AGL > 200m do k=1,n1 z_agl = zt(k)*rtgt(i,j) if(z_agl > 2000. .and. z_agl <= 20000.) then tracer(k,i,j) = 1. endif enddo enddo enddo end select end if else ! Just initialize to zero do j=1,n3 do i=1,n2 do k=1,n1 tracer(k,i,j) = ZERO_f4 if (ng <= 3 .and. tID == 3) tracer(k,i,j) = 1. enddo enddo enddo ! select case(tID) ! Operating on which tracer? ! !case(1) ! ! do j=1,n3 ! do i=1,n2 ! do k=1,n1 ! tracer(k,i,j) = 0.9532_f4 ! CO2 gas fraction ! enddo ! enddo ! enddo ! ! end select endif end subroutine tracer_init