Changeset 159 for codes/icosagcm


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

Towards Lagrangian vertical coordinate (not there yet)

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

Legend:

Unmodified
Added
Removed

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

    r157 r159  
    3131  END SUBROUTINE init_caldyn 
    3232   
    33   SUBROUTINE caldyn(write_out,f_phis, f_ps, f_theta_rhodz, f_u, f_q, & 
     33  SUBROUTINE caldyn(write_out,f_phis, f_ps, f_mass, f_theta_rhodz, f_u, f_q, & 
    3434       f_hflux, f_wflux, f_dps, f_dtheta_rhodz, f_du) 
    3535  USE icosa 
    … …  
    4040  TYPE(t_field),POINTER :: f_phis(:) 
    4141  TYPE(t_field),POINTER :: f_ps(:) 
     42  TYPE(t_field),POINTER :: f_mass(:) 
    4243  TYPE(t_field),POINTER :: f_theta_rhodz(:) 
    4344  TYPE(t_field),POINTER :: f_u(:) 
    … …  
    5152    SELECT CASE (TRIM(caldyn_type)) 
    5253      CASE('gcm') 
    53         CALL caldyn_gcm(write_out,f_phis, f_ps, f_theta_rhodz, f_u, f_q, & 
     54        CALL caldyn_gcm(write_out,f_phis, f_ps, f_mass, f_theta_rhodz, f_u, f_q, & 
    5455             f_hflux, f_wflux, f_dps, f_dtheta_rhodz, f_du) 
    5556      CASE('adv') 

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

    r157 r159  
    55 
    66  INTEGER, PARAMETER :: energy=1, enstrophy=2 
    7   TYPE(t_field),POINTER :: f_out_u(:), f_rhodz(:), f_qu(:) 
    8   REAL(rstd),POINTER :: out_u(:,:), p(:,:), rhodz(:,:), qu(:,:) 
     7  TYPE(t_field),POINTER :: f_out_u(:), f_qu(:) 
     8  REAL(rstd),POINTER :: out_u(:,:), p(:,:), qu(:,:) 
    99 
    1010  TYPE(t_field),POINTER :: f_buf_i(:), f_buf_ulon(:), f_buf_ulat(:), f_buf_u3d(:) 
    … …  
    1818  TYPE(t_field),POINTER :: f_wwuu(:) 
    1919    
    20   INTEGER :: caldyn_hydrostat, caldyn_conserv 
     20  INTEGER :: caldyn_conserv 
    2121   
    2222  TYPE(t_message) :: req_ps, req_theta_rhodz, req_u, req_qu 
    2323 
    24   PUBLIC init_caldyn, caldyn_BC, caldyn, write_output_fields,req_ps 
     24  PUBLIC init_caldyn, caldyn_BC, caldyn, write_output_fields,req_ps, & 
     25       caldyn_eta, eta_mass, eta_lag 
    2526 
    2627CONTAINS 
    … …  
    4142       caldyn_conserv=enstrophy 
    4243    CASE DEFAULT 
    43        IF (is_mpi_root) PRINT *,'Bad selector for variable caldyn_conserv : <', TRIM(def),'> options are <energy>, <enstrophy>' 
     44       IF (is_mpi_root) PRINT *,'Bad selector for variable caldyn_conserv : <', & 
     45            TRIM(def),'> options are <energy>, <enstrophy>' 
    4446       STOP 
    4547    END SELECT 
    … …  
    5557 
    5658    CALL allocate_field(f_out_u,field_u,type_real,llm)  
    57     CALL allocate_field(f_rhodz,field_t,type_real,llm)  
    5859    CALL allocate_field(f_qu,field_u,type_real,llm)  
    5960   
    60     CALL allocate_field(f_buf_i,field_t,type_real,llm) 
    61     CALL allocate_field(f_buf_p,field_t,type_real,llm+1)  
    62     CALL allocate_field(f_buf_u3d,field_t,type_real,3,llm)  ! 3D vel at cell centers 
     61    CALL allocate_field(f_buf_i,   field_t,type_real,llm) 
     62    CALL allocate_field(f_buf_p,   field_t,type_real,llm+1)  
     63    CALL allocate_field(f_buf_u3d, field_t,type_real,3,llm)  ! 3D vel at cell centers 
    6364    CALL allocate_field(f_buf_ulon,field_t,type_real,llm) 
    6465    CALL allocate_field(f_buf_ulat,field_t,type_real,llm) 
    65     CALL allocate_field(f_buf_v,field_z,type_real,llm) 
    66     CALL allocate_field(f_buf_s,field_t,type_real) 
    67  
    68     CALL allocate_field(f_theta,field_t,type_real,llm)     ! potential temperature 
    69     CALL allocate_field(f_pk,field_t,type_real,llm) 
    70     CALL allocate_field(f_geopot,field_t,type_real,llm+1)  ! geopotential 
    71     CALL allocate_field(f_divm,field_t,type_real,llm)      ! mass flux divergence 
    72     CALL allocate_field(f_wwuu,field_u,type_real,llm+1)  
     66    CALL allocate_field(f_buf_v,   field_z,type_real,llm) 
     67    CALL allocate_field(f_buf_s,   field_t,type_real) 
     68 
     69    CALL allocate_field(f_theta, field_t,type_real,llm,  name='theta')   ! potential temperature 
     70    CALL allocate_field(f_pk,    field_t,type_real,llm,  name='pk') 
     71    CALL allocate_field(f_geopot,field_t,type_real,llm+1,name='geopot')  ! geopotential 
     72    CALL allocate_field(f_divm,  field_t,type_real,llm,  name='divm')    ! mass flux divergence 
     73    CALL allocate_field(f_wwuu,  field_u,type_real,llm+1,name='wwuu') 
    7374 
    7475  END SUBROUTINE allocate_caldyn 
    … …  
    119120  END SUBROUTINE caldyn_BC 
    120121    
    121   SUBROUTINE caldyn(write_out,f_phis, f_ps, f_theta_rhodz, f_u, f_q, & 
     122  SUBROUTINE caldyn(write_out,f_phis, f_ps, f_mass, f_theta_rhodz, f_u, f_q, & 
    122123       f_hflux, f_wflux, f_dps, f_dtheta_rhodz, f_du) 
    123124    USE icosa 
    … …  
    132133    TYPE(t_field),POINTER :: f_phis(:) 
    133134    TYPE(t_field),POINTER :: f_ps(:) 
     135    TYPE(t_field),POINTER :: f_mass(:) 
    134136    TYPE(t_field),POINTER :: f_theta_rhodz(:) 
    135137    TYPE(t_field),POINTER :: f_u(:) 
    … …  
    140142    TYPE(t_field),POINTER :: f_du(:) 
    141143     
    142     REAL(rstd),POINTER :: phis(:), ps(:), dps(:) 
    143     REAL(rstd),POINTER :: theta_rhodz(:,:), dtheta_rhodz(:,:) 
     144    REAL(rstd),POINTER :: ps(:), dps(:) 
     145    REAL(rstd),POINTER :: mass(:,:), theta_rhodz(:,:), dtheta_rhodz(:,:) 
    144146    REAL(rstd),POINTER :: u(:,:), du(:,:), hflux(:,:), wflux(:,:) 
    145     REAL(rstd),POINTER :: rhodz(:,:), qu(:,:) 
     147    REAL(rstd),POINTER :: qu(:,:) 
    146148 
    147149! temporary shared variable 
    … …  
    178180          u=f_u(ind) 
    179181          theta_rhodz = f_theta_rhodz(ind) 
    180           rhodz=f_rhodz(ind) 
     182          mass=f_mass(ind) 
    181183          theta = f_theta(ind) 
    182184          qu=f_qu(ind) 
    183           CALL compute_pvort(ps,u,theta_rhodz, rhodz,theta,qu) 
     185          CALL compute_pvort(ps,u,theta_rhodz, mass,theta,qu) 
    184186       ENDDO 
    185187 
    … …  
    192194          u=f_u(ind) 
    193195          theta_rhodz=f_theta_rhodz(ind) 
    194           rhodz=f_rhodz(ind) 
     196          mass=f_mass(ind) 
    195197          theta = f_theta(ind) 
    196198          qu=f_qu(ind) 
    197           phis=f_phis(ind) 
    198199          pk = f_pk(ind) 
    199200          geopot = f_geopot(ind)   
    200           CALL compute_geopot(ps,phis,rhodz,theta, pk,geopot) 
     201          CALL compute_geopot(ps,mass,theta, pk,geopot) 
    201202          hflux=f_hflux(ind) 
    202203          divm = f_divm(ind) 
    203204          dtheta_rhodz=f_dtheta_rhodz(ind) 
    204205          du=f_du(ind) 
    205           CALL compute_caldyn_horiz(u,rhodz,qu,theta,pk,geopot, hflux,divm,dtheta_rhodz,du) 
     206          CALL compute_caldyn_horiz(u,mass,qu,theta,pk,geopot, hflux,divm,dtheta_rhodz,du) 
    206207          wflux=f_wflux(ind) 
    207208          wwuu=f_wwuu(ind) 
    208209          dps=f_dps(ind) 
    209           CALL compute_caldyn_vert(u,theta,rhodz,divm, wflux,wwuu, dps, dtheta_rhodz, du) 
     210          CALL compute_caldyn_vert(u,theta,mass,divm, wflux,wwuu, dps, dtheta_rhodz, du) 
    210211       ENDDO        
    211212        
    … …  
    217218          u=f_u(ind) 
    218219          theta_rhodz=f_theta_rhodz(ind) 
    219           rhodz=f_rhodz(ind) 
     220          mass=f_mass(ind) 
    220221          theta = f_theta(ind) 
    221222          qu=f_qu(ind) 
    222           CALL compute_pvort(ps,u,theta_rhodz, rhodz,theta,qu) 
    223           phis=f_phis(ind) 
     223          CALL compute_pvort(ps,u,theta_rhodz, mass,theta,qu) 
    224224          pk = f_pk(ind) 
    225225          geopot = f_geopot(ind)   
    226           CALL compute_geopot(ps,phis,rhodz,theta, pk,geopot) 
     226          CALL compute_geopot(ps,mass,theta, pk,geopot) 
    227227          hflux=f_hflux(ind) 
    228228          divm = f_divm(ind) 
    229229          dtheta_rhodz=f_dtheta_rhodz(ind) 
    230230          du=f_du(ind) 
    231           CALL compute_caldyn_horiz(u,rhodz,qu,theta,pk,geopot, hflux,divm,dtheta_rhodz,du) 
     231          CALL compute_caldyn_horiz(u,mass,qu,theta,pk,geopot, hflux,divm,dtheta_rhodz,du) 
    232232          wflux=f_wflux(ind) 
    233233          wwuu=f_wwuu(ind) 
    234234          dps=f_dps(ind) 
    235           CALL compute_caldyn_vert(u,theta,rhodz,divm, wflux,wwuu, dps, dtheta_rhodz, du) 
     235          CALL compute_caldyn_vert(u,theta,mass,divm, wflux,wwuu, dps, dtheta_rhodz, du) 
    236236       ENDDO 
    237237        
    … …  
    248248!       CALL write_output_fields(f_ps, f_phis, f_dps, f_u, f_theta_rhodz, f_q, & 
    249249!            f_buf_i, f_buf_v, f_buf_u3d, f_buf_ulon, f_buf_ulat, f_buf_s, f_buf_p) 
     250 
    250251       CALL writefield("ps",f_ps) 
    251252       CALL writefield("dps",f_dps) 
     253       CALL writefield("mass",f_mass) 
    252254       CALL writefield("vort",f_qu) 
    253255       CALL writefield("theta",f_theta) 
    … …  
    344346  END SUBROUTINE compute_pvort 
    345347   
    346   SUBROUTINE compute_geopot(ps,phis,rhodz,theta,pk,geopot)  
     348  SUBROUTINE compute_geopot(ps,rhodz,theta,pk,geopot)  
    347349  USE icosa 
    348350  USE disvert_mod 
    … …  
    351353  USE omp_para 
    352354  IMPLICIT NONE 
    353     REAL(rstd),INTENT(IN)  :: phis(iim*jjm) 
    354     REAL(rstd),INTENT(IN)  :: ps(iim*jjm) 
    355     REAL(rstd),INTENT(IN)  :: rhodz(iim*jjm,llm) 
    356     REAL(rstd),INTENT(IN)  :: theta(iim*jjm,llm)  ! potential temperature 
    357     REAL(rstd),INTENT(INOUT) :: pk(iim*jjm,llm)   ! Exner function 
    358     REAL(rstd),INTENT(INOUT) :: geopot(iim*jjm,llm+1)    ! geopotential 
     355    REAL(rstd),INTENT(INOUT) :: ps(iim*jjm) 
     356    REAL(rstd),INTENT(IN)    :: rhodz(iim*jjm,llm) 
     357    REAL(rstd),INTENT(IN)    :: theta(iim*jjm,llm)    ! potential temperature 
     358    REAL(rstd),INTENT(INOUT) :: pk(iim*jjm,llm)       ! Exner function 
     359    REAL(rstd),INTENT(INOUT) :: geopot(iim*jjm,llm+1) ! geopotential 
    359360     
    360361    INTEGER :: i,j,ij,l 
    … …  
    363364    CALL trace_start("compute_geopot") 
    364365 
     366    IF(caldyn_eta==eta_mass) THEN 
     367     
    365368!!! Compute exner function and geopotential 
    366     DO l = 1,llm 
    367 !$OMP DO SCHEDULE(STATIC)  
     369       DO l = 1,llm 
     370          !$OMP DO SCHEDULE(STATIC)  
     371          DO j=jj_begin-1,jj_end+1 
     372             DO i=ii_begin-1,ii_end+1 
     373                ij=(j-1)*iim+i 
     374                p_ik = ptop + mass_ak(l) + mass_bk(l)*ps(ij) ! FIXME : leave ps for the moment ; change ps to Ms later 
     375                !         p_ik = ptop + g*(mass_ak(l)+ mass_bk(l)*ps(i,j)) 
     376                exner_ik = cpp * (p_ik/preff) ** kappa 
     377                pk(ij,l) = exner_ik 
     378                ! specific volume v = kappa*theta*pi/p = dphi/g/rhodz 
     379                geopot(ij,l+1) = geopot(ij,l) + (g*kappa)*rhodz(ij,l)*theta(ij,l)*exner_ik/p_ik  
     380             ENDDO 
     381          ENDDO 
     382       ENDDO 
     383 
     384    ELSE 
     385 
     386       ! We are using a Lagrangian vertical coordinate 
     387       ! Pressure must be computed first top-down (temporarily stored in pk) 
     388       ! Then Exner pressure and geopotential are computed bottom-up 
     389       ! Notice that the computation below should work also when caldyn_eta=eta_mass 
     390 
     391       ! uppermost layer 
    368392       DO j=jj_begin-1,jj_end+1 
    369393          DO i=ii_begin-1,ii_end+1 
    370394             ij=(j-1)*iim+i 
    371              p_ik = ptop + mass_ak(l) + mass_bk(l)*ps(ij) ! FIXME : leave ps for the moment ; change ps to Ms later 
    372              !         p_ik = ptop + g*(mass_ak(l)+ mass_bk(l)*ps(i,j)) 
    373              exner_ik = cpp * (p_ik/preff) ** kappa 
    374              pk(ij,l) = exner_ik 
    375              ! specific volume v = kappa*theta*pi/p = dphi/g/rhodz 
    376              geopot(ij,l+1) = geopot(ij,l) + (g*kappa)*rhodz(ij,l)*theta(ij,l)*exner_ik/p_ik  
     395             pk(ij,llm) = ptop + (.5*g)*rhodz(ij,llm) 
     396          END DO 
     397       END DO 
     398       ! other layers 
     399       DO l = llm-1, 1, -1 
     400          !$OMP DO SCHEDULE(STATIC)  
     401          DO j=jj_begin-1,jj_end+1 
     402             DO i=ii_begin-1,ii_end+1 
     403                ij=(j-1)*iim+i 
     404                pk(ij,l) = pk(ij,l+1) + (.5*g)*(rhodz(ij,l)+rhodz(ij,l+1)) 
     405             END DO 
     406          END DO 
     407       END DO 
     408       ! surface pressure (for diagnostics) 
     409       DO j=jj_begin-1,jj_end+1 
     410          DO i=ii_begin-1,ii_end+1 
     411             ij=(j-1)*iim+i 
     412             ps(ij) = pk(ij,1) + (.5*g)*rhodz(ij,1) 
     413          END DO 
     414       END DO 
     415 
     416       DO l = 1,llm 
     417          !$OMP DO SCHEDULE(STATIC)  
     418          DO j=jj_begin-1,jj_end+1 
     419             DO i=ii_begin-1,ii_end+1 
     420                ij=(j-1)*iim+i 
     421                p_ik = pk(ij,l) 
     422                exner_ik = cpp * (p_ik/preff) ** kappa 
     423                ! specific volume v = kappa*theta*pi/p = dphi/g/rhodz 
     424                geopot(ij,l+1) = geopot(ij,l) + (g*kappa)*rhodz(ij,l)*theta(ij,l)*exner_ik/p_ik  
     425                pk(ij,l) = exner_ik 
     426             ENDDO 
    377427          ENDDO 
    378428       ENDDO 
    379     ENDDO 
     429 
     430    END IF 
    380431 
    381432  CALL trace_end("compute_geopot") 
    … …  
    798849     
    799850    ! Temperature 
    800     CALL theta_rhodz2temperature(f_ps,f_theta_rhodz,f_buf_i) ; 
    801      
     851!    CALL theta_rhodz2temperature(f_ps,f_theta_rhodz,f_buf_i) ; ! FIXME 
     852      
    802853    CALL getin('physics',physics_type) 
    803854    IF (TRIM(physics_type)=='dcmip') THEN 
    … …  
    828879    ENDIF 
    829880     
    830     ! geopotential 
     881    ! geopotential ! FIXME 
    831882    CALL thetarhodz2geopot(f_ps,f_phis,f_theta_rhodz, f_buf_s,f_buf_p,f_buf1_i,f_buf2_i,f_buf_i) 
    832883    CALL writefield("p",f_buf_p) 
    833     CALL writefield("phi",f_buf_i) 
     884    CALL writefield("phi",f_geopot)   ! geopotential 
    834885    CALL writefield("theta",f_buf1_i) ! potential temperature 
    835     CALL writefield("pk",f_buf2_i) ! Exner pressure 
    836    
     886    CALL writefield("pk",f_buf2_i)    ! Exner pressure 
    837887   
    838888  END SUBROUTINE write_output_fields 

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

    r156 r159  
    55  REAL(rstd), SAVE, POINTER :: presnivs(:) 
    66  REAL(rstd), SAVE, POINTER :: mass_al(:), mass_bl(:), mass_ak(:), mass_bk(:), mass_dak(:), mass_dbk(:) 
     7 
    78  REAL(rstd) :: ptop ! pressure at top of atmosphere l=llm+1 
     9  ! vertical coordinate : Lagrangian or mass-based 
     10  INTEGER, SAVE :: caldyn_eta 
     11  INTEGER, PARAMETER :: eta_mass=1, eta_lag=2 
    812 
    913CONTAINS 
    … …  
    9195  END SUBROUTINE init_disvert   
    9296   
     97  SUBROUTINE compute_rhodz(comp, ps, rhodz) 
     98    USE icosa 
     99    USE omp_para 
     100    LOGICAL, INTENT(IN) :: comp ! .TRUE. to compute, .FALSE. to check 
     101    REAL(rstd), INTENT(IN) :: ps(iim*jjm) 
     102    REAL(rstd), INTENT(INOUT) :: rhodz(iim*jjm,llm) 
     103    REAL(rstd) :: m, err 
     104    INTEGER :: l,i,j,ij,dd 
     105    err=0. 
     106 
     107    dd=0 
     108 
     109    DO l = ll_begin, ll_end 
     110       DO j=jj_begin-dd,jj_end+dd 
     111          DO i=ii_begin-dd,ii_end+dd 
     112             ij=(j-1)*iim+i 
     113             m = ( ap(l) - ap(l+1) + (bp(l)-bp(l+1))*ps(ij) )/g  
     114             IF(comp) THEN 
     115                rhodz(ij,l) = m 
     116             ELSE 
     117                err = MAX(err,abs(m-rhodz(ij,l))) 
     118             END IF 
     119          ENDDO 
     120       ENDDO 
     121    ENDDO 
     122 
     123    IF(.NOT. comp) THEN 
     124       IF(err>1e-10) THEN 
     125          PRINT *, 'Discrepancy between ps and rhodz detected', err 
     126          STOP 
     127       ELSE 
     128!          PRINT *, 'No discrepancy between ps and rhodz detected' 
     129       END IF 
     130    END IF 
     131 
     132  END SUBROUTINE compute_rhodz 
     133 
     134  
    93135  SUBROUTINE write_apbp 
    94136  USE icosa 

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

    r149 r159  
    44CONTAINS 
    55   
    6   SUBROUTINE etat0(f_ps,f_phis,f_theta_rhodz,f_u, f_q) 
     6  SUBROUTINE etat0(f_ps,f_mass,f_phis,f_theta_rhodz,f_u, f_q) 
    77    USE icosa 
     8    USE disvert_mod 
    89    USE etat0_jablonowsky06_mod, ONLY : etat0_jablonowsky06=>etat0 
    910    USE etat0_academic_mod, ONLY : etat0_academic=>etat0   
    … …  
    1920    IMPLICIT NONE 
    2021    TYPE(t_field),POINTER :: f_ps(:) 
     22    TYPE(t_field),POINTER :: f_mass(:) 
    2123    TYPE(t_field),POINTER :: f_phis(:) 
    2224    TYPE(t_field),POINTER :: f_theta_rhodz(:) 
    2325    TYPE(t_field),POINTER :: f_u(:) 
    2426    TYPE(t_field),POINTER :: f_q(:) 
    25      
     27    REAL(rstd),POINTER :: ps(:), mass(:,:) 
     28    LOGICAL :: init_mass 
     29    INTEGER :: ind,i,j,ij,l 
     30 
     31    ! most etat0 routines set ps and not mass 
     32    ! in that case and if caldyn_eta == eta_lag 
     33    ! the initial distribution of mass is taken to be the same 
     34    ! as what the mass coordinate would dictate 
     35    ! however if etat0_XXX defines mass then the flag init_mass must be set to .FALSE. 
     36    ! otherwise mass will be overwritten  
     37    init_mass = (caldyn_eta == eta_lag) 
     38 
    2639    etat0_type='jablonowsky06' 
    2740    CALL getin("etat0",etat0_type) 
    … …  
    5669       STOP 
    5770    END SELECT 
    58      
     71 
     72    IF(init_mass) THEN ! initialize mass distribution using ps 
     73       !$OMP BARRIER 
     74       DO ind=1,ndomain 
     75          CALL swap_dimensions(ind) 
     76          CALL swap_geometry(ind) 
     77          mass=f_mass(ind); ps=f_ps(ind) 
     78          CALL compute_rhodz(.TRUE., ps, mass) 
     79       END DO 
     80    END IF 
     81 
    5982  END SUBROUTINE etat0 
    6083          
    6184END MODULE etat0_mod 
     85 

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

    r138 r159  
    6363         field(ind)%name = name 
    6464      ELSE 
    65          field(ind)%name = '(unkown)' 
     65         field(ind)%name = '(undefined)' 
    6666      END IF 
    6767 
    … …  
    232232     
    233233    IF (field%ndim/=2 .OR. field%data_type/=type_real) THEN 
    234        PRINT *, 'get_val_r2d : bad pointer assignation with ' // TRIM(field%name)  
     234       PRINT *, 'get_val_r2d : bad pointer assignment with ' // TRIM(field%name)  
    235235       STOP 
    236236    END IF 
    … …  
    244244     
    245245    IF (field%ndim/=3 .OR. field%data_type/=type_real) THEN 
    246        PRINT *, 'get_val_r3d : bad pointer assignation with ' // TRIM(field%name)  
     246       PRINT *, 'get_val_r3d : bad pointer assignment with ' // TRIM(field%name)  
    247247       STOP 
     248!       CALL ABORT 
    248249    END IF 
    249250    field_pt=>field%rval3d 
    … …  
    256257     
    257258    IF (field%ndim/=4 .OR. field%data_type/=type_real) THEN 
    258        PRINT *, 'get_val_r4d : bad pointer assignation with ' // TRIM(field%name) 
     259       PRINT *, 'get_val_r4d : bad pointer assignment with ' // TRIM(field%name) 
    259260       STOP 
    260261    END IF 
    … …  
    268269    TYPE(t_field),INTENT(IN) :: field 
    269270     
    270     IF (field%ndim/=2 .OR. field%data_type/=type_integer) STOP 'get_val_i2d : bad pointer assignation'         
     271    IF (field%ndim/=2 .OR. field%data_type/=type_integer) STOP 'get_val_i2d : bad pointer assignment'         
    271272    field_pt=>field%ival2d 
    272273  END SUBROUTINE  getval_i2d 
    … …  
    277278    TYPE(t_field),INTENT(IN) :: field 
    278279     
    279     IF (field%ndim/=3 .OR. field%data_type/=type_integer) STOP 'get_val_i3d : bad pointer assignation'         
     280    IF (field%ndim/=3 .OR. field%data_type/=type_integer) STOP 'get_val_i3d : bad pointer assignment'         
    280281    field_pt=>field%ival3d 
    281282  END SUBROUTINE  getval_i3d 
    … …  
    286287    TYPE(t_field),INTENT(IN) :: field 
    287288     
    288     IF (field%ndim/=4 .OR. field%data_type/=type_integer) STOP 'get_val_i4d : bad pointer assignation'         
     289    IF (field%ndim/=4 .OR. field%data_type/=type_integer) STOP 'get_val_i4d : bad pointer assignment'         
    289290    field_pt=>field%ival4d 
    290291  END SUBROUTINE  getval_i4d 
    … …  
    295296    TYPE(t_field),INTENT(IN) :: field 
    296297     
    297     IF (field%ndim/=2 .OR. field%data_type/=type_logical) STOP 'get_val_l2d : bad pointer assignation'         
     298    IF (field%ndim/=2 .OR. field%data_type/=type_logical) STOP 'get_val_l2d : bad pointer assignment'         
    298299    field_pt=>field%lval2d 
    299300  END SUBROUTINE  getval_l2d 
    … …  
    304305    TYPE(t_field),INTENT(IN) :: field 
    305306     
    306     IF (field%ndim/=3 .OR. field%data_type/=type_logical) STOP 'get_val_l3d : bad pointer assignation'         
     307    IF (field%ndim/=3 .OR. field%data_type/=type_logical) STOP 'get_val_l3d : bad pointer assignment'         
    307308    field_pt=>field%lval3d 
    308309  END SUBROUTINE  getval_l3d 
    … …  
    313314    TYPE(t_field),INTENT(IN) :: field 
    314315     
    315     IF (field%ndim/=4 .OR. field%data_type/=type_logical) STOP 'get_val_l4d : bad pointer assignation'         
     316    IF (field%ndim/=4 .OR. field%data_type/=type_logical) STOP 'get_val_l4d : bad pointer assignment'         
    316317    field_pt=>field%lval4d 
    317318  END SUBROUTINE  getval_l4d     

  • 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.