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