module tuv_bndry implicit none ! ! Output: Lower boundary conditions for t,u,v: ! ubnd,vbnd equivalent to un(itc),vn(itc) at the LB ! real,allocatable,dimension(:,:) :: ! (lon0:lon1,lat0:lat1) | tbnd, ! tn lower boundary (not in use as of Mar, 2004) | ubnd, ! un lower boundary (used in duv) | vbnd ! vn lower boundary (used in duv) contains !----------------------------------------------------------------------- subroutine tuvbnd(tn,tn_nm,un,un_nm,vn,vn_nm,w,z,barm,cp, | hdt,hdu,hdv,cool_imp,cool_exp,ulbc,ulbc_nm,vlbc,vlbc_nm, | lev0,lev1,lon0,lon1,lat0,lat1) ! ! Calculate lower boundaries of u,v ! stored in ubnd,vbnd ! equivalent to LB for un(itc),vn(itc) ! needs LB fields: ! ulbc,vlbc,ulbc_nm,vlbc_nm ! which are equivalent to ! un(itp),vn(itp),un_nm(itp),vn_nm(itp) ! use params_module,only: dz,nlonp4 use cons_module,only: shapiro,dtx2inv,expz,gask,racs,grav,re_inv, | dtr,cor,tanphi=>tn use init_module,only: glat ! ! 3/9/06 btf: use raykk from gwsource instead of rayk+eqfric. ! raykk from mgw is rayk + gswm rfricu (see sub gwsource in mgw.F) ! use mgw,only: raykk ! rayleigh friction (lev,lon,lat) use qrj_module,only: qtotal ! ! Routine to add fields to secondary histories: use addfld_module,only: addfld ! ! Input args: integer,intent(in) :: lev0,lev1,lon0,lon1,lat0,lat1 real,dimension(lev0:lev1,lon0-2:lon1+2,lat0-2:lat1+2),intent(in):: | tn, ! temperature (f(:,nj+nt )) not used | tn_nm, ! temperature, time n-1 (f(:,nj+ntnm)) not used | un, ! zonal wind (f(:,nj+nu)) | un_nm, ! zonal wind, time n-1 (f(:,nj+nunm)) | vn, ! meridional wind (f(:,nj+nv)) | vn_nm, ! meridional wind, time n-1 (f(:,nj+nvnm)) | w, ! omega (for vertical velocity) (f(:,nj+nw)) | z, ! geopotential height (f(:,nj+nz)) | barm, ! mean mol weight (f(:,nj+nms)) | cp, ! specific heat (f(:,ncp)) | hdt, ! horizontal diffusion of T (f(:,nqdh)) | hdu, ! horizontal diffusion of U (f(:,nflh)) | hdv, ! horizontal diffusion of V (f(:,nfph)) | cool_exp, ! explicit cooling (radcool.F) (f(:,nwte)) | cool_imp ! implicit cooling (radcool.F) (f(:,nwti)) real,dimension(lon0-2:lon1+2,lat0-2:lat1+2),intent(in):: | ulbc, ! LB zonal wind | ulbc_nm, ! LB zonal wind, time n-1 | vlbc, ! LB meridional wind | vlbc_nm ! LB meridional wind, time n-1 ! ! Local: integer :: k,i,lat,lonbeg,lonend integer ::i0,i1,nk,nkm1,ni ! for addfsech real,dimension(lon0:lon1) :: | advect, ! horizontal advection (t) | advecu, ! horizontal advection (u) | advecv, ! horizontal advection (v) | cpmbar, ! cp*barm | shapiro2, ! t shapiro stage 2 | ushapiro2, ! u shapiro stage 2 | vshapiro2, ! v shapiro stage 2 | qterm1, ! 1./(2*dt)*t(n-1)-advect + q/cp | qterm2, ! 1./(2*dt)+cool_imp+w*r/(cp*mbar) | heatlbc, ! total heating at lower boundary | a11,a12,a21,a22, ! final terms for ubnd,vbnd | deta ! a11*a22-a12*a21 real :: shapiro1(lon0-2:lon1+2) ! t shapiro stage 1 real :: ushapiro1(lon0-2:lon1+2) ! u shapiro stage 1 real :: vshapiro1(lon0-2:lon1+2) ! v shapiro stage 1 real,dimension(1,lon0:lon1) :: | zl,zp ! derivatives of z at bottom boundary real :: rlat,eqfric ! for rayleigh friction real,dimension(lev0:lev1,lon0:lon1) :: ! diag only (addfsech calls) | advect_diag,heatlbc_diag,qterm1_diag, | qterm2_diag,smooth1,smooth2,tbnd_diag, | usmooth1,vsmooth1,usmooth2,vsmooth2 real,dimension(lon0:lon1,lat0:lat1) :: ulbc_diag,vlbc_diag, | ulbcnm_diag,vlbcnm_diag,advecu_diag,advecv_diag, | zl_diag,zp_diag,z_diag,advecub_diag,advecvb_diag, | usmo_diag,vsmo_diag,usmo2_diag,vsmo2_diag,advecu2_diag, | advecv2_diag,a11_diag,a22_diag,a21_diag,a12_diag,delta_diag, | vbnd_diag,ubnd_diag real,parameter :: hlat = 40./57.295 ! for rayleigh friction ! ! For addfsech calls: i0 = lon0 ; i1 = lon1 ; nk = lev1-lev0+1 ; nkm1 = nk-1 ni = i1-i0+1 ! lonbeg = lon0 if (lon0==1) lonbeg = 3 lonend = lon1 if (lon1==nlonp4) lonend = lon1-2 ! ! Latitude scan: do lat=lat0,lat1 ! ! Equatorial rayleigh friction (tunable) (s.a. hlat above): rlat = glat(lat)*dtr eqfric = 2.E-5*exp(-(rlat/hlat)**2) ! eqfric = 1.e-20 ! ! Derive boundary conditions for U,V. ! ! Horizontal advection of u,v at lower boundary. ! (u,v bottom boundary in ulbc,vlbc): ! call advecl(ulbc(:,lat-2:lat+2),advecu,lon0,lon1,lat) call advecl(vlbc(:,lat-2:lat+2),advecv,lon0,lon1,lat) ! ! Save to secondary history: advecu_diag(:,lat)= advecu(:) advecv_diag(:,lat)= advecv(:) ! ! Longitudinal and latitudinal derivatives of geopotential bottom boundary: call dldp(z(lev0:lev0,:,lat-2:lat+2),zl,zp,lev0,lev0,lon0,lon1, | lat) ! ! zl = G/(RAD*COS(PHI))*DZ/DLAMDA ! zp = G/RA*DZ/DPHI ! (formerly sub glpl in tgcm24) ! do i=lonbeg,lonend zl(1,i) = zl(1,i)*grav*racs(lat) zp(1,i) = zp(1,i)*grav*re_inv enddo ! i=lonbeg,lonend ! ! 5/19/04 btf: May need to do periodic points on zl,zp here. ! Try setting to zero for now to avoid FP exception: ! if (lon0==1) then zl(1,1:2) = 0. zp(1,1:2) = 0. endif if (lon1==nlonp4) then zl(1,nlonp4-1:nlonp4) = 0. zp(1,nlonp4-1:nlonp4) = 0. endif ! ! Save to secondary history: ! zl_diag(:,lat) = zl(1,:) ! zp_diag(:,lat) = zp(1,:) ! z_diag(:,lat) = z(lev0,:,lat) ! call addfld('ZL_LBC ',' ',' ',zl_diag, ! | 'lev',lev0,lev1,'lon',i0,i1,lat) ! call addfld('ZP_LBC ',' ',' ',zp_diag, ! | 'lev',lev0,lev1,'lon',i0,i1,lat) ! ! Add in vertical advection (bottom boundary is in ulbc, vlbc): do i=lon0,lon1 advecu(i) = advecu(i)+2.*w(lev0,i,lat)* | (un(lev0,i,lat)-ulbc(i,lat))/dz advecv(i) = advecv(i)+2.*w(lev0,i,lat)* | (vn(lev0,i,lat)-vlbc(i,lat))/dz ! ! Add horizontal derivatives of geopotential: advecu(i) = advecu(i)+zl(1,i) advecv(i) = advecv(i)+zp(1,i) ! ! Fourth order horizontal diffusion: ! hdu,v: zonal,meridional diffusion from hdif3 (hdif.F) ! gwu,v: gravity wave tendencies from mgw (mgw.F). ! (gw may be turned off, check calc_gw) ! ! advecu(i) = advecu(i)- ! | 1.5*(hdu(lev0 ,i,lat)+gwu(lev0, i,lat))+ ! | 0.5*(hdu(lev0+1,i,lat)+gwu(lev0+1,i,lat)) ! advecv(i) = advecv(i)- ! | 1.5*(hdv(lev0 ,i,lat)+gwv(lev0, i,lat))+ ! | 0.5*(hdv(lev0+1,i,lat)+gwv(lev0+1,i,lat)) advecu(i) = advecu(i)- | 1.5*(hdu(lev0 ,i,lat))+ | 0.5*(hdu(lev0+1,i,lat)) advecv(i) = advecv(i)- | 1.5*(hdv(lev0 ,i,lat))+ | 0.5*(hdv(lev0+1,i,lat)) enddo ! i=lon0,lon1 ! Save to secondary history: advecub_diag(:,lat) = advecu(:) advecvb_diag(:,lat) = advecv(:) ! ! Shapiro smoother for u,v, stage 1 (lower boundary of u,v_nm is in top slot): do i=lon0-2,lon1+2 ushapiro1(i) = ulbc_nm(i,lat)-shapiro* ! s6 | (ulbc_nm(i,lat+2)+ulbc_nm(i,lat-2)- | 4.*(ulbc_nm(i,lat+1)+ulbc_nm(i,lat-1))+ | 6.* ulbc_nm(i,lat)) vshapiro1(i) = vlbc_nm(i,lat)-shapiro* ! s7 | (vlbc_nm(i,lat+2)+vlbc_nm(i,lat-2)- | 4.*(vlbc_nm(i,lat+1)+vlbc_nm(i,lat-1))+ | 6.* vlbc_nm(i,lat)) enddo ! i=lon0-2,lon1+2 ! Save to secondary history: usmo_diag(:,lat) = ushapiro1(:) vsmo_diag(:,lat) = vshapiro1(:) ! ! Shapiro smoother for u,v, stage 2: ushapiro2 = 0. ! whole array vshapiro2 = 0. ! whole array do i=lonbeg,lonend ushapiro2(i)= ushapiro1(i)-shapiro* ! s8 | (ushapiro1(i+2)+ushapiro1(i-2)- | 4.*(ushapiro1(i+1)+ushapiro1(i-1))+6.*ushapiro1(i)) vshapiro2(i)= vshapiro1(i)-shapiro* ! s9 | (vshapiro1(i+2)+vshapiro1(i-2)- | 4.*(vshapiro1(i+1)+vshapiro1(i-1))+6.*vshapiro1(i)) enddo ! i=lonbeg,lonend ! Save to secondary history: usmo2_diag(:,lat) = ushapiro2(:) vsmo2_diag(:,lat) = vshapiro2(:) ! tgcm24 adds periodic points to ushapiro2,vspairo2 at this point. ! ! advecu = U(N-1)/(2.*DT)-advecu ! advecv = V(N-1)/(2.*DT)-advecv ! do i=lon0,lon1 advecu(i) = ushapiro2(i)*dtx2inv-advecu(i) ! s1 advecv(i) = vshapiro2(i)*dtx2inv-advecv(i) ! s2 enddo ! i=lon0,lon1 ! Save to secondary history: advecu2_diag(:,lat) = advecu(:) advecv2_diag(:,lat) = advecv(:) ! ! Final terms for ubnd,vbnd (bottom boundary in top slot): ! 3/9/06 btf: use raykk from sub gwsource (raykk=rayk+gswm_rayfric) ! do i=lon0,lon1 ! a11(i) = dtx2inv+(rayk(lev0)+eqfric)**1.5/ ! s3 ! | (rayk(lev0+1)+eqfric)**0.5 ! a12(i) = -(cor(lat)+ulbc(i,lat)*tanphi(lat)*re_inv) ! s4 ! a21(i) = -a12(i) ! s5 ! a22(i) = dtx2inv+(rayk(lev0)+eqfric)**1.5/ ! s6 ! | (rayk(lev0+1)+eqfric)**0.5 a11(i) = dtx2inv+(raykk(lev0,i,lat)+eqfric)**1.5/ ! s3 | (raykk(lev0+1,i,lat)+eqfric)**0.5 a12(i) = -(cor(lat)+ulbc(i,lat)*tanphi(lat)*re_inv) ! s4 a21(i) = -a12(i) ! s5 a22(i) = dtx2inv+(raykk(lev0,i,lat)+eqfric)**1.5/ ! s6 | (raykk(lev0+1,i,lat)+eqfric)**0.5 ! ! deta = DET(A) = A11*A22-A12*A21 deta(i) = a11(i)*a22(i)-a12(i)*a21(i) ! s7 ! ! Solve for ubnd,vbnd: ubnd(i,lat)= ( a22(i)*advecu(i)-a12(i)*advecv(i))/deta(i) vbnd(i,lat)= (-a21(i)*advecu(i)+a11(i)*advecv(i))/deta(i) ! ulbc_diag(i,lat) = ulbc(i,lat) vlbc_diag(i,lat) = vlbc(i,lat) ulbcnm_diag(i,lat) = ulbc_nm(i,lat) vlbcnm_diag(i,lat) = vlbc_nm(i,lat) ! a11_diag(i,lat) = a11(i) ! a22_diag(i,lat) = a22(i) ! a12_diag(i,lat) = a12(i) ! a21_diag(i,lat) = a21(i) ! delta_diag(i,lat) = deta(i) enddo ! i=lon0,lon1 ! ! call addfld('A11 ',' ',' ',a11_diag, ! | 'lon',lon0,lon1,'lat',lat0,lat1,0) ! call addfld('A12 ',' ',' ',a12_diag, ! | 'lon',lon0,lon1,'lat',lat0,lat1,0) ! call addfld('A21 ',' ',' ',a21_diag, ! | 'lon',lon0,lon1,'lat',lat0,lat1,0) ! call addfld('A22 ',' ',' ',a22_diag, ! | 'lon',lon0,lon1,'lat',lat0,lat1,0) ! ! Save u,v lower boundaries to secondary history: ! call addfld('UBND ',' ',' ',ubnd_diag, ! | 'lev',lev0,lev1,'lon',i0,i1,lat) ! call addfld('VBND ',' ',' ',vbnd_diag, ! | 'lev',lev0,lev1,'lon',i0,i1,lat) ! ! End latitude scan: enddo ! lat=lat0,lat1 ! call addfld('UBND_TUV','UBND_TUV','m/s', ! | ubnd,'lon',lon0,lon1,'lat',lat0,lat1,0) ! call addfld('VBND_TUV','VBND_TUV','m/s', ! | vbnd,'lon',lon0,lon1,'lat',lat0,lat1,0) ! call addfld('ULBC_TUV','ULBC_TUV','m/s', ! | ulbc_diag,'lon',lon0,lon1,'lat',lat0,lat1,0) ! call addfld('VLBC_TUV','VLBC_TUV','m/s', ! | vlbc_diag,'lon',lon0,lon1,'lat',lat0,lat1,0) ! call addfld('ULBCNM_TUV','ULBCNM_TUV','m/s', ! | ulbcnm_diag,'lon',lon0,lon1,'lat',lat0,lat1,0) ! call addfld('VLBCNM_TUV','VLBCNM_TUV','m/s', ! | vlbcnm_diag,'lon',lon0,lon1,'lat',lat0,lat1,0) ! call addfld('ADVU_TUV','ADVU_TUV',' ', ! | advecu_diag,'lon',lon0,lon1,'lat',lat0,lat1,0) ! call addfld('ADVV_TUV','ADVV_TUV',' ', ! | advecv_diag,'lon',lon0,lon1,'lat',lat0,lat1,0) ! call addfld('ZP_TUV','ZP_TUV',' ', ! | zp_diag,'lon',lon0,lon1,'lat',lat0,lat1,0) ! call addfld('ZL_TUV','ZL_TUV',' ', ! | zl_diag,'lon',lon0,lon1,'lat',lat0,lat1,0) ! call addfld('Z_TUV','Z_TUV',' ', ! | z_diag,'lon',lon0,lon1,'lat',lat0,lat1,0) ! ! call addfld('ADUB_TUV','ADUB_TUV',' ', ! | advecub_diag,'lon',lon0,lon1,'lat',lat0,lat1,0) ! call addfld('ADVB_TUV','ADVB_TUV',' ', ! | advecvb_diag,'lon',lon0,lon1,'lat',lat0,lat1,0) ! call addfld('USMO_TUV','USMO_TUV',' ', ! | usmo_diag ,'lon',lon0,lon1,'lat',lat0,lat1,0) ! call addfld('VSMO_TUV','VSMO_TUV',' ', ! | vsmo_diag,'lon',lon0,lon1,'lat',lat0,lat1,0) ! call addfld('USMO2_TUV','USMO2_TUV',' ', ! | usmo2_diag,'lon',lon0,lon1,'lat',lat0,lat1,0) ! call addfld('VSMO2_TUV','VSMO2_TUV',' ', ! | vsmo2_diag,'lon',lon0,lon1,'lat',lat0,lat1,0) ! call addfld('ADU2_TUV','ADU2_TUV',' ', ! | advecu2_diag,'lon',lon0,lon1,'lat',lat0,lat1,0) ! call addfld('ADV2_TUV','ADV2_TUV',' ', ! | advecv2_diag,'lon',lon0,lon1,'lat',lat0,lat1,0) ! call addfld('A11_TUV','A11_TUV',' ', ! | a11_diag,'lon',lon0,lon1,'lat',lat0,lat1,0) ! call addfld('A22_TUV','A22_TUV',' ', ! | a22_diag,'lon',lon0,lon1,'lat',lat0,lat1,0) ! call addfld('A12_TUV','A12_TUV',' ', ! | a12_diag,'lon',lon0,lon1,'lat',lat0,lat1,0) ! call addfld('A21_TUV','A21_TUV',' ', ! | a21_diag,'lon',lon0,lon1,'lat',lat0,lat1,0) ! call addfld('DELTA_TUV','DELTA_TUV',' ', ! | delta_diag,'lon',lon0,lon1,'lat',lat0,lat1,0) ! end subroutine tuvbnd !----------------------------------------------------------------------- subroutine dldp(f,xl,xp,lev0,lev1,lon0,lon1,lat) ! ! Derivatives of f are returned in xl,xp. Bottom boundary only. ! use cons_module,only: dlamda_2div3, dlamda_1div12, | dphi_2div3, dphi_1div12 use params_module,only: nlonp4 ! ! Args: implicit none integer,intent(in) :: lev0,lev1,lon0,lon1,lat real,dimension(lev0:lev0,lon0-2:lon1+2,lat-2:lat+2), | intent(in) :: f real,dimension(lev0:lev0,lon0:lon1),intent(out) :: xl,xp ! ! Local: integer :: i,lonbeg,lonend ! lonbeg = lon0 if (lon0==1) lonbeg = 3 lonend = lon1 if (lon1==nlonp4) lonend = lon1-2 ! do i=lonbeg,lonend ! ! Longitudinal derivative: xl(lev0,i) = dlamda_2div3 *(f(lev0,i+1,lat)-f(lev0,i-1,lat))- | dlamda_1div12*(f(lev0,i+2,lat)-f(lev0,i-2,lat)) ! ! Latitudinal derivative: xp(lev0,i) = dphi_2div3 *(f(lev0,i,lat+1)-f(lev0,i,lat-1))- | dphi_1div12*(f(lev0,i,lat+2)-f(lev0,i,lat-2)) ! write(6,"('dldp: lat=',i3,' i=',i3,' xl=',/,(6e12.4))") ! | lat,i,xl(:,i) enddo ! i=lonbeg,lonend end subroutine dldp !----------------------------------------------------------------------- ! subroutine xferdiag(fin,idin, fout,idout1,idout2) ! ! Transfer 1d array fin to 2d fout. First dimension of fout (usually zp) ! will be redundant. ! This is used to transfer (lon) arrays to diagnostic (lev,lon) arrays ! for addfsech. ! ! 8/8/06 btf: This routine is only necessary if one needs to post-process ! secondary history fields with tgcmproc_f90. Currently, ! sub addfld will write global 2d fields, and tgcmproc_idl ! or ncl will read and display them. ! ! integer,intent(in) :: idin,idout1,idout2 ! real,intent(in) :: fin(idin) ! real,intent(out) :: fout(idout1,idout2) ! ! integer :: k ! do k=1,idout1 ! fout(k,:) = fin(:) ! enddo ! end subroutine xferdiag !----------------------------------------------------------------------- subroutine alloc_tuvbnd(lon0,lon1,lat0,lat1) ! ! Args: integer,intent(in) :: lon0,lon1,lat0,lat1 ! ! Local: integer :: istat ! ! Allocate lower boundaries for t,u,v: ! write(6,*) 'alloc_tuvbnd:',lon0,lon1,lat0,lat1 allocate(tbnd(lon0:lon1,lat0:lat1),stat=istat) if (istat /= 0) write(6,"(/,'>>> alloc_tuvbnd: error allocating', | ' tbnd: stat=',i3)") istat ! allocate(ubnd(lon0:lon1,lat0:lat1),stat=istat) if (istat /= 0) write(6,"(/,'>>> alloc_tuvbnd: error allocating', | ' ubnd: stat=',i3)") istat ! allocate(vbnd(lon0:lon1,lat0:lat1),stat=istat) if (istat /= 0) write(6,"(/,'>>> alloc_tuvbnd: error allocating', | ' vbnd: stat=',i3)") istat end subroutine alloc_tuvbnd !----------------------------------------------------------------------- end module tuv_bndry