Changeset 196 for codes/icosagcm/trunk/src/physics.f90

Timestamp:

07/06/14 09:10:13 (10 years ago)

Author:

dubos

Message:

First draft of generic dynamics-physics interface - works with DCMIP5.1

File:

• codes/icosagcm/trunk/src/physics.f90 (modified) (4 diffs)

codes/icosagcm/trunk/src/physics.f90

@@ -1 +1 @@
 MODULE physics_mod
+
+  USE field_mod
+
+  PRIVATE
+
+  INTEGER, PARAMETER :: phys_none=0, phys_HS94=1, phys_DCMIP=2
+
+  INTEGER :: nb_extra_physics_2D, nb_extra_physics_3D, phys_type
+  TYPE(t_field),POINTER :: f_extra_physics_2D(:), f_extra_physics_3D(:)
 
   CHARACTER(LEN=255) :: physics_type="automatic"
 !$OMP THREADPRIVATE(physics_type)
 
+  PUBLIC :: physics, init_physics
 CONTAINS
 
@@ -10 +20 @@
     USE etat0_mod
     USE icosa
-    USE physics_dcmip_mod,init_physics_dcmip=>init_physics
-!    USE physics_dry_mod
+    USE physics_interface_mod
+    USE physics_dcmip_mod, ONLY : init_physics_dcmip=>init_physics
     IMPLICIT NONE
 
@@ -23 +33 @@
        SELECT CASE(TRIM(etat0_type))
        CASE('held_suarez')
+          phys_type = phys_HS94
        CASE DEFAULT
           IF(is_mpi_root) PRINT*,"NO PHYSICAL PACKAGE USED"
+          phys_type = phys_none
        END SELECT
-
 
     CASE ('dcmip')
        CALL init_physics_dcmip
-
-    CASE ('dry')
-!       CALL init_physics_dry
-
+       nb_extra_physics_2D=1 ! precl
+       nb_extra_physics_3D=0
+       phys_type = phys_DCMIP
     CASE DEFAULT
-       PRINT*, 'init_physics : Bad selector for variable physics <',TRIM(physics_type), &
-            '> options are <automatic>, <dcmip>, <dry>'
+       IF(is_mpi_root) PRINT*, 'init_physics : Bad selector for variable physics <',&
+            TRIM(physics_type), '> options are <automatic>, <dcmip>, <dry>'
        STOP
     END SELECT
 
+    IF(is_mpi_root) THEN
+       PRINT *, 'phys_type = ',phys_type
+       PRINT *, 'nb_extra_physics_2D = ', nb_extra_physics_2D
+       PRINT *, 'nb_extra_physics_3D = ', nb_extra_physics_3D
+    END IF
+    IF(nb_extra_physics_2D>0) CALL allocate_field(f_extra_physics_2D,field_t,type_real,nb_extra_physics_2D)
+    IF(nb_extra_physics_3D>0) CALL allocate_field(f_extra_physics_3D,field_t,type_real,llm,nb_extra_physics_3D)
+
   END SUBROUTINE init_physics
 
   SUBROUTINE physics(it,jD_cur,jH_cur,f_phis, f_ps, f_theta_rhodz, f_ue, f_q)
     USE icosa
-!    USE physics_dry_mod
-    USE physics_dcmip_mod, physics_dcmip=>physics
-    USE etat0_mod
+    USE physics_interface_mod
+ !   USE etat0_mod
     USE etat0_heldsz_mod
     IMPLICIT NONE
@@ -56 +73 @@
     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(:,:,:)
+
     LOGICAL:: firstcall,lastcall
-
-    SELECT CASE(TRIM(physics_type))
-    CASE ('automatic')
-
-       SELECT CASE(TRIM(etat0_type))
-       CASE('held_suarez')
-          !     CALL transfert_request(f_ps,req_i1)
-          !     CALL transfert_request(f_theta_rhodz,req_i1)
-          !     CALL transfert_request(f_ue,req_e1_vect)
-          CALL held_suarez(f_ps,f_theta_rhodz,f_ue) 
-       CASE DEFAULT
-!          PRINT*,"NO PHYSICAL PACAKAGE USED" ! FIXME MPI
-       END SELECT
-
-    CASE ('dcmip')
-       CALL physics_dcmip(it, f_phis, f_ps, f_theta_rhodz, f_ue, f_q)
-
-    CASE ('dry')
-!       CALL physics_dry(it,jD_cur,jH_cur,f_phis, f_ps, f_theta_rhodz, f_ue, f_q)
-
+    INTEGER :: ind
+    TYPE(physics_inout) :: args
+
+    SELECT CASE(phys_type)
+    CASE (phys_none)
+       ! No physics, do nothing
+    CASE(phys_HS94)
+       CALL held_suarez(f_ps,f_theta_rhodz,f_ue) 
     CASE DEFAULT
