Changeset 213 for codes/icosagcm/trunk


Ignore:
Timestamp:
07/15/14 18:15:39 (10 years ago)
Author:
dubos
Message:

New dyn/phys interface - halo points not passed to physics any more (cleanup follows)

Location:
codes/icosagcm/trunk/src
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • codes/icosagcm/trunk/src/physics.f90

    r196 r213  
    77  INTEGER, PARAMETER :: phys_none=0, phys_HS94=1, phys_DCMIP=2 
    88 
    9   INTEGER :: nb_extra_physics_2D, nb_extra_physics_3D, phys_type 
     9  INTEGER :: phys_type 
    1010  TYPE(t_field),POINTER :: f_extra_physics_2D(:), f_extra_physics_3D(:) 
     11  TYPE(t_field),POINTER :: f_dulon(:), f_dulat(:) 
    1112 
    1213  CHARACTER(LEN=255) :: physics_type="automatic" 
    … …  
    1415 
    1516  PUBLIC :: physics, init_physics 
     17 
    1618CONTAINS 
    1719 
    … …  
    2123    USE icosa 
    2224    USE physics_interface_mod 
    23     USE physics_dcmip_mod, ONLY : init_physics_dcmip=>init_physics 
    24     IMPLICIT NONE 
    25  
     25    USE physics_dcmip_mod, ONLY : & 
     26         init_physics_dcmip=>init_physics, init_physics_dcmip_new=>init_physics_new 
     27    IMPLICIT NONE 
     28 
     29    CALL allocate_field(f_dulon,field_t,type_real,llm, name='dulon') 
     30    CALL allocate_field(f_dulat,field_t,type_real,llm, name='dulat') 
     31    CALL init_pack_before ! Compute physics_inout%ngrid and offsets used by pack/unpack 
     32     
    2633    physics_type='automatic' 
    2734    CALL getin("physics",physics_type) 
    … …  
    4047 
    4148    CASE ('dcmip') 
    42        CALL init_physics_dcmip 
    43        nb_extra_physics_2D=1 ! precl 
    44        nb_extra_physics_3D=0 
     49       CALL init_physics_dcmip_new 
     50!       CALL init_physics_dcmip 
    4551       phys_type = phys_DCMIP 
    4652    CASE DEFAULT 
    … …  
    5561       PRINT *, 'nb_extra_physics_3D = ', nb_extra_physics_3D 
    5662    END IF 
    57     IF(nb_extra_physics_2D>0) CALL allocate_field(f_extra_physics_2D,field_t,type_real,nb_extra_physics_2D) 
    58     IF(nb_extra_physics_3D>0) CALL allocate_field(f_extra_physics_3D,field_t,type_real,llm,nb_extra_physics_3D) 
    59  
     63 
     64    IF(.FALSE.) THEN ! draft interface 
     65       IF(nb_extra_physics_2D>0) CALL allocate_field(f_extra_physics_2D,field_t,type_real,nb_extra_physics_2D) 
     66       IF(nb_extra_physics_3D>0) CALL allocate_field(f_extra_physics_3D,field_t,type_real,llm,nb_extra_physics_3D) 
     67    ELSE 
     68       CALL init_pack_after ! Defines Ai, lon, lat in physics_inout 
     69    END IF 
    6070  END SUBROUTINE init_physics 
    6171 
    … …  
    6373    USE icosa 
    6474    USE physics_interface_mod 
    65  !   USE etat0_mod 
     75    USE physics_dcmip_mod, ONLY : write_physics_dcmip => write_physics 
    6676    USE etat0_heldsz_mod 
    6777    IMPLICIT NONE 
    … …  
    8191    LOGICAL:: firstcall,lastcall 
    8292    INTEGER :: ind 
    83     TYPE(physics_inout) :: args 
     93    TYPE(t_physics_inout) :: args 
     94 
     95    IF(MOD(it+1,itau_physics)==0) THEN 
     96     
     97       SELECT CASE(phys_type) 
     98       CASE (phys_none) 
     99          ! No physics, do nothing 
     100       CASE(phys_HS94) 
     101          CALL held_suarez(f_ps,f_theta_rhodz,f_ue)  
     102       CASE DEFAULT 
     103          IF(.FALSE.) THEN ! draft interface 
     104             args%dt_phys = dt*itau_physics 
     105             DO ind=1,ndomain 
     106                IF (.NOT. assigned_domain(ind)) CYCLE 
     107                CALL swap_dimensions(ind) 
     108                CALL swap_geometry(ind) 
     109                 
     110                phis=f_phis(ind) 
     111                ps=f_ps(ind) 
     112                theta_rhodz=f_theta_rhodz(ind) 
     113                ue=f_ue(ind) 
     114                q=f_q(ind) 
     115                 
     116                IF(nb_extra_physics_2D>0) args%extra_2D=f_extra_physics_2D(ind) 
     117                IF(nb_extra_physics_3D>0) args%extra_3D=f_extra_physics_3D(ind) 
     118                CALL physics_column(args, phis, ps, theta_rhodz, ue, q) 
     119             ENDDO 
     120              
     121             IF (mod(it,itau_out)==0 ) THEN 
     122                IF(nb_extra_physics_2D>0) CALL writefield("extra_physics_2D",f_extra_physics_2D) 
     123                IF(nb_extra_physics_3D>0) CALL writefield("extra_physics_3D",f_extra_physics_3D) 
     124             ENDIF 
     125          ELSE ! new interface 
     126             physics_inout%dt_phys = dt*itau_physics 
     127             CALL physics_column_new(it, f_phis, f_ps, f_theta_rhodz, f_ue, f_q) 
     128          END IF 
     129       END SELECT 
     130 
     131       CALL transfert_request(f_theta_rhodz,req_i0) 
     132       CALL transfert_request(f_ue,req_e0_vect) 
     133       CALL transfert_request(f_q,req_i0) 
     134    END IF 
     135 
     136    IF (mod(it,itau_out)==0 ) THEN 
     137       SELECT CASE(phys_type) 
     138       CASE (phys_DCMIP) 
     139          CALL write_physics_dcmip 
     140       END SELECT 
     141    END IF 
     142     
     143  END SUBROUTINE physics 
     144 
     145!--------------------------------- New interface -------------------------------------- 
     146 
     147  SUBROUTINE physics_column_new(it, f_phis, f_ps, f_theta_rhodz, f_ue, f_q) 
     148    USE icosa 
     149    USE physics_interface_mod 
     150    USE physics_dcmip_mod, ONLY : full_physics_dcmip => full_physics 
     151    IMPLICIT NONE 
     152    TYPE(t_field),POINTER :: f_phis(:) 
     153    TYPE(t_field),POINTER :: f_ps(:) 
     154    TYPE(t_field),POINTER :: f_theta_rhodz(:) 
     155    TYPE(t_field),POINTER :: f_ue(:) 
     156    TYPE(t_field),POINTER :: f_q(:) 
     157    REAL(rstd),POINTER :: phis(:) 
     158    REAL(rstd),POINTER :: ps(:) 
     159    REAL(rstd),POINTER :: theta_rhodz(:,:) 
     160    REAL(rstd),POINTER :: ue(:,:) 
     161    REAL(rstd),POINTER :: dulon(:,:) 
     162    REAL(rstd),POINTER :: dulat(:,:) 
     163    REAL(rstd),POINTER :: q(:,:,:) 
     164    INTEGER :: it, ind 
     165 
     166    DO ind=1,ndomain 
     167       IF (.NOT. assigned_domain(ind)) CYCLE 
     168       CALL swap_dimensions(ind) 
     169       CALL swap_geometry(ind) 
     170       phis=f_phis(ind) 
     171       ps=f_ps(ind) 
     172       theta_rhodz=f_theta_rhodz(ind) 
     173       ue=f_ue(ind) 
     174       q=f_q(ind) 
     175       CALL pack_physics(pack_info(ind), phis, ps, theta_rhodz, ue, q) 
     176    END DO 
    84177 
    85178    SELECT CASE(phys_type) 
    86     CASE (phys_none) 
    87        ! No physics, do nothing 
    88     CASE(phys_HS94) 
    89        CALL held_suarez(f_ps,f_theta_rhodz,f_ue)  
    90     CASE DEFAULT 
    91        CALL transfert_request(f_ps,req_i1) 
    92        CALL transfert_request(f_theta_rhodz,req_i1) 
    93        CALL transfert_request(f_ue,req_e1_vect) 
    94        CALL transfert_request(f_q,req_i1) 
    95         
    96        args%dt_phys = dt 
    97        DO ind=1,ndomain 
    98           IF (.NOT. assigned_domain(ind)) CYCLE 
    99           CALL swap_dimensions(ind) 
    100           CALL swap_geometry(ind) 
    101        
    102           phis=f_phis(ind) 
    103           ps=f_ps(ind) 
    104           theta_rhodz=f_theta_rhodz(ind) 
    105           ue=f_ue(ind) 
    106           q=f_q(ind) 
    107  
    108           IF(nb_extra_physics_2D>0) args%extra_2D=f_extra_physics_2D(ind) 
    109           IF(nb_extra_physics_3D>0) args%extra_3D=f_extra_physics_3D(ind) 
    110           CALL physics_column(args, phis, ps, theta_rhodz, ue, q) 
    111        ENDDO 
    112  
    113        IF (mod(it,itau_out)==0 ) THEN 
    114           IF(nb_extra_physics_2D>0) CALL writefield("extra_physics_2D",f_extra_physics_2D) 
    115           IF(nb_extra_physics_3D>0) CALL writefield("extra_physics_3D",f_extra_physics_3D) 
    116        ENDIF 
     179    CASE (phys_DCMIP) 
     180       CALL full_physics_dcmip 
    117181    END SELECT 
    118182 
    119   END SUBROUTINE physics 
    120  
    121   SUBROUTINE physics_column(args, phis, ps, theta_rhodz, ue, q) 
     183    DO ind=1,ndomain 
     184       IF (.NOT. assigned_domain(ind)) CYCLE 
     185       CALL swap_dimensions(ind) 
     186       CALL swap_geometry(ind) 
     187       ps=f_ps(ind) 
     188       theta_rhodz=f_theta_rhodz(ind) 
     189       q=f_q(ind) 
     190       dulon=f_dulon(ind) 
     191       dulat=f_dulat(ind) 
     192       CALL unpack_physics(pack_info(ind), ps, theta_rhodz, q, dulon, dulat) 
     193    END DO 
     194 
     195    ! Transfer dulon, dulat 
     196    CALL transfert_request(f_dulon,req_i0) 
     197    CALL transfert_request(f_dulat,req_i0) 
     198 
     199    DO ind=1,ndomain 
     200       IF (.NOT. assigned_domain(ind)) CYCLE 
     201       CALL swap_dimensions(ind) 
     202       CALL swap_geometry(ind) 
     203       ue=f_ue(ind) 
     204       dulon=f_dulon(ind) 
     205       dulat=f_dulat(ind) 
     206       CALL compute_update_velocity(dulon, dulat, ue) 
     207    END DO 
     208 
     209  END SUBROUTINE physics_column_new 
     210 
     211  SUBROUTINE pack_physics(info, phis, ps, theta_rhodz, ue, q) 
    122212    USE icosa 
    123213    USE wind_mod 
    … …  
    125215    USE theta2theta_rhodz_mod 
    126216    USE physics_interface_mod 
     217    IMPLICIT NONE 
     218    TYPE(t_pack_info) :: info 
     219    REAL(rstd) :: phis(iim*jjm) 
     220    REAL(rstd) :: ps(iim*jjm) 
     221    REAL(rstd) :: theta_rhodz(iim*jjm,llm) 
     222    REAL(rstd) :: ue(3*iim*jjm,llm) 
     223    REAL(rstd) :: q(iim*jjm,llm,nqtot) 
     224 
     225    REAL(rstd) :: p(iim*jjm,llm+1) 
     226    REAL(rstd) :: Temp(iim*jjm,llm) 
     227    REAL(rstd) :: uc(iim*jjm,3,llm) 
     228    REAL(rstd) :: ulon(iim*jjm,llm) 
     229    REAL(rstd) :: ulat(iim*jjm,llm) 
     230 
     231    CALL compute_pression(ps,p,0) 
     232    CALL compute_theta_rhodz2temperature(ps,theta_rhodz,Temp,0) 
     233    CALL compute_wind_centered(ue,uc) 
     234    CALL compute_wind_centered_lonlat_compound(uc, ulon, ulat) 
     235 
     236    CALL pack_domain(info, phis, physics_inout%phis) 
     237    CALL pack_domain(info, p, physics_inout%p) 
     238    CALL pack_domain(info, Temp, physics_inout%Temp) 
     239    CALL pack_domain(info, ulon, physics_inout%ulon) 
     240    CALL pack_domain(info, ulat, physics_inout%ulat) 
     241    CALL pack_domain(info, q, physics_inout%q) 
     242  END SUBROUTINE pack_physics 
     243 
     244  SUBROUTINE unpack_physics(info, ps,theta_rhodz, q, dulon, dulat) 
     245    USE icosa 
     246    USE physics_interface_mod 
     247    USE theta2theta_rhodz_mod 
     248    IMPLICIT NONE 
     249    TYPE(t_pack_info) :: info 
     250    REAL(rstd) :: ps(iim*jjm) 
     251    REAL(rstd) :: theta_rhodz(iim*jjm,llm) 
     252    REAL(rstd) :: Temp(iim*jjm,llm) 
     253    REAL(rstd) :: q(iim*jjm,llm,nqtot) 
     254    REAL(rstd) :: dulon(iim*jjm,llm) 
     255    REAL(rstd) :: dulat(iim*jjm,llm) 
     256 
     257    REAL(rstd) :: dq(iim*jjm,llm,nqtot) 
     258    REAL(rstd) :: dTemp(iim*jjm,llm) 
     259    CALL unpack_domain(info, dulon, physics_inout%dulon) 
     260    CALL unpack_domain(info, dulat, physics_inout%dulat) 
     261    CALL unpack_domain(info, dq, physics_inout%dq) 
     262    CALL unpack_domain(info, Temp, physics_inout%Temp) 
     263    CALL unpack_domain(info, dTemp, physics_inout%dTemp) 
     264    q = q + physics_inout%dt_phys * dq 
     265    Temp = Temp + physics_inout%dt_phys * dTemp 
     266    CALL compute_temperature2theta_rhodz(ps,Temp,theta_rhodz,0) 
     267  END SUBROUTINE unpack_physics 
     268 
     269  SUBROUTINE compute_update_velocity(dulon, dulat, ue) 
     270    USE icosa 
     271    USE physics_interface_mod 
     272    USE wind_mod 
     273    IMPLICIT NONE 
     274    REAL(rstd) :: dulon(iim*jjm,llm) 
     275    REAL(rstd) :: dulat(iim*jjm,llm) 
     276    REAL(rstd) :: ue(3*iim*jjm,llm) 
     277    REAL(rstd) :: duc(iim*jjm,3,llm) 
     278    REAL(rstd) :: dt2, due 
     279    INTEGER :: i,j,ij,l 
     280    ! Reconstruct wind tendencies at edges and add to normal wind 
     281    CALL compute_wind_centered_from_lonlat_compound(dulon,dulat,duc) 
     282    dt2=.5*physics_inout%dt_phys 
     283    DO l=1,llm 
     284      DO j=jj_begin,jj_end 
     285        DO i=ii_begin,ii_end 
     286          ij=(j-1)*iim+i 
     287          due = sum( (duc(ij,:,l) + duc(ij+t_right,:,l))*ep_e(ij+u_right,:) ) 
     288          ue(ij+u_right,l) = ue(ij+u_right,l) + dt2*due 
     289 
     290          due = sum( (duc(ij,:,l) + duc(ij+t_lup,:,l))*ep_e(ij+u_lup,:) ) 
     291          ue(ij+u_lup,l)=ue(ij+u_lup,l) + dt2*due 
     292 
     293          due = sum( (duc(ij,:,l) + duc(ij+t_ldown,:,l))*ep_e(ij+u_ldown,:) ) 
     294          ue(ij+u_ldown,l)=ue(ij+u_ldown,l) + dt2*due 
     295        ENDDO 
     296      ENDDO 
     297    ENDDO 
     298  END SUBROUTINE compute_update_velocity 
     299 
     300!--------------------------------- Draft interface -------------------------------------- 
     301 
     302  SUBROUTINE physics_column(args, phis, ps, theta_rhodz, ue, q) 
     303    USE icosa 
     304    USE wind_mod 
     305    USE pression_mod 
     306    USE theta2theta_rhodz_mod 
     307    USE physics_interface_mod 
    127308    USE physics_dcmip_mod 
    128309    IMPLICIT NONE 
    129     TYPE(physics_inout) :: args     
     310    TYPE(t_physics_inout) :: args     
    130311    REAL(rstd) :: phis(iim*jjm) 
    131312    REAL(rstd) :: ps(iim*jjm) 
    … …  
    190371        DO i=ii_begin,ii_end 
    191372          ij=(j-1)*iim+i 
    192           due = sum((uc(ij,:,l) + uc(ij+t_right,:,l))*ep_e(ij+u_right,:) ) 
     373          due = sum( (uc(ij,:,l) + uc(ij+t_right,:,l))*ep_e(ij+u_right,:) ) 
    193374          ue(ij+u_right,l) = ue(ij+u_right,l) + dt2*due 
    194375 

  • codes/icosagcm/trunk/src/physics_dcmip.f90

    r196 r213  
    77 
    88  TYPE(t_field),POINTER :: f_out_i(:) 
    9   TYPE(t_field),POINTER :: f_precl(:) 
    109  REAL(rstd),POINTER :: out_i(:,:) 
    1110 
    12   PUBLIC :: compute_phys_wrap, init_physics 
     11  TYPE(t_field),POINTER  :: f_precl(:) 
     12  REAL(rstd),ALLOCATABLE :: precl_packed(:) 
     13 
     14  PUBLIC :: compute_phys_wrap, init_physics, & 
     15       init_physics_new, full_physics, write_physics 
    1316 
    1417CONTAINS 
    1518 
     19!-------------------------- New interface ---------------------- 
     20 
     21  SUBROUTINE init_physics_new 
     22    USE physics_interface_mod 
     23    IMPLICIT NONE 
     24    INTEGER :: ngrid 
     25    testcase=1 ! OK for 4.2 (moist baroclinic instability) 
     26    CALL getin("dcmip_physics",testcase) 
     27    ngrid = physics_inout%ngrid 
     28    ! Input 
     29    ALLOCATE(physics_inout%Ai(ngrid)) 
     30    ALLOCATE(physics_inout%lon(ngrid)) 
     31    ALLOCATE(physics_inout%lat(ngrid)) 
     32    ALLOCATE(physics_inout%phis(ngrid)) 
     33    ALLOCATE(physics_inout%p(ngrid,llm+1)) 
     34    ALLOCATE(physics_inout%Temp(ngrid,llm)) 
     35    ALLOCATE(physics_inout%ulon(ngrid,llm)) 
     36    ALLOCATE(physics_inout%ulat(ngrid,llm)) 
     37    ALLOCATE(physics_inout%q(ngrid,llm,nqtot)) 
     38    ! Output (tendencies) 
     39    ALLOCATE(physics_inout%dTemp(ngrid,llm)) 
     40    ALLOCATE(physics_inout%dulon(ngrid,llm)) 
     41    ALLOCATE(physics_inout%dulat(ngrid,llm)) 
     42    ALLOCATE(physics_inout%dq(ngrid,llm,nqtot)) 
     43    ! Physics-specific data 
     44    ALLOCATE(precl_packed(ngrid)) 
     45    CALL allocate_field(f_precl, field_t,type_real) 
     46 
     47    PRINT *, 'init_physics_new', SIZE(physics_inout%Ai) 
     48  END SUBROUTINE init_physics_new 
     49 
     50  SUBROUTINE full_physics 
     51    USE physics_interface_mod 
     52    CALL compute_physics(physics_inout%ngrid, physics_inout%dt_phys, & 
     53         physics_inout%lat, physics_inout%p, physics_inout%Temp, &  
     54         physics_inout%ulon, physics_inout%ulat, physics_inout%q(:,:,1), & 
     55         physics_inout%dTemp, physics_inout%dulon, physics_inout%dulat, & 
     56         physics_inout%dq(:,:,1), precl_packed) 
     57  END SUBROUTINE full_physics 
     58 
     59  SUBROUTINE write_physics 
     60    USE output_field_mod 
     61    USE physics_interface_mod 
     62    CALL unpack_field(f_precl, precl_packed) 
     63    CALL output_field("precl",f_precl) 
     64  END SUBROUTINE write_physics 
     65 
     66!------------------------ Draft interface ----------------------- 
     67 
    1668  SUBROUTINE init_physics 
     69    USE physics_interface_mod 
    1770    IMPLICIT NONE 
    1871    testcase=1 ! OK for 4.2 (moist baroclinic instability) 
    1972    CALL getin("dcmip_physics",testcase) 
     73    nb_extra_physics_2D=1 ! precl 
     74    nb_extra_physics_3D=0 
    2075  END SUBROUTINE init_physics 
    2176 
    2277  SUBROUTINE compute_phys_wrap(args) 
    2378    USE physics_interface_mod 
    24     TYPE(physics_inout) :: args 
     79    TYPE(t_physics_inout) :: args 
    2580    CALL compute_physics(args%ngrid, args%dt_phys, args%lat, & 
    2681         args%p, args%Temp, args%ulon, args%ulat, args%q(:,:,1), & 
    2782         args%dTemp, args%dulon, args%dulat, args%dq(:,:,1), args%extra_2D(:,1)) 
    2883  END SUBROUTINE compute_phys_wrap 
     84 
     85!------------------ Interface-independent wrapper --------------------------- 
    2986 
    3087  SUBROUTINE compute_physics(ngrid,dt_phys,lat, p,Temp,u,v,q, dTemp,du,dv,dq, precl) 
    … …  
    54111    REAL(rstd) :: qfi(ngrid,llm) 
    55112 
    56     INTEGER :: i,j,l,ll,ij 
     113    INTEGER :: l,ll,ij 
    57114    REAL(rstd) :: dt_phys, inv_dt 
    58115 
    … …  
    61118 
    62119    DO l=1,llm+1 
    63       DO j=jj_begin,jj_end 
    64         DO i=ii_begin,ii_end 
    65           ij=(j-1)*iim+i 
     120      DO ij=1,ngrid 
    66121          pint(ij,l)=p(ij,llm+2-l) 
    67         ENDDO 
    68122      ENDDO 
    69123    ENDDO 
    … …  
    71125    DO l=1,llm 
    72126       ll=llm+1-l 
    73        DO j=jj_begin,jj_end 
    74           DO i=ii_begin,ii_end 
    75              ij=(j-1)*iim+i 
    76              pmid(ij,l)=0.5*(pint(ij,l)+pint(ij,l+1)) ! Pressure inside layers 
    77              pdel(ij,l)=pint(ij,l+1)-pint(ij,l)       ! Pressure difference between two layers 
    78              ufi(ij,l)=u(ij,ll) 
    79              vfi(ij,l)=v(ij,ll) 
    80              qfi(ij,l)=q(ij,ll) 
    81              Tfi(ij,l)=Temp(ij,ll)/(1+0.608*qfi(ij,l)) 
    82           ENDDO 
     127       DO ij=1,ngrid 
     128          pmid(ij,l)=0.5*(pint(ij,l)+pint(ij,l+1)) ! Pressure inside layers 
     129          pdel(ij,l)=pint(ij,l+1)-pint(ij,l)       ! Pressure difference between two layers 
     130          ufi(ij,l)=u(ij,ll) 
     131          vfi(ij,l)=v(ij,ll) 
     132          qfi(ij,l)=q(ij,ll) 
     133          Tfi(ij,l)=Temp(ij,ll)/(1+0.608*qfi(ij,l))   
    83134       ENDDO 
    84135    ENDDO 
    … …  
    90141    DO l=1,llm 
    91142       ll=llm+1-l 
    92        DO j=jj_begin,jj_end 
    93           DO i=ii_begin,ii_end 
    94              ij=(j-1)*iim+i 
    95              dTemp(ij,l) = inv_dt * ( Tfi(ij,ll)*(1+0.608*qfi(ij,ll)) - Temp(ij,l) ) 
    96              du(ij,l) = inv_dt * (ufi(ij,ll) - u(ij,l)) 
    97              dv(ij,l) = inv_dt * (vfi(ij,ll) - v(ij,l)) 
    98              dq(ij,l)  = inv_dt * (qfi(ij,ll) - q(ij,l)) 
    99           ENDDO 
     143       DO ij=1,ngrid 
     144          dTemp(ij,l) = inv_dt * ( Tfi(ij,ll)*(1+0.608*qfi(ij,ll)) - Temp(ij,l) ) 
     145          du(ij,l) = inv_dt * (ufi(ij,ll) - u(ij,l)) 
     146          dv(ij,l) = inv_dt * (vfi(ij,ll) - v(ij,l)) 
     147          dq(ij,l)  = inv_dt * (qfi(ij,ll) - q(ij,l)) 
    100148       ENDDO 
    101149    ENDDO 

  • codes/icosagcm/trunk/src/physics_interface.f90

    r196 r213  
    55  PRIVATE 
    66 
    7   TYPE :: physics_inout 
    8      ! Input 
     7  INTEGER :: nb_extra_physics_2D, nb_extra_physics_3D 
     8 
     9  TYPE t_physics_inout 
     10     ! Input, time-independent 
    911     INTEGER :: ngrid 
    1012     REAL(rstd) :: dt_phys 
    11      REAL(rstd), DIMENSION(:), POINTER :: lon, lat 
     13     REAL(rstd), DIMENSION(:), POINTER :: Ai, lon, lat, phis 
     14     ! Input, time-dependent 
    1215     REAL(rstd), DIMENSION(:,:), POINTER :: p, Temp, ulon, ulat 
    1316     REAL(rstd), DIMENSION(:,:,:), POINTER :: q 
    … …  
    1821     REAL(rstd), DIMENSION(:,:), POINTER :: extra_2D 
    1922     REAL(rstd), DIMENSION(:,:,:), POINTER :: extra_3D 
    20   END TYPE physics_inout 
    21  
    22   PUBLIC :: physics_inout 
     23  END TYPE t_physics_inout 
     24 
     25!------------------------ (new interface) -------------------------- 
     26! physics_inout is used to exchange information with physics 
     27! Field ngrid is initialized by physics.f90/init_physics. Its other fields 
     28! must be defined by XX/init_physics (where XX =  e.g. physics_dcmip.f90)  
     29! by either pointing to internal data of the physics package 
     30! or by a specific allocation 
     31! size : (ngrid), (ngrid,llm) except p(ngrid,llm+1), (ngrid,llm,nqtot)  
     32 
     33  TYPE(t_physics_inout), SAVE :: physics_inout 
     34  
     35!------------------------ (new interface) -------------------------- 
     36! pack_info contains indices used by pack/unpack routines 
     37! to pack together the data of all the domains managed by the MPI process 
     38! It is initialized by physics.f90/init_physics  
     39  TYPE t_pack_info 
     40     INTEGER :: ngrid, & ! number of non-halo points in that domain 
     41          nseg ! number of segments (contigous parts) in that domain 
     42     ! size and start of each segment : ij domain index, k packed index 
     43     INTEGER, ALLOCATABLE :: n(:), ij(:), k(:) 
     44  END TYPE t_pack_info 
     45 
     46  TYPE(t_pack_info), ALLOCATABLE, SAVE :: pack_info(:) 
     47   
     48  INTERFACE pack_field 
     49     MODULE PROCEDURE pack_2D 
     50     MODULE PROCEDURE pack_3D 
     51     MODULE PROCEDURE pack_4D 
     52  END INTERFACE pack_field 
     53 
     54  INTERFACE unpack_field 
     55     MODULE PROCEDURE unpack_2D 
     56     MODULE PROCEDURE unpack_3D 
     57     MODULE PROCEDURE unpack_4D 
     58  END INTERFACE unpack_field 
     59 
     60  INTERFACE pack_domain 
     61     MODULE PROCEDURE pack_domain_2D 
     62     MODULE PROCEDURE pack_domain_3D 
     63     MODULE PROCEDURE pack_domain_4D 
     64  END INTERFACE pack_domain 
     65 
     66  INTERFACE unpack_domain 
     67     MODULE PROCEDURE unpack_domain_2D 
     68     MODULE PROCEDURE unpack_domain_3D 
     69     MODULE PROCEDURE unpack_domain_4D 
     70  END INTERFACE unpack_domain 
     71 
     72  PUBLIC :: nb_extra_physics_2D, nb_extra_physics_3D, & 
     73       t_physics_inout, physics_inout, & 
     74       t_pack_info, pack_info, init_pack_before, init_pack_after, & 
     75       pack_domain, pack_field, unpack_domain, unpack_field, & 
     76       garbage_3D 
     77 
     78CONTAINS 
     79   
     80    SUBROUTINE init_pack_before 
     81    USE icosa 
     82    IMPLICIT NONE 
     83    INTEGER :: ind, offset 
     84    !$OMP MASTER 
     85    offset=0 
     86    ALLOCATE(pack_info(ndomain)) 
     87    DO ind=1,ndomain 
     88       CALL swap_dimensions(ind) 
     89       CALL swap_geometry(ind) 
     90       CALL count_segments(domain(ind)%own, pack_info(ind)) 
     91       pack_info(ind)%k = pack_info(ind)%k + offset 
     92       offset = offset + pack_info(ind)%ngrid 
     93    END DO 
     94    physics_inout%ngrid = offset 
     95    !$OMP END MASTER 
     96    !$OMP BARRIER 
     97  END SUBROUTINE init_pack_before 
     98 
     99  SUBROUTINE count_segments(own, info) 
     100    USE icosa 
     101    IMPLICIT NONE 
     102    LOGICAL, DIMENSION(:,:) :: own 
     103    TYPE(t_pack_info) :: info 
     104 
     105    INTEGER, DIMENSION(jjm) :: n 
     106    INTEGER :: ngrid, nseg, i, j, jj, k 
     107    INTEGER, PARAMETER :: method=4 
     108    SELECT CASE(method) 
     109    CASE(1) ! Copy all points, including halo (works) 
     110       info%nseg=1 
     111       info%ngrid=iim*jjm 
     112       ALLOCATE(info%n(1)) 
     113       ALLOCATE(info%ij(1)) 
     114       ALLOCATE(info%k(1)) 
     115       info%n(1)=iim*jjm 
     116       info%ij(1)=1 
     117       info%k(1)=1 
     118    CASE(2) ! Copy all points, including halo, one at a time (works, slow ?) 
     119       info%nseg=iim*jjm 
     120       info%ngrid=iim*jjm 
     121       ALLOCATE(info%n(iim*jjm)) 
     122       ALLOCATE(info%ij(iim*jjm)) 
     123       ALLOCATE(info%k(iim*jjm)) 
     124       DO jj=1,iim*jjm 
     125          info%n(jj) =1 
     126          info%ij(jj)=jj 
     127          info%k(jj) =jj 
     128       END DO 
     129    CASE(3) ! Copy non-halo points only, one at a time (works, slow ?) 
     130       n=COUNT(own,1) 
     131       ngrid=SUM(n) 
     132       info%ngrid=ngrid 
     133       info%nseg=ngrid 
     134       ALLOCATE(info%n(ngrid)) 
     135       ALLOCATE(info%ij(ngrid)) 
     136       ALLOCATE(info%k(ngrid)) 
     137       jj=1 
     138       DO j=1,jjm 
     139          DO i=1,iim 
     140             IF(own(i,j)) THEN 
     141                info%n(jj)=1 
     142                info%k(jj)=jj 
     143                info%ij(jj) = iim*(j-1)+i 
     144                jj=jj+1 
     145             END IF 
     146          END DO 
     147       END DO 
     148        
     149    CASE DEFAULT ! Copy non-halo points only, as contiguous segments (works) 
     150       n=0 
     151       n=COUNT(own,1) 
     152       ngrid=SUM(n) 
     153       info%ngrid=ngrid 
     154       nseg=COUNT(n>0) 
     155       info%nseg=nseg 
     156       ALLOCATE(info%n(nseg)) 
     157       ALLOCATE(info%ij(nseg)) 
     158       ALLOCATE(info%k(nseg)) 
     159 
     160       jj=1 
     161       k=1 
     162       DO j=1,jjm 
     163          IF(n(j)>0) THEN 
     164             ! find first .TRUE. value in own(:,j) 
     165             DO i=1,iim 
     166                IF(own(i,j)) THEN 
     167                   info%n(jj)=n(j) 
     168                   info%k(jj)=k 
     169                   info%ij(jj) = iim*(j-1)+i 
     170                   IF(COUNT(own(i:i+n(j)-1,j)) /= n(j)) STOP 
     171                   EXIT 
     172                END IF 
     173             END DO 
     174             k = k + n(j) 
     175             jj=jj+1 
     176          END IF 
     177       END DO 
     178 
     179       IF(k-1/=ngrid) THEN 
     180          PRINT *, 'Total number of grid points inconsistent', k-1, ngrid 
     181          STOP 
     182       END IF 
     183       IF(jj-1/=nseg) THEN 
     184          PRINT *, 'Number of segments inconsistent', jj-1, nseg 
     185          STOP 
     186       END IF 
     187 
     188    END SELECT 
     189     
     190    PRINT *, 'count_segments', info%nseg, info%ngrid, SUM(info%n), COUNT(own), iim*jjm 
     191  END SUBROUTINE count_segments 
     192 
     193  SUBROUTINE init_pack_after 
     194    USE icosa 
     195    IMPLICIT NONE 
     196    INTEGER :: ind, offset 
     197    DO ind=1,ndomain 
     198       IF (.NOT. assigned_domain(ind)) CYCLE 
     199       CALL swap_dimensions(ind) 
     200       CALL swap_geometry(ind) 
     201       CALL pack_domain_2D(pack_info(ind), Ai, physics_inout%Ai) 
     202       CALL pack_lonlat(pack_info(ind)) 
     203    END DO 
     204  END SUBROUTINE init_pack_after 
     205 
     206  SUBROUTINE pack_lonlat(info) 
     207    USE icosa 
     208    IMPLICIT NONE 
     209    TYPE(t_pack_info) :: info 
     210    REAL(rstd) :: lon(iim*jjm), lat(iim*jjm) 
     211    INTEGER :: ij 
     212    DO ij=1,iim*jjm 
     213       CALL xyz2lonlat(xyz_i(ij,:),lon(ij),lat(ij))  
     214    ENDDO 
     215    CALL pack_domain_2D(info, lon, physics_inout%lon) 
     216    CALL pack_domain_2D(info, lat, physics_inout%lat) 
     217  END SUBROUTINE pack_lonlat 
     218 
     219!-------------------------------- Pack / Unpack 2D --------------------------- 
     220 
     221  SUBROUTINE pack_2D(f_2D, packed) 
     222    USE icosa 
     223    IMPLICIT NONE 
     224    TYPE(t_field),POINTER :: f_2D(:) 
     225    REAL(rstd)            :: packed(:) 
     226    REAL(rstd), POINTER   :: loc(:) 
     227    INTEGER :: ind 
     228    DO ind=1,ndomain 
     229       IF (.NOT. assigned_domain(ind)) CYCLE 
     230       loc = f_2D(ind) 
     231       CALL pack_domain_2D(pack_info(ind), loc, packed) 
     232    END DO 
     233  END SUBROUTINE pack_2D 
     234 
     235  SUBROUTINE unpack_2D(f_2D, packed)  
     236    USE icosa 
     237    IMPLICIT NONE 
     238    TYPE(t_field),POINTER :: f_2D(:) 
     239    REAL(rstd)            :: packed(:) 
     240    REAL(rstd), POINTER   :: loc(:) 
     241    INTEGER :: ind 
     242    DO ind=1,ndomain 
     243       IF (.NOT. assigned_domain(ind)) CYCLE 
     244       loc = f_2D(ind) 
     245       CALL unpack_domain_2D(pack_info(ind), loc, packed) 
     246    END DO 
     247  END SUBROUTINE unpack_2D 
     248 
     249  SUBROUTINE pack_domain_2D(info, loc, glob) 
     250    USE icosa 
     251    IMPLICIT NONE 
     252    TYPE(t_pack_info) :: info 
     253    REAL(rstd), DIMENSION(:) :: loc, glob 
     254    INTEGER :: jj,n,k,ij 
     255    DO jj=1, info%nseg 
     256       n = info%n(jj)-1 
     257       ij = info%ij(jj) 
     258       k = info%k(jj) 
     259       glob(k:k+n) = loc(ij:ij+n) 
     260    END DO 
     261  END SUBROUTINE pack_domain_2D 
     262 
     263  SUBROUTINE unpack_domain_2D(info, loc, glob) 
     264    IMPLICIT NONE 
     265    TYPE(t_pack_info) :: info 
     266    REAL(rstd), DIMENSION(:) :: loc, glob 
     267    INTEGER :: jj,n,k,ij 
     268    DO jj=1, info%nseg 
     269       n = info%n(jj)-1 
     270       ij = info%ij(jj) 
     271       k = info%k(jj) 
     272       loc(ij:ij+n) = glob(k:k+n) 
     273    END DO 
     274  END SUBROUTINE unpack_domain_2D 
     275 
     276!-------------------------------- Pack / Unpack 3D --------------------------- 
     277 
     278  SUBROUTINE pack_3D(f_3D, packed) 
     279    USE icosa 
     280    IMPLICIT NONE 
     281    TYPE(t_field),POINTER :: f_3D(:) 
     282    REAL(rstd)            :: packed(:,:) 
     283    REAL(rstd), POINTER   :: loc(:,:) 
     284    INTEGER :: ind 
     285    DO ind=1,ndomain 
     286       IF (.NOT. assigned_domain(ind)) CYCLE 
     287       loc = f_3D(ind) 
     288       CALL pack_domain_3D(pack_info(ind), loc, packed) 
     289    END DO 
     290  END SUBROUTINE pack_3D 
     291 
     292  SUBROUTINE unpack_3D(f_3D, packed)  
     293    USE icosa 
     294    IMPLICIT NONE 
     295    TYPE(t_field),POINTER :: f_3D(:) 
     296    REAL(rstd)            :: packed(:,:) 
     297    REAL(rstd), POINTER   :: loc(:,:) 
     298    INTEGER :: ind 
     299    DO ind=1,ndomain 
     300       IF (.NOT. assigned_domain(ind)) CYCLE 
     301       loc = f_3D(ind) 
     302       CALL unpack_domain_3D(pack_info(ind), loc, packed) 
     303    END DO 
     304  END SUBROUTINE unpack_3D 
     305 
     306  SUBROUTINE pack_domain_3D(info, loc, glob) 
     307    IMPLICIT NONE 
     308    TYPE(t_pack_info) :: info 
     309    REAL(rstd), DIMENSION(:,:) :: loc, glob 
     310    INTEGER :: jj,n,k,ij 
     311    DO jj=1, info%nseg 
     312       n = info%n(jj)-1 
     313       ij = info%ij(jj) 
     314       k = info%k(jj) 
     315       glob(k:k+n,:) = loc(ij:ij+n,:) 
     316    END DO 
     317  END SUBROUTINE pack_domain_3D 
     318 
     319  SUBROUTINE unpack_domain_3D(info, loc, glob) 
     320    IMPLICIT NONE 
     321    TYPE(t_pack_info) :: info 
     322    REAL(rstd), DIMENSION(:,:) :: loc, glob 
     323    INTEGER :: jj,n,k,ij 
     324    DO jj=1, info%nseg 
     325       n = info%n(jj)-1 
     326       ij = info%ij(jj) 
     327       k = info%k(jj) 
     328       loc(ij:ij+n,:) = glob(k:k+n,:) 
     329    END DO     
     330  END SUBROUTINE unpack_domain_3D 
     331 
     332  SUBROUTINE garbage_3D(loc,own) 
     333    USE icosa 
     334    IMPLICIT NONE 
     335    LOGICAL :: own(iim,jjm) 
     336    REAL(rstd) :: loc(iim*jjm,llm) 
     337    INTEGER :: i,j,ij 
     338    ! write garbage in non-owned points 
     339    DO j=1,jjm 
     340       DO i=1,iim 
     341          IF(.NOT.own(i,j)) THEN 
     342             ij=iim*(j-1)+i 
     343             loc(ij,:)=-1e30 
     344          END IF 
     345       END DO 
     346    END DO    
     347  END SUBROUTINE garbage_3D 
     348 
     349!-------------------------------- Pack / Unpack 4D --------------------------- 
     350 
     351  SUBROUTINE pack_4D(f_4D, packed) 
     352    USE icosa 
     353    IMPLICIT NONE 
     354    TYPE(t_field),POINTER :: f_4D(:) 
     355    REAL(rstd)            :: packed(:,:,:) 
     356    REAL(rstd), POINTER   :: loc(:,:,:) 
     357    INTEGER :: ind 
     358    DO ind=1,ndomain 
     359       IF (.NOT. assigned_domain(ind)) CYCLE 
     360       loc = f_4D(ind) 
     361       CALL pack_domain_4D(pack_info(ind), loc, packed) 
     362    END DO 
     363  END SUBROUTINE pack_4D 
     364 
     365  SUBROUTINE unpack_4D(f_4D, packed)  
     366    USE icosa 
     367    IMPLICIT NONE 
     368    TYPE(t_field),POINTER :: f_4D(:) 
     369    REAL(rstd)            :: packed(:,:,:) 
     370    REAL(rstd), POINTER   :: loc(:,:,:) 
     371    INTEGER :: ind 
     372    DO ind=1,ndomain 
     373       IF (.NOT. assigned_domain(ind)) CYCLE 
     374       loc = f_4D(ind) 
     375       CALL unpack_domain_4D(pack_info(ind), loc, packed) 
     376    END DO 
     377  END SUBROUTINE unpack_4D 
     378 
     379  SUBROUTINE pack_domain_4D(info, loc, glob) 
     380    IMPLICIT NONE 
     381    TYPE(t_pack_info) :: info 
     382    REAL(rstd), DIMENSION(:,:,:) :: loc, glob 
     383    INTEGER :: jj,n,k,ij 
     384    DO jj=1, info%nseg 
     385       n = info%n(jj)-1 
     386       ij = info%ij(jj) 
     387       k = info%k(jj) 
     388       glob(k:k+n,:,:) = loc(ij:ij+n,:,:) 
     389    END DO 
     390  END SUBROUTINE pack_domain_4D 
     391 
     392  SUBROUTINE unpack_domain_4D(info, loc, glob) 
     393    IMPLICIT NONE 
     394    TYPE(t_pack_info) :: info 
     395    REAL(rstd), DIMENSION(:,:,:) :: loc, glob 
     396    INTEGER :: jj,n,k,ij 
     397    DO jj=1, info%nseg 
     398       n = info%n(jj)-1 
     399       ij = info%ij(jj) 
     400       k = info%k(jj) 
     401       loc(ij:ij+n,:,:) = glob(k:k+n,:,:) 
     402    END DO 
     403  END SUBROUTINE unpack_domain_4D 
    23404 
    24405END MODULE physics_interface_mod 
Note: See TracChangeset for help on using the changeset viewer.