! $Id: physiq.F 1565 2011-08-31 12:53:29Z jghattas $ MODULE physiq_mod IMPLICIT NONE CONTAINS SUBROUTINE physiq(nlon,nlev, & debut,lafin,pdtphys, & paprs,pplay,pphi,pphis,presnivs, & u,v,t,qx, & flxmass_w, & d_u, d_v, d_t, d_qx, d_ps) USE dimphy, ONLY: klon,klev USE infotrac_phy, ONLY: nqtot USE lmdz_geometry, ONLY: latitude USE comcstphy, ONLY: rg USE iophy, ONLY: histbeg_phy,histwrite_phy USE ioipsl, ONLY: getin,histvert,histdef,histend,ymds2ju USE lmdz_phys_para, ONLY: jj_nb USE phys_state_var_mod, ONLY: phys_state_var_init USE lmdz_grid_phy, ONLY: nbp_lon,nbp_lat USE lmdz_xios, ONLY: xios_update_calendar, using_xios USE lmdz_wxios, ONLY: wxios_add_vaxis, wxios_set_cal, wxios_closedef USE iophy, ONLY: histwrite_phy IMPLICIT none ! Routine argument: INTEGER,INTENT(IN) :: nlon ! number of atmospheric colums INTEGER,INTENT(IN) :: nlev ! number of vertical levels (should be =klev) logical,INTENT(IN) :: debut ! signals first CALL to physics logical,INTENT(IN) :: lafin ! signals last CALL to physics REAL,INTENT(IN) :: pdtphys ! physics time step (s) REAL,INTENT(IN) :: paprs(klon,klev+1) ! interlayer pressure (Pa) REAL,INTENT(IN) :: pplay(klon,klev) ! mid-layer pressure (Pa) REAL,INTENT(IN) :: pphi(klon,klev) ! geopotential at mid-layer REAL,INTENT(IN) :: pphis(klon) ! surface geopotential REAL,INTENT(IN) :: presnivs(klev) ! pseudo-pressure (Pa) of mid-layers REAL,INTENT(IN) :: u(klon,klev) ! eastward zonal wind (m/s) REAL,INTENT(IN) :: v(klon,klev) ! northward meridional wind (m/s) REAL,INTENT(IN) :: t(klon,klev) ! temperature (K) REAL,INTENT(IN) :: qx(klon,klev,nqtot) ! tracers (.../kg_air) REAL,INTENT(IN) :: flxmass_w(klon,klev) ! vertical mass flux REAL,INTENT(OUT) :: d_u(klon,klev) ! physics tendency on u (m/s/s) REAL,INTENT(OUT) :: d_v(klon,klev) ! physics tendency on v (m/s/s) REAL,INTENT(OUT) :: d_t(klon,klev) ! physics tendency on t (K/s) REAL,INTENT(OUT) :: d_qx(klon,klev,nqtot) ! physics tendency on tracers REAL,INTENT(OUT) :: d_ps(klon) ! physics tendency on surface pressure INTEGER,save :: itau=0 ! counter to count number of calls to physics !$OMP THREADPRIVATE(itau) REAL :: temp_newton(klon,klev) INTEGER :: k logical, save :: first=.TRUE. !$OMP THREADPRIVATE(first) ! For I/Os INTEGER :: itau0 REAL :: zjulian REAL :: dtime INTEGER :: nhori ! horizontal coordinate ID INTEGER,save :: nid_hist ! output file ID !$OMP THREADPRIVATE(nid_hist) INTEGER :: zvertid ! vertical coordinate ID INTEGER,save :: iwrite_phys ! output every iwrite_phys physics step !$OMP THREADPRIVATE(iwrite_phys) INTEGER,save :: iwrite_phys_omp ! intermediate variable to read iwrite_phys ! (must be shared by all threads) REAL :: t_ops ! frequency of the IOIPSL operations (eg average over...) REAL :: t_wrt ! frequency of the IOIPSL outputs ! initializations IF (debut) then ! Things to do only for the first CALL to physics ! load initial conditions for physics (including the grid) CALL phys_state_var_init() ! some initializations, required before calling phyetat0 CALL phyetat0("startphy.nc") ! Initialize outputs: itau0=0 !$OMP MASTER iwrite_phys_omp=1 !default: output every physics timestep ! NB: getin() is not threadsafe; only one thread should CALL it. CALL getin("iwrite_phys",iwrite_phys_omp) !$OMP END MASTER !$OMP BARRIER iwrite_phys=iwrite_phys_omp t_ops=pdtphys*iwrite_phys ! frequency of the IOIPSL operation t_wrt=pdtphys*iwrite_phys ! frequency of the outputs in the file ! compute zjulian for annee0=1979 and month=1 dayref=1 and hour=0.0 !CALL ymds2ju(annee0, month, dayref, hour, zjulian) CALL ymds2ju(1979, 1, 1, 0.0, zjulian) dtime=pdtphys #ifndef CPP_IOIPSL_NO_OUTPUT ! Initialize IOIPSL output file CALL histbeg_phy("histins.nc",itau0,zjulian,dtime,nhori,nid_hist) #endif !$OMP MASTER #ifndef CPP_IOIPSL_NO_OUTPUT ! IOIPSL ! define vertical coordinate CALL histvert(nid_hist,"presnivs","Vertical levels","Pa",klev, & presnivs,zvertid,'down') ! define variables which will be written in "histins.nc" file CALL histdef(nid_hist,'temperature','Atmospheric temperature','K', & nbp_lon,jj_nb,nhori,klev,1,klev,zvertid,32, & 'inst(X)',t_ops,t_wrt) CALL histdef(nid_hist,'u','Eastward Zonal Wind','m/s', & nbp_lon,jj_nb,nhori,klev,1,klev,zvertid,32, & 'inst(X)',t_ops,t_wrt) CALL histdef(nid_hist,'v','Northward Meridional Wind','m/s', & nbp_lon,jj_nb,nhori,klev,1,klev,zvertid,32, & 'inst(X)',t_ops,t_wrt) CALL histdef(nid_hist,'ps','Surface Pressure','Pa', & nbp_lon,jj_nb,nhori,1,1,1,zvertid,32, & 'inst(X)',t_ops,t_wrt) ! end definition sequence CALL histend(nid_hist) #endif IF (using_xios) THEN !XIOS ! Declare available vertical axes to be used in output files: CALL wxios_add_vaxis("presnivs", klev, presnivs) ! Declare calendar and time step CALL wxios_set_cal(dtime,"earth_360d",1,1,1,0.0,1,1,1,0.0) !Finalize the context: CALL wxios_closedef() ENDIF !$OMP END MASTER !$OMP BARRIER END IF ! of if (debut) ! increment local time counter itau itau=itau+1 ! set all tendencies to zero d_u(1:klon,1:klev)=0. d_v(1:klon,1:klev)=0. d_t(1:klon,1:klev)=0. d_qx(1:klon,1:klev,1:nqtot)=0. d_ps(1:klon)=0. ! compute tendencies to return to the dynamics: ! "friction" on the first layer d_u(1:klon,1)=-u(1:klon,1)/86400. d_v(1:klon,1)=-v(1:klon,1)/86400. ! newtonian relaxation towards temp_newton() DO k=1,klev temp_newton(1:klon,k)=280.+cos(latitude(1:klon))*40.-pphi(1:klon,k)/rg*6.e-3 d_t(1:klon,k)=(temp_newton(1:klon,k)-t(1:klon,k))/1.e5 END DO PRINT*,'PHYDEV: itau=',itau ! write some outputs: ! IOIPSL #ifndef CPP_IOIPSL_NO_OUTPUT IF (modulo(itau,iwrite_phys)==0) THEN CALL histwrite_phy(nid_hist,.FALSE.,"temperature",itau,t) CALL histwrite_phy(nid_hist,.FALSE.,"u",itau,u) CALL histwrite_phy(nid_hist,.FALSE.,"v",itau,v) CALL histwrite_phy(nid_hist,.FALSE.,"ps",itau,paprs(:,1)) END IF #endif !XIOS IF (using_xios) THEN !$OMP MASTER !Increment XIOS time CALL xios_update_calendar(itau) !$OMP END MASTER !$OMP BARRIER !Send fields to XIOS: (NB these fields must also be defined as ! in iodef.xml to be correctly used CALL histwrite_phy("temperature",t) CALL histwrite_phy("temp_newton",temp_newton) CALL histwrite_phy("u",u) CALL histwrite_phy("v",v) CALL histwrite_phy("ps",paprs(:,1)) ENDIF ! if lastcall, then it is time to write "restartphy.nc" file IF (lafin) THEN CALL phyredem("restartphy.nc") END IF END SUBROUTINE physiq END MODULE physiq_mod