#define MODHS 1 subroutine tphysidl(ztodt, state, tend) !----------------------------------------------------------------------- ! ! Purpose: ! algorithm 1: Held/Suarez IDEALIZED physics ! algorithm 2: Held/Suarez IDEALIZED physics (Williamson modified stratosphere ! algorithm 3: Held/Suarez IDEALIZED physics (Lin/Williamson modified strato/meso-sphere ! algorithm 4: Boer/Denis IDEALIZED physics ! ! Author: J. Olson ! !----------------------------------------------------------------------- use shr_kind_mod, only: r8 => shr_kind_r8 use ppgrid, only: pcols, pver,begchunk,endchunk use phys_grid, only: get_rlat_all_p, get_rlon_all_p,get_lat_all_p,get_lon_all_p,gather_chunk_to_field,scatter_field_to_chunk use spmd_utils, only: masterproc use physics_types, only: physics_state, physics_tend, physics_ptend, & physics_ptend_init, physics_update use physconst, only: gravit, cappa, rair, cpair use abortutils, only: endrun use ref_pres, only: pref_mid_norm, psurf_ref use cam_history, only: outfld,addfld,phys_decomp use cam_logfile, only: iulog use time_manager, only: get_nstep,is_first_step,is_first_restart_step use check_energy, only: check_energy_chng implicit none ! ! Input arguments ! real(r8), intent(in) :: ztodt ! Two times model timestep (2 delta-t) type(physics_state), intent(inout) :: state ! ! Output arguments ! type(physics_tend ), intent(inout) :: tend ! !---------------------------Local workspace----------------------------- ! type(physics_ptend) :: ptend integer :: lchnk ! chunk identifier integer :: ncol ! number of atmospheric columns integer :: nstep ! timestep number real(r8) clat(pcols) ! latitudes(radians) for columns real(r8) clon(pcols) ! wykang lon for columns integer ilat(pcols) ! wykang latitudes index for columns integer ilon(pcols) ! wykang longitude index for columns real(r8) pmid(pcols,pver) ! mid-point pressure integer i,k ! Longitude, level indices integer linenum integer pre real(r8) k_top_damp ! wanying added damp to model top real(r8) lapsewy ! lapse rate integer free ! no relaxation at all integer ndim integer tforcing ! have tforcing? integer tforcingtype ! additive gaussian or sine integer FT ! temperature forcing integer FUV !wind forcing integer flag ! idealized experiment type integer clima ! use climatology integer relaxwind ! if relax wind fields to model? integer damp ! damping flag integer topo ! if there is topography integer tightrelax integer checkread integer solverlx integer nowind integer fricuv2clima integer moving character*32:: ffname character*32:: model character*32:: season character*32:: Forclima character*64:: IDEALRLX character*32:: toponame character*64:: forcmodel character*32:: modified character*32:: FTUVname real(r8) tmp ! temporary real(r8),save:: tau_t(pcols,pver,55) real(r8),save:: tau_uv(pcols,pver,55) ! timescale for u restoring real(r8) l1 real(r8) l2 real(r8) ump ! timescale, u friction real(r8) kf ! 1./efolding_time for wind dissipation real(r8) ka ! 1./efolding_time for temperature diss. real(r8) kaa ! 1./efolding_time for temperature diss. real(r8) ks ! 1./efolding_time for temperature diss. real(r8) kv ! 1./efolding_time (normalized) for wind real(r8) kt ! 1./efolding_time for temperature diss. real(r8) alp integer read_in_tmp(pver) real(r8),save :: TRLX(pcols,pver,55) real(r8) TERR(pcols,pver) real(r8) UERR(pcols,pver) real(r8) VERR(pcols,pver) real(r8) heatforce(pcols,pver) real(r8),save :: t_rlx(pcols,pver,55) !"radiative equilibrium" T real(r8),save :: u_rlx(pcols,pver,55) real(r8),save :: v_rlx(pcols,pver,55) real(r8),save :: t_forc(pcols,pver,55) real(r8),save :: u_forc(pcols,pver,55) real(r8),save :: v_forc(pcols,pver,55) real(r8) Amp !Amplitude of forcing real(r8) lon1 real(r8) Fperiod real(r8) lonshift real(r8) lonL real(r8) lonR real(r8) latS real(r8) latN real(r8) Frclat real(r8) Frclon real(r8) Frchigh real(r8) Frclow real(r8) Frclev real(r8) sigmalat real(r8) sigmalon real(r8) Frcthick real(r8) TimeVary real(r8) step0 logical BLK logical Blklon real(r8) Blklon1 real(r8) Blklon2 logical Blklat real(r8) Blklat1 real(r8) Blklat2 logical Blklev real(r8) Blklev1 real(r8) Blklev2 real(r8) Blkloncenter real(r8) Blklonsigma real(r8) Blklatcenter real(r8) Blklatsigma real(r8) Blklevcenter real(r8) Blklevsigma real(r8) Ttend_blk real(r8) Utend_blk real(r8) Vtend_blk real(r8) weight real(r8) trefc ! used in calc of "radiative equilibrium" T real(r8) cossq(pcols) ! coslat**2 real(r8) cossqsq(pcols) ! coslat**4 real(r8) sinsq(pcols) ! sinlat**2 real(r8) onemsig ! 1. - sigma_reference real(r8) efoldf ! efolding time for wind dissipation real(r8) efolda ! efolding time for T dissipation real(r8) efoldaa ! efolding time for T dissipation real(r8) efolds ! efolding time for T dissipation real(r8) efold_strat ! efolding time for T dissipation in Strat real(r8) efold_meso ! efolding time for T dissipation in Meso real(r8) efoldv ! efolding time for wind dissipation real(r8) p_infint ! effective top of model real(r8) constw ! constant real(r8) lapsew ! lapse rate real(r8) p0strat ! threshold pressure real(r8) phi0 ! threshold latitude real(r8) dphi0 ! del-latitude real(r8) a0 ! coefficient real(r8) aeq ! 100 mb real(r8) apole ! 2 mb real(r8) pi ! 3.14159... real(r8) coslat(pcols) ! cosine(latitude) real(r8) acoslat ! abs(acos(coslat)) real(r8) constc ! constant real(r8) lapsec ! lapse rate real(r8) lapse ! lapse rate real(r8) h0 ! scale height (7 km) real(r8) sigmab ! threshold sigma level real(r8) pressmb ! model pressure in mb real(r8) t00 ! minimum reference temperature integer idlflag ! Flag to choose which idealized physics real(r8) :: zero(pcols) logical :: firststep logical :: firstrststep logical :: lonlim logical :: latlim logical :: Hall !CONTROL !----------------------------------------------------------------------- !CONTROL idlflag =1 free =0 pi = 4._r8*atan(1._r8) !no relaxation to anything tforcing=1 !1:add forcing 0:no tforcing tforcingtype=1 !0:additive gaussian >=1:sine wave number. Only work for tforcing=1,2. If !this is active the sigmalon is neglicted Amp=1._r8/86400._r8 !forcing strength in unit of K/40days. Only work for tforcing=1,2 Frclat=0._r8 Frclon=135._r8 Frclow=900._r8*100._r8 Frchigh=100._r8*100._r8 sigmalat=5._r8 sigmalon=20._r8 Frcthick=sqrt(Frclow/Frchigh) Frclev=sqrt(Frchigh*Frclow) moving=1 ! 1: if add moving forcing. Only work for tforcing=1,2 Fperiod=40._r8 ! Period of moving force. lonlim=.true. lonshift=120._r8/180._r8*pi ! only work for sine heating lonL=60._r8 lonR=180._r8 if (.NOT.lonlim) then lonL=0._r8 lonR=360._r8 endif !limit forcing to some longitude range latlim=.false. latN=10._r8 latS=-10._r8 if (.NOT.latlim) then latN=90._r8 latS=-90._r8 endif TimeVary=0._r8 !the time for forcing to increase from 0 to Amp, unit in year. 0._r8 turn off !negative value means standing periodical wave, using period of Fperiod. step0=0._r8*365._r8*48._r8 !used only when TimeVary!=0 damp=1 !1:have damping on the top !2:add stronger top_damp which also depends on lat !3:change tau_t directly k_top_damp = 1._r8/(30._r8*86400._r8) !damping strength !!!!!!!!!!!!! !Following only used in idlflag==1 !!!!!!!!!!!!! !!!general tightrelax=0 !only in idlflag=1 ! turn on to reduce relaxation time scale to 1day.Better !turn off! lapsewy = 6.5_r8/cpair*1000._r8 !only in idlflag=1 flag=3 !only in idlflag=1 !rest_atm: 1 solidbody:2 readfile:3 topo=1 !only in idlflag=1 !have topography or not relaxwind=0 !only in idlflag=1 ! relax not only temperature but also wind fields !!!iteration solverlx=0 !only in idlflag=1 ! switch to calculate relaxation temperature field, currently no lonavg is ! taken no matter ndim=2 or 3 alp=1000._r8 !only in idlflag=1 ! strength of relax in updating relax fields Hall=.false. ! output T_ U_ V_ now !!!these used for reading background and forcing from files model="waccmctl" season="Jan" !!forcing to maintain climatology (Hall) FT=1 FUV=1 !0:no FUV/FT 1: FUV/FT from file(only in idlflag=1) forcmodel="IDEAL_topo_wa1jan_noEfix_damp30d" !"topo_spdjf_noEfix_damp30d" FTUVname="waccm_Jan35" !model_season,"notopo_waccmdjf2d" Forclima="" !"":extra forced by -\bar{\Psi_1-\Psi_0}/dt !"clima":extra forced by -(\bar{\Psi_1}-\bar{\Psi_0}) !only in idlflag=1 FT=1 !!background clima=1 !1:use GCM climatology (only valid when flag==3) 0:use RLX result IDEALRLX="" !Case name of RLX run for this model. Only when clima==0 nowind=1 !have not relaxed wind in RLX. Only when clima==0 ndim=3 !only in idlflag=1 modified="35" ! "modified1_damp30noEfix","modified2_damp30noEfix","modified1","modified2","_x4CO2","_x4CO2s" ! set to modifed when using spcamctl to change the relaxation temperature ! compensate the difference in checkingforc_daily1_0 ! if not set to "" ! for x4CO2 experiment: "_x4CO2s"(symmetric JJA SH),"_x4CO2longavgs"(JJA SH ! for NH, DJF SH for SH, lonavg) !!check checkread=1 !only in idlflag=1 ! whether to check read file. Need to change the chk !output's name!! !!friction uv to what fricuv2clima=1 !set to 1 when using step0-IC as forcing. Only when FT=FUV=1 !!block the wave BLK=.False. Blklon=.TRUE. Blklon1=90._r8/180._r8*pi Blklon2=180._r8/180._r8*pi Blkloncenter=(Blklon1+Blklon2)/2._r8 Blklonsigma=(Blklon2-Blklon1)/4._r8 Blklat=.TRUE. Blklat1=10._r8/180._r8*pi Blklat2=35._r8/180._r8*pi Blklatcenter=(Blklat1+Blklat2)/2._r8 Blklatsigma=(Blklat2-Blklat1)/4._r8 Blklev=.TRUE. Blklev1=1000._r8 Blklev2=20000._r8 Blklevcenter=(log(Blklev1) + log(Blklev2))/2._r8 Blklevsigma=(log(Blklev2) - log(Blklev1))/4._r8 !------------------------------------------------------------------------ if (topo==1) then toponame="topo_" else toponame="" endif lchnk = state%lchnk ncol = state%ncol nstep = get_nstep() zero = 0._r8 firststep=is_first_step() firstrststep= is_first_restart_step() firststep=firststep.OR.firstrststep call get_rlat_all_p(lchnk, ncol, clat) call get_rlon_all_p(lchnk, ncol, clon) call get_lat_all_p(lchnk, ncol, ilat) call get_lon_all_p(lchnk, ncol, ilon) ! get lat and lon indices for columns in this chunck. do i=1,ncol coslat (i) = cos(clat(i)) sinsq (i) = sin(clat(i))*sin(clat(i)) cossq (i) = coslat(i)*coslat(i) cossqsq(i) = cossq (i)*cossq (i) enddo do k = 1, pver do i = 1, ncol pmid(i,k) = state%pmid(i,k) enddo enddo ! initialize individual parameterization tendencies call physics_ptend_init(ptend, state%psetcols, 'tphysidl', ls=.true., lu=.true., lv=.true.) !write(iulog,*) 'pver',pver,'ncol',ncol,'pcols',pcols if (free ==1) then do k=1,pver do i=1,pcols ptend%u(i,k) = 0._r8 ptend%v(i,k) = 0._r8 ptend%s(i,k) = 0._r8 enddo enddo elseif (idlflag == 1 ) then ! !----------------------------------------------------------------------- ! ! Held/Suarez IDEALIZED physics algorithm: ! ! Held, I. M., and M. J. Suarez, 1994: A proposal for the ! intercomparison of the dynamical cores of atmospheric general ! circulation models. ! Bulletin of the Amer. Meteor. Soc., vol. 75, pp. 1825-1830. ! !----------------------------------------------------------------------- ! ! Compute idealized radiative heating rates (as dry static energy) ! efoldf = 1._r8 efolda = 40._r8 !default efolds = 4._r8 !default if (tightrelax==1) then efolda = 1._r8 efolds = 0.1_r8 endif sigmab = 0.7_r8 t00 = 200._r8 ! onemsig = 1._r8 - sigmab ! ka = 1._r8/(86400._r8*efolda) ks = 1._r8/(86400._r8*efolds) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !add fields !call addfld('TERR ','K ',pver,'I','Temperature difference from target',phys_decomp) !call addfld('RLXCT ','/s ',pver,'I','Relax strength for T',phys_decomp) !if (firststep) then !initialize !tau_t=0._r8 !TRLX=0._r8 !TERR=0._r8 !!write(iulog,*) 'Have timestep 0!!!!!!!!!!===wykang' !endif if(firststep) then do k=1,pver if (pref_mid_norm(k) > sigmab) then do i=1,ncol kt = ka + (ks - ka)*cossqsq(i)*(pref_mid_norm(k) - sigmab)/onemsig #ifdef MODHS tau_t(i,k,lchnk-begchunk+1) = kt/(1._r8+ ztodt*kt) !implicit integration for T, using T^(n+1) instead of T^(n-1) on the rhs of T !eq. #else tau_t(i,k,lchnk-begchunk+1) = kt #endif if (flag==0) then trefc = 315._r8 - 60._r8*sinsq(i) t_rlx(i,k,lchnk-begchunk+1) = (trefc - 10._r8*cossq(i)*log((pmid(i,k)/psurf_ref)))*(pmid(i,k)/psurf_ref)**cappa t_rlx(i,k,lchnk-begchunk+1) = max(t00,t_rlx(i,k,lchnk-begchunk+1)) elseif (flag==2) then !superrotation !trefc=300._r8 - 48.7_r8*sinsq(i)*(pmid(i,k)/psurf_ref)**4.3_r8 trefc=300._r8 - 120._r8*sinsq(i)*(pmid(i,k)/psurf_ref)**4.3_r8 t_rlx(i,k,lchnk-begchunk+1)=trefc +lapsewy*10._r8/log(10._r8)*log((pmid(i,k)/psurf_ref)) elseif (flag==1) then !rest atmosphere trefc=285._r8 t_rlx(i,k,lchnk-begchunk+1)=trefc+lapsewy*10._r8/log(10._r8)*log((pmid(i,k)/psurf_ref)) endif enddo else do i=1,ncol #ifdef MODHS tau_t(i,k,lchnk-begchunk+1) = ka/(1._r8+ ztodt*ka) #else tau_t(i,k,lchnk-begchunk+1) = ka #endif if (damp==3) then if (k<4) then tau_t(i,k,lchnk-begchunk+1)=(ka*(pmid(i,k)-pmid(i,1))+k_top_damp*(pmid(i,4)-pmid(i,k)))/(pmid(i,4)-pmid(i,1)) #ifdef MODHS tau_t(i,k,lchnk-begchunk+1) = tau_t(i,k,lchnk-begchunk+1)/(1._r8+ztodt*tau_t(i,k,lchnk-begchunk+1)) #endif endif endif if (flag==0) then trefc = 315._r8 - 60._r8*sinsq(i) t_rlx(i,k,lchnk-begchunk+1) = (trefc - 10._r8*cossq(i)*log((pmid(i,k)/psurf_ref)))*(pmid(i,k)/psurf_ref)**cappa t_rlx(i,k,lchnk-begchunk+1) = max(t00,t_rlx(i,k,lchnk-begchunk+1)) elseif (flag==2) then !superrotation !trefc=300._r8 - 48.7_r8*sinsq(i)*(pmid(i,k)/psurf_ref)**4.3_r8 trefc=300._r8 - 120._r8*sinsq(i)*(pmid(i,k)/psurf_ref)**4.3_r8 t_rlx(i,k,lchnk-begchunk+1)=trefc+lapsewy*10._r8/log(10._r8)*log((pmid(i,k)/psurf_ref)) elseif (flag==1) then !rest atmosphere trefc=285._r8 t_rlx(i,k,lchnk-begchunk+1)=trefc+lapsewy*10._r8/log(10._r8)*log((pmid(i,k)/psurf_ref)) endif enddo endif enddo !readin relaxation temperature file if realistic atmosphere (flag==3) if (flag==3) then do i = 1,ncol !open(51,file='/glade/p/ucsu0031/wykang/output/spcam4xco2/spcam4xco2_TDJF.txt',status='old',form='formatted',ACCESS='DIRECT',RECL=240) if (topo==0) then if (ndim==3) then if (clima==1) then open(51,file='/glade/p/ucsu0031/wykang/output/'//trim(model)//'/'//trim(model)//'_vinth_T'//trim(season)//trim(modified)//'.txt',status='old',form='formatted') else open(51,file='/glade/scratch/wykang/'//trim(IDEALRLX)//'/atm/hist/TRLX_'//trim(model)//trim(season)//trim(modified)//'.txt',status='old',form='formatted') endif elseif (ndim==2) then if (clima==1) then open(51,file='/glade/p/ucsu0031/wykang/output/'//trim(model)//'/'//trim(model)//'_vinth_T'//trim(season)//trim(modified)//'_lonavg.txt',status='old',form='formatted') else open(51,file='/glade/scratch/wykang/'//trim(IDEALRLX)//'/atm/hist/TRLX_'//trim(model)//trim(season)//trim(modified)//'.txt',status='old',form='formatted') endif endif else !if there is topography,there should be zonal asymmetry,ndim=3 if (clima==1) then if (ndim==3) then open(51,file='/glade/p/ucsu0031/wykang/output/'//trim(model)//'/'//trim(model)//'_T'//trim(season)//trim(modified)//'.txt',status='old',form='formatted') !!!PLEASE_SET1!!! elseif (ndim==2) then open(51,file='/glade/p/ucsu0031/wykang/output/'//trim(model)//'/'//trim(model)//'_T'//trim(season)//trim(modified)//'_lonavg.txt',status='old',form='formatted') endif else open(51,file='/glade/p/ucsu0031/wykang/'//trim(IDEALRLX)//'/atm/hist/TRLX_'//trim(model)//trim(season)//trim(modified)//'.txt',status='old',form='formatted') endif endif ! must use exactly the same vertical resolution to run. All files are ! interpolated to WACCM level ! (waccm ones: level0=10 nlev=56) ! (spcam ones: level0=0 nlev=30) linenum=(ilon(i)-1)*96+ilat(i) !write(iulog,*) 'wykang==============ilon:',ilon(i),'ilat:',ilat(i) if (linenum/=1) then do pre = 1,linenum-1 read(unit=51,fmt="(35i8)") read_in_tmp !write(iulog,"(35i8)") read_in_tmp enddo endif !read(unit=51,fmt="(35i8)",REC=linenum) read_in_tmp read(unit=51,fmt="(35i8)") read_in_tmp !write(iulog,fmt=*) 'yeah...ilon',ilon(i),'ilat',ilat(i),' data:',read_in_tmp ! the number before i should be the same as nlev!~ t_rlx(i,1:pver,lchnk-begchunk+1) = read_in_tmp(1:pver) / 10000._r8 close(51) enddo endif !flag==3 !============================================================== !If use Hall's Method calculate forcing to maintain background !============================================================== if (FT==1) then do i=1,ncol if (ndim==3) then open(54,file='/glade/p/ucsu0031/wykang/'//trim(forcmodel)//'/atm/hist/FT_'//trim(FTUVname)//trim(Forclima)//'.txt',status='old',form='formatted')!!!PLEASE_SET2!!! write(iulog,*) 'wykang read','/glade/p/ucsu0031/wykang/',trim(forcmodel),'/atm/hist/FT_',trim(FTUVname),trim(Forclima),'.txt' elseif (ndim==2) then open(54,file='/glade/p/ucsu0031/wykang/'//trim(forcmodel)//'/atm/hist/FT_'//trim(FTUVname)//trim(Forclima)//'_lonavg.txt',status='old',form='formatted') endif linenum=(ilon(i)-1)*96+ilat(i) if (linenum/=1) then do pre = 1,linenum-1 read(unit=54,fmt="(35i8)") read_in_tmp enddo endif read(unit=54,fmt="(35i8)") read_in_tmp t_forc(i,1:pver,lchnk-begchunk+1) = read_in_tmp(1:pver) /1000000000._r8 close(54) enddo else t_forc(i,1:pver,lchnk-begchunk+1) = 0._r8 endif !FT from Hall's method endif !nstep==0 !do i=1,ncol ! if (abs(t_rlx(i,pver,lchnk-begchunk+1))<1._r8) then ! write(iulog,*) 'wykang===000 step:',nstep,'lchnk:',lchnk,'i:',i ! endif !enddo !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !Add Forcing !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! if (tforcing==1) then !Forcing @ Frclat do k=1,pver do i=1,ncol if (pmid(i,k) .ge. Frchigh .and. pmid(i,k) .le. Frclow .and. mod(clon(i)+lonshift,2*pi) .ge. lonL/180._r8*pi .and. mod(clon(i)+lonshift,2*pi) .le. lonR/180._r8*pi .and.clat(i) .ge. latS/180._r8*pi .and. clat(i) .le. latN/180._r8*pi) then if (tforcingtype==0) then ! gaussian structure if (moving==1) then ! add moving forcing lon1=(MOD(clon(i)*180._r8/pi-(lonR-lonL)/Fperiod/48._r8*nstep,lonR-lonL)+(lonR-lonL)/2._r8+Frclon)*pi/180._r8 else lon1=clon(i) endif heatforce(i,k)=Amp*cos((log(psurf_ref/pmid(i,k))-log(psurf_ref/Frclev))/log(Frcthick)*pi/2._r8) if (sigmalon .NE. 0._r8) then heatforce(i,k)=heatforce(i,k)*exp(-(lon1-Frclon/180._r8*pi)**2._r8/(2._r8*(sigmalon*pi/180._r8)**2._r8)) endif if (sigmalat .NE. 0._r8) then heatforce(i,k)=heatforce(i,k)*exp(-(asin(sqrt(sinsq(i)))-Frclat*pi/180._r8)**2._r8/(2._r8*(sigmalat*pi/180._r8)**2._r8)) endif else ! sine structure if (moving==1) then ! add moving forcing lon1=MOD((clon(i)+lonshift)*180._r8/pi-360._r8/Fperiod/48._r8*nstep,360._r8)*pi/180._r8 else lon1=clon(i) endif heatforce(i,k)=Amp*cos((lon1-Frclon/180._r8*pi)*real(tforcingtype,8))*cos((log(psurf_ref/pmid(i,k))-log(psurf_ref/Frclev))/log(Frcthick)*pi/2._r8) if (sigmalat.NE.0._r8) then heatforce(i,k)=heatforce(i,k)*exp(-(asin(sqrt(sinsq(i)))-Frclat*pi/180._r8)**2._r8/(2._r8*(sigmalat*pi/180._r8)**2._r8)) endif endif endif enddo enddo endif if (tforcing .NE. 0 .AND. lonlim) then do i=1,ncol do k=1,pver heatforce(i,k)=heatforce(i,k)*(tanh((mod(clon(i)+lonshift,2*pi)-lonL/180._r8*pi)/((lonR-lonL)/180._r8*pi)*10._r8)*tanh((lonR/180._r8*pi-mod(clon(i)+lonshift,2*pi))/((lonR-lonL)/180._r8*pi)*10._r8)+1._r8)/2._r8 enddo enddo endif if (tforcing .NE. 0 .AND. latlim) then do i=1,ncol do k=1,pver heatforce(i,k)=heatforce(i,k)*(tanh((clat(i)-latS/180._r8*pi)/((latN-latS)/180._r8*pi)*10._r8)*tanh((latN/180._r8*pi-clat(i))/((latN-latS)/180._r8*pi)*10._r8)+1._r8)/2._r8 enddo enddo endif !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !check !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! if (checkread==1) then if (firststep) then open(54,file='/glade/p/ucsu0031/wykang/'//trim(forcmodel)//'/atm/hist/FT_'//trim(FTUVname)//trim(Forclima)//'check.txt',Access='append',form='formatted') do i=1,ncol l1=clon(i)*180._r8/3.141593_r8 l2=clat(i)*180._r8/3.141593_r8 linenum=(ilon(i)-1)*96+ilat(i) write(unit=54,fmt="(3i5,2f8.1,f10.2)") linenum,ilon(i),ilat(i),l1,l2,t_forc(i,1,lchnk-begchunk+1)*1000000000._r8 enddo close(54) endif if (nstep==8) then open(54,file='/glade/p/ucsu0031/wykang/'//trim(forcmodel)//'/atm/hist/FT_'//trim(FTUVname)//trim(Forclima)//'check8.txt',Access='append',form='formatted') do i=1,ncol l1=clon(i)*180._r8/3.141593_r8 l2=clat(i)*180._r8/3.141593_r8 linenum=(ilon(i)-1)*96+ilat(i) write(unit=54,fmt="(3i5,2f8.1,f10.1)") linenum,ilon(i),ilat(i),l1,l2,t_forc(i,pver,lchnk-begchunk+1)*1000000000._r8 enddo close(54) endif endif !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !Calculate tendency !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! if (solverlx==1.AND.firststep) then do k=1,pver do i=1,ncol TRLX(i,k,lchnk-begchunk+1)=t_rlx(i,k,lchnk-begchunk+1) enddo enddo endif !do k=1,pver !do i=1,ncol !ptend%s(i,k)=0._r8 !enddo !enddo if (solverlx==1) then do k=1,pver do i=1,ncol ptend%s(i,k) = (TRLX(i,k,lchnk-begchunk+1) -state%t(i,k))*tau_t(i,k,lchnk-begchunk+1)*cpair TERR(i,k)=state%t(i,k)-t_rlx(i,k,lchnk-begchunk+1) if (nstep.NE.0) then alp=1000._r8*(nstep**(0.3_r8)) endif TRLX(i,k,lchnk-begchunk+1)=TRLX(i,k,lchnk-begchunk+1)-TERR(i,k)/alp !this method does not take zonal average for 2d relaxation fields, to do !zonal !average, need to gather TRLX from chunks, take the average and scatter !it back !to chunks. enddo enddo else !solverlx!=1 do k=1,pver do i=1,ncol ptend%s(i,k) = (t_rlx(i,k,lchnk-begchunk+1) -state%t(i,k))*tau_t(i,k,lchnk-begchunk+1)*cpair TERR(i,k)=state%t(i,k)-t_rlx(i,k,lchnk-begchunk+1) enddo enddo endif if (FT==1) then do k=1,pver do i=1,ncol ptend%s(i,k) = ptend%s(i,k) +t_forc(i,k,lchnk-begchunk+1)*cpair enddo enddo endif if (tforcing==1 ) then !Extra heating or cooling if (TimeVary.GT.0._r8) then do k=1,pver do i=1,ncol ptend%s(i,k) = ptend%s(i,k)+heatforce(i,k)*cpair*(1-ABS((nstep-step0)/48._r8/365._r8-TimeVary)/TimeVary) enddo enddo else if (TimeVary.EQ.0._r8) then do k=1,pver do i=1,ncol ptend%s(i,k) = ptend%s(i,k) +heatforce(i,k)*cpair enddo enddo else do k=1,pver do i=1,ncol ptend%s(i,k) = ptend%s(i,k)+heatforce(i,k)*cpair*sin((nstep)/48._r8/Fperiod*2*pi) enddo enddo endif endif if (BLK) then do k=1,pver do i=1,ncol if (((clat(i).GE.Blklat1.and.clat(i).LE.Blklat2).OR.(.NOT.Blklat)).AND.& ((clon(i).GE.Blklon1.and.clon(i).LE.Blklon2).OR.(.NOT.Blklon)).AND.& ((pmid(i,k).GE.Blklev1.and.pmid(i,k).LE.Blklev2).OR.(.NOT.Blklev))) then Ttend_blk=(t_rlx(i,k,lchnk-begchunk+1)-state%t(i,k))/(ztodt/2._r8)/2._r8 *cpair weight=1._r8 !if (Blklon) then ! weight=weight*tanh((Blklon2-clon(i))/(Blklon2-Blklon1)*10._r8)*tanh((clon(i)-Blklon1)/(Blklon2-Blklon1)*10._r8) !endif !if (Blklat) then ! weight=weight*tanh((Blklat2-clat(i))/(Blklat2-Blklat1)*10._r8)*tanh((clat(i)-Blklat1)/(Blklat2-Blklat1)*10._r8) !endif !if (Blklev) then ! weight=weight*tanh((log(Blklev2)-log(pmid(i,k)))/(log(Blklev2)-log(Blklev1))*10._r8)*tanh((log(pmid(i,k))-log(Blklev1))/(log(Blklev2)-log(Blklev1))*10._r8) !endif if (Blklon) then weight=weight*exp(-(clon(i)-Blkloncenter)**2._r8/(2._r8*(Blklonsigma)**2)) endif if (Blklat) then weight=weight*exp(-(clat(i)-Blklatcenter)**2._r8/(2._r8*(Blklatsigma)**2)) endif if (Blklev) then weight=weight*exp(-(log(pmid(i,k))-Blklevcenter)**2._r8/(2._r8*(Blklevsigma)**2)) endif ptend%s(i,k)=(1-weight)*ptend%s(i,k)+weight*Ttend_blk endif enddo enddo endif !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! initialize uv tendency !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! do k=1,pver do i=1,pcols ptend%u(i,k) = 0._r8 ptend%v(i,k) = 0._r8 enddo enddo !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !relax wind field to model output too !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! if (firststep) then do i=1,ncol do k=1,pver tau_uv(i,k,lchnk-begchunk+1)=tau_t(i,k,lchnk-begchunk+1) !keep the same relaxation time scale as temperature enddo enddo endif if (relaxwind==1 .OR. FUV==1 .OR. solverlx==1 .OR. flag==3) then if (firststep) then do i = 1,ncol if (topo==0) then if (ndim==3) then if (clima==1 .or. nowind==1) then open(52,file='/glade/p/ucsu0031/wykang/output/'//trim(model)//'/'//trim(model)//'_vinth_U'//trim(season)//trim(modified)//'.txt',status='old',form='formatted') open(53,file='/glade/p/ucsu0031/wykang/output/'//trim(model)//'/'//trim(model)//'_vinth_V'//trim(season)//trim(modified)//'.txt',status='old',form='formatted') endif elseif (ndim==2) then if (clima==1 .or. nowind==1) then open(52,file='/glade/p/ucsu0031/wykang/output/'//trim(model)//'/'//trim(model)//'_vinth_U'//trim(season)//trim(modified)//'_lonavg.txt',status='old',form='formatted') open(53,file='/glade/p/ucsu0031/wykang/output/'//trim(model)//'/'//trim(model)//'_vinth_V'//trim(season)//trim(modified)//'_lonavg.txt',status='old',form='formatted') endif endif else !if there is topography...(not working on this right now) if (clima==1 .or. nowind==1) then if (ndim==3) then open(52,file='/glade/p/ucsu0031/wykang/output/'//trim(model)//'/'//trim(model)//'_U'//trim(season)//trim(modified)//'.txt',status='old',form='formatted')!!!PLEASE_SET1!!! open(53,file='/glade/p/ucsu0031/wykang/output/'//trim(model)//'/'//trim(model)//'_V'//trim(season)//trim(modified)//'.txt',status='old',form='formatted')!!!PLEASE_SET1!!! elseif (ndim==2) then open(52,file='/glade/p/ucsu0031/wykang/output/'//trim(model)//'/'//trim(model)//'_U'//trim(season)//trim(modified)//'_lonavg.txt',status='old',form='formatted') open(53,file='/glade/p/ucsu0031/wykang/output/'//trim(model)//'/'//trim(model)//'_V'//trim(season)//trim(modified)//'_lonavg.txt',status='old',form='formatted') endif endif endif linenum=(ilon(i)-1)*96+ilat(i) if (linenum/=1) then do pre = 1,linenum-1 read(unit=52,fmt="(35i8)") read_in_tmp read(unit=53,fmt="(35i8)") read_in_tmp enddo endif read(unit=52,fmt="(35i8)") read_in_tmp u_rlx(i,1:pver,lchnk-begchunk+1) = read_in_tmp(1:pver) / 10000._r8 read(unit=53,fmt="(35i8)") read_in_tmp v_rlx(i,1:pver,lchnk-begchunk+1) = read_in_tmp(1:pver) / 10000._r8 close(52) close(53) enddo endif !nstep==1 if (relaxwind.EQ.1) then do k=1,pver do i=1,ncol ptend%u(i,k) = (u_rlx(i,k,lchnk-begchunk+1) - state%u(i,k))*tau_uv(i,k,lchnk-begchunk+1) ptend%v(i,k) = (v_rlx(i,k,lchnk-begchunk+1) - state%v(i,k))*tau_uv(i,k,lchnk-begchunk+1) UERR(i,k)=state%u(i,k)-u_rlx(i,k,lchnk-begchunk+1) VERR(i,k)=state%v(i,k)-v_rlx(i,k,lchnk-begchunk+1) enddo enddo endif endif !relaxwind==1 .or. flag==3 .or. FUV==1 ! !surface friction ! kf = 1._r8/(86400._r8*efoldf) ! do k=1,pver if (pref_mid_norm(k) > sigmab) then kv = kf*(pref_mid_norm(k) - sigmab)/onemsig #ifdef MODHS ump = -kv/(1._r8+ ztodt*kv) #else ump = -kv #endif do i=1,ncol if (fricuv2clima==1) then ptend%u(i,k) = ptend%u(i,k)+ump*(state%u(i,k)-u_rlx(i,k,lchnk-begchunk+1)) ptend%v(i,k) = ptend%v(i,k)+ump*(state%v(i,k)-v_rlx(i,k,lchnk-begchunk+1)) else ptend%u(i,k) = ptend%u(i,k)+ump*state%u(i,k) ptend%v(i,k) = ptend%v(i,k)+ump*state%v(i,k) endif enddo endif enddo ! ! wanying: adding damping at the top: ! if (damp>0) then do k=1,pver if (pmid(i,k).lt.11._r8) then do i=1,ncol if (fricuv2clima==1) then ptend%u(i,k) = ptend%u(i,k)-1._r8/30._r8/86400._r8*(state%u(i,k)-u_rlx(i,k,lchnk-begchunk+1)) ptend%v(i,k) = ptend%v(i,k)-1._r8/30._r8/86400._r8*(state%v(i,k)-v_rlx(i,k,lchnk-begchunk+1)) else ptend%u(i,k) = ptend%u(i,k)-1._r8/30._r8/86400._r8*state%u(i,k) ptend%v(i,k) = ptend%v(i,k)-1._r8/30._r8/86400._r8*state%v(i,k) endif enddo endif enddo endif if (damp==2) then do k=1,3 do i=1,ncol weight=400._r8/pmid(i,k)*cos(clat(i)) if (fricuv2clima==1) then ptend%s(i,k) = (1-weight)*ptend%s(i,k)-weight*k_top_damp*(state%t(i,k)-t_rlx(i,k,lchnk-begchunk+1))*cpair ptend%u(i,k) = (1-weight)*ptend%u(i,k)-weight*k_top_damp*(state%u(i,k)-u_rlx(i,k,lchnk-begchunk+1)) ptend%v(i,k) = (1-weight)*ptend%v(i,k)-weight*k_top_damp*(state%v(i,k)-v_rlx(i,k,lchnk-begchunk+1)) else ptend%u(i,k) = (1-weight)*ptend%u(i,k)-weight*k_top_damp*state%u(i,k) ptend%v(i,k) = (1-weight)*ptend%v(i,k)-weight*k_top_damp*state%v(i,k) endif enddo enddo endif if (FUV==1) then if (firststep) then do i=1,ncol if (ndim==3) then open(55,file='/glade/p/ucsu0031/wykang/'//trim(forcmodel)//'/atm/hist/FU_'//trim(FTUVname)//trim(Forclima)//'.txt',status='old',form='formatted')!!!PLEASE_SET2!!! open(56,file='/glade/p/ucsu0031/wykang/'//trim(forcmodel)//'/atm/hist/FV_'//trim(FTUVname)//trim(Forclima)//'.txt',status='old',form='formatted')!!!PLEASE_SET2!!! elseif (ndim==2) then open(55,file='/glade/p/ucsu0031/wykang/'//trim(forcmodel)//'/atm/hist/FU_'//trim(FTUVname)//trim(Forclima)//'_lonavg.txt',status='old',form='formatted') open(56,file='/glade/p/ucsu0031/wykang/'//trim(forcmodel)//'/atm/hist/FV_'//trim(FTUVname)//trim(Forclima)//'_lonavg.txt',status='old',form='formatted') endif linenum=(ilon(i)-1)*96+ilat(i) if (linenum/=1) then do pre = 1,linenum-1 read(unit=55,fmt="(35i8)") read_in_tmp read(unit=56,fmt="(35i8)") read_in_tmp enddo endif read(unit=55,fmt="(35i8)") read_in_tmp u_forc(i,1:pver,lchnk-begchunk+1) = read_in_tmp(1:pver) / 1000000000._r8 read(unit=56,fmt="(35i8)") read_in_tmp v_forc(i,1:pver,lchnk-begchunk+1) = read_in_tmp(1:pver) / 1000000000._r8 close(55) close(56) enddo endif !first step initiate force uv do i=1,ncol do k=1,pver ptend%u(i,k) = ptend%u(i,k)+u_forc(i,k,lchnk-begchunk+1) ptend%v(i,k) = ptend%v(i,k)+v_forc(i,k,lchnk-begchunk+1) enddo enddo else ! if FUV==1 u_forc(i,1:pver,lchnk-begchunk+1) =0._r8 v_forc(i,1:pver,lchnk-begchunk+1) =0._r8 endif ! if FUV==1 add wind forcing from file if (BLK) then do k=1,pver do i=1,ncol if (((clat(i).GE.Blklat1.and.clat(i).LE.Blklat2).OR.(.NOT.Blklat)).AND.& ((clon(i).GE.Blklon1.and.clon(i).LE.Blklon2).OR.(.NOT.Blklon)).AND.& ((pmid(i,k).GE.Blklev1.and.pmid(i,k).LE.Blklev2).OR.(.NOT.Blklev))) then Utend_blk=(u_rlx(i,k,lchnk-begchunk+1) -state%u(i,k))/(ztodt/2._r8)/2._r8 Vtend_blk=(v_rlx(i,k,lchnk-begchunk+1) -state%v(i,k))/(ztodt/2._r8)/2._r8 weight=1._r8 !if (Blklon) then ! weight=weight*tanh((Blklon2-clon(i))/(Blklon2-Blklon1)*10._r8)*tanh((clon(i)-Blklon1)/(Blklon2-Blklon1)*10._r8) !endif !if (Blklat) then ! weight=weight*tanh((Blklat2-clat(i))/(Blklat2-Blklat1)*10._r8)*tanh((clat(i)-Blklat1)/(Blklat2-Blklat1)*10._r8) !endif !if (Blklev) then ! weight=weight*tanh((log(Blklev2)-log(pmid(i,k)))/(log(Blklev2)-log(Blklev1))*10._r8)*tanh((log(pmid(i,k))-log(Blklev1))/(log(Blklev2)-log(Blklev1))*10._r8) !endif if (Blklon) then weight=weight*exp(-(clon(i)-Blkloncenter)**2._r8/(2._r8*(Blklonsigma)**2)) endif if (Blklat) then weight=weight*exp(-(clat(i)-Blklatcenter)**2._r8/(2._r8*(Blklatsigma)**2)) endif if (Blklev) then weight=weight*exp(-(log(pmid(i,k))-Blklevcenter)**2._r8/(2._r8*(Blklevsigma)**2)) endif ptend%u(i,k)=(1-weight)*ptend%u(i,k)+weight*Utend_blk ptend%v(i,k)=(1-weight)*ptend%v(i,k)+weight*Vtend_blk endif enddo enddo endif if (solverlx==1) then call outfld('TRLX ', TRLX(:,:,lchnk-begchunk+1), pcols, lchnk ) endif!only relax temperature at this point call outfld('TERR ', TERR, pcols, lchnk ) call outfld('UERR ', UERR, pcols, lchnk ) call outfld('VERR ', VERR, pcols, lchnk ) call outfld('TFORC ', t_forc(:,:,lchnk-begchunk+1), pcols, lchnk ) call outfld('UFORC ', u_forc(:,:,lchnk-begchunk+1), pcols, lchnk ) call outfld('VFORC ', v_forc(:,:,lchnk-begchunk+1), pcols, lchnk ) call outfld('UTARGET ', u_rlx(:,:,lchnk-begchunk+1),pcols,lchnk) call outfld('VTARGET ', v_rlx(:,:,lchnk-begchunk+1),pcols,lchnk) call outfld('RLXCT ', tau_t(:,:,lchnk-begchunk+1), pcols, lchnk) call outfld('TTARGET ', t_rlx(:,:,lchnk-begchunk+1),pcols,lchnk) !if (Hall) then !call outfld('T_ ',state%t,pcols,lchnk) !call outfld('U_ ',state%u,pcols,lchnk) !call outfld('V_ ',state%v,pcols,lchnk) !call outfld('P_ ',state%pmid,pcols,lchnk) !endif elseif (idlflag == 2 ) then ! !----------------------------------------------------------------------- ! ! Modified Held/Suarez IDEALIZED physics algorithm ! (modified with Williamson stratosphere): ! ! Williamson, D. L., J. G. Olson and B. A. Boville, 1998: A comparison ! of semi--Lagrangian and Eulerian tropical climate simulations. ! Mon. Wea. Rev., vol 126, pp. 1001-1012. ! !----------------------------------------------------------------------- ! ! Compute idealized radiative heating rates (as dry static energy) ! efoldf = 1._r8 efolda = 40._r8 efoldaa = 40._r8 efolds = 4._r8 sigmab = 0.7_r8 t00 = 200._r8 ! onemsig = 1._r8 - sigmab ! ka = 1._r8/(86400._r8*efolda) kaa = 1._r8/(86400._r8*efoldaa) ks = 1._r8/(86400._r8*efolds) ! pi = 4._r8*atan(1._r8) phi0 = 60._r8*pi/180._r8 dphi0 = 15._r8*pi/180._r8 a0 = 2.65_r8/dphi0 aeq = 10000._r8 apole = 200._r8 lapsew = -3.345e-03_r8 constw = rair*lapsew/gravit lapsec = 2.00e-03_r8 constc = rair*lapsec/gravit if (firststep) then do k=1,pver if (pref_mid_norm(k) > sigmab) then do i=1,ncol kt = ka + (ks - ka)*cossqsq(i)*(pref_mid_norm(k) - sigmab)/onemsig acoslat = abs(acos(coslat(i))) p0strat = aeq - (aeq - apole)*0.5_r8*(1._r8 + tanh(a0*(acoslat - phi0))) tau_t(i,k,lchnk-begchunk+1) = kt/(1._r8+ ztodt*kt) trefc = 315._r8 - 60._r8*sinsq(i) t_rlx(i,k,lchnk-begchunk+1) = (trefc - 10._r8*cossq(i)*log((pmid(i,k)/psurf_ref)))*(pmid(i,k)/psurf_ref)** & cappa t_rlx(i,k,lchnk-begchunk+1) = max(t00,t_rlx(i,k,lchnk-begchunk+1)) if (pmid(i,k) < 10000._r8) then t_rlx(i,k,lchnk-begchunk+1) = t00*((pmid(i,k)/10000._r8))**constc tau_t(i,k,lchnk-begchunk+1) = kaa/(1._r8+ ztodt*kaa) endif if (pmid(i,k) < p0strat) then t_rlx(i,k,lchnk-begchunk+1) = t_rlx(i,k,lchnk-begchunk+1) + t00*( ((pmid(i,k)/p0strat))**constw - 1._r8 ) tau_t(i,k,lchnk-begchunk+1) = kaa/(1._r8+ ztodt*kaa) endif enddo else do i=1,ncol acoslat = abs(acos(coslat(i))) p0strat = aeq - (aeq - apole)*0.5_r8*(1._r8 + tanh(a0*(acoslat - phi0))) tau_t(i,k,lchnk-begchunk+1)= ka/(1._r8+ ztodt*ka) trefc = 315._r8 - 60._r8*sinsq(i) t_rlx(i,k,lchnk-begchunk+1) = (trefc - 10._r8*cossq(i)*log((pmid(i,k)/psurf_ref)))*(pmid(i,k)/psurf_ref)** & cappa t_rlx(i,k,lchnk-begchunk+1) = max(t00,t_rlx(i,k,lchnk-begchunk+1)) if (pmid(i,k) < 10000._r8) then t_rlx(i,k,lchnk-begchunk+1) = t00*((pmid(i,k)/10000._r8))**constc tau_t(i,k,lchnk-begchunk+1)= kaa/(1._r8+ ztodt*kaa) endif if (pmid(i,k) < p0strat) then t_rlx(i,k,lchnk-begchunk+1) = t_rlx(i,k,lchnk-begchunk+1) + t00*( ((pmid(i,k)/p0strat))**constw - 1._r8 ) tau_t(i,k,lchnk-begchunk+1)= kaa/(1._r8+ ztodt*kaa) endif enddo endif enddo endif ! firststep T_rlx set finished ! !Add Forcing ! if (tforcing==1) then !Forcing @ Frclat pi = 4._r8*atan(1._r8) do k=1,pver do i=1,ncol if (pmid(i,k) .ge. Frchigh .and. pmid(i,k) .le. Frclow .and. clon(i) .ge. lonL/180._r8*pi .and.clon(i) .le. lonR/180._r8*pi .and.clat(i) .ge. latS/180._r8*pi .and. clat(i) .le. latN/180._r8*pi) then if (tforcingtype==0) then if (moving==1) then ! add moving forcing lon1=(MOD(clon(i)*180._r8/pi-(lonR-lonL)/Fperiod/48._r8*nstep,lonR-lonL)+(lonR-lonL)/2._r8+Frclon)*pi/180._r8 else lon1=clon(i) endif heatforce(i,k)=Amp*cos((log(psurf_ref/pmid(i,k))-log(psurf_ref/Frclev))/log(Frcthick)*pi/2._r8) if (sigmalon .NE. 0._r8) then heatforce(i,k)=heatforce(i,k)*exp(-(lon1-Frclon/180._r8*pi)**2._r8/(2._r8*(sigmalon*pi/180._r8)**2._r8)) endif if (sigmalat .NE. 0._r8) then heatforce(i,k)=heatforce(i,k)*exp(-(asin(sqrt(sinsq(i)))-Frclat*pi/180._r8)**2._r8/(2._r8*(sigmalat*pi/180._r8)**2._r8)) endif else if (moving==1) then ! add moving forcing lon1=MOD(clon(i)*180._r8/pi-(lonR-lonL)/Fperiod/48._r8*nstep,360._r8)*pi/180._r8 else lon1=clon(i) endif heatforce(i,k)=Amp*cos((lon1-Frclon/180._r8*pi)*real(tforcingtype,8))*cos((log(psurf_ref/pmid(i,k))-log(psurf_ref/Frclev))/log(Frcthick)*pi/2._r8) if (sigmalat.NE.0._r8) then heatforce(i,k)=heatforce(i,k)*exp(-(asin(sqrt(sinsq(i)))-Frclat*pi/180._r8)**2._r8/(2._r8*(sigmalat*pi/180._r8)**2._r8)) endif endif endif enddo enddo endif if (tforcing .NE. 0 .AND. lonlim) then do i=1,ncol do k=1,pver heatforce(i,k)=heatforce(i,k)*(tanh((clon(i)-lonL/180._r8*pi)/((lonR-lonL)/180._r8*pi)*10._r8)*tanh((lonR/180._r8*pi-clon(i))/((lonR-lonL)/180._r8*pi)*10._r8)+1._r8)/2._r8 enddo enddo endif if (tforcing .NE. 0 .AND. latlim) then do i=1,ncol do k=1,pver heatforce(i,k)=heatforce(i,k)*(tanh((clat(i)-latS/180._r8*pi)/((latN-latS)/180._r8*pi)*10._r8)*tanh((latN/180._r8*pi-clat(i))/((latN-latS)/180._r8*pi)*10._r8)+1._r8)/2._r8 enddo enddo endif !Forcing end call outfld('TTARGET ', t_rlx(:,:,lchnk-begchunk+1), pcols, lchnk) call outfld('RLXCT ', tau_t(:,:,lchnk-begchunk+1), pcols, lchnk) ! !relax to equilibrium temperature field ! do k=1,pver do i=1,ncol ptend%s(i,k) = (t_rlx(i,k,lchnk-begchunk+1) - state%t(i,k))*tau_t(i,k,lchnk-begchunk+1)*cpair enddo enddo if (tforcing==1 .or. tforcing==2) then !Extra heating or cooling if (TimeVary.NE.0._r8) then do k=1,pver do i=1,ncol ptend%s(i,k) =ptend%s(i,k)+heatforce(i,k)*cpair*(1-ABS((nstep-step0)/48._r8/365._r8-TimeVary)/TimeVary) enddo enddo else do k=1,pver do i=1,ncol ptend%s(i,k) = ptend%s(i,k) +heatforce(i,k)*cpair enddo enddo endif endif TERR(i,k)=-t_rlx(i,k,lchnk-begchunk+1)+state%t(i,k) call outfld('TERR ', TERR, pcols, lchnk) ! ! Add diffusion near the surface for the wind fields ! do k=1,pver do i=1,pcols ptend%u(i,k) = 0._r8 ptend%v(i,k) = 0._r8 enddo enddo kf = 1._r8/(86400._r8*efoldf) ! do k=1,pver if (pref_mid_norm(k) > sigmab) then kv = kf*(pref_mid_norm(k) - sigmab)/onemsig ump = -kv/(1._r8+ ztodt*kv) do i=1,ncol ptend%u(i,k) = ump*state%u(i,k) ptend%v(i,k) = ump*state%v(i,k) enddo endif enddo !!!!!TOP DAMPING if (damp==1) then do k=1,4 do i=1,ncol if (fricuv2clima==1) then ptend%u(i,k) =ptend%u(i,k)-k_top_damp*(state%u(i,k)-u_rlx(i,k,lchnk-begchunk+1)) ptend%v(i,k) =ptend%v(i,k)-k_top_damp*(state%v(i,k)-v_rlx(i,k,lchnk-begchunk+1)) else ptend%u(i,k) = ptend%u(i,k)-k_top_damp*state%u(i,k) ptend%v(i,k) = ptend%v(i,k)-k_top_damp*state%v(i,k) endif enddo enddo endif elseif (idlflag == 3 ) then ! !----------------------------------------------------------------------- ! ! Held/Suarez IDEALIZED physics algorithm: ! (modified with Lin/Williamson stratosphere/mesosphere): ! ! Held, I. M., and M. J. Suarez, 1994: A proposal for the ! intercomparison of the dynamical cores of atmospheric general ! circulation models. ! Bulletin of the Amer. Meteor. Soc., vol. 75, pp. 1825-1830. ! !----------------------------------------------------------------------- ! ! Compute idealized radiative heating rates (as dry static energy) ! efoldf = 1._r8 efolda = 40._r8 efolds = 4._r8 efold_strat = 40._r8 efold_meso = 10._r8 efoldv = 0.5_r8 sigmab = 0.7_r8 lapse = 0.00225_r8 h0 = 7000._r8 t00 = 200._r8 p_infint = 0.01_r8 ! onemsig = 1._r8 - sigmab ! ka = 1._r8/(86400._r8*efolda) ks = 1._r8/(86400._r8*efolds) ! do k=1,pver if (pref_mid_norm(k) > sigmab) then do i=1,ncol kt = ka + (ks - ka)*cossqsq(i)*(pref_mid_norm(k) - sigmab)/onemsig tmp = kt/(1._r8+ ztodt*kt) trefc = 315._r8 - 60._r8*sinsq(i) t_rlx(i,k,lchnk-begchunk+1) = (trefc - 10._r8*cossq(i)*log((pmid(i,k)/psurf_ref)))*(pmid(i,k)/psurf_ref)**cappa t_rlx(i,k,lchnk-begchunk+1) = max(t00,t_rlx(i,k,lchnk-begchunk+1)) ptend%s(i,k) = (t_rlx(i,k,lchnk-begchunk+1) - state%t(i,k))*tmp*cpair enddo else do i=1,ncol tmp = ka/(1._r8+ ztodt*ka) pressmb = pmid(i,k)*0.01_r8 trefc = 315._r8 - 60._r8*sinsq(i) t_rlx(i,k,lchnk-begchunk+1) = (trefc - 10._r8*cossq(i)*log((pmid(i,k)/psurf_ref)))*(pmid(i,k)/psurf_ref)** & cappa t_rlx(i,k,lchnk-begchunk+1) = max(t00,t_rlx(i,k,lchnk-begchunk+1)) if (pressmb <= 100._r8 .and. pressmb > 1._r8) then t_rlx(i,k,lchnk-begchunk+1) = t00 + lapse*h0*coslat(i)*log(100._r8/pressmb) tmp = (efold_strat-efold_meso)*log(pressmb)/log(100._r8) tmp = efold_meso + tmp tmp = 1._r8/(86400._r8*tmp) tmp = tmp/(1._r8+ ztodt*tmp) endif if (pressmb <= 1._r8 .and. pressmb > 0.01_r8) then t_rlx(i,k,lchnk-begchunk+1) = t00 + lapse*h0*coslat(i)*log(100._r8*pressmb) tmp = 1._r8/(86400._r8*efold_meso) tmp = tmp/(1._r8+ ztodt*tmp) endif if (pressmb <= 0.01_r8) then tmp = 1._r8/(86400._r8*efold_meso) tmp = tmp/(1._r8+ ztodt*tmp) endif ptend%s(i,k) = (t_rlx(i,k,lchnk-begchunk+1) - state%t(i,k))*tmp*cpair enddo endif enddo ! ! Add diffusion near the surface for the wind fields ! do k=1,pver do i=1,pcols ptend%u(i,k) = 0._r8 ptend%v(i,k) = 0._r8 enddo enddo kf = 1._r8/(86400._r8*efoldf) ! do k=1,pver if (pref_mid_norm(k) > sigmab) then kv = kf*(pref_mid_norm(k) - sigmab)/onemsig ump = -kv/(1._r8+ ztodt*kv) do i=1,ncol ptend%u(i,k) = ump*state%u(i,k) ptend%v(i,k) = ump*state%v(i,k) enddo else do i=1,ncol pressmb = pmid(i,k)*0.01_r8 if (pressmb <= 100._r8) then kv = 1._r8/(86400._r8*efoldv) ump = 1._r8 + tanh(1.5_r8*log10(p_infint/pressmb)) kv = kv*ump ump = -kv/(1._r8+ ztodt*kv) ptend%u(i,k) = ump*state%u(i,k) ptend%v(i,k) = ump*state%v(i,k) endif enddo endif enddo else write(iulog,*) 'TPHYSIDL: flag for choosing desired type of idealized ', & 'physics ("idlflag") is set incorrectly.' write(iulog,*) 'The valid options are 1, 2, or 3.' write(iulog,*) 'idlflag is currently set to: ',idlflag call endrun('tphysidl: invalid option') endif ! update the state and total physics tendency call physics_update(state, ptend, ztodt, tend) ! Can't turn on conservation error messages unless the appropriate heat ! surface flux is computed and supplied as an argument to ! check_energy_chng to account for how the ideal physics tforcings are ! changing the total energy. call check_energy_chng(state, tend, "tphysidl", nstep, ztodt, zero, zero, zero, zero) call outfld('QRS', tend%dtdt, pcols, lchnk) end subroutine tphysidl