Changeset 4211 for dynamico_lmdz


Ignore:
Timestamp:
Jan 7, 2020, 12:08:48 PM (5 years ago)
Author:
dubos
Message:

simple_physics : cleanup physiq_mod

File:
1 edited

Legend:

Unmodified
Added
Removed

  • dynamico_lmdz/simple_physics/phyparam/param/physiq_mod.F90

    r4178 r4211  
    11! $Id: physiq.F 1565 2011-08-31 12:53:29Z jghattas $
    22MODULE 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
    519
    620CONTAINS
    721
    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   
    25158#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
    134160    ! Declare available vertical axes to be used in output files:   
    135161    CALL wxios_add_vaxis("presnivs", klev, presnivs)
    136 
     162   
    137163    ! Declare calendar and time step
    138164    CALL wxios_set_cal(dtime,"earth_360d",1,1,1,0.0,1,1,1,0.0)
    … …  
    141167    CALL wxios_closedef()
    142168#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
    191185#ifndef CPP_IOIPSL_NO_OUTPUT
    192 if (modulo(itau,iwrite_phys)==0) then
    193   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 endif
    198 #endif
    199 
    200 !XIOS
     186    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
    201195#ifdef CPP_XIOS
    202 !$OMP MASTER
     196    !$OMP MASTER
    203197    !Increment XIOS time
    204198    CALL xios_update_calendar(itau)
    205 !$OMP END MASTER
    206 !$OMP BARRIER
    207 
     199    !$OMP END MASTER
     200    !$OMP BARRIER
     201   
    208202    !Send fields to XIOS: (NB these fields must also be defined as
    209203    ! <field id="..." /> in iodef.xml to be correctly used
    … …  
    214208    CALL histwrite_phy("ps",paprs(:,1))
    215209#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 
    224213END MODULE physiq_mod
Note: See TracChangeset for help on using the changeset viewer.