Index: /codes/icosagcm/trunk/src/physics.f90
===================================================================
--- /codes/icosagcm/trunk/src/physics.f90	(revision 81)
+++ /codes/icosagcm/trunk/src/physics.f90	(revision 81)
@@ -0,0 +1,54 @@
+MODULE physics_mod
+
+  CHARACTER(LEN=255) :: physics_type="none"
+
+
+CONTAINS
+
+  SUBROUTINE init_physics(dt)
+  USE icosa
+  USE physics_dcmip_mod, init_physics_dcmip=>init_physics
+  IMPLICIT NONE
+    REAL(rstd),INTENT(IN) :: dt
+    
+    CALL getin("physics",physics_type)
+    
+    SELECT CASE(TRIM(physics_type))
+      CASE ('none')
+    
+      CASE ('dcmip')
+        CALL init_physics_dcmip(dt)
+      
+      CASE DEFAULT
+         PRINT*, 'Bad selector for variable physics <',physics_type, &
+              '> options are <none>, <dcmip>,'
+         STOP
+    END SELECT
+    
+  END SUBROUTINE init_physics
+
+  SUBROUTINE physics(f_phis, f_ps, f_theta_rhodz, f_ue, f_q)
+  USE icosa
+  USE physics_dcmip_mod, physics_dcmip=>physics
+  IMPLICIT NONE
+    TYPE(t_field),POINTER :: f_phis(:)
+    TYPE(t_field),POINTER :: f_ps(:)
+    TYPE(t_field),POINTER :: f_theta_rhodz(:)
+    TYPE(t_field),POINTER :: f_ue(:)
+    TYPE(t_field),POINTER :: f_q(:)
+    
+    SELECT CASE(TRIM(physics_type))
+      CASE ('none')
+    
+      CASE ('dcmip')
+        CALL physics_dcmip(f_phis, f_ps, f_theta_rhodz, f_ue, f_q)
+      
+      CASE DEFAULT
+         PRINT*, 'Bad selector for variable physics <',physics_type, &
+              '> options are <none>, <dcmip>,'
+         STOP
+    END SELECT
+    
+  END SUBROUTINE physics
+
+END MODULE physics_mod
Index: /codes/icosagcm/trunk/src/physics_dcmip.f90
===================================================================
--- /codes/icosagcm/trunk/src/physics_dcmip.f90	(revision 81)
+++ /codes/icosagcm/trunk/src/physics_dcmip.f90	(revision 81)
@@ -0,0 +1,624 @@
+MODULE physics_dcmip_mod
+  USE ICOSA
+  PRIVATE
+  
+  REAL(rstd),SAVE :: dt
+  INTEGER,SAVE :: testcase
+  PUBLIC init_physics, physics
+  TYPE(t_field),POINTER :: f_out_i(:)
+    REAL(rstd),POINTER :: out_i(:,:)
+
+CONTAINS
+
+  SUBROUTINE init_physics(dt_in)
+  IMPLICIT NONE
+    REAL(rstd),INTENT(IN) :: dt_in
+
+    dt=dt_in
+    testcase=1
+    CALL getin("dcmip_physics",testcase)
+    CALL allocate_field(f_out_i,field_t,type_real,llm)
+        
+  END SUBROUTINE init_physics
+
+
+  SUBROUTINE physics( f_phis, f_ps, f_theta_rhodz, f_ue, f_q)
+  USE icosa
+  IMPLICIT NONE
+    TYPE(t_field),POINTER :: f_phis(:)
+    TYPE(t_field),POINTER :: f_ps(:)
+    TYPE(t_field),POINTER :: f_theta_rhodz(:)
+    TYPE(t_field),POINTER :: f_ue(:)
+    TYPE(t_field),POINTER :: f_q(:)
+  
+    REAL(rstd),POINTER :: phis(:)
+    REAL(rstd),POINTER :: ps(:)
+    REAL(rstd),POINTER :: theta_rhodz(:,:)
+    REAL(rstd),POINTER :: ue(:,:)
+    REAL(rstd),POINTER :: q(:,:,:)
+    INTEGER :: ind
+
+    DO ind=1,ndomain
+      CALL swap_dimensions(ind)
+      CALL swap_geometry(ind)
+      
+      phis=f_phis(ind)
+      ps=f_ps(ind)
+      theta_rhodz=f_theta_rhodz(ind)
+      ue=f_ue(ind)
+      q=f_q(ind)
+      out_i=f_out_i(ind)
+    
+      CALL compute_physics( phis, ps, theta_rhodz, ue, q(:,:,3))
+  
+    ENDDO
+!    CALL writefield("out_i",f_out_i)
+
+  END SUBROUTINE physics
+  
+  SUBROUTINE compute_physics(phis, ps, theta_rhodz, ue, q)
+  USE icosa
+  USE pression_mod
+  USE exner_mod
+  USE theta2theta_rhodz_mod
+  USE geopotential_mod
+  USE wind_mod
+  IMPLICIT NONE
+    REAL(rstd) :: phis(iim*jjm)
+    REAL(rstd) :: ps(iim*jjm)
+    REAL(rstd) :: theta_rhodz(iim*jjm,llm)
+    REAL(rstd) :: ue(3*iim*jjm,llm)
+    REAL(rstd) :: q(iim*jjm,llm)
+
+    REAL(rstd) :: p(iim*jjm,llm+1)
+    REAL(rstd) :: pks(iim*jjm)
+    REAL(rstd) :: pk(iim*jjm,llm)
+    REAL(rstd) :: phi(iim*jjm,llm)
+    REAL(rstd) :: T(iim*jjm,llm)
+    REAL(rstd) :: theta(iim*jjm,llm)
+
+   REAL(rstd) :: uc(iim*jjm,3,llm)
+   REAL(rstd) :: u(iim*jjm,llm)
+   REAL(rstd) :: v(iim*jjm,llm)
+    REAL(rstd) :: ufi(iim*jjm,llm)
+    REAL(rstd) :: vfi(iim*jjm,llm)
+    REAL(rstd) :: lat(iim*jjm)
+    REAL(rstd) :: lon
+    REAL(rstd) :: pmid(iim*jjm,llm)
+    REAL(rstd) :: pdel(iim*jjm,llm)
+    REAL(rstd) :: precl(iim*jjm)
+    INTEGER :: i,j,l,ij
+     
+    PRINT *,'Entering in DCMIP physics'    
+    CALL compute_pression(ps,p,0)
+    CALL compute_exner(ps,p,pks,pk,0)
+    CALL compute_theta_rhodz2theta(ps,theta_rhodz,theta,0)
+    CALL compute_geopotential(phis,pks,pk,theta,phi,0)
+    CALL compute_theta_rhodz2temperature(ps,theta_rhodz,T,0)
+   CALL compute_wind_centered(ue,uc)
+   CALL compute_wind_centered_lonlat_compound(uc, u, v)
+
+    DO j=jj_begin,jj_end
+      DO i=ii_begin,ii_end
+        ij=(j-1)*iim+i
+        CALL xyz2lonlat(xyz_i(ij,:),lon,lat(ij)) 
+      ENDDO
+    ENDDO
+   
+    DO l=1,llm
+      DO j=jj_begin,jj_end
+        DO i=ii_begin,ii_end
+          ij=(j-1)*iim+i
+          pmid(ij,l)=0.5*(p(ij,l)+p(ij,l+1)) 
+        ENDDO
+      ENDDO
+    ENDDO
+    
+    DO l=1,llm
+      DO j=jj_begin,jj_end
+        DO i=ii_begin,ii_end
+          ij=(j-1)*iim+i
+          pdel(ij,l)=p(ij,l+1)-p(ij,l)
+        ENDDO
+      ENDDO
+    ENDDO       
+       
+    
+    ufi=u
+    vfi=v    
+    
+    CALL simple_physics(iim*jjm, llm, dt, lat, t, q , ufi, vfi, pmid, p, pdel, 1/pdel, ps, precl, testcase) 
+    
+        
+    ufi=ufi-u
+    vfi=vfi-v
+
+    DO l=1,llm
+      DO j=jj_begin,jj_end
+        DO i=ii_begin,ii_end
+          ij=(j-1)*iim+i
+          uc(ij,:,l)=ufi(ij,l)*elon_i(ij,:)+vfi(ij,l)*elat_i(ij,:)
+        ENDDO
+      ENDDO
+    ENDDO
+
+    out_i=ufi
+    
+    DO l=1,llm
+      DO j=jj_begin,jj_end
+        DO i=ii_begin,ii_end
+          ij=(j-1)*iim+i
+          ue(ij+u_right,l)=ue(ij+u_right,l)+sum( 0.5*(uc(ij,:,l) + uc(ij+t_right,:,l))*ep_e(ij+u_right,:) )
+          ue(ij+u_lup,l)=ue(ij+u_lup,l)+sum( 0.5*(uc(ij,:,l) + uc(ij+t_lup,:,l))*ep_e(ij+u_lup,:) )
+          ue(ij+u_ldown,l)=ue(ij+u_ldown,l)+sum( 0.5*(uc(ij,:,l) + uc(ij+t_ldown,:,l))*ep_e(ij+u_ldown,:) )
+        ENDDO
+      ENDDO
+    ENDDO
+   
+    CALL compute_temperature2theta_rhodz(ps,T,theta_rhodz,0)
+
+
+  END SUBROUTINE compute_physics
+   
+
+subroutine simple_physics (pcols, pver, dtime, lat, t, q, u, v, pmid, pint, pdel, rpdel, ps, precl, test)
+!----------------------------------------------------------------------- 
+! 
+! Purpose: Simple Physics Package
+!
+! Author: K. A. Reed (University of Michigan, kareed@umich.edu)
+!         version 5 
+!         July/8/2012
+!
+!  Change log:
+!  v2: removal of some NCAR CAM-specific 'use' associations
+!  v3: corrected precl(i) computation, the precipitation rate is now computed via a vertical integral, the previous single-level computation in v2 was a bug
+!  v3: corrected dtdt(i,1) computation, the term '-(i,1)' was missing the temperature variable: '-t(i,1)'
+!  v4: modified and enhanced parameter list to make the routine truly standalone, the number of columns and vertical levels have been added: pcols, pver
+!  v4: 'ncol' has been removed, 'pcols' is used instead
+!  v5: the sea surface temperature (SST) field Tsurf is now an array, the SST now depends on the latitude
+!  v5: addition of the latitude array 'lat' and the flag 'test' in the parameter list
+!      if test = 0: constant SST is used, correct setting for the tropical cyclone test case 5-1
+!      if test = 1: newly added latitude-dependent SST is used, correct setting for the moist baroclinic wave test with simple-physics (test 4-3)
+! 
+! Description: Includes large-scale precipitation, surface fluxes and
+!              boundary-leyer mixing. The processes are time-split
+!              in that order. A partially implicit formulation is
+!              used to foster numerical stability.
+!              The routine assumes that the model levels are ordered
+!              in a top-down approach, e.g. level 1 denotes the uppermost
+!              full model level.
+!
+!              This routine is based on an implementation which was
+!              developed for the NCAR Community Atmosphere Model (CAM).
+!              Adjustments for other models will be necessary.
+!
+!              The routine provides both updates of the state variables
+!              u, v, T, q (these are local copies of u,v,T,q within this physics
+!              routine) and also collects their time tendencies.
+!              The latter might be used to couple the physics and dynamics
+!              in a process-split way. For a time-split coupling, the final
+!              state should be given to the dynamical core for the next time step.
+! Test:      0 = Reed and Jablonowski (2011) tropical cyclone test case (test 5-1)
+!            1 = Moist baroclinic instability test (test 4-3)
+!
+!
+! Reference: Reed, K. A. and C. Jablonowski (2012), Idealized tropical cyclone 
+!            simulations of intermediate complexity: A test case for AGCMs, 
+!            J. Adv. Model. Earth Syst., Vol. 4, M04001, doi:10.1029/2011MS000099
+!-----------------------------------------------------------------------
+  ! use physics_types     , only: physics_dme_adjust   ! This is for CESM/CAM
+  ! use cam_diagnostics,    only: diag_phys_writeout   ! This is for CESM/CAM
+
+   implicit none
+
+   integer, parameter :: r8 = selected_real_kind(12)
+
+!
+! Input arguments - MODEL DEPENDENT
+!
+   integer, intent(in)  :: pcols        ! Set number of atmospheric columns       
+   integer, intent(in)  :: pver         ! Set number of model levels
+   real(r8), intent(in) :: dtime        ! Set model physics timestep
+   real(r8), intent(in) :: lat(pcols)   ! Latitude 
+   integer, intent(in)  :: test         ! Test number
+   
+!
+! Input/Output arguments 
+!
+!  pcols is the maximum number of vertical columns per 'chunk' of atmosphere
+!
+   real(r8), intent(inout) :: t(pcols,pver)      ! Temperature at full-model level (K)
+   real(r8), intent(inout) :: q(pcols,pver)      ! Specific Humidity at full-model level (kg/kg)
+   real(r8), intent(inout) :: u(pcols,pver)      ! Zonal wind at full-model level (m/s)
+   real(r8), intent(inout) :: v(pcols,pver)      ! Meridional wind at full-model level (m/s)
+   real(r8), intent(inout) :: pmid(pcols,pver)   ! Pressure is full-model level (Pa)
+   real(r8), intent(inout) :: pint(pcols,pver+1) ! Pressure at model interfaces (Pa)
+   real(r8), intent(inout) :: pdel(pcols,pver)   ! Layer thickness (Pa)
+   real(r8), intent(in) :: rpdel(pcols,pver)  ! Reciprocal of layer thickness (1/Pa)
+   real(r8), intent(inout) :: ps(pcols)          ! Surface Pressue (Pa)
+
+!
+! Output arguments 
+!
+   real(r8), intent(out) :: precl(pcols)         ! Precipitation rate (m_water / s)
+
+!
+!---------------------------Local workspace-----------------------------
+!
+
+! Integers for loops
+
+   integer  i,k                         ! Longitude, level indices
+
+! Physical Constants - Many of these may be model dependent
+
+   real(r8) gravit                      ! Gravity
+   real(r8) rair                        ! Gas constant for dry air
+   real(r8) cpair                       ! Specific heat of dry air 
+   real(r8) latvap                      ! Latent heat of vaporization
+   real(r8) rh2o                        ! Gas constant for water vapor
+   real(r8) epsilo                      ! Ratio of gas constant for dry air to that for vapor
+   real(r8) zvir                        ! Constant for virtual temp. calc. =(rh2o/rair) - 1
+   real(r8) a                           ! Reference Earth's Radius (m)
+   real(r8) omega                       ! Reference rotation rate of the Earth (s^-1)
+   real(r8) pi                          ! pi
+
+! Simple Physics Specific Constants 
+
+!++++++++                     
+   real(r8) Tsurf(pcols)                ! Sea Surface Temperature (constant for tropical cyclone)
+!++++++++                                 Tsurf needs to be dependent on latitude for the
+                                        ! moist baroclinic wave test 4-3 with simple-physics, adjust
+
+   real(r8) SST_tc                      ! Sea Surface Temperature for tropical cyclone test
+   real(r8) T0                          ! Control temp for calculation of qsat
+   real(r8) e0                          ! Saturation vapor pressure at T0 for calculation of qsat
+   real(r8) rhow                        ! Density of Liquid Water
+   real(r8) p0                          ! Constant for calculation of potential temperature
+   real(r8) Cd0                         ! Constant for calculating Cd from Smith and Vogl 2008
+   real(r8) Cd1                         ! Constant for calculating Cd from Smith and Vogl 2008
+   real(r8) Cm                          ! Constant for calculating Cd from Smith and Vogl 2008
+   real(r8) v20                         ! Threshold wind speed for calculating Cd from Smith and Vogl 2008
+   real(r8) C                           ! Drag coefficient for sensible heat and evaporation
+   real(r8) T00                         ! Horizontal mean T at surface for moist baro test
+   real(r8) u0                          ! Zonal wind constant for moist baro test
+   real(r8) latw                        ! halfwidth for  for baro test
+   real(r8) eta0                        ! Center of jets (hybrid) for baro test
+   real(r8) etav                        ! Auxiliary variable for baro test
+   real(r8) q0                          ! Maximum specific humidity for baro test
+
+! Physics Tendency Arrays
+  real(r8) dtdt(pcols,pver)             ! Temperature tendency 
+  real(r8) dqdt(pcols,pver)             ! Specific humidity tendency
+  real(r8) dudt(pcols,pver)             ! Zonal wind tendency
+  real(r8) dvdt(pcols,pver)             ! Meridional wind tendency
+
+! Temporary variables for tendency calculations
+
+   real(r8) tmp                         ! Temporary
+   real(r8) qsat                        ! Saturation vapor pressure
+   real(r8) qsats                       ! Saturation vapor pressure of SST
+
+! Variables for Boundary Layer Calculation
+
+   real(r8) wind(pcols)                 ! Magnitude of Wind
+   real(r8) Cd(pcols)                   ! Drag coefficient for momentum
+   real(r8) Km(pcols,pver+1)            ! Eddy diffusivity for boundary layer calculations 
+   real(r8) Ke(pcols,pver+1)            ! Eddy diffusivity for boundary layer calculations
+   real(r8) rho                         ! Density at lower/upper interface
+   real(r8) za(pcols)                   ! Heights at midpoints of first model level
+   real(r8) dlnpint                     ! Used for calculation of heights
+   real(r8) pbltop                      ! Top of boundary layer
+   real(r8) pblconst                    ! Constant for the calculation of the decay of diffusivity 
+   real(r8) CA(pcols,pver)              ! Matrix Coefficents for PBL Scheme 
+   real(r8) CC(pcols,pver)              ! Matrix Coefficents for PBL Scheme 
+   real(r8) CE(pcols,pver+1)            ! Matrix Coefficents for PBL Scheme
+   real(r8) CAm(pcols,pver)             ! Matrix Coefficents for PBL Scheme
+   real(r8) CCm(pcols,pver)             ! Matrix Coefficents for PBL Scheme
+   real(r8) CEm(pcols,pver+1)           ! Matrix Coefficents for PBL Scheme
+   real(r8) CFu(pcols,pver+1)           ! Matrix Coefficents for PBL Scheme
+   real(r8) CFv(pcols,pver+1)           ! Matrix Coefficents for PBL Scheme
+   real(r8) CFt(pcols,pver+1)           ! Matrix Coefficents for PBL Scheme
+   real(r8) CFq(pcols,pver+1)           ! Matrix Coefficents for PBL Scheme
+
+
+! Variable for Dry Mass Adjustment, this dry air adjustment is necessary to
+! conserve the mass of the dry air
+
+   real(r8) qini(pcols,pver)            ! Initial specific humidity
+
+!===============================================================================
+!
+! Physical Constants - MAY BE MODEL DEPENDENT 
+!
+!===============================================================================
+   gravit = 9.80616_r8                   ! Gravity (9.80616 m/s^2)
+   rair   = 287.0_r8                     ! Gas constant for dry air: 287 J/(kg K)
+   cpair  = 1.0045e3_r8                  ! Specific heat of dry air: here we use 1004.5 J/(kg K)
+   latvap = 2.5e6_r8                     ! Latent heat of vaporization (J/kg)
+   rh2o   = 461.5_r8                     ! Gas constant for water vapor: 461.5 J/(kg K)
+   epsilo = rair/rh2o                    ! Ratio of gas constant for dry air to that for vapor
+   zvir   = (rh2o/rair) - 1._r8          ! Constant for virtual temp. calc. =(rh2o/rair) - 1 is approx. 0.608
+   a      = 6371220.0_r8                 ! Reference Earth's Radius (m)
+   omega  = 7.29212d-5                   ! Reference rotation rate of the Earth (s^-1)
+   pi     = 4._r8*atan(1._r8)            ! pi
+
+!===============================================================================
+!
+! Local Constants for Simple Physics
+!
+!===============================================================================
+      C        = 0.0011_r8      ! From Smith and Vogl 2008
+      SST_tc   = 302.15_r8      ! Constant Value for SST for tropical cyclone test
+      T0       = 273.16_r8      ! control temp for calculation of qsat
+      e0       = 610.78_r8      ! saturation vapor pressure at T0 for calculation of qsat
+      rhow     = 1000.0_r8      ! Density of Liquid Water 
+      Cd0      = 0.0007_r8      ! Constant for Cd calc. Smith and Vogl 2008
+      Cd1      = 0.000065_r8    ! Constant for Cd calc. Smith and Vogl 2008
+      Cm       = 0.002_r8       ! Constant for Cd calc. Smith and Vogl 2008
+      v20      = 20.0_r8        ! Threshold wind speed for calculating Cd from Smith and Vogl 2008
+      p0       = 100000.0_r8    ! Constant for potential temp calculation
+      pbltop   = 85000._r8      ! Top of boundary layer
+      pblconst = 10000._r8      ! Constant for the calculation of the decay of diffusivity
+      T00      = 288.0_r8         ! Horizontal mean T at surface for moist baro test
+      u0       = 35.0_r8          ! Zonal wind constant for moist baro test
+      latw     = 2.0_r8*pi/9.0_r8 ! Halfwidth for  for baro test
+      eta0     = 0.252_r8         ! Center of jets (hybrid) for baro test
+      etav     = (1._r8-eta0)*0.5_r8*pi ! Auxiliary variable for baro test
+      q0       = 0.021            ! Maximum specific humidity for baro test
+
+!===============================================================================
+!
+! Definition of local arrays
+!
+!===============================================================================
+!
+! Calculate hydrostatic height za of the lowest model level
+!
+     do i=1,pcols 
+        dlnpint = log(ps(i)) - log(pint(i,pver))                 ! ps(i) is identical to pint(i,pver+1), note: this is the correct sign (corrects typo in JAMES paper)
+        za(i) = rair/gravit*t(i,pver)*(1._r8+zvir*q(i,pver))*0.5_r8*dlnpint
+     end do
+!
+! Set Initial Specific Humidity - For dry mass adjustment at the end
+!
+     qini(:pcols,:pver) = q(:pcols,:pver)
+!--------------------------------------------------------------
+! Set Sea Surface Temperature (constant for tropical cyclone)
+! Tsurf needs to be dependent on latitude for the
+! moist baroclinic wave test 4-3 with simple-physics 
+!--------------------------------------------------------------
+     if (test .eq. 1) then     ! moist baroclinic wave with simple-physics
+        do i=1,pcols
+           Tsurf(i) = (T00 + pi*u0/rair * 1.5_r8 * sin(etav) * (cos(etav))**0.5_r8 *                 &
+                     ((-2._r8*(sin(lat(i)))**6 * ((cos(lat(i)))**2 + 1._r8/3._r8) + 10._r8/63._r8)* &
+                     u0 * (cos(etav))**1.5_r8  +                                                    &
+                     (8._r8/5._r8*(cos(lat(i)))**3 * ((sin(lat(i)))**2 + 2._r8/3._r8) - pi/4._r8)*a*omega*0.5_r8 ))/ &
+                     (1._r8+zvir*q0*exp(-(lat(i)/latw)**4))
+
+        end do
+     else
+        do i=1,pcols          ! constant SST for the tropical cyclone test case
+           Tsurf(i) = SST_tc
+        end do
+     end if
+
+!===============================================================================
+!
+! Set initial physics time tendencies and precipitation field to zero
+!
+!===============================================================================
+     dtdt(:pcols,:pver)  = 0._r8            ! initialize temperature tendency with zero
+     dqdt(:pcols,:pver)  = 0._r8            ! initialize specific humidity tendency with zero
+     dudt(:pcols,:pver)  = 0._r8            ! initialize zonal wind tendency with zero
+     dvdt(:pcols,:pver)  = 0._r8            ! initialize meridional wind tendency with zero
+     precl(:pcols) = 0._r8                  ! initialize precipitation rate with zero
+
+!===============================================================================
+!
+! Large-Scale Condensation and Precipitation Rate
+!
+!===============================================================================
+!
+! Calculate Tendencies
+!
+      do k=1,pver
+         do i=1,pcols
+            qsat = epsilo*e0/pmid(i,k)*exp(-latvap/rh2o*((1._r8/t(i,k))-1._r8/T0))  ! saturation specific humidity
+            if (q(i,k) > qsat) then                                                 ! saturated?
+               tmp  = 1._r8/dtime*(q(i,k)-qsat)/(1._r8+(latvap/cpair)*(epsilo*latvap*qsat/(rair*t(i,k)**2)))
+               dtdt(i,k) = dtdt(i,k)+latvap/cpair*tmp
+               dqdt(i,k) = dqdt(i,k)-tmp
+               precl(i) = precl(i) + tmp*pdel(i,k)/(gravit*rhow)                    ! precipitation rate, computed via a vertical integral
+                                                                                    ! corrected in version 1.3
+            end if
+         end do
+      end do
+!
+! Update moisture and temperature fields from Large-Scale Precipitation Scheme
+!
+      do k=1,pver
+         do i=1,pcols
+            t(i,k) =  t(i,k) + dtdt(i,k)*dtime    ! update the state variables T and q
+            q(i,k) =  q(i,k) + dqdt(i,k)*dtime
+         end do
+      end do
+
+!===============================================================================
+! Send variables to history file - THIS PROCESS WILL BE MODEL SPECIFIC
+!
+! note: The variables, as done in many GCMs, are written to the history file
+!       after the moist physics process.  This ensures that the moisture fields
+!       are somewhat in equilibrium.
+!===============================================================================
+  !  call diag_phys_writeout(state)   ! This is for CESM/CAM
+
+!===============================================================================
+!
+! Turbulent mixing coefficients for the PBL mixing of horizontal momentum,
+! sensible heat and latent heat
+!
+! We are using Simplified Ekman theory to compute the diffusion coefficients 
+! Kx for the boundary-layer mixing. The Kx values are calculated at each time step
+! and in each column.
+!
+!===============================================================================
+!
+! Compute magnitude of the wind and drag coeffcients for turbulence scheme:
+! they depend on the conditions at the lowest model level and stay constant
+! up to the 850 hPa level. Above this level the coefficients are decreased
+! and tapered to zero. At the 700 hPa level the strength of the K coefficients
+! is about 10% of the maximum strength. 
+!
+     do i=1,pcols
+        wind(i) = sqrt(u(i,pver)**2+v(i,pver)**2)    ! wind magnitude at the lowest level
+     end do
+     do i=1,pcols
+        Ke(i,pver+1) = C*wind(i)*za(i)
+        if( wind(i) .lt. v20) then
+           Cd(i) = Cd0+Cd1*wind(i) 
+           Km(i,pver+1) = Cd(i)*wind(i)*za(i)
+        else
+           Cd(i) = Cm
+           Km(i,pver+1) = Cm*wind(i)*za(i)
+        endif
+     end do
+
+      do k=1,pver
+         do i=1,pcols
+            if( pint(i,k) .ge. pbltop) then
+               Km(i,k) = Km(i,pver+1)                 ! constant Km below 850 hPa level
+               Ke(i,k) = Ke(i,pver+1)                 ! constant Ke below 850 hPa level
+            else
+               Km(i,k) = Km(i,pver+1)*exp(-(pbltop-pint(i,k))**2/(pblconst)**2)  ! Km tapered to 0
+               Ke(i,k) = Ke(i,pver+1)*exp(-(pbltop-pint(i,k))**2/(pblconst)**2)  ! Ke tapered to 0
+            end if 
+         end do
+      end do     
+
+
+!===============================================================================
+! Update the state variables u, v, t, q with the surface fluxes at the
+! lowest model level, this is done with an implicit approach
+! see Reed and Jablonowski (JAMES, 2012)
+!
+! Sea Surface Temperature Tsurf is constant for tropical cyclone test 5-1
+! Tsurf needs to be dependent on latitude for the
+! moist baroclinic wave test 4-3 with simple-physics, adjust
+!===============================================================================
+
+     do i=1,pcols
+        qsats = epsilo*e0/ps(i)*exp(-latvap/rh2o*((1._r8/Tsurf(i))-1._r8/T0))  ! saturation specific humidity at the surface
+        dudt(i,pver) = dudt(i,pver) + (u(i,pver) &
+                            /(1._r8+Cd(i)*wind(i)*dtime/za(i))-u(i,pver))/dtime
+        dvdt(i,pver) = dvdt(i,pver) + (v(i,pver) &
+                            /(1._r8+Cd(i)*wind(i)*dtime/za(i))-v(i,pver))/dtime
+        u(i,pver)   = u(i,pver)/(1._r8+Cd(i)*wind(i)*dtime/za(i))
+        v(i,pver)   = v(i,pver)/(1._r8+Cd(i)*wind(i)*dtime/za(i))
+        dtdt(i,pver) = dtdt(i,pver) +((t(i,pver)+C*wind(i)*Tsurf(i)*dtime/za(i)) &
+                            /(1._r8+C*wind(i)*dtime/za(i))-t(i,pver))/dtime 
+        t(i,pver)   = (t(i,pver)+C*wind(i)*Tsurf(i)*dtime/za(i)) &
+                            /(1._r8+C*wind(i)*dtime/za(i))  
+        dqdt(i,pver) = dqdt(i,pver) +((q(i,pver)+C*wind(i)*qsats*dtime/za(i)) &
+                            /(1._r8+C*wind(i)*dtime/za(i))-q(i,pver))/dtime
+        q(i,pver) = (q(i,pver)+C*wind(i)*qsats*dtime/za(i))/(1._r8+C*wind(i)*dtime/za(i))
+     end do
+!===============================================================================
+
+
+!===============================================================================
+! Boundary layer mixing, see Reed and Jablonowski (JAMES, 2012)
+!===============================================================================
+! Calculate Diagonal Variables for Implicit PBL Scheme
+!
+      do k=1,pver-1
+         do i=1,pcols
+            rho = (pint(i,k+1)/(rair*(t(i,k+1)+t(i,k))/2.0_r8))
+            CAm(i,k)   = rpdel(i,k)*dtime*gravit*gravit*Km(i,k+1)*rho*rho   &
+                         /(pmid(i,k+1)-pmid(i,k))    
+            CCm(i,k+1) = rpdel(i,k+1)*dtime*gravit*gravit*Km(i,k+1)*rho*rho &
+                         /(pmid(i,k+1)-pmid(i,k))
+            CA(i,k)    = rpdel(i,k)*dtime*gravit*gravit*Ke(i,k+1)*rho*rho   &
+                         /(pmid(i,k+1)-pmid(i,k))
+            CC(i,k+1)  = rpdel(i,k+1)*dtime*gravit*gravit*Ke(i,k+1)*rho*rho &
+                         /(pmid(i,k+1)-pmid(i,k))
+         end do
+      end do
+      do i=1,pcols
+         CAm(i,pver) = 0._r8
+         CCm(i,1) = 0._r8
+         CEm(i,pver+1) = 0._r8
+         CA(i,pver) = 0._r8
+         CC(i,1) = 0._r8
+         CE(i,pver+1) = 0._r8
+         CFu(i,pver+1) = 0._r8
+         CFv(i,pver+1) = 0._r8
+         CFt(i,pver+1) = 0._r8
+         CFq(i,pver+1) = 0._r8 
+      end do
+      do i=1,pcols
+         do k=pver,1,-1
+            CE(i,k)  = CC(i,k)/(1._r8+CA(i,k)+CC(i,k)-CA(i,k)*CE(i,k+1)) 
+            CEm(i,k) = CCm(i,k)/(1._r8+CAm(i,k)+CCm(i,k)-CAm(i,k)*CEm(i,k+1))
+            CFu(i,k) = (u(i,k)+CAm(i,k)*CFu(i,k+1)) &
+                       /(1._r8+CAm(i,k)+CCm(i,k)-CAm(i,k)*CEm(i,k+1))
+            CFv(i,k) = (v(i,k)+CAm(i,k)*CFv(i,k+1)) &
+                       /(1._r8+CAm(i,k)+CCm(i,k)-CAm(i,k)*CEm(i,k+1))
+            CFt(i,k) = ((p0/pmid(i,k))**(rair/cpair)*t(i,k)+CA(i,k)*CFt(i,k+1)) &
+                       /(1._r8+CA(i,k)+CC(i,k)-CA(i,k)*CE(i,k+1)) 
+            CFq(i,k) = (q(i,k)+CA(i,k)*CFq(i,k+1)) &
+                       /(1._r8+CA(i,k)+CC(i,k)-CA(i,k)*CE(i,k+1))
+        end do
+      end do
+
+!
+! Calculate the updated temperature, specific humidity and horizontal wind
+!
+! First we need to calculate the updates at the top model level
+!
+      do i=1,pcols
+            dudt(i,1)  = dudt(i,1)+(CFu(i,1)-u(i,1))/dtime
+            dvdt(i,1)  = dvdt(i,1)+(CFv(i,1)-v(i,1))/dtime
+            u(i,1)    = CFu(i,1)
+            v(i,1)    = CFv(i,1)
+            dtdt(i,1)  = dtdt(i,1)+(CFt(i,1)*(pmid(i,1)/p0)**(rair/cpair)-t(i,1))/dtime  ! corrected in version 1.3
+            t(i,1)    = CFt(i,1)*(pmid(i,1)/p0)**(rair/cpair)
+            dqdt(i,1)  = dqdt(i,1)+(CFq(i,1)-q(i,1))/dtime
+            q(i,1)  = CFq(i,1)
+      end do
+!
+! Loop over the remaining level
+!
+      do i=1,pcols
+         do k=2,pver
+            dudt(i,k)  = dudt(i,k)+(CEm(i,k)*u(i,k-1)+CFu(i,k)-u(i,k))/dtime
+            dvdt(i,k)  = dvdt(i,k)+(CEm(i,k)*v(i,k-1)+CFv(i,k)-v(i,k))/dtime
+            u(i,k)    = CEm(i,k)*u(i,k-1)+CFu(i,k) 
+            v(i,k)    = CEm(i,k)*v(i,k-1)+CFv(i,k)
+            dtdt(i,k)  = dtdt(i,k)+((CE(i,k)*t(i,k-1) &
+                              *(p0/pmid(i,k-1))**(rair/cpair)+CFt(i,k)) &
+                              *(pmid(i,k)/p0)**(rair/cpair)-t(i,k))/dtime 
+            t(i,k)    = (CE(i,k)*t(i,k-1)*(p0/pmid(i,k-1))**(rair/cpair)+CFt(i,k)) &
+                              *(pmid(i,k)/p0)**(rair/cpair)
+            dqdt(i,k)  = dqdt(i,k)+(CE(i,k)*q(i,k-1)+CFq(i,k)-q(i,k))/dtime
+            q(i,k)  = CE(i,k)*q(i,k-1)+CFq(i,k)
+         end do
+      end do
+
+!===============================================================================
+!
+! Dry Mass Adjustment - THIS PROCESS WILL BE MODEL SPECIFIC 
+!
+! note: Care needs to be taken to ensure that the model conserves the total
+!       dry air mass. Add your own routine here.
+!===============================================================================
+  !  call physics_dme_adjust(state, tend, qini, dtime)   ! This is for CESM/CAM
+
+   return
+end subroutine simple_physics 
+
+
+
+
+END MODULE physics_dcmip_mod
Index: /codes/icosagcm/trunk/src/timeloop_gcm.f90
===================================================================
--- /codes/icosagcm/trunk/src/timeloop_gcm.f90	(revision 80)
+++ /codes/icosagcm/trunk/src/timeloop_gcm.f90	(revision 81)
@@ -10,4 +10,5 @@
   USE guided_mod
   USE advect_tracer_mod
+  USE physics_mod
   
   IMPLICIT NONE
@@ -98,11 +99,15 @@
   CALL init_guided(dt)
   CALL init_advect_tracer(dt)
+  CALL init_physics(dt)
   
   CALL etat0(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
+  CALL transfert_request(f_q,req_i1) 
 
   DO it=0,itaumax
+
     PRINT *,"It No :",It,"   t :",dt*It
     IF (mod(it,itau_out)==0 ) THEN
       CALL writefield("q",f_q)
+      CALL update_time_counter(dt*it)
     ENDIF
     
@@ -110,5 +115,5 @@
     CALL caldyn(it,f_phis,f_ps,f_theta_rhodz,f_u, f_dps, f_dtheta_rhodz, f_du)
     CALL advect_tracer(f_ps,f_u,f_q)
-
+    CALL physics(f_phis, f_ps, f_theta_rhodz, f_u, f_q)
     
     SELECT CASE (TRIM(scheme))
@@ -132,5 +137,5 @@
         
     END SELECT
-
+    
  ENDDO
   
