Set up emission routine as physics plugin

.gitignore

@@ -2,5 +2,9 @@ ECMWF*.tar
 dynamico
 XIOS
 Grib/*
+xios_client*
+fort.*
+
 *~
 *.log
+\#*\#

DYNAMICO_emission/Makefile

+
+camelot : F90=mpif90
+camelot : F90FLAGS=-I ../dynamico/inc -i4 -r8 -auto -align all -assume realloc_lhs
+camelot : NCFLAGS=-L/opt/netcdf4/4.4.1.1-parallel/ifort/lib -lnetcdff -Wl,-rpath -Wl,/opt/netcdf4/4.4.1.1-parallel/ifort/lib -lnetcdf
+camelot : ICOFLAGS=-L ../dynamico/lib -licosa -L${MKLROOT}/lib/intel64 -lmkl_intel_lp64 -lmkl_sequential -lmkl_core -lm
+camelot : XIOSFLAGS=-L ../XIOS/lib -licosa -lxios -lstdc++
+camelot : LDFLAGS=$(ICOFLAGS) $(XIOSFLAGS) $(NCFLAGS)
+camelot : all
+
+all :
+	. ../dynamico/arch.env && module list && F90="$(F90)" F90FLAGS="$(F90FLAGS)" LDFLAGS="$(LDFLAGS)" make a.out
+
+clean :
+	@rm -rf obj include bin *~ */*~
+	@mkdir obj include bin
+
+
+%.so : $(OBJECTS)
+	$(F90) -shared $^ -o $@ 
+
+SRC = $(basename $(notdir $@))
+
+a.out : obj/driver.o
+	$(F90) obj/driver.o obj/emission.o $(LDFLAGS)
+	mv a.out bin/DYNAMICO_emission.exe
+
+obj/driver.o : obj/emission.o
+
+obj/%.o: src/%.F90
+	$(F90) -I include -fPIC $(F90FLAGS) -c $<
+	mv $(SRC).o obj
+	mv $(SRC).mod include || true
+

DYNAMICO_emission/src/driver.F90

+PROGRAM driver
+  USE icosa_init_mod, ONLY : icosa_init
+  USE physics_mod, ONLY : init_physics_plugin, physics_plugin
+  USE emission, ONLY : init_physics, physics
+  
+  init_physics_plugin => init_physics
+  physics_plugin => physics
+  
+  CALL icosa_init
+
+END PROGRAM driver
+
+SUBROUTINE initialize_external_physics
+END SUBROUTINE initialize_external_physics
+
+SUBROUTINE external_physics
+END SUBROUTINE external_physics

DYNAMICO_emission/src/emission.F90

