Changeset 4236 for dynamico_lmdz/simple_physics/phyparam

Timestamp:

Jan 22, 2020, 9:46:46 PM (5 years ago)

Author:

dubos

Message:

simple_physics : some Python bindings

Location:

dynamico_lmdz/simple_physics/phyparam

Files:

• Makefile (modified) (1 diff)
• bash (added)
• bash/beautify.sh (moved) (moved from dynamico_lmdz/simple_physics/bash/beautify.sh)
• bash/concatenate_all_code.sh (moved) (moved from dynamico_lmdz/simple_physics/bash/concatenate_all_code.sh)
• bash/extract_cython.sh (added)
• bash/makedeps.sh (moved) (moved from dynamico_lmdz/simple_physics/bash/makedeps.sh)
• dynphy_lonlat/iniphysiq_mod.F90 (modified) (1 diff)
• physics/error_mod.F90 (modified) (1 diff)
• physics/iniphyparam_mod.F90 (modified) (2 diffs)
• physics/logging.F90 (modified) (7 diffs)
• physics/phyparam_mod.F90 (modified) (5 diffs)
• physics/read_param_mod.F90 (modified) (1 diff)
• physics/restart.F90 (added)
• physics/writefield_mod.F90 (modified) (1 diff)
• python (added)
• python/Extension.help (added)
• python/lib (added)
• python/plugins.pyx (added)
• python/setup.py (added)

dynamico_lmdz/simple_physics/phyparam/Makefile

@@ -16 +16 @@
 all : obj $(OBJECTS) lib/libphyparam.so
 
+py : clean all
+        bash/extract_cython.sh
+        cd python ; rm -rf build *.c ; python setup.py build_ext --inplace
+        cd python ; python -c 'import phyparam as phys ; phys.setup(100.); phys.alloc(30,100) ; phys.coldstart(30, 100.)'
 obj :
         @rm -rf obj include lib xml tmp *~ */*~
         @mkdir obj include lib xml tmp
 
-clean :
+clean : nice
         @rm -rf obj include lib xml tmp *~ */*~
         @mkdir obj include lib xml tmp
-        @../bash/concatenate_all_code.sh $(NAMES) > tmp/all_code
-        @../bash/makedeps.sh tmp/all_code >  Makefile.inc
+        @bash/concatenate_all_code.sh $(NAMES) > tmp/all_code
+        @bash/makedeps.sh tmp/all_code >  Makefile.inc
 
 nice:
-        ../bash/beautify.sh physics/*.F90
+        bash/beautify.sh physics/*.F90
 
 %.so : $(OBJECTS)

dynamico_lmdz/simple_physics/phyparam/dynphy_lonlat/iniphysiq_mod.F90

@@ -87 +87 @@
     CALL setup_phyparam
 
-    CALL iniphyparam(klon_omp,nlayer, &
-         punjours, &
-         pdayref,ptimestep, &
-         prad,pg,pr,pcpp)
+    CALL iniphyparam(ptimestep, punjours, prad, pg, pr, pcpp)
     
     !$OMP END PARALLEL

dynamico_lmdz/simple_physics/phyparam/physics/error_mod.F90

@@ -11 +11 @@
   INTERFACE check_NaN
      MODULE PROCEDURE check_NaN1, check_NaN2
-  END INTERFACE check_NaN
+  END INTERFACE
 
   PUBLIC :: check_NaN

dynamico_lmdz/simple_physics/phyparam/physics/iniphyparam_mod.F90

@@ -10 +10 @@
 CONTAINS
 
-  SUBROUTINE iniphyparam(ngrid,nlayer, &
-       &           punjours,              &
-       &           pdayref,ptimestep,     &
-       &           prad,pg,pr,pcpp)
-    USE callkeys
+  SUBROUTINE read_params(ptimestep) BIND(C, name='phyparam_setup')
+    !$cython header void phyparam_setup(double);
+    !$cython wrapper def setup(timestep) : phy.phyparam_setup(timestep)
+    USE read_param_mod
     USE phys_const, ONLY : planet_rad,g,r,cpp,rcp,dtphys,unjours,mugaz
+    USE astronomy
     USE planet, ONLY : coefir, coefvis
-    USE astronomy
     USE turbulence, ONLY : lmixmin, emin_turb
     USE surface
-    USE read_param_mod
+    USE callkeys
+    REAL, INTENT(IN), VALUE :: ptimestep
 
-    INTEGER, INTENT(IN) :: &
-         ngrid, &       ! Size of the horizontal grid
-         nlayer, &      ! Number of vertical layers.
-         pdayref        ! Day of reference for the simulation
-    REAL, INTENT(IN)  :: ptimestep, prad, pg, pr, pcpp, punjours
-
-    CALL read_param('unjours', 86400., unjours,'unjours')
-    CALL read_param('planet_rad',prad,planet_rad,'planet_rad')
-    CALL read_param('g',9.8           ,g,'g')
-    CALL read_param('cpp',1004.       ,cpp,'cpp')
-    CALL read_param('mugaz',28.       ,mugaz,'mugaz')
+    CALL read_param('planet_rad',6.4e6 ,planet_rad,'planet_rad')
+    CALL read_param('g',9.8            ,g,'g')
+    CALL read_param('cpp',1004.        ,cpp,'cpp')
+    CALL read_param('mugaz',28.        ,mugaz,'mugaz')
     r=perfect_gas_const/mugaz
     rcp=r/cpp
 
-    CALL read_param('year_day',360.   ,year_day,'year_day')
-    CALL read_param('periheli',150.   ,periheli,'periheli')
-    CALL read_param('aphelie',150.    ,aphelie,'aphelie')
-    CALL read_param('peri_day',0.     ,peri_day,'peri_day')
-    CALL read_param('obliquit',23.    ,obliquit,'obliquit')
-    CALL read_param('Cd_mer',.01      ,Cd_mer,'Cd_mer')
-    CALL read_param('Cd_ter',.01      ,Cd_ter,'Cd_ter')
-    CALL read_param('I_mer',3000.    ,I_mer,'I_mer')
-    CALL read_param('I_ter',3000.    ,I_ter,'I_ter')
-    CALL read_param('alb_ter',.112    ,alb_ter,'alb_ter')
-    CALL read_param('alb_mer',.112    ,alb_mer,'alb_mer')
-    CALL read_param('emi_mer',1.      ,emi_mer,'emi_mer')
-    CALL read_param('emi_mer',1.      ,emi_mer,'emi_mer')
+    CALL read_param('unjours', 86400.,  unjours,'unjours')
+    CALL read_param('year_day',360.    ,year_day,'year_day')
+    CALL read_param('periheli',150.    ,periheli,'periheli')
+    CALL read_param('aphelie',150.     ,aphelie,'aphelie')
+    CALL read_param('peri_day',0.      ,peri_day,'peri_day')
+    CALL read_param('obliquit',23.     ,obliquit,'obliquit')
+
+    CALL read_param('Cd_mer',.01       ,Cd_mer,'Cd_mer')
+    CALL read_param('Cd_ter',.01       ,Cd_ter,'Cd_ter')
+    CALL read_param('I_mer',3000.      ,I_mer,'I_mer')
+    CALL read_param('I_ter',3000.      ,I_ter,'I_ter')
+    CALL read_param('alb_ter',.112     ,alb_ter,'alb_ter')
+    CALL read_param('alb_mer',.112     ,alb_mer,'alb_mer')
+    CALL read_param('emi_mer',1.       ,emi_mer,'emi_mer')
+    CALL read_param('emi_mer',1.       ,emi_mer,'emi_mer')
     CALL read_param('emin_turb',1.e-16 ,emin_turb,'emin_turb')
-    CALL read_param('lmixmin',100.    ,lmixmin,'lmixmin')
-    CALL read_param('coefvis',.99     ,coefvis,'coefvis')
-    CALL read_param('coefir',.08      ,coefir,'coefir')
+    CALL read_param('lmixmin',100.     ,lmixmin,'lmixmin')
 
-    CALL read_param('callrad',.true.,callrad,'appel rayonnement')
-    CALL read_param('calldifv',.true.,calldifv,'appel difv')
-    CALL read_param('calladj',.true.,calladj,'appel adj')
-    CALL read_param('callsoil',.true.,callsoil,'appel soil')
-    CALL read_param('season',.true.,season,'with seasonal cycle')
-    CALL read_param('diurnal',.false.,diurnal,'with diurnal cycle')
-    CALL read_param('lverbose',.true.,lverbose,'lverbose')
-    CALL read_param('period_sort',1.,period_sort,'period sorties en jour')
+    CALL read_param('coefvis',.99      ,coefvis,'coefvis')
+    CALL read_param('coefir',.08       ,coefir,'coefir')
+
+    CALL read_param('callrad',  .true.,  callrad,   'appel rayonnement')
+    CALL read_param('calldifv', .true.,  calldifv,  'appel difv')
+    CALL read_param('calladj',  .true.,  calladj,   'appel adj')
+    CALL read_param('callsoil', .true.,  callsoil,  'appel soil')
+    CALL read_param('season',   .true.,  season,    'with seasonal cycle')
+    CALL read_param('diurnal',  .false., diurnal,   'with diurnal cycle')
+    CALL read_param('lverbose', .true.,  lverbose,  'lverbose')
+    CALL read_param('period_sort', 1., period_sort, 'period sorties en jour')
+
+    !   choice of the frequency of the computation of radiations
+    IF(diurnal) THEN
+       iradia=NINT(unjours/(20.*ptimestep))
+    ELSE
+       iradia=NINT(unjours/(4.*ptimestep))
+    ENDIF
+    iradia=1
+
+    dtphys=ptimestep
+  END SUBROUTINE read_params
+
+  SUBROUTINE iniphyparam(ptimestep, punjours, prad, pg, pr, pcpp)
+    USE phys_const, ONLY : planet_rad,g,r,cpp,rcp,dtphys,unjours
+    USE callkeys
+    REAL, INTENT(IN)  :: ptimestep, punjours, prad, pg, pr, pcpp
+
+    CALL read_params(ptimestep)
 
     CALL check_mismatch('unjours', punjours, unjours)
@@ -78 +92 @@
     WRITELOG(*,*) ' Cycle diurne ',diurnal
 
-    !   choice of the frequency of the computation of radiations
-    IF(diurnal) THEN
-       iradia=NINT(punjours/(20.*ptimestep))
-    ELSE
-       iradia=NINT(punjours/(4.*ptimestep))
-    ENDIF
-    iradia=1
-    WRITELOG(*,*) 'unjours',punjours
+    WRITELOG(*,*) 'unjours',unjours
     WRITELOG(*,*) 'The radiative transfer is computed each ', &
          &   iradia,' physical time-step or each ', &
          &   iradia*ptimestep,' seconds'
-
-    dtphys=ptimestep
 
     LOG_INFO('iniphyparam')

dynamico_lmdz/simple_physics/phyparam/physics/logging.F90

@@ -3 +3 @@
   ! see also use_logging.h
   ! macro LOGBUF accumulates log output into logging_buffer
-
   IMPLICIT NONE
   SAVE
@@ -9 +8 @@
 
   INTERFACE  ! Explicit interfaces for plugins
-
      ! Plugin that typically prints all lines in the loggin buffer 'buf' and prepends tags (log level, timestamp, ...)
-     SUBROUTINE plugin(lev, tag, buf)
-       INTEGER, INTENT(IN) :: lev
-       CHARACTER(*), INTENT(IN) :: tag, buf(:)
+     SUBROUTINE plugin(lev, taglen, tag, buflen, bufsize, buf) BIND(C)
+       USE, INTRINSIC :: iso_c_binding, ONLY : c_char, c_null_char, c_int
+       INTEGER(c_int), INTENT(IN), VALUE :: lev, taglen, buflen, bufsize
+       CHARACTER(KIND=c_char), INTENT(IN) :: tag(taglen), buf(buflen, bufsize)
      END SUBROUTINE plugin
 
      ! Plugin that writes into string 'line' information about the gridpoint of index 'index'
-     SUBROUTINE plugin_log_gridpoint(index, line)
-       INTEGER, INTENT(IN) :: index ! index of gridpoint
-       CHARACTER(*), INTENT(OUT) :: line
+     SUBROUTINE plugin_log_gridpoint(index, line_len, line) BIND(C)
+       USE, INTRINSIC :: iso_c_binding, ONLY : c_char, c_null_char, c_int
+       INTEGER(c_int), INTENT(IN), VALUE :: index, line_len ! index of gridpoint, LEN(line)
+       CHARACTER(KIND=c_char), INTENT(OUT) :: line(line_len)
      END SUBROUTINE plugin_log_gridpoint
 
@@ -30 +30 @@
 #ifndef XCODEML
   ! Note F2003/F2008: XCodeML cannot parse procedure pointers
-  PROCEDURE(plugin), POINTER, PUBLIC :: flush_plugin => NULL()
-  PROCEDURE(plugin_log_gridpoint), POINTER, PUBLIC :: log_gridpoint_plugin => NULL()
+  PUBLIC :: flush_plugin, log_gridpoint_plugin
+  PROCEDURE(plugin), POINTER :: flush_plugin => NULL()
+  PROCEDURE(plugin_log_gridpoint), POINTER :: log_gridpoint_plugin => NULL()
 #endif
 
@@ -49 +50 @@
 CONTAINS
 
+  SUBROUTINE set_plugins(flush_plugin_c) BIND(C, name='phyparam_set_plugins_logging')
+    !$cython header void phyparam_set_plugins_logging(void *);
+    USE, INTRINSIC :: ISO_C_BINDING
+    TYPE(C_FUNPTR), INTENT(IN), VALUE :: flush_plugin_c
+    CALL C_F_PROCPOINTER(flush_plugin_c, flush_plugin)
+  END SUBROUTINE set_plugins
+
   SUBROUTINE missing_plugin(name, mod)
     CHARACTER(*), INTENT(IN) :: name, mod
@@ -56 +64 @@
     CALL flush_log(log_level_warn, 'missing_plugin')
   END SUBROUTINE missing_plugin
+
+  !-------------------------------------------------------------------------------------------------
 
   SUBROUTINE flush_log(lev,tag)
@@ -66 +76 @@
        WRITE(*,*) '        see logging.F90'
     END IF
-    IF(logging_lineno>0 .AND. lev<=max_log_level) CALL flush_plugin(lev, TRIM(tag), logging_buf(1:logging_lineno))
+    IF(logging_lineno>0 .AND. lev<=max_log_level) &
+         CALL flush_plugin(lev, LEN(tag), TRIM(tag), linesize, logging_lineno, logging_buf)
     logging_lineno=0
 #endif
   END SUBROUTINE flush_log
 
-  SUBROUTINE default_flush_plugin(lev, tag, buf)
-    INTEGER, INTENT(IN) :: lev
-    CHARACTER(*), INTENT(IN) :: tag, buf(:)
+  SUBROUTINE default_flush_plugin(lev, taglen, tag, buflen, bufsize, buf)
+    USE, INTRINSIC :: iso_c_binding, ONLY : c_char, c_null_char, c_int
+    INTEGER(c_int), INTENT(IN), VALUE :: lev, taglen, buflen, bufsize
+    CHARACTER(KIND=c_char), INTENT(IN) :: tag(taglen), buf(buflen, bufsize)
     CHARACTER(100) :: prefix
+    CHARACTER(buflen+1) :: line
     INTEGER :: i
     WRITE(prefix,*) '[', dbtag(lev), ' ', tag, ']'
-    DO i=1, SIZE(buf)
-       WRITE(*,*) TRIM(prefix), TRIM(buf(i))
+    DO i=1, bufsize
+       WRITE(line,*) buf(:,i)
+       WRITE(*,*) TRIM(prefix), TRIM(line)
     END DO
   END SUBROUTINE default_flush_plugin
+
+  !-------------------------------------------------------------------------------------------------
 
   SUBROUTINE log_gridpoint(index)
@@ -91 +107 @@
        WRITE(*,*) '        see logging.F90'
     END IF
-    CALL log_gridpoint_plugin(index, logging_buf(logging_lineno))
+    CALL log_gridpoint_plugin(index, linesize, logging_buf(logging_lineno))
 #endif
   END SUBROUTINE log_gridpoint
 
-  SUBROUTINE default_log_gridpoint(index, line)
-    INTEGER, INTENT(IN) :: index ! index of gridpoint
-    CHARACTER(*), INTENT(OUT) :: line
+  SUBROUTINE default_log_gridpoint(index, line_len, line)
+    USE, INTRINSIC :: iso_c_binding, ONLY : c_char, c_null_char, c_int
+    INTEGER(c_int), INTENT(IN), VALUE :: index, line_len ! index of gridpoint, LEN(line)
+    CHARACTER(KIND=c_char), INTENT(OUT) :: line(line_len)
     line=' '
   END SUBROUTINE default_log_gridpoint

dynamico_lmdz/simple_physics/phyparam/physics/phyparam_mod.F90

@@ -34 +34 @@
        &            pplev,pplay,pphi,           &
        &            pu,pv,pt,                   &
-       &            pdu,pdv,pdt,pdpsrf)
+       &            pdu,pdv,pdt,pdpsrf) BIND(C, name='phyparam_phyparam')
     USE phys_const, ONLY : g, rcp, r, unjours
     USE surface,    ONLY : soil
@@ -45 +45 @@
     !   20 parameters GCM for planetary atmospheres.
     !   It includes:
-    !   raditive transfer (long and shortwave) for CO2 and dust.
+    !   radiative transfer (long and shortwave) for CO2 and dust.
     !   vertical turbulent mixing
     !   convective adjsutment
+    !   heat diffusion in the soil
     !
     !   author: Frederic Hourdin 15 / 10 /93
     !=======================================================================
 
-    INTEGER, INTENT(IN) ::      &
+    INTEGER, INTENT(IN), VALUE :: &
          ngrid,                 & ! Size of the horizontal grid.
          nlayer                   ! Number of vertical layers.
-    LOGICAL, INTENT(IN) ::      &
+    LOGICAL, INTENT(IN), VALUE  :: &
          firstcall,             & ! True at the first call
          lastcall                 ! True at the last call
-    REAL, INTENT(IN)    ::      &
+    REAL, INTENT(IN), VALUE     ::      &
          rjourvrai,             & ! Number of days counted from the North. Spring equinox
          gmtime,                & ! time of the day in seconds
-         ptimestep,             & ! timestep (s)
+         ptimestep                ! timestep (s)
+    REAL, INTENT(IN) :: &
          pplev(ngrid,nlayer+1), & ! Pressure at interfaces between layers (pa)
          pplay(ngrid,nlayer),   & ! Pressure at the middle of the layers (Pa)
@@ -431 +433 @@
   END SUBROUTINE radiative_tendencies
 
-  SUBROUTINE alloc(ngrid, nlayer)
+  SUBROUTINE alloc(ngrid, nlayer) BIND(C, name='phyparam_alloc')
+    !$cython header void phyparam_alloc(int, int);
+    !$cython wrapper def alloc(ngrid, nlayer) : phy.phyparam_alloc(ngrid, nlayer)
     USE astronomy, ONLY : iniorbit
     USE surface, ONLY : zc,zd
-    INTEGER, INTENT(IN) :: ngrid, nlayer
+    INTEGER, INTENT(IN), VALUE :: ngrid, nlayer
     LOGICAL, PARAMETER :: firstcall=.TRUE.
     ! allocate arrays for internal state
@@ -448 +452 @@
   END SUBROUTINE alloc
 
-  SUBROUTINE precompute
+  SUBROUTINE precompute() BIND(C, name='phyparam_precompute')
+    !$cython header void phyparam_precompute();
+    !$cython wrapper def precompute() : phy.phyparam_precompute()
     USE surface
     ! precompute time-independent arrays
@@ -458 +464 @@
   END SUBROUTINE precompute
 
-  SUBROUTINE coldstart(ngrid, ptimestep)
+  SUBROUTINE coldstart(ngrid, ptimestep) BIND(C, name='phyparam_coldstart')
+    !$cython header void phyparam_coldstart(int, double);
+    !$cython wrapper def coldstart (ngrid, timestep): phy.phyparam_coldstart(ngrid, timestep)
     ! create internal state to start a run without a restart file
     USE surface, ONLY : soil
-    INTEGER, INTENT(IN) :: ngrid
-    REAL, INTENT(iN)    :: ptimestep
+    INTEGER, INTENT(IN), VALUE :: ngrid
+    REAL, INTENT(IN),    VALUE :: ptimestep
     tsurf(:)   = tsoil_init
     tsoil(:,:) = tsoil_init

dynamico_lmdz/simple_physics/phyparam/physics/read_param_mod.F90

@@ -43 +43 @@
   INTERFACE read_param
      PROCEDURE read_paramr, read_parami, read_paramb
-  END INTERFACE read_param
+  END INTERFACE
 
   PUBLIC :: read_param

dynamico_lmdz/simple_physics/phyparam/physics/writefield_mod.F90

@@ -24 +24 @@
   INTERFACE writefield
      PROCEDURE writefield1, writefield2
-  END INTERFACE writefield
+  END INTERFACE
 
   PUBLIC :: writefield