Ignore:
Timestamp:
06/25/13 09:32:16 (11 years ago)
Author:
dubos
Message:

Towards Lagrangian vertical coordinate (not there yet)

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r158 r159  
    1111  TYPE(t_message) :: req_ps0, req_theta_rhodz0, req_u0, req_q0 
    1212 
    13   TYPE(t_field),POINTER :: f_phis(:) 
    1413  TYPE(t_field),POINTER :: f_q(:) 
    15   TYPE(t_field),POINTER :: f_dtheta(:), f_rhodz(:) 
    16   TYPE(t_field),POINTER :: f_ps(:),f_psm1(:), f_psm2(:) 
    17   TYPE(t_field),POINTER :: f_u(:),f_um1(:),f_um2(:) 
    18   TYPE(t_field),POINTER :: f_theta_rhodz(:),f_theta_rhodzm1(:),f_theta_rhodzm2(:) 
    19   TYPE(t_field),POINTER :: f_dps(:),f_dpsm1(:), f_dpsm2(:) 
    20   TYPE(t_field),POINTER :: f_du(:),f_dum1(:),f_dum2(:) 
    21   TYPE(t_field),POINTER :: f_dtheta_rhodz(:),f_dtheta_rhodzm1(:),f_dtheta_rhodzm2(:) 
     14  TYPE(t_field),POINTER :: f_rhodz(:), f_mass(:), f_massm1(:), f_massm2(:), f_dmass(:) 
     15  TYPE(t_field),POINTER :: f_phis(:), f_ps(:),f_psm1(:), f_psm2(:), f_dps(:) 
     16  TYPE(t_field),POINTER :: f_u(:),f_um1(:),f_um2(:), f_du(:) 
     17  TYPE(t_field),POINTER :: f_theta_rhodz(:),f_theta_rhodzm1(:),f_theta_rhodzm2(:), f_dtheta_rhodz(:) 
    2218  TYPE(t_field),POINTER :: f_hflux(:), f_wflux(:), f_hfluxt(:), f_wfluxt(:)   
    2319 
    … …  
    3430  USE caldyn_mod 
    3531  USE etat0_mod 
     32  USE disvert_mod 
    3633  USE guided_mod 
    3734  USE advect_tracer_mod 
    … …  
    4542  IMPLICIT NONE 
    4643 
    47     CHARACTER(len=255) :: scheme_name 
    48  
    49     CALL allocate_field(f_dps,field_t,type_real) 
    50     CALL allocate_field(f_du,field_u,type_real,llm) 
    51     CALL allocate_field(f_dtheta_rhodz,field_t,type_real,llm) 
    52 ! Model state at current time step (RK/MLF/Euler) 
    53     CALL allocate_field(f_phis,field_t,type_real) 
    54     CALL allocate_field(f_ps,field_t,type_real) 
    55     CALL allocate_field(f_u,field_u,type_real,llm) 
    56     CALL allocate_field(f_theta_rhodz,field_t,type_real,llm) 
    57 ! Model state at previous time step (RK/MLF) 
    58     CALL allocate_field(f_psm1,field_t,type_real) 
    59     CALL allocate_field(f_um1,field_u,type_real,llm) 
    60     CALL allocate_field(f_theta_rhodzm1,field_t,type_real,llm) 
    61 ! Tracers 
    62     CALL allocate_field(f_q,field_t,type_real,llm,nqtot) 
    63     CALL allocate_field(f_rhodz,field_t,type_real,llm) 
    64 ! Mass fluxes 
    65     CALL allocate_field(f_hflux,field_u,type_real,llm)    ! instantaneous mass fluxes 
    66     CALL allocate_field(f_wflux,field_t,type_real,llm+1)  ! computed by caldyn 
    67     CALL allocate_field(f_hfluxt,field_u,type_real,llm)   ! mass "fluxes" accumulated in time 
    68     CALL allocate_field(f_wfluxt,field_t,type_real,llm+1) 
    69  
     44    CHARACTER(len=255) :: def 
    7045 
    7146!---------------------------------------------------- 
    … …  
    12297 
    12398 
    124  
    125  
    126  
    127     scheme_name='runge_kutta' 
    128     CALL getin('scheme',scheme_name) 
    129  
    130     SELECT CASE (TRIM(scheme_name)) 
     99    def='eta_mass' 
     100    CALL getin('caldyn_eta',def) 
     101    SELECT CASE(TRIM(def)) 
     102    CASE('eta_mass') 
     103       caldyn_eta=eta_mass 
     104    CASE('eta_lag') 
     105       caldyn_eta=eta_lag 
     106    CASE DEFAULT 
     107       IF (is_mpi_root) PRINT *,'Bad selector for variable caldyn_eta : <', TRIM(def),'> options are <eta_mass>, <eta_lag>' 
     108       STOP 
     109    END SELECT 
     110    IF (is_mpi_root) PRINT *, 'caldyn_eta=',TRIM(def), caldyn_eta, eta_lag, eta_mass 
     111 
     112! Time-independant orography 
     113    CALL allocate_field(f_phis,field_t,type_real,name='phis') 
     114! Trends 
     115    CALL allocate_field(f_du,field_u,type_real,llm,name='du') 
     116    CALL allocate_field(f_dtheta_rhodz,field_t,type_real,llm,name='dtheta_rhodz') 
     117! Model state at current time step (RK/MLF/Euler) 
     118    CALL allocate_field(f_ps,field_t,type_real, name='ps') 
     119    CALL allocate_field(f_mass,field_t,type_real,llm,name='mass') 
     120    CALL allocate_field(f_u,field_u,type_real,llm,name='u') 
     121    CALL allocate_field(f_theta_rhodz,field_t,type_real,llm,name='theta_rhodz') 
     122! Model state at previous time step (RK/MLF) 
     123    CALL allocate_field(f_um1,field_u,type_real,llm,name='um1') 
     124    CALL allocate_field(f_theta_rhodzm1,field_t,type_real,llm,name='theta_rhodzm1') 
     125! Tracers 
     126    CALL allocate_field(f_q,field_t,type_real,llm,nqtot) 
     127    CALL allocate_field(f_rhodz,field_t,type_real,llm,name='rhodz') 
     128! Mass fluxes 
     129    CALL allocate_field(f_hflux,field_u,type_real,llm)    ! instantaneous mass fluxes 
     130    CALL allocate_field(f_hfluxt,field_u,type_real,llm)   ! mass "fluxes" accumulated in time 
     131    CALL allocate_field(f_wflux,field_t,type_real,llm+1)  ! vertical mass fluxes 
     132 
     133    IF(caldyn_eta == eta_mass) THEN ! eta = mass coordinate (default) 
     134       CALL allocate_field(f_dps,field_t,type_real,name='dps') 
     135       CALL allocate_field(f_psm1,field_t,type_real,name='psm1') 
     136       CALL allocate_field(f_wfluxt,field_t,type_real,llm+1,name='wfluxt') 
     137       ! the following are unused but must point to something 
     138       f_massm1 => f_mass 
     139       f_dmass => f_mass 
     140    ELSE ! eta = Lagrangian vertical coordinate 
     141       CALL allocate_field(f_mass,field_t,type_real,llm) 
     142       CALL allocate_field(f_massm1,field_t,type_real,llm) 
     143       CALL allocate_field(f_dmass,field_t,type_real,llm) 
     144       ! the following are unused but must point to something 
     145       f_wfluxt => f_wflux 
     146       f_dps  => f_phis 
     147       f_psm1 => f_phis 
     148    END IF 
     149 
     150 
     151    def='runge_kutta' 
     152    CALL getin('scheme',def) 
     153 
     154    SELECT CASE (TRIM(def)) 
    131155      CASE('euler') 
    132156        scheme=euler 
    … …  
    141165        nb_stage=matsuno_period+1 
    142166     ! Model state 2 time steps ago (MLF) 
    143         CALL allocate_field(f_psm2,field_t,type_real) 
    144167        CALL allocate_field(f_theta_rhodzm2,field_t,type_real,llm) 
    145168        CALL allocate_field(f_um2,field_u,type_real,llm) 
     169        IF(caldyn_eta == eta_mass) THEN ! eta = mass coordinate (default) 
     170           CALL allocate_field(f_psm2,field_t,type_real) 
     171           ! the following are unused but must point to something 
     172           f_massm2 => f_mass 
     173        ELSE ! eta = Lagrangian vertical coordinate 
     174           CALL allocate_field(f_massm2,field_t,type_real,llm) 
     175           ! the following are unused but must point to something 
     176           f_psm2 => f_phis 
     177        END IF 
     178 
    146179    CASE default 
    147        PRINT*,'Bad selector for variable scheme : <', TRIM(scheme_name),             & 
     180       PRINT*,'Bad selector for variable scheme : <', TRIM(def),             & 
    148181            ' > options are <euler>, <runge_kutta>, <leapfrog_matsuno>' 
    149182       STOP 
    … …  
    157190    CALL init_check_conserve 
    158191    CALL init_physics 
    159     CALL etat0(f_ps,f_phis,f_theta_rhodz,f_u, f_q) 
     192    CALL etat0(f_ps,f_mass,f_phis,f_theta_rhodz,f_u, f_q) 
    160193 
    161194    CALL transfert_request(f_phis,req_i0)  
    … …  
    173206  USE icosa 
    174207  USE dissip_gcm_mod 
     208  USE disvert_mod 
    175209  USE caldyn_mod 
     210  USE caldyn_gcm_mod, ONLY : req_ps 
    176211  USE etat0_mod 
    177212  USE guided_mod 
    … …  
    184219  USE check_conserve_mod 
    185220  IMPLICIT NONE   
    186     REAL(rstd),POINTER :: phis(:) 
    187221    REAL(rstd),POINTER :: q(:,:,:) 
    188     REAL(rstd),POINTER :: ps(:) ,psm1(:), psm2(:) 
    189     REAL(rstd),POINTER :: u(:,:) , um1(:,:), um2(:,:) 
    190     REAL(rstd),POINTER :: rhodz(:,:), theta_rhodz(:,:) , theta_rhodzm1(:,:), theta_rhodzm2(:,:) 
    191     REAL(rstd),POINTER :: dps(:), dpsm1(:), dpsm2(:) 
    192     REAL(rstd),POINTER :: du(:,:), dum1(:,:), dum2(:,:) 
    193     REAL(rstd),POINTER :: dtheta_rhodz(:,:),dtheta_rhodzm1(:,:),dtheta_rhodzm2(:,:) 
     222    REAL(rstd),POINTER :: phis(:), ps(:) ,psm1(:), psm2(:), dps(:) 
     223    REAL(rstd),POINTER :: u(:,:) , um1(:,:), um2(:,:), du(:,:) 
     224    REAL(rstd),POINTER :: rhodz(:,:), mass(:,:), massm1(:,:), massm2(:,:), dmass(:,:) 
     225    REAL(rstd),POINTER :: theta_rhodz(:,:) , theta_rhodzm1(:,:), theta_rhodzm2(:,:), dtheta_rhodz(:,:) 
    194226    REAL(rstd),POINTER :: hflux(:,:),wflux(:,:),hfluxt(:,:),wfluxt(:,:) 
    195227 
    196228    INTEGER :: ind 
    197229    INTEGER :: it,i,j,n,  stage 
    198     CHARACTER(len=255) :: scheme_name 
    199230    LOGICAL :: fluxt_zero(ndomain) ! set to .TRUE. to start accumulating fluxes in time 
    200231    LOGICAL, PARAMETER :: check=.FALSE. 
    201232 
    202     CALL caldyn_BC(f_phis, f_wflux) 
     233    CALL caldyn_BC(f_phis, f_wflux) ! set constant values in first/last interfaces 
    203234 
    204235!$OMP BARRIER 
    … …  
    234265    DO stage=1,nb_stage 
    235266       CALL caldyn((stage==1) .AND. (MOD(it,itau_out)==0), & 
    236             f_phis,f_ps,f_theta_rhodz,f_u, f_q, & 
     267            f_phis,f_ps,f_mass,f_theta_rhodz,f_u, f_q, & 
    237268            f_hflux, f_wflux, f_dps, f_dtheta_rhodz, f_du) 
    238269       SELECT CASE (scheme) 
    … …  
    333364 
    334365    SUBROUTINE RK_scheme(stage) 
    335       USE caldyn_gcm_mod 
    336366      IMPLICIT NONE 
    337367      INTEGER :: ind, stage 
    … …  
    344374      tau = dt*coef(stage) 
    345375 
    346       DO ind=1,ndomain 
    347          CALL swap_dimensions(ind) 
    348          CALL swap_geometry(ind) 
    349          ps=f_ps(ind)    
    350          psm1=f_psm1(ind)  
    351          dps=f_dps(ind)  
    352           
    353          IF (stage==1) THEN ! first stage : save model state in XXm1 
    354            IF (omp_first) THEN 
    355              DO j=jj_begin,jj_end 
    356                DO i=ii_begin,ii_end 
    357                  ij=(j-1)*iim+i 
    358                  psm1(ij)=ps(ij) 
     376      ! if mass coordinate, deal with ps first on one core 
     377      IF(caldyn_eta==eta_mass) THEN 
     378         IF (omp_first) THEN 
     379            DO ind=1,ndomain 
     380               CALL swap_dimensions(ind) 
     381               CALL swap_geometry(ind) 
     382               ps=f_ps(ind)    
     383               psm1=f_psm1(ind)  
     384               dps=f_dps(ind)  
     385                
     386               IF (stage==1) THEN ! first stage : save model state in XXm1 
     387                  DO j=jj_begin,jj_end 
     388                     DO i=ii_begin,ii_end 
     389                        ij=(j-1)*iim+i 
     390                        psm1(ij)=ps(ij) 
     391                     ENDDO 
     392                  ENDDO 
     393               ENDIF 
     394             
     395               ! updates are of the form x1 := x0 + tau*f(x1) 
     396               DO j=jj_begin,jj_end 
     397                  DO i=ii_begin,ii_end 
     398                     ij=(j-1)*iim+i 
     399                     ps(ij)=psm1(ij)+tau*dps(ij) 
     400                  ENDDO 
    359401               ENDDO 
    360              ENDDO 
    361            ENDIF 
    362          END IF 
    363           
    364          ! updates are of the form x1 := x0 + tau*f(x1) 
    365            IF (omp_first) THEN 
    366              DO j=jj_begin,jj_end 
    367                DO i=ii_begin,ii_end 
    368                  ij=(j-1)*iim+i 
    369                  ps(ij)=psm1(ij)+tau*dps(ij) 
    370                ENDDO 
    371              ENDDO 
    372            ENDIF 
    373         
    374        ENDDO 
    375  
    376        CALL send_message(f_ps,req_ps)       
    377  
    378  
    379        
     402            ENDDO 
     403         ENDIF 
     404    
     405         CALL send_message(f_ps,req_ps) 
     406      END IF 
     407 
     408      ! now deal with other prognostic variables 
    380409      DO ind=1,ndomain 
    381410         CALL swap_dimensions(ind) 
    … …  
    386415          
    387416         IF (stage==1) THEN ! first stage : save model state in XXm1 
    388                
    389          DO l=ll_begin,ll_end 
     417             
     418           DO l=ll_begin,ll_end 
    390419             DO j=jj_begin,jj_end 
    391420               DO i=ii_begin,ii_end 
    … …  
    398427             ENDDO 
    399428           ENDDO 
     429 
     430           IF(caldyn_eta==eta_lag) THEN ! mass = additional prognostic variable 
     431              DO l=ll_begin,ll_end 
     432                 DO j=jj_begin,jj_end 
     433                    DO i=ii_begin,ii_end 
     434                       ij=(j-1)*iim+i 
     435                       massm1(ij,l)=mass(ij,l) 
     436                    ENDDO 
     437                 ENDDO 
     438              ENDDO 
     439           END IF 
    400440 
    401441         END IF 
    … …  
    413453           ENDDO 
    414454         ENDDO 
    415  
     455         IF(caldyn_eta==eta_lag) THEN ! mass = additional prognostic variable 
     456            DO l=ll_begin,ll_end 
     457               DO j=jj_begin,jj_end 
     458                  DO i=ii_begin,ii_end 
     459                     ij=(j-1)*iim+i 
     460                     mass(ij,l)=massm1(ij,l)+tau*dmass(ij,l) 
     461                  ENDDO 
     462               ENDDO 
     463            ENDDO 
     464         END IF 
     465          
    416466         IF(stage==nb_stage) THEN ! accumulate mass fluxes at last stage 
    417467            hflux=f_hflux(ind);     hfluxt=f_hfluxt(ind) 
    … …  
    476526  END SUBROUTINE timeloop     
    477527 
    478   SUBROUTINE compute_rhodz(comp, ps, rhodz) 
    479     USE icosa 
    480     USE disvert_mod 
    481     USE omp_para 
    482     LOGICAL, INTENT(IN) :: comp ! .TRUE. to compute, .FALSE. to check 
    483     REAL(rstd), INTENT(IN) :: ps(iim*jjm) 
    484     REAL(rstd), INTENT(INOUT) :: rhodz(iim*jjm,llm) 
    485     REAL(rstd) :: m, err 
    486     INTEGER :: l,i,j,ij,dd 
    487     err=0. 
    488  
    489     dd=0 
    490  
    491     DO l = ll_begin, ll_end 
    492        DO j=jj_begin-dd,jj_end+dd 
    493           DO i=ii_begin-dd,ii_end+dd 
    494              ij=(j-1)*iim+i 
    495              m = ( ap(l) - ap(l+1) + (bp(l)-bp(l+1))*ps(ij) )/g  
    496              IF(comp) THEN 
    497                 rhodz(ij,l) = m 
    498              ELSE 
    499                 err = MAX(err,abs(m-rhodz(ij,l))) 
    500              END IF 
    501           ENDDO 
    502        ENDDO 
    503     ENDDO 
    504  
    505     IF(.NOT. comp) THEN 
    506        IF(err>1e-10) THEN 
    507           PRINT *, 'Discrepancy between ps and rhodz detected', err 
    508           STOP 
    509        ELSE 
    510 !          PRINT *, 'No discrepancy between ps and rhodz detected' 
    511        END IF 
    512     END IF 
    513  
    514   END SUBROUTINE compute_rhodz 
    515  
    516   SUBROUTINE accumulate_fluxes(hflux,wflux, hfluxt,wfluxt, tau,fluxt_zero) 
     528 SUBROUTINE accumulate_fluxes(hflux,wflux, hfluxt,wfluxt, tau,fluxt_zero) 
    517529    USE icosa 
    518530    USE omp_para 
     531    USE disvert_mod 
     532    IMPLICIT NONE 
    519533    REAL(rstd), INTENT(IN)    :: hflux(3*iim*jjm,llm), wflux(iim*jjm,llm+1) 
    520534    REAL(rstd), INTENT(INOUT) :: hfluxt(3*iim*jjm,llm), wfluxt(iim*jjm,llm+1) 
    … …  
    538552       ENDDO 
    539553 
    540        DO l=ll_begin,ll_endp1 
    541          DO j=jj_begin,jj_end 
    542            DO i=ii_begin,ii_end 
    543              ij=(j-1)*iim+i 
    544              wfluxt(ij,l) = tau*wflux(ij,l) 
    545            ENDDO 
    546          ENDDO 
    547        ENDDO 
     554       IF(caldyn_eta==eta_mass) THEN ! no need for vertical fluxes if eta_lag 
     555          DO l=ll_begin,ll_endp1 
     556             DO j=jj_begin,jj_end 
     557                DO i=ii_begin,ii_end 
     558                   ij=(j-1)*iim+i 
     559                   wfluxt(ij,l) = tau*wflux(ij,l) 
     560                ENDDO 
     561             ENDDO 
     562          ENDDO 
     563       END IF 
    548564    ELSE 
    549565 
    … …  
    559575       ENDDO 
    560576 
    561        DO l=ll_begin,ll_endp1 
    562          DO j=jj_begin,jj_end 
    563            DO i=ii_begin,ii_end 
    564              ij=(j-1)*iim+i 
    565              wfluxt(ij,l) = wfluxt(ij,l)+tau*wflux(ij,l) 
    566            ENDDO 
    567          ENDDO 
    568        ENDDO 
    569  
    570     END IF 
     577       IF(caldyn_eta==eta_mass) THEN ! no need for vertical fluxes if eta_lag 
     578          DO l=ll_begin,ll_endp1 
     579             DO j=jj_begin,jj_end 
     580                DO i=ii_begin,ii_end 
     581                   ij=(j-1)*iim+i 
     582                   wfluxt(ij,l) = wfluxt(ij,l)+tau*wflux(ij,l) 
     583                ENDDO 
     584             ENDDO 
     585          ENDDO 
     586       END IF 
     587 
     588   END IF 
    571589 
    572590  END SUBROUTINE accumulate_fluxes 
Note: See TracChangeset for help on using the changeset viewer.