Changeset 4211 for dynamico_lmdz/simple_physics/phyparam
- Timestamp:
- Jan 7, 2020, 12:08:48 PM (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
dynamico_lmdz/simple_physics/phyparam/param/physiq_mod.F90
r4178 r4211 1 1 ! $Id: physiq.F 1565 2011-08-31 12:53:29Z jghattas $ 2 2 MODULE physiq_mod 3 4 IMPLICIT NONE 3 USE dimphy, only : klon,klev 4 5 IMPLICIT NONE 6 PRIVATE 7 SAVE 8 9 REAL, PARAMETER :: oneday = 86400. ! hard-coded 10 11 INTEGER :: itau=0 ! counter to count number of calls to physics 12 !$OMP THREADPRIVATE(itau) 13 INTEGER :: nid_hist ! NetCDF ID of output file 14 !$OMP THREADPRIVATE(nid_hist) 15 16 INTEGER :: iwrite_phys ! output every iwrite_phys physics step 17 18 PUBLIC :: physiq 5 19 6 20 CONTAINS 7 21 8 SUBROUTINE physiq (nlon,nlev, & 9 & debut,lafin,pdtphys, & 10 & paprs,pplay,pphi,pphis,presnivs, & 11 & u,v,t,qx, & 12 & flxmass_w, & 13 & d_u, d_v, d_t, d_qx, d_ps) 14 15 USE dimphy, only : klon,klev 16 USE infotrac_phy, only : nqtot 17 USE geometry_mod, only : latitude 18 USE comcstphy, only : rg 19 USE iophy, only : histbeg_phy,histwrite_phy 20 USE ioipsl, only : getin,histvert,histdef,histend,ymds2ju 21 USE mod_phys_lmdz_para, only : jj_nb 22 USE phys_state_var_mod, only : phys_state_var_init 23 USE mod_grid_phy_lmdz, ONLY: nbp_lon,nbp_lat 24 22 SUBROUTINE physiq (nlon,nlev, & 23 & debut,lafin,pdtphys, & 24 & paprs,pplay,pphi,pphis,presnivs, & 25 & u,v,t,qx, & 26 & flxmass_w, & 27 & d_u, d_v, d_t, d_qx, d_ps) 28 29 USE infotrac_phy, only : nqtot 30 ! 31 ! Routine argument: 32 ! 33 integer,intent(in) :: nlon ! number of atmospheric colums 34 integer,intent(in) :: nlev ! number of vertical levels (should be =klev) 35 logical,intent(in) :: debut ! signals first call to physics 36 logical,intent(in) :: lafin ! signals last call to physics 37 real,intent(in) :: pdtphys ! physics time step (s) 38 real,intent(in) :: paprs(klon,klev+1) ! interlayer pressure (Pa) 39 real,intent(in) :: pplay(klon,klev) ! mid-layer pressure (Pa) 40 real,intent(in) :: pphi(klon,klev) ! geopotential at mid-layer 41 real,intent(in) :: pphis(klon) ! surface geopotential 42 real,intent(in) :: presnivs(klev) ! pseudo-pressure (Pa) of mid-layers 43 real,intent(in) :: u(klon,klev) ! eastward zonal wind (m/s) 44 real,intent(in) :: v(klon,klev) ! northward meridional wind (m/s) 45 real,intent(in) :: t(klon,klev) ! temperature (K) 46 real,intent(in) :: qx(klon,klev,nqtot) ! tracers (.../kg_air) 47 real,intent(in) :: flxmass_w(klon,klev) ! vertical mass flux 48 real,intent(out) :: d_u(klon,klev) ! physics tendency on u (m/s/s) 49 real,intent(out) :: d_v(klon,klev) ! physics tendency on v (m/s/s) 50 real,intent(out) :: d_t(klon,klev) ! physics tendency on t (K/s) 51 real,intent(out) :: d_qx(klon,klev,nqtot) ! physics tendency on tracers 52 real,intent(out) :: d_ps(klon) ! physics tendency on surface pressure 53 54 real :: temp_newton(klon,klev) 55 integer :: k 56 logical, save :: first=.true. 57 !$OMP THREADPRIVATE(first) 58 59 ! For I/Os 60 REAL,SAVE :: rjourvrai=0. 61 REAL,SAVE :: gmtime=0. 62 !$OMP THREADPRIVATE(rjourvrai,gmtime) 63 64 65 ! initializations 66 if (debut) CALL init_physiq_mod(pdtphys, presnivs) ! Things to do only for the first call to physics 67 68 ! increment local time counter itau 69 itau=itau+1 70 gmtime=gmtime+pdtphys/oneday 71 IF (gmtime>=1.) THEN 72 gmtime=gmtime-1. 73 rjourvrai=rjourvrai+1. 74 ENDIF 75 76 CALL phyparam(klon,klev,nqtot, & 77 debut,lafin, & 78 rjourvrai,gmtime,pdtphys, & 79 paprs,pplay,pphi,pphis,presnivs, & 80 u,v,t,qx, & 81 flxmass_w, & 82 d_u,d_v,d_t,d_qx,d_ps) 83 84 print*,'PHYDEV: itau=',itau 85 86 CALL write_hist(paprs, t, u, v) ! write some outputs 87 88 ! if lastcall, then it is time to write "restartphy.nc" file 89 IF (lafin) CALL phyredem("restartphy.nc") 90 91 END SUBROUTINE physiq 92 93 SUBROUTINE init_physiq_mod(dtime, presnivs) 94 USE iophy, only : histbeg_phy 95 USE phys_state_var_mod, only : phys_state_var_init 96 USE ioipsl, only : getin,histvert,histdef,histend,ymds2ju 97 USE mod_phys_lmdz_para, only : jj_nb 98 USE mod_grid_phy_lmdz, ONLY: nbp_lon,nbp_lat 99 REAL, INTENT(IN) :: dtime 100 REAL, INTENT(IN) :: presnivs(klev) ! pseudo-pressure (Pa) of mid-layers 101 102 ! For I/Os 103 INTEGER, PARAMETER :: itau0=0, annee0=1979, month=1, dayref=1 104 REAL, PARAMETER :: hour=0.0 105 INTEGER :: nid_hori ! NetCDF ID of horizontal coordinate 106 INTEGER :: nid_vert ! NetCDF ID of vertical coordinate 107 REAL :: t_ops ! frequency of the IOIPSL operations (eg average over...) 108 REAL :: t_wrt ! frequency of the IOIPSL outputs 109 REAL :: zjulian 110 111 ! load initial conditions for physics (including the grid) 112 call phys_state_var_init() ! some initializations, required before calling phyetat0 113 call phyetat0("startphy.nc") 114 115 ! Initialize outputs: 116 117 ! NB: getin() is not threadsafe; only one thread should call it. 118 !$OMP MASTER 119 iwrite_phys=1 !default: output every physics timestep 120 call getin("iwrite_phys",iwrite_phys) 121 !$OMP END MASTER 122 !$OMP BARRIER 123 124 t_ops=dtime*iwrite_phys ! frequency of the IOIPSL operation 125 t_wrt=dtime*iwrite_phys ! frequency of the outputs in the file 126 127 ! compute zjulian for annee0=1979 and month=1 dayref=1 and hour=0.0 128 CALL ymds2ju(annee0, month, dayref, hour, zjulian) 129 #ifndef CPP_IOIPSL_NO_OUTPUT 130 ! Initialize IOIPSL output file 131 call histbeg_phy("histins.nc",itau0,zjulian,dtime,nid_hori,nid_hist) 132 #endif 133 134 !$OMP MASTER 135 136 #ifndef CPP_IOIPSL_NO_OUTPUT 137 ! IOIPSL 138 ! define vertical coordinate 139 call histvert(nid_hist,"presnivs","Vertical levels","Pa",klev, & 140 presnivs,nid_vert,'down') 141 ! define variables which will be written in "histins.nc" file 142 call histdef(nid_hist,'temperature','Atmospheric temperature','K', & 143 nbp_lon,jj_nb,nid_hori,klev,1,klev,nid_vert,32, & 144 'inst(X)',t_ops,t_wrt) 145 call histdef(nid_hist,'u','Eastward Zonal Wind','m/s', & 146 nbp_lon,jj_nb,nid_hori,klev,1,klev,nid_vert,32, & 147 'inst(X)',t_ops,t_wrt) 148 call histdef(nid_hist,'v','Northward Meridional Wind','m/s', & 149 nbp_lon,jj_nb,nid_hori,klev,1,klev,nid_vert,32, & 150 'inst(X)',t_ops,t_wrt) 151 call histdef(nid_hist,'ps','Surface Pressure','Pa', & 152 nbp_lon,jj_nb,nid_hori,1,1,1,nid_vert,32, & 153 'inst(X)',t_ops,t_wrt) 154 ! end definition sequence 155 call histend(nid_hist) 156 #endif 157 25 158 #ifdef CPP_XIOS 26 USE xios, ONLY: xios_update_calendar 27 USE wxios, only: wxios_add_vaxis, wxios_set_cal, wxios_closedef 28 USE iophy, ONLY: histwrite_phy 29 #endif 30 31 IMPLICIT none 32 ! 33 ! Routine argument: 34 ! 35 integer,intent(in) :: nlon ! number of atmospheric colums 36 integer,intent(in) :: nlev ! number of vertical levels (should be =klev) 37 logical,intent(in) :: debut ! signals first call to physics 38 logical,intent(in) :: lafin ! signals last call to physics 39 real,intent(in) :: pdtphys ! physics time step (s) 40 real,intent(in) :: paprs(klon,klev+1) ! interlayer pressure (Pa) 41 real,intent(in) :: pplay(klon,klev) ! mid-layer pressure (Pa) 42 real,intent(in) :: pphi(klon,klev) ! geopotential at mid-layer 43 real,intent(in) :: pphis(klon) ! surface geopotential 44 real,intent(in) :: presnivs(klev) ! pseudo-pressure (Pa) of mid-layers 45 real,intent(in) :: u(klon,klev) ! eastward zonal wind (m/s) 46 real,intent(in) :: v(klon,klev) ! northward meridional wind (m/s) 47 real,intent(in) :: t(klon,klev) ! temperature (K) 48 real,intent(in) :: qx(klon,klev,nqtot) ! tracers (.../kg_air) 49 real,intent(in) :: flxmass_w(klon,klev) ! vertical mass flux 50 real,intent(out) :: d_u(klon,klev) ! physics tendency on u (m/s/s) 51 real,intent(out) :: d_v(klon,klev) ! physics tendency on v (m/s/s) 52 real,intent(out) :: d_t(klon,klev) ! physics tendency on t (K/s) 53 real,intent(out) :: d_qx(klon,klev,nqtot) ! physics tendency on tracers 54 real,intent(out) :: d_ps(klon) ! physics tendency on surface pressure 55 56 integer,save :: itau=0 ! counter to count number of calls to physics 57 !$OMP THREADPRIVATE(itau) 58 real :: temp_newton(klon,klev) 59 integer :: k 60 logical, save :: first=.true. 61 !$OMP THREADPRIVATE(first) 62 63 ! For I/Os 64 integer :: itau0 65 real :: zjulian 66 real :: dtime 67 integer :: nhori ! horizontal coordinate ID 68 integer,save :: nid_hist ! output file ID 69 !$OMP THREADPRIVATE(nid_hist) 70 integer :: zvertid ! vertical coordinate ID 71 integer,save :: iwrite_phys ! output every iwrite_phys physics step 72 !$OMP THREADPRIVATE(iwrite_phys) 73 integer,save :: iwrite_phys_omp ! intermediate variable to read iwrite_phys 74 ! (must be shared by all threads) 75 real :: t_ops ! frequency of the IOIPSL operations (eg average over...) 76 real :: t_wrt ! frequency of the IOIPSL outputs 77 REAL,SAVE :: rjourvrai=0. 78 REAL,SAVE :: gmtime=0. 79 !$OMP THREADPRIVATE(rjourvrai,gmtime) 80 81 82 ! initializations 83 if (debut) then ! Things to do only for the first call to physics 84 ! load initial conditions for physics (including the grid) 85 call phys_state_var_init() ! some initializations, required before calling phyetat0 86 call phyetat0("startphy.nc") 87 88 ! Initialize outputs: 89 itau0=0 90 !$OMP MASTER 91 iwrite_phys_omp=1 !default: output every physics timestep 92 ! NB: getin() is not threadsafe; only one thread should call it. 93 call getin("iwrite_phys",iwrite_phys_omp) 94 !$OMP END MASTER 95 !$OMP BARRIER 96 iwrite_phys=iwrite_phys_omp 97 t_ops=pdtphys*iwrite_phys ! frequency of the IOIPSL operation 98 t_wrt=pdtphys*iwrite_phys ! frequency of the outputs in the file 99 ! compute zjulian for annee0=1979 and month=1 dayref=1 and hour=0.0 100 !CALL ymds2ju(annee0, month, dayref, hour, zjulian) 101 call ymds2ju(1979, 1, 1, 0.0, zjulian) 102 dtime=pdtphys 103 #ifndef CPP_IOIPSL_NO_OUTPUT 104 ! Initialize IOIPSL output file 105 call histbeg_phy("histins.nc",itau0,zjulian,dtime,nhori,nid_hist) 106 #endif 107 108 !$OMP MASTER 109 110 #ifndef CPP_IOIPSL_NO_OUTPUT 111 ! IOIPSL 112 ! define vertical coordinate 113 call histvert(nid_hist,"presnivs","Vertical levels","Pa",klev, & 114 presnivs,zvertid,'down') 115 ! define variables which will be written in "histins.nc" file 116 call histdef(nid_hist,'temperature','Atmospheric temperature','K', & 117 nbp_lon,jj_nb,nhori,klev,1,klev,zvertid,32, & 118 'inst(X)',t_ops,t_wrt) 119 call histdef(nid_hist,'u','Eastward Zonal Wind','m/s', & 120 nbp_lon,jj_nb,nhori,klev,1,klev,zvertid,32, & 121 'inst(X)',t_ops,t_wrt) 122 call histdef(nid_hist,'v','Northward Meridional Wind','m/s', & 123 nbp_lon,jj_nb,nhori,klev,1,klev,zvertid,32, & 124 'inst(X)',t_ops,t_wrt) 125 call histdef(nid_hist,'ps','Surface Pressure','Pa', & 126 nbp_lon,jj_nb,nhori,1,1,1,zvertid,32, & 127 'inst(X)',t_ops,t_wrt) 128 ! end definition sequence 129 call histend(nid_hist) 130 #endif 131 132 #ifdef CPP_XIOS 133 !XIOS 159 !XIOS 134 160 ! Declare available vertical axes to be used in output files: 135 161 CALL wxios_add_vaxis("presnivs", klev, presnivs) 136 162 137 163 ! Declare calendar and time step 138 164 CALL wxios_set_cal(dtime,"earth_360d",1,1,1,0.0,1,1,1,0.0) … … 141 167 CALL wxios_closedef() 142 168 #endif 143 !$OMP END MASTER 144 !$OMP BARRIER 145 endif ! of if (debut) 146 147 ! increment local time counter itau 148 itau=itau+1 149 gmtime=gmtime+pdtphys/86400. 150 IF (gmtime>=1.) THEN 151 gmtime=gmtime-1. 152 rjourvrai=rjourvrai+1. 153 ENDIF 154 155 IF ( 1 == 0 ) THEN 156 ! set all tendencies to zero 157 d_u(1:klon,1:klev)=0. 158 d_v(1:klon,1:klev)=0. 159 d_t(1:klon,1:klev)=0. 160 d_qx(1:klon,1:klev,1:nqtot)=0. 161 d_ps(1:klon)=0. 162 163 ! compute tendencies to return to the dynamics: 164 ! "friction" on the first layer 165 d_u(1:klon,1)=-u(1:klon,1)/86400. 166 d_v(1:klon,1)=-v(1:klon,1)/86400. 167 ! newtonian relaxation towards temp_newton() 168 do k=1,klev 169 temp_newton(1:klon,k)=280.+cos(latitude(1:klon))*40.-pphi(1:klon,k)/rg*6.e-3 170 d_t(1:klon,k)=(temp_newton(1:klon,k)-t(1:klon,k))/1.e5 171 enddo 172 173 ELSE 174 175 CALL phyparam(klon,klev,nqtot, & 176 debut,lafin, & 177 rjourvrai,gmtime,pdtphys, & 178 paprs,pplay,pphi,pphis,presnivs, & 179 u,v,t,qx, & 180 flxmass_w, & 181 d_u,d_v,d_t,d_qx,d_ps) 182 183 184 ENDIF 185 186 187 print*,'PHYDEV: itau=',itau 188 189 ! write some outputs: 190 ! IOIPSL 169 !$OMP END MASTER 170 !$OMP BARRIER 171 END SUBROUTINE init_physiq_mod 172 173 SUBROUTINE write_hist(paprs, t, u, v) 174 USE iophy, only : histwrite_phy 175 #ifdef CPP_XIOS 176 USE xios, ONLY: xios_update_calendar 177 USE wxios, only: wxios_add_vaxis, wxios_set_cal, wxios_closedef 178 USE iophy, ONLY: histwrite_phy 179 #endif 180 REAL, INTENT(IN) :: paprs(klon,klev+1), & ! interlayer pressure (Pa) 181 & t(klon,klev), & ! temperature (K) 182 & u(klon,klev), & ! eastward zonal wind (m/s) 183 & v(klon,klev) ! northward meridional wind (m/s) 184 ! output using IOIPSL 191 185 #ifndef CPP_IOIPSL_NO_OUTPUT 192 if (modulo(itau,iwrite_phys)==0) then193 call histwrite_phy(nid_hist,.false.,"temperature",itau,t)194 call histwrite_phy(nid_hist,.false.,"u",itau,u)195 call histwrite_phy(nid_hist,.false.,"v",itau,v)196 call histwrite_phy(nid_hist,.false.,"ps",itau,paprs(:,1))197 endif198 #endif 199 200 !XIOS186 if (modulo(itau,iwrite_phys)==0) then 187 call histwrite_phy(nid_hist,.false.,"temperature",itau,t) 188 call histwrite_phy(nid_hist,.false.,"u",itau,u) 189 call histwrite_phy(nid_hist,.false.,"v",itau,v) 190 call histwrite_phy(nid_hist,.false.,"ps",itau,paprs(:,1)) 191 endif 192 #endif 193 194 ! output using XIOS 201 195 #ifdef CPP_XIOS 202 !$OMP MASTER196 !$OMP MASTER 203 197 !Increment XIOS time 204 198 CALL xios_update_calendar(itau) 205 !$OMP END MASTER206 !$OMP BARRIER207 199 !$OMP END MASTER 200 !$OMP BARRIER 201 208 202 !Send fields to XIOS: (NB these fields must also be defined as 209 203 ! <field id="..." /> in iodef.xml to be correctly used … … 214 208 CALL histwrite_phy("ps",paprs(:,1)) 215 209 #endif 216 217 ! if lastcall, then it is time to write "restartphy.nc" file 218 if (lafin) then 219 call phyredem("restartphy.nc") 220 endif 221 222 end subroutine physiq 223 210 211 END SUBROUTINE write_hist 212 224 213 END MODULE physiq_mod
Note: See TracChangeset
for help on using the changeset viewer.
