| 1 | subroutine newtrelax(ngrid,nlayer,mu0,sinlat,popsk,temp,pplay,pplev,dtrad,firstcall) |
|---|
| 2 | |
|---|
| 3 | use comcstfi_mod, only: rcp, pi |
|---|
| 4 | use callkeys_mod, only: tau_relax |
|---|
| 5 | implicit none |
|---|
| 6 | |
|---|
| 7 | #include "netcdf.inc" |
|---|
| 8 | |
|---|
| 9 | !================================================================== |
|---|
| 10 | ! |
|---|
| 11 | ! Purpose |
|---|
| 12 | ! ------- |
|---|
| 13 | ! Alternative Newtonian radiative transfer scheme. |
|---|
| 14 | ! |
|---|
| 15 | ! Authors |
|---|
| 16 | ! ------- |
|---|
| 17 | ! R. Wordsworth (2010) |
|---|
| 18 | ! |
|---|
| 19 | !================================================================== |
|---|
| 20 | |
|---|
| 21 | |
|---|
| 22 | ! Input |
|---|
| 23 | integer,intent(in) :: ngrid, nlayer |
|---|
| 24 | logical,intent(in) :: firstcall |
|---|
| 25 | real,intent(in) :: mu0(ngrid) ! cosine of sun incident angle |
|---|
| 26 | real,intent(in) :: sinlat(ngrid) ! sine of latitude |
|---|
| 27 | real,intent(in) :: temp(ngrid,nlayer) ! temperature at each layer (K) |
|---|
| 28 | real,intent(in) :: pplay(ngrid,nlayer) ! pressure at each layer (Pa) |
|---|
| 29 | real,intent(in) :: pplev(ngrid,nlayer+1) ! pressure at each level (Pa) |
|---|
| 30 | real,intent(in) :: popsk(ngrid,nlayer) ! pot. T to T converter |
|---|
| 31 | |
|---|
| 32 | ! Output |
|---|
| 33 | real,intent(out) :: dtrad(ngrid,nlayer) |
|---|
| 34 | |
|---|
| 35 | ! Internal |
|---|
| 36 | real Trelax_V, Trelax_H |
|---|
| 37 | real,allocatable,dimension(:,:),save :: Trelax |
|---|
| 38 | !$OMP THREADPRIVATE(Trelax) |
|---|
| 39 | |
|---|
| 40 | real T_trop ! relaxation temperature at tropopause (K) |
|---|
| 41 | real T_surf ! relaxation temperature at surface (K) |
|---|
| 42 | real dT_EP ! Equator-Pole relaxation temperature difference (K) |
|---|
| 43 | |
|---|
| 44 | real sig, f_sig, sig_trop |
|---|
| 45 | integer l,ig |
|---|
| 46 | |
|---|
| 47 | |
|---|
| 48 | logical tidallocked |
|---|
| 49 | parameter (tidallocked = .true.) |
|---|
| 50 | |
|---|
| 51 | ! Setup relaxation temperature |
|---|
| 52 | if(firstcall) then |
|---|
| 53 | |
|---|
| 54 | ALLOCATE(Trelax(ngrid,nlayer)) |
|---|
| 55 | |
|---|
| 56 | print*,'-----------------------------------------------------' |
|---|
| 57 | print*,'| ATTENTION: You are using a Newtonian cooling scheme' |
|---|
| 58 | print*,'| for the radiative transfer. This means that ALL' |
|---|
| 59 | print*,'| other physics subroutines must be switched off.' |
|---|
| 60 | print*,'-----------------------------------------------------' |
|---|
| 61 | |
|---|
| 62 | if(tidallocked)then |
|---|
| 63 | do ig=1,ngrid |
|---|
| 64 | |
|---|
| 65 | T_surf = 126. + 239.*mu0(ig) |
|---|
| 66 | T_trop = 140. + 52.*mu0(ig) |
|---|
| 67 | do l=1,nlayer |
|---|
| 68 | |
|---|
| 69 | if(mu0(ig).le.0.0)then ! night side |
|---|
| 70 | Trelax(ig,l)=0.0 |
|---|
| 71 | else ! day side |
|---|
| 72 | Trelax(ig,l) = T_surf*popsk(ig,l) |
|---|
| 73 | if (Trelax(ig,l).lt.T_trop) Trelax(ig,l) = T_trop |
|---|
| 74 | endif |
|---|
| 75 | |
|---|
| 76 | enddo |
|---|
| 77 | enddo |
|---|
| 78 | |
|---|
| 79 | else |
|---|
| 80 | |
|---|
| 81 | T_trop = 200. |
|---|
| 82 | T_surf = 288. |
|---|
| 83 | dT_EP = 70. |
|---|
| 84 | |
|---|
| 85 | sig_trop=(T_trop/T_surf)**(1./rcp) |
|---|
| 86 | |
|---|
| 87 | do l=1,nlayer |
|---|
| 88 | do ig=1,ngrid |
|---|
| 89 | |
|---|
| 90 | ! vertically varying component |
|---|
| 91 | Trelax_V = T_surf*popsk(ig,l) |
|---|
| 92 | if (Trelax_V.lt.T_trop) Trelax_V = T_trop |
|---|
| 93 | |
|---|
| 94 | ! horizontally varying component |
|---|
| 95 | sig = pplay(ig,l)/pplev(ig,1) |
|---|
| 96 | if(sig.ge.sig_trop)then |
|---|
| 97 | f_sig=sin((pi/2)*((sig-sig_trop)/(1-sig_trop))) |
|---|
| 98 | else |
|---|
| 99 | f_sig=0.0 |
|---|
| 100 | endif |
|---|
| 101 | Trelax_H = -f_sig*dT_EP*(sinlat(ig)**2 - 1./3.) |
|---|
| 102 | |
|---|
| 103 | Trelax(ig,l) = Trelax_V + Trelax_H |
|---|
| 104 | |
|---|
| 105 | enddo |
|---|
| 106 | enddo |
|---|
| 107 | |
|---|
| 108 | endif |
|---|
| 109 | |
|---|
| 110 | endif |
|---|
| 111 | |
|---|
| 112 | ! Calculate radiative forcing |
|---|
| 113 | do l=1,nlayer |
|---|
| 114 | do ig=1,ngrid |
|---|
| 115 | dtrad(ig,l) = -(temp(ig,l) - Trelax(ig,l)) / tau_relax |
|---|
| 116 | if(temp(ig,l).gt.500.)then ! Trelax(ig,l))then |
|---|
| 117 | print*,'ig=',ig |
|---|
| 118 | print*,'l=',l |
|---|
| 119 | print*,'temp=',temp(ig,l) |
|---|
| 120 | print*,'Trelax=',Trelax(ig,l) |
|---|
| 121 | endif |
|---|
| 122 | enddo |
|---|
| 123 | enddo |
|---|
| 124 | |
|---|
| 125 | call writediagfi(ngrid,'Tref','rad forc temp','K',3,Trelax) |
|---|
| 126 | !call writediagfi(ngrid,'ThetaZ','stellar zenith angle','deg',2,mu0) |
|---|
| 127 | |
|---|
| 128 | end subroutine newtrelax |
|---|