+MODULE emission
+  
+  ! FCM gets confused when external modules are USEd at module level
+  ! => USE statements to DYNAMICO modules go into subroutines
+  
+  USE icosa, ONLY : t_field
+  
+  IMPLICIT NONE
+  PRIVATE
+  SAVE
+  
+  REAL, PARAMETER :: oneday = 86400. ! hard-coded                                                                                                                                                         
+  INTEGER, PARAMETER :: log_unit = 15
+  
+  !  TYPE(t_field),POINTER :: f_write2d(:), f_write_llm(:), f_write_llmp1(:)
+  
+  PUBLIC :: init_physics, physics
+  
+CONTAINS
+  
+  SUBROUTINE init_physics
+    ! DYNAMICO
+    USE mpipara,  ONLY : is_mpi_master
+    USE icosa,    ONLY : llm, rstd
+    
+    USE icosa,     ONLY : g, radius, cpp, kappa 
+    USE getin_mod, ONLY : getin
+    USE physics_interface_mod, ONLY : inout => physics_inout
+    REAL(rstd) :: unjours
+
+    WRITE(*,*) 'init_emission called'
+    
+    unjours = 86400.
+    CALL getin('unjours', unjours)
+    
+  END SUBROUTINE init_physics
+  
+  SUBROUTINE physics
+    USE mpipara,  ONLY : is_mpi_master
+    USE icosa,    ONLY : llm
+    USE physics_interface_mod, ONLY : inout => physics_inout
+    REAL :: dps(inout%ngrid), play(inout%ngrid, llm), pphi(inout%ngrid, llm)
+    REAL :: timestep, time, jourvrai, gmtime
+    INTEGER :: l
+    IF(is_mpi_master) WRITE(log_unit,*) 'emission called', SHAPE(inout%p), SHAPE(inout%pk)      
+    
+    timestep = inout%dt_phys
+    time = timestep * inout%it
+    gmtime = time/oneday
+    jourvrai = FLOOR(gmtime)
+    gmtime   = gmtime - jourvrai
+    
+    ! compute pressure and geopotential at full levels
+    CALL compute_play(inout%ngrid, llm, inout%p, play)
+    CALL compute_play(inout%ngrid, llm, inout%geopot, pphi)
+    
+    ! substract surface geopotential
+    DO l=1,llm
+       pphi(:,l) = pphi(:,l) - inout%geopot(:,1)
+    END DO
+    
+    !      ! go
+    !      CALL emission(inout%ngrid,llm,                       &
+    !           &        firstcall,lastcall,                    &
+    !           &        jourvrai, gmtime, timestep,            &
+    !           &        inout%p, play, pphi,                   &
+    !           &        inout%ulon,  inout%ulat,  inout%temp,  &
+    !           &        inout%dulon, inout%dulat, inout%dtemp, dps)
+    
+  END SUBROUTINE physics
+  
+  SUBROUTINE compute_play(ngrid, llm, plev, play)
+    INTEGER, INTENT(IN) :: ngrid, llm
+    REAL, INTENT(IN)    :: plev(ngrid, llm+1) ! pressure at interfaces (half-levels)
+    REAL, INTENT(OUT)   :: play(ngrid, llm)   ! pressure in layers (full levels)
+    INTEGER :: ij, l
+    DO l = 1,llm
+       DO ij = 1,ngrid
+          play(ij,l) = .5*(plev(ij,l)+plev(ij,l+1))
+       END DO
+    END DO
+  END SUBROUTINE compute_play
+  
+  SUBROUTINE plugin_writefield1(name,longname,unit, var)
+    USE physics_interface_mod, ONLY : unpack_field, inout => physics_inout
+    USE output_field_mod, ONLY : output_field
+    CHARACTER(*), INTENT(IN) :: name, longname, unit
+    REAL, INTENT(IN)         :: var(:)
+    WRITE(*,*) TRIM(name), ' : ', TRIM(longname), SHAPE(var), inout%it
+    WRITE(*,*) TRIM(name), ' : ', MINVAL(var), MAXVAL(var)
+!    CALL unpack_field(f_write2d, var)
+!    CALL output_field('phyparam_'//TRIM(name), f_write2d)
+  END SUBROUTINE plugin_writefield1
+  
+  SUBROUTINE plugin_writefield2(name,longname,unit, var)
+    USE physics_interface_mod, ONLY : unpack_field, inout => physics_inout
+    USE output_field_mod, ONLY : output_field
+    USE icosa,    ONLY : llm
+    CHARACTER(*), INTENT(IN) :: name, longname, unit
+    REAL, INTENT(IN)         :: var(:,:)
+    INTEGER :: nlev
+    WRITE(*,*) TRIM(name), ' : ', TRIM(longname), SHAPE(var), inout%it
+    WRITE(*,*) TRIM(name), ' : ', MINVAL(var), MAXVAL(var)
+    nlev = SIZE(var, 2)
+    IF(nlev==llm) THEN
+!       CALL unpack_field(f_write_llm, var)
+!       CALL output_field('phyparam_'//TRIM(name), f_write_llm)
+    ELSEIF(nlev==llm+1) THEN
+!       CALL unpack_field(f_write_llmp1, var)
+!       CALL output_field('phyparam_'//TRIM(name), f_write_llmp1)
+    END IF
+  END SUBROUTINE plugin_writefield2
+  
+END MODULE emission

RUN/gcm.exe

-../dynamico/bin/icosa_gcm.exe
\ No newline at end of file
+../DYNAMICO_emission/bin/DYNAMICO_emission.exe
\ No newline at end of file

RUN/run.def

@@ -44,6 +44,8 @@ nudging_lat_start=20
 nudging_lat_end=70
 
 #-------------- Physics -------------
+physics=plugin
+itau_physics=1
 #---------------- Run ---------------
 run_length=86400
 # run_length=777600