-       PRINT*, 'Bad selector for variable physics <',TRIM(physics_type), &
-            '> options are <automatic>, <dcmip>, <dry>'
-       STOP
+       CALL transfert_request(f_ps,req_i1)
+       CALL transfert_request(f_theta_rhodz,req_i1)
+       CALL transfert_request(f_ue,req_e1_vect)
+       CALL transfert_request(f_q,req_i1)
+       
+       args%dt_phys = dt
+       DO ind=1,ndomain
+          IF (.NOT. assigned_domain(ind)) CYCLE
+          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)
+
+          IF(nb_extra_physics_2D>0) args%extra_2D=f_extra_physics_2D(ind)
+          IF(nb_extra_physics_3D>0) args%extra_3D=f_extra_physics_3D(ind)
+          CALL physics_column(args, phis, ps, theta_rhodz, ue, q)
+       ENDDO
+
+       IF (mod(it,itau_out)==0 ) THEN
+          IF(nb_extra_physics_2D>0) CALL writefield("extra_physics_2D",f_extra_physics_2D)
+          IF(nb_extra_physics_3D>0) CALL writefield("extra_physics_3D",f_extra_physics_3D)
+       ENDIF
     END SELECT
 
   END SUBROUTINE physics
 
+  SUBROUTINE physics_column(args, phis, ps, theta_rhodz, ue, q)
+    USE icosa
+    USE wind_mod
+    USE pression_mod
+    USE theta2theta_rhodz_mod
+    USE physics_interface_mod
+    USE physics_dcmip_mod
+    IMPLICIT NONE
+    TYPE(physics_inout) :: args    
+    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), TARGET :: q(iim*jjm,llm,nqtot)
+    ! local arrays
+    REAL(rstd), TARGET :: lat(iim*jjm)
+    REAL(rstd), TARGET :: lon(iim*jjm)
+    REAL(rstd), TARGET :: p(iim*jjm,llm+1)
+    REAL(rstd), TARGET :: Temp(iim*jjm,llm)
+    REAL(rstd), TARGET :: ulon(iim*jjm,llm)
+    REAL(rstd), TARGET :: ulat(iim*jjm,llm)
+    REAL(rstd), TARGET :: dTemp(iim*jjm,llm)
+    REAL(rstd), TARGET :: dulon(iim*jjm,llm)
+    REAL(rstd), TARGET :: dulat(iim*jjm,llm)
+    REAL(rstd), TARGET :: dq(iim*jjm,llm,nqtot)
+    REAL(rstd) :: uc(iim*jjm,3,llm)  ! 3D velocity at cell centers
+
+    INTEGER :: i,j,ij,l
+    REAL(rstd) :: due, dt2
+
+    DO j=jj_begin,jj_end
+      DO i=ii_begin,ii_end
+        ij=(j-1)*iim+i
+        CALL xyz2lonlat(xyz_i(ij,:),lon(ij),lat(ij)) 
+      ENDDO
+    ENDDO
+
+    ! Reconstruct wind vector at hexagons
+    CALL compute_pression(ps,p,0)
+    CALL compute_theta_rhodz2temperature(ps,theta_rhodz,Temp,0)
+    CALL compute_wind_centered(ue,uc)
+    CALL compute_wind_centered_lonlat_compound(uc, ulon, ulat)
+    args%ngrid = iim*jjm
+    args%lon => lon
+    args%lat => lat
+    args%p => p
+    args%Temp => Temp
+    args%ulon => ulon
+    args%ulat => ulat
+    args%q => q
+    args%dTemp => dTemp
+    args%dulon => dulon
+    args%dulat => dulat
+    args%dq => dq
+
+    SELECT CASE(phys_type)
+    CASE (phys_DCMIP)
+       CALL compute_phys_wrap(args)
+    END SELECT
+
+    q = q + args%dt_phys * dq
+    Temp = Temp + args%dt_phys * dTemp
+    CALL compute_temperature2theta_rhodz(ps,Temp,theta_rhodz,0)
+    
+    ! Reconstruct wind tendencies at edges and add
+    CALL compute_wind_centered_from_lonlat_compound(dulon,dulat,uc)
+    dt2=.5*args%dt_phys
+    DO l=1,llm
+      DO j=jj_begin,jj_end
+        DO i=ii_begin,ii_end
+          ij=(j-1)*iim+i
+          due = sum((uc(ij,:,l) + uc(ij+t_right,:,l))*ep_e(ij+u_right,:) )
+          ue(ij+u_right,l) = ue(ij+u_right,l) + dt2*due
+
+          due = sum( (uc(ij,:,l) + uc(ij+t_lup,:,l))*ep_e(ij+u_lup,:) )
+          ue(ij+u_lup,l)=ue(ij+u_lup,l) + dt2*due
+
+          due = sum( (uc(ij,:,l) + uc(ij+t_ldown,:,l))*ep_e(ij+u_ldown,:) )
+          ue(ij+u_ldown,l)=ue(ij+u_ldown,l) + dt2*due
+        ENDDO
+      ENDDO
+    ENDDO
+
+  END SUBROUTINE physics_column
+
 END MODULE physics_mod