| 1 | module newton_cooling_hotJ |
|---|
| 2 | |
|---|
| 3 | !========================================================================================== |
|---|
| 4 | ! Purpose |
|---|
| 5 | ! ------- |
|---|
| 6 | ! Compute a Newtonian cooling scheme for Hot Jupiters |
|---|
| 7 | ! for scenario 1 of the MOCHA intercomparaison project |
|---|
| 8 | ! (add citation to protocol paper here when it's live). |
|---|
| 9 | ! Check this paper's equations (1) and (4):https://iopscience.iop.org/article/10.3847/PSJ/ac9dfe/pdf |
|---|
| 10 | ! |
|---|
| 11 | ! We aim at having a generic code but you never know, it might need improving at some point. |
|---|
| 12 | ! The current (at time of writing) newtrelax.F90 routine is hardcoded for telluric temperate planets and untested. |
|---|
| 13 | ! Thus, we don't use it and use this one instead. |
|---|
| 14 | ! |
|---|
| 15 | ! Authors |
|---|
| 16 | ! ------- |
|---|
| 17 | ! Lucas Teinturier (2024) |
|---|
| 18 | ! |
|---|
| 19 | !========================================================================================== |
|---|
| 20 | implicit none |
|---|
| 21 | |
|---|
| 22 | ! Module variables |
|---|
| 23 | real, allocatable, save :: T0(:) |
|---|
| 24 | real, allocatable, save :: tau_relax(:) |
|---|
| 25 | real, allocatable, save :: delta_Teq(:) |
|---|
| 26 | real, allocatable, save :: Trelax(:,:) |
|---|
| 27 | character(100),save :: planetary_suffix |
|---|
| 28 | !$OMP THREADPRIVATE(planetary_suffix, Trelax,tau_relax,T0,delta_Teq) |
|---|
| 29 | |
|---|
| 30 | contains |
|---|
| 31 | |
|---|
| 32 | subroutine newtcool_MOCHA(ngrid,nlayer,coslon,coslat,temp,pplay,firstcall,lastcall,dtrad) |
|---|
| 33 | |
|---|
| 34 | ! use callkeys_mod, only: planetary_suffix ! this is to know which profiles to load for the T0, the delta Teq and the tau_rad |
|---|
| 35 | use mod_phys_lmdz_para, only : is_master, bcast ! for OpenMP stuff |
|---|
| 36 | implicit none |
|---|
| 37 | |
|---|
| 38 | ! Inputs |
|---|
| 39 | integer, intent(in) :: ngrid,nlayer |
|---|
| 40 | logical, intent(in) :: firstcall ! is it the first call of physiq_mod ? |
|---|
| 41 | logical, intent(in) :: lastcall !is it the last call of physiq_mod ? |
|---|
| 42 | real, intent(in) :: coslon(ngrid) !cosine of the longitude |
|---|
| 43 | real, intent(in) :: coslat(ngrid) ! cosine of the latitude |
|---|
| 44 | real, intent(in) :: temp(ngrid,nlayer) ! Temperature at each layer (K) |
|---|
| 45 | real, intent(in) :: pplay(ngrid,nlayer) ! Pressure mid-layers (Pa) |
|---|
| 46 | |
|---|
| 47 | ! Output |
|---|
| 48 | real, intent(out) :: dtrad(ngrid,nlayer) ! Tendency on temperature dT/dt (K/s) |
|---|
| 49 | |
|---|
| 50 | !! Internal variable |
|---|
| 51 | integer ig,l |
|---|
| 52 | character(100) :: filename |
|---|
| 53 | |
|---|
| 54 | if (firstcall) then |
|---|
| 55 | ! Allocation of the dynamical arrays |
|---|
| 56 | allocate(T0(nlayer)) |
|---|
| 57 | allocate(tau_relax(nlayer)) |
|---|
| 58 | allocate(delta_Teq(nlayer)) |
|---|
| 59 | allocate(Trelax(ngrid,nlayer)) |
|---|
| 60 | |
|---|
| 61 | if (is_master) then |
|---|
| 62 | print*,'-----------------------------------------------------' |
|---|
| 63 | print*,'| ATTENTION: You are using a Newtonian cooling scheme' |
|---|
| 64 | print*,'| for the radiative transfer. This means that ALL' |
|---|
| 65 | print*,'| other physics subroutines must be switched off.' |
|---|
| 66 | print*,'| Check that you have the required files in the ' |
|---|
| 67 | print*,'| simulation directory !' |
|---|
| 68 | print*,'-----------------------------------------------------' |
|---|
| 69 | print*,"the planetary suffix is ",planetary_suffix |
|---|
| 70 | |
|---|
| 71 | !! We load the data using the subroutine load_input |
|---|
| 72 | |
|---|
| 73 | ! Loading T0 |
|---|
| 74 | filename = trim(planetary_suffix) // "T0.dat" |
|---|
| 75 | ! print*,"filename = ",filename |
|---|
| 76 | call read_input(nlayer,filename,T0) |
|---|
| 77 | print*,"I successfully read",filename |
|---|
| 78 | |
|---|
| 79 | ! Loading tau_relax |
|---|
| 80 | filename = trim(planetary_suffix) // "tau_relax.dat" |
|---|
| 81 | call read_input(nlayer,filename,tau_relax) |
|---|
| 82 | print*,"I successfully read",filename |
|---|
| 83 | |
|---|
| 84 | ! Loading delta_Teq |
|---|
| 85 | filename = trim(planetary_suffix) // "delta_Teq.dat" |
|---|
| 86 | call read_input(nlayer,filename,delta_Teq) |
|---|
| 87 | print*,"I successfully read",filename |
|---|
| 88 | |
|---|
| 89 | endif ! of is_master |
|---|
| 90 | |
|---|
| 91 | ! Broadcast tau_relax and Trelax to everyone |
|---|
| 92 | call bcast(tau_relax) |
|---|
| 93 | call bcast(Trelax) |
|---|
| 94 | call bcast(T0) |
|---|
| 95 | call bcast(delta_Teq) |
|---|
| 96 | |
|---|
| 97 | ! now initialising Trelax depending on day or night side |
|---|
| 98 | do l=1,nlayer |
|---|
| 99 | do ig=1,ngrid |
|---|
| 100 | ! if we're on the day-side (the sub-stellar point is at lon =0, dayside is where the coslon >=0) |
|---|
| 101 | if (coslon(ig) .ge. 0) then |
|---|
| 102 | Trelax(ig,l) = T0(l)+delta_Teq(l)*(ABS(coslon(ig)*coslat(ig))-0.5) |
|---|
| 103 | else !we're on the night-side |
|---|
| 104 | Trelax(ig,l) = T0(l)-0.5*delta_Teq(l) |
|---|
| 105 | endif |
|---|
| 106 | enddo !ig=1,ngrid |
|---|
| 107 | enddo ! l=1,nlayer |
|---|
| 108 | |
|---|
| 109 | ! deallocate T0 and delta_Teq, we don't need them anymore |
|---|
| 110 | if (allocated(T0)) deallocate(T0) |
|---|
| 111 | if (allocated(delta_Teq)) deallocate(delta_Teq) |
|---|
| 112 | endif ! of firstcall |
|---|
| 113 | |
|---|
| 114 | ! call writediagfi(ngrid,"Trelax","Relaxation temperature ","K",3,Trelax) |
|---|
| 115 | ! Calculation of the radiative forcing |
|---|
| 116 | do l=1,nlayer |
|---|
| 117 | do ig=1,ngrid |
|---|
| 118 | if (pplay(ig,l) .le. 1.0e6) then |
|---|
| 119 | ! if pressure is lower than 10 bar |
|---|
| 120 | dtrad(ig,l) = (Trelax(ig,l)-temp(ig,l))/tau_relax(l) |
|---|
| 121 | else |
|---|
| 122 | ! Deeper than 10 bar, no relaxation, dtrad = 0 |
|---|
| 123 | dtrad(ig,l) = 0. |
|---|
| 124 | endif !(pplay(ig,l) .le. 1.e6) |
|---|
| 125 | enddo !ig =1,ngrid |
|---|
| 126 | enddo !l = 1,nlayer |
|---|
| 127 | |
|---|
| 128 | if (lastcall) then |
|---|
| 129 | deallocate(tau_relax) |
|---|
| 130 | deallocate(Trelax) |
|---|
| 131 | endif |
|---|
| 132 | |
|---|
| 133 | end subroutine newtcool_MOCHA |
|---|
| 134 | |
|---|
| 135 | subroutine read_input(nlayer,filename, field) |
|---|
| 136 | |
|---|
| 137 | !======================================== |
|---|
| 138 | ! Purpose |
|---|
| 139 | ! ------- |
|---|
| 140 | ! Read the input file for this module |
|---|
| 141 | ! Each file starts with an integer that should |
|---|
| 142 | ! be equal to nlayer (stops if that's not true) |
|---|
| 143 | ! |
|---|
| 144 | ! Author |
|---|
| 145 | ! ------ |
|---|
| 146 | ! Lucas Teinturier(2024) |
|---|
| 147 | ! |
|---|
| 148 | !======================================== |
|---|
| 149 | |
|---|
| 150 | implicit none |
|---|
| 151 | |
|---|
| 152 | ! Inputs |
|---|
| 153 | integer,intent(in) :: nlayer |
|---|
| 154 | character(100),intent(in) :: filename |
|---|
| 155 | |
|---|
| 156 | ! Output |
|---|
| 157 | real, intent(out) :: field(nlayer) |
|---|
| 158 | |
|---|
| 159 | !! Internal variables |
|---|
| 160 | integer ierr, nline, ii |
|---|
| 161 | |
|---|
| 162 | ! Opening the file |
|---|
| 163 | open(401,form='formatted',status='old',file=trim(filename) ,iostat=ierr) |
|---|
| 164 | if (ierr /=0) then |
|---|
| 165 | print*,"Problem in newton_cooling_hotJ.F90" |
|---|
| 166 | print*,"I have an issue opening file ",trim(filename) |
|---|
| 167 | stop |
|---|
| 168 | endif |
|---|
| 169 | ! Checking that we have the right number of atmospheric layers |
|---|
| 170 | read(401,*) nline |
|---|
| 171 | if (nline /= nlayer) then |
|---|
| 172 | print*,"Error, you're not using the right # of atmospheric layers in ",trim(filename) |
|---|
| 173 | stop |
|---|
| 174 | endif |
|---|
| 175 | ! Now reading the content of the file |
|---|
| 176 | do ii = 1,nline |
|---|
| 177 | read(401,*) field(ii) |
|---|
| 178 | enddo |
|---|
| 179 | close(401) |
|---|
| 180 | end subroutine read_input |
|---|
| 181 | |
|---|
| 182 | end module newton_cooling_hotJ |
|---|
