Changeset 157


Ignore:
Timestamp:
06/04/13 11:29:14 (11 years ago)
Author:
dubos
Message:

caldyn cleanup, tested with JW06 & MPI (no OpenMP)

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

Legend:

Unmodified
Added
Removed

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

    r151 r157  
    11MODULE caldyn_mod 
    22  USE icosa 
     3  USE caldyn_gcm_mod, ONLY : caldyn_BC 
    34  PRIVATE 
    45  CHARACTER(LEN=255),SAVE :: caldyn_type 
    56   
    6   PUBLIC init_caldyn, caldyn 
     7  PUBLIC init_caldyn, caldyn, caldyn_BC 
    78   
    89CONTAINS 

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

    r156 r157  
    55 
    66  INTEGER, PARAMETER :: energy=1, enstrophy=2 
    7   TYPE(t_field),POINTER :: f_out_u(:), f_p(:), f_rhodz(:), f_qu(:) 
     7  TYPE(t_field),POINTER :: f_out_u(:), f_rhodz(:), f_qu(:) 
    88  REAL(rstd),POINTER :: out_u(:,:), p(:,:), rhodz(:,:), qu(:,:) 
    99 
    … …  
    1414  TYPE(t_field),POINTER :: f_theta(:) 
    1515  TYPE(t_field),POINTER :: f_pk(:) 
    16 !  TYPE(t_field),POINTER :: f_pks(:) 
    17   TYPE(t_field),POINTER :: f_phi(:) 
    1816  TYPE(t_field),POINTER :: f_geopot(:) 
    1917  TYPE(t_field),POINTER :: f_divm(:) 
    … …  
    2422  TYPE(t_message) :: req_ps, req_theta_rhodz, req_u, req_qu 
    2523 
    26   PUBLIC init_caldyn, caldyn, write_output_fields,req_ps 
     24  PUBLIC init_caldyn, caldyn_BC, caldyn, write_output_fields,req_ps 
    2725 
    2826CONTAINS 
    … …  
    4846    IF (is_mpi_root) PRINT *, 'caldyn_conserv=',def 
    4947 
    50     def='direct' 
    51     CALL getin('caldyn_exner',def) 
    52     SELECT CASE(TRIM(def)) 
    53     CASE('lmdz') 
    54        caldyn_exner=lmdz 
    55     CASE('direct') 
    56        caldyn_exner=direct 
    57     CASE DEFAULT 
    58        IF (is_mpi_root) PRINT*,'Bad selector for variable caldyn_exner : <', TRIM(def),'> options are <lmdz>, <direct>' 
    59        STOP 
    60     END SELECT 
    61  
    62     def='direct' 
    63     CALL getin('caldyn_hydrostat',def) 
    64     SELECT CASE(TRIM(def)) 
    65     CASE('lmdz') 
    66        caldyn_hydrostat=lmdz 
    67     CASE('direct') 
    68        caldyn_hydrostat=direct 
    69     CASE DEFAULT 
    70        IF (is_mpi_root) PRINT*,'Bad selector for variable caldyn_hydrostat : <', TRIM(def),'> options are <lmdz>, <direct>' 
    71        STOP 
    72     END SELECT 
    73      
    7448    CALL allocate_caldyn 
    7549   
    … …  
    8155 
    8256    CALL allocate_field(f_out_u,field_u,type_real,llm)  
    83     CALL allocate_field(f_p,field_t,type_real,llm+1) 
    8457    CALL allocate_field(f_rhodz,field_t,type_real,llm)  
    8558    CALL allocate_field(f_qu,field_u,type_real,llm)  
    … …  
    9366    CALL allocate_field(f_buf_s,field_t,type_real) 
    9467 
    95     CALL allocate_field(f_theta,field_t,type_real,llm)  ! potential temperature 
     68    CALL allocate_field(f_theta,field_t,type_real,llm)     ! potential temperature 
    9669    CALL allocate_field(f_pk,field_t,type_real,llm) 
    97 !    CALL allocate_field(f_pks,field_t,type_real) ! Exner function 
    98 !    CALL allocate_field(f_phi,field_t,type_real,llm)     ! geopotential 
    99     CALL allocate_field(f_geopot,field_t,type_real,llm+1)  ! geopotential (new) 
    100     CALL allocate_field(f_divm,field_t,type_real,llm)  ! mass flux divergence 
     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 
    10172    CALL allocate_field(f_wwuu,field_u,type_real,llm+1)  
    10273 
    10374  END SUBROUTINE allocate_caldyn 
     75 
     76  SUBROUTINE caldyn_BC(f_phis, f_wflux) 
     77    USE icosa 
     78    USE mpipara 
     79    USE omp_para 
     80    IMPLICIT NONE 
     81    TYPE(t_field),POINTER :: f_phis(:) 
     82    TYPE(t_field),POINTER :: f_wflux(:) 
     83    REAL(rstd),POINTER  :: phis(:) 
     84    REAL(rstd),POINTER  :: wflux(:,:) 
     85    REAL(rstd),POINTER  :: geopot(:,:) 
     86    REAL(rstd),POINTER  :: wwuu(:,:) 
     87 
     88    INTEGER :: ind,i,j,ij,l 
     89 
     90    IF (omp_first) THEN 
     91       DO ind=1,ndomain 
     92          CALL swap_dimensions(ind) 
     93          CALL swap_geometry(ind) 
     94          geopot=f_geopot(ind) 
     95          phis=f_phis(ind) 
     96          wflux=f_wflux(ind) 
     97          wwuu=f_wwuu(ind) 
     98           
     99          DO j=jj_begin-1,jj_end+1 
     100             DO i=ii_begin-1,ii_end+1 
     101                ij=(j-1)*iim+i 
     102                ! lower BCs : geopot=phis, wflux=0, wwuu=0 
     103                geopot(ij,1) = phis(ij) 
     104                wflux(ij,1) = 0. 
     105                wwuu(ij+u_right,1)=0    
     106                wwuu(ij+u_lup,1)=0    
     107                wwuu(ij+u_ldown,1)=0 
     108                ! top BCs : wflux=0, wwuu=0 
     109                wflux(ij,llm+1)  = 0. 
     110                wwuu(ij+u_right,llm+1)=0    
     111                wwuu(ij+u_lup,llm+1)=0    
     112                wwuu(ij+u_ldown,llm+1)=0 
     113             ENDDO 
     114          ENDDO 
     115       END DO 
     116    ENDIF 
     117 
     118    !$OMP BARRIER 
     119  END SUBROUTINE caldyn_BC 
    104120    
    105121  SUBROUTINE caldyn(write_out,f_phis, f_ps, f_theta_rhodz, f_u, f_q, & 
    … …  
    127143    REAL(rstd),POINTER :: theta_rhodz(:,:), dtheta_rhodz(:,:) 
    128144    REAL(rstd),POINTER :: u(:,:), du(:,:), hflux(:,:), wflux(:,:) 
    129     REAL(rstd),POINTER :: p(:,:), rhodz(:,:), qu(:,:) 
     145    REAL(rstd),POINTER :: rhodz(:,:), qu(:,:) 
    130146 
    131147! temporary shared variable 
    132148    REAL(rstd),POINTER  :: theta(:,:)   
    133     REAL(rstd),POINTER  :: pk(:,:) !, pks(:) 
    134     REAL(rstd),POINTER  :: phi(:,:)     
     149    REAL(rstd),POINTER  :: pk(:,:) 
    135150    REAL(rstd),POINTER  :: geopot(:,:) 
    136151    REAL(rstd),POINTER  :: divm(:,:)  
    137152    REAL(rstd),POINTER  :: wwuu(:,:) 
    138      
    139      
    140     INTEGER :: ind,ij 
     153         
     154    INTEGER :: ind 
    141155    LOGICAL,SAVE :: first=.TRUE. 
    142156!$OMP THREADPRIVATE(first) 
    143  
    144157     
    145158    IF (first) THEN 
    … …  
    163176          CALL swap_geometry(ind) 
    164177          ps=f_ps(ind) 
     178          u=f_u(ind) 
     179          theta_rhodz = f_theta_rhodz(ind) 
    165180          rhodz=f_rhodz(ind) 
    166           p=f_p(ind) 
     181          theta = f_theta(ind) 
    167182          qu=f_qu(ind) 
    168           u=f_u(ind) 
    169  
    170           CALL compute_pvort(ps, u, p,rhodz,qu) 
    171  
     183          CALL compute_pvort(ps,u,theta_rhodz, rhodz,theta,qu) 
    172184       ENDDO 
    173185 
    … …  
    177189          CALL swap_dimensions(ind) 
    178190          CALL swap_geometry(ind) 
     191          ps=f_ps(ind) 
     192          u=f_u(ind) 
     193          theta_rhodz=f_theta_rhodz(ind) 
     194          rhodz=f_rhodz(ind) 
     195          theta = f_theta(ind) 
     196          qu=f_qu(ind) 
    179197          phis=f_phis(ind) 
     198          pk = f_pk(ind) 
     199          geopot = f_geopot(ind)   
     200          CALL compute_geopot(ps,phis,rhodz,theta, pk,geopot) 
    180201          hflux=f_hflux(ind) 
     202          divm = f_divm(ind) 
     203          dtheta_rhodz=f_dtheta_rhodz(ind) 
     204          du=f_du(ind) 
     205          CALL compute_caldyn_horiz(u,rhodz,qu,theta,pk,geopot, hflux,divm,dtheta_rhodz,du) 
    181206          wflux=f_wflux(ind) 
    182           ps=f_ps(ind) 
     207          wwuu=f_wwuu(ind) 
    183208          dps=f_dps(ind) 
    184           theta_rhodz=f_theta_rhodz(ind) 
    185           dtheta_rhodz=f_dtheta_rhodz(ind) 
    186           rhodz=f_rhodz(ind) 
    187           p=f_p(ind) 
    188           qu=f_qu(ind) 
    189           u=f_u(ind) 
    190           du=f_du(ind) 
    191           out_u=f_out_u(ind) 
    192           theta = f_theta(ind) 
    193           pk = f_pk(ind) 
    194 !          pks = f_pks(ind) 
    195 !          phi = f_phi(ind) 
    196           geopot = f_geopot(ind)   
    197           divm = f_divm(ind) 
    198           wwuu=f_wwuu(ind) 
    199  
    200           CALL compute_caldyn(ps, u, p,rhodz,qu, phis, theta_rhodz, hflux, wflux, dps, dtheta_rhodz, du, & 
    201                               theta,pk, geopot, divm, wwuu) 
    202 !                              theta,pk, pks, phi, geopot, divm, wwuu) 
     209          CALL compute_caldyn_vert(u,theta,rhodz,divm, wflux,wwuu, dps, dtheta_rhodz, du) 
    203210       ENDDO        
    204211        
    … …  
    207214          CALL swap_dimensions(ind) 
    208215          CALL swap_geometry(ind) 
     216          ps=f_ps(ind) 
     217          u=f_u(ind) 
     218          theta_rhodz=f_theta_rhodz(ind) 
     219          rhodz=f_rhodz(ind) 
     220          theta = f_theta(ind) 
     221          qu=f_qu(ind) 
     222          CALL compute_pvort(ps,u,theta_rhodz, rhodz,theta,qu) 
    209223          phis=f_phis(ind) 
    210           phi=f_phi(ind) 
     224          pk = f_pk(ind) 
    211225          geopot = f_geopot(ind)   
    212           ps=f_ps(ind) 
     226          CALL compute_geopot(ps,phis,rhodz,theta, pk,geopot) 
     227          hflux=f_hflux(ind) 
     228          divm = f_divm(ind) 
     229          dtheta_rhodz=f_dtheta_rhodz(ind) 
     230          du=f_du(ind) 
     231          CALL compute_caldyn_horiz(u,rhodz,qu,theta,pk,geopot, hflux,divm,dtheta_rhodz,du) 
     232          wflux=f_wflux(ind) 
     233          wwuu=f_wwuu(ind) 
    213234          dps=f_dps(ind) 
    214           hflux=f_hflux(ind) 
    215           wflux=f_wflux(ind) 
    216           theta_rhodz=f_theta_rhodz(ind) 
    217           dtheta_rhodz=f_dtheta_rhodz(ind) 
    218           rhodz=f_rhodz(ind) 
    219           p=f_p(ind) 
    220           qu=f_qu(ind) 
    221           u=f_u(ind) 
    222           du=f_du(ind) 
    223           out_u=f_out_u(ind)  
    224           wwuu=f_wwuu(ind) 
    225           CALL compute_pvort(ps, u, p,rhodz,qu) 
    226           CALL compute_caldyn(ps, u, p,rhodz,qu, phis, theta_rhodz, hflux, wflux, dps, dtheta_rhodz, du,   & 
    227                               theta,pk, phi, divm, wwuu) 
    228 !                              theta,pk, pks, phi, geopot, divm, wwuu) 
     235          CALL compute_caldyn_vert(u,theta,rhodz,divm, wflux,wwuu, dps, dtheta_rhodz, du) 
    229236       ENDDO 
    230237        
    … …  
    254261END SUBROUTINE caldyn 
    255262     
    256 SUBROUTINE compute_pvort(ps, u, p,rhodz,qu) 
     263SUBROUTINE compute_pvort(ps,u,theta_rhodz, rhodz,theta,qu) 
    257264  USE icosa 
    258265  USE disvert_mod 
    … …  
    263270  REAL(rstd),INTENT(IN)  :: u(iim*3*jjm,llm) 
    264271  REAL(rstd),INTENT(IN)  :: ps(iim*jjm) 
    265   REAL(rstd),INTENT(OUT) :: p(iim*jjm,llm+1) 
     272  REAL(rstd),INTENT(IN)  :: theta_rhodz(iim*jjm,llm) 
    266273  REAL(rstd),INTENT(OUT) :: rhodz(iim*jjm,llm) 
     274  REAL(rstd),INTENT(OUT) :: theta(iim*jjm,llm) 
    267275  REAL(rstd),INTENT(OUT) :: qu(iim*3*jjm,llm) 
    268276   
    … …  
    274282  CALL trace_start("compute_pvort")   
    275283 
    276   CALL wait_message(req_ps)    
    277 !!! Compute pressure 
    278    DO    l    = ll_begin, ll_endp1 
    279     CALL test_message(req_u)  
    280  
    281      DO j=jj_begin-1,jj_end+1 
    282         DO i=ii_begin-1,ii_end+1 
    283            ij=(j-1)*iim+i 
    284            p(ij,l) = ap(l) + bp(l) * ps(ij) 
    285         ENDDO 
    286      ENDDO 
    287   ENDDO 
    288  
    289 !$OMP BARRIER   
    290  
    291 !!! Compute mass 
     284  CALL wait_message(req_ps)   
     285  CALL wait_message(req_theta_rhodz)  
     286 
     287!!! Compute mass & theta 
    292288   DO l = ll_begin,ll_end 
    293289     CALL test_message(req_u)  
    … …  
    295291        DO i=ii_begin-1,ii_end+1 
    296292           ij=(j-1)*iim+i 
    297            rhodz(ij,l) = ( p(ij,l) - p(ij,l+1) )/g 
     293           rhodz(ij,l) = ( mass_dak(l)+ps(ij)*mass_dbk(l) )/g 
     294           theta(ij,l) = theta_rhodz(ij,l)/rhodz(ij,l) 
    298295        ENDDO 
    299296     ENDDO 
    … …  
    304301!!! Compute shallow-water potential vorticity 
    305302  DO l = ll_begin,ll_end 
    306     CALL test_message(req_theta_rhodz)  
    307303 
    308304    DO j=jj_begin-1,jj_end+1 
    … …  
    347343                                                        
    348344  END SUBROUTINE compute_pvort 
    349  
    350    
    351   SUBROUTINE compute_caldyn(ps, u, p,rhodz,qu, phis, theta_rhodz, hflux, wflux, dps, dtheta_rhodz, du,   & 
    352                             theta,pk, geopot, divm, wwuu)  
    353 !                            theta,pk, pks, phi, geopot, divm, wwuu)  
     345   
     346  SUBROUTINE compute_geopot(ps,phis,rhodz,theta,pk,geopot)  
    354347  USE icosa 
    355348  USE disvert_mod 
    … …  
    359352  IMPLICIT NONE 
    360353    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 
     359     
     360    INTEGER :: i,j,ij,l 
     361    REAL(rstd) :: p_ik, exner_ik 
     362 
     363    CALL trace_start("compute_geopot") 
     364 
     365!!! Compute exner function and geopotential 
     366    DO l = 1,llm 
     367!$OMP DO SCHEDULE(STATIC)  
     368       DO j=jj_begin-1,jj_end+1 
     369          DO i=ii_begin-1,ii_end+1 
     370             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  
     377          ENDDO 
     378       ENDDO 
     379    ENDDO 
     380 
     381  CALL trace_end("compute_geopot") 
     382 
     383  END SUBROUTINE compute_geopot 
     384 
     385  SUBROUTINE compute_caldyn_horiz(u,rhodz,qu,theta,pk,geopot, hflux,divm, dtheta_rhodz, du) 
     386  USE icosa 
     387  USE disvert_mod 
     388  USE exner_mod 
     389  USE trace 
     390  USE omp_para 
     391  IMPLICIT NONE 
    361392    REAL(rstd),INTENT(IN)  :: u(iim*3*jjm,llm) 
    362     REAL(rstd),INTENT(IN)  :: theta_rhodz(iim*jjm,llm) 
    363     REAL(rstd),INTENT(IN)  :: ps(iim*jjm) 
    364     REAL(rstd),INTENT(IN)  :: p(iim*jjm,llm+1) 
    365393    REAL(rstd),INTENT(IN)  :: rhodz(iim*jjm,llm) 
    366394    REAL(rstd),INTENT(IN)  :: qu(iim*3*jjm,llm) 
    367  
    368     REAL(rstd),INTENT(OUT) :: du(iim*3*jjm,llm), hflux(iim*3*jjm,llm) ! hflux in kg/s 
     395    REAL(rstd),INTENT(IN)  :: theta(iim*jjm,llm)  ! potential temperature 
     396    REAL(rstd),INTENT(IN)  :: pk(iim*jjm,llm) ! Exner function 
     397    REAL(rstd),INTENT(IN)  :: geopot(iim*jjm,llm+1)    ! geopotential 
     398 
     399    REAL(rstd),INTENT(OUT) :: hflux(iim*3*jjm,llm) ! hflux in kg/s 
     400    REAL(rstd),INTENT(OUT) :: divm(iim*jjm,llm)  ! mass flux divergence 
    369401    REAL(rstd),INTENT(OUT) :: dtheta_rhodz(iim*jjm,llm) 
    370     REAL(rstd),INTENT(OUT) :: dps(iim*jjm) 
    371     REAL(rstd),INTENT(OUT) :: wflux(iim*jjm,llm+1) ! vertical mass flux (kg/m2/s) 
    372  
    373 ! temporary variable 
    374     REAL(rstd),INTENT(INOUT) :: theta(iim*jjm,llm)  ! potential temperature 
    375     REAL(rstd),INTENT(INOUT) :: pk(iim*jjm,llm) !, pks(iim*jjm) ! Exner function 
    376 !    REAL(rstd),INTENT(INOUT) :: phi(iim*jjm,llm)    ! geopotential 
    377     REAL(rstd),INTENT(INOUT) :: geopot(iim*jjm,llm+1)    ! geopotential 
    378     REAL(rstd),INTENT(INOUT) :: divm(iim*jjm,llm)  ! mass flux divergence 
    379     REAL(rstd),INTENT(INOUT) :: wwuu(iim*3*jjm,llm+1) 
    380      
     402    REAL(rstd),INTENT(OUT) :: du(iim*3*jjm,llm) 
     403 
    381404    REAL(rstd) :: Ftheta(3*iim*jjm,llm) ! theta flux 
    382     REAL(rstd) :: berni(iim*jjm,llm)  ! Bernouilli function 
     405    REAL(rstd) :: berni(iim*jjm,llm)  ! Bernoulli function 
    383406     
    384407    INTEGER :: i,j,ij,l 
    385     REAL(rstd) :: ww,uu, p_ik, exner_ik 
    386  
    387     CALL trace_start("compute_caldyn") 
     408    REAL(rstd) :: ww,uu 
     409 
     410    CALL trace_start("compute_caldyn_horiz") 
    388411 
    389412    CALL wait_message(req_theta_rhodz)  
    390  
    391 !!! Compute theta 
    392    DO l = ll_begin, ll_end 
    393       IF (caldyn_conserv==energy) CALL test_message(req_qu)  
    394       DO j=jj_begin-1,jj_end+1 
    395          DO i=ii_begin-1,ii_end+1 
    396             ij=(j-1)*iim+i 
    397             theta(ij,l) = theta_rhodz(ij,l)/rhodz(ij,l) 
    398          ENDDO 
    399       ENDDO 
    400    ENDDO 
    401413 
    402414  DO l = ll_begin, ll_end 
    … …  
    438450      ENDDO 
    439451    ENDDO 
    440   ENDDO 
    441  
    442 !$OMP BARRIER    
    443 !!! cumulate mass flux divergence from top to bottom 
    444   DO  l = llm-1, 1, -1 
    445     IF (caldyn_conserv==energy) CALL test_message(req_qu)  
    446 !$OMP DO SCHEDULE(STATIC)  
    447     DO j=jj_begin,jj_end 
    448       DO i=ii_begin,ii_end 
    449         ij=(j-1)*iim+i 
    450         divm(ij,l) = divm(ij,l) + divm(ij,l+1) 
    451       ENDDO 
    452     ENDDO 
    453   ENDDO 
    454  
    455 ! IMPLICIT FLUSH on divm 
    456 !!!!!!!!!!!!!!!!!!!!!!!!! 
    457    
    458 !!! Compute vertical mass flux 
    459 !  DO l = 2,llm 
    460   DO l=ll_beginp1,ll_end 
    461     IF (caldyn_conserv==energy) CALL test_message(req_qu)  
    462     DO j=jj_begin,jj_end 
    463       DO i=ii_begin,ii_end 
    464         ij=(j-1)*iim+i 
    465         ! w = int(z,ztop,div(flux)dz) + B(eta)dps/dt 
    466         ! => w>0 for upward transport 
    467         wflux( ij, l ) = divm( ij, l ) - bp(l) * divm( ij, 1 ) 
    468       ENDDO 
    469     ENDDO 
    470   ENDDO 
    471  
    472 ! compute dps, set vertical mass flux at the surface to 0 
    473   IF (omp_first) THEN 
    474     DO j=jj_begin,jj_end 
    475       DO i=ii_begin,ii_end 
    476         ij=(j-1)*iim+i 
    477         wflux(ij,1)  = 0. 
    478         ! dps/dt = -int(div flux)dz 
    479         dps(ij)=-divm(ij,1) * g  
    480       ENDDO 
    481     ENDDO 
    482   ENDIF 
    483  
    484  
    485   IF (omp_last) THEN 
    486     DO j=jj_begin,jj_end 
    487       DO i=ii_begin,ii_end 
    488         ij=(j-1)*iim+i 
    489         wflux(ij,llm+1)  = 0. 
    490       ENDDO 
    491     ENDDO 
    492   ENDIF 
     452 
     453 END DO 
     454 
    493455!!! Compute potential vorticity (Coriolis) contribution to du 
    494  
     456  
    495457    SELECT CASE(caldyn_conserv) 
    496458    CASE(energy) ! energy-conserving TRiSK 
    … …  
    596558    END SELECT 
    597559   
    598 !!!! Compute Exner function 
    599 !!    CALL compute_exner(ps,p,pks,pk,1)  
    600 !! replaced in source 
    601 !    IF (omp_first) THEN     
    602 !       DO j=jj_begin-1,jj_end+1 
    603 !          DO i=ii_begin-1,ii_end+1 
    604 !             ij=(j-1)*iim+i 
    605 !             pks(ij) = cpp * ( ps(ij)/preff ) ** kappa 
    606 !          ENDDO 
    607 !       ENDDO 
    608 !     ENDIF 
    609 !      
    610 !       ! 3D : pk 
    611 !     DO l = ll_begin, ll_end 
    612 !        DO j=jj_begin-1,jj_end+1 
    613 !           DO i=ii_begin-1,ii_end+1 
    614 !              ij=(j-1)*iim+i 
    615 !              pk(ij,l) = cpp * ((.5/preff)*(p(ij,l)+p(ij,l+1))) ** kappa 
    616 !           ENDDO 
    617 !        ENDDO 
    618 !     ENDDO 
    619 ! 
    620 !!---> flush pk,pks, theta 
    621 !!$OMP BARRIER 
    622 ! 
    623 !!! Compute geopotential (old) 
    624 ! 
    625 !  ! for first layer 
    626 !  IF (omp_first) THEN 
    627 !    DO j=jj_begin-1,jj_end+1 
    628 !      DO i=ii_begin-1,ii_end+1 
    629 !        ij=(j-1)*iim+i 
    630 !        phi( ij,1 ) = phis( ij ) + theta(ij,1) * ( pks(ij) - pk(ij,1) ) 
    631 !      ENDDO 
    632 !    ENDDO 
    633 !  ENDIF 
    634 !!!-> implicit flush on phi(:,1) 
    635 !           
    636 !  ! for other layers 
    637 !   DO l = ll_beginp1, ll_end 
    638 !     DO j=jj_begin-1,jj_end+1 
    639 !       DO i=ii_begin-1,ii_end+1 
    640 !         ij=(j-1)*iim+i 
    641 !         phi(ij,l) =  0.5 * ( theta(ij,l)  + theta(ij,l-1) )  &  
    642 !                          * (  pk(ij,l-1) -  pk(ij,l)    ) 
    643 !       ENDDO 
    644 !     ENDDO 
    645 !   ENDDO 
    646 ! 
    647 !!$OMP BARRIER 
    648 !   DO l = 2, llm 
    649 !!$OMP DO 
    650 !     DO j=jj_begin-1,jj_end+1 
    651 !       DO i=ii_begin-1,ii_end+1 
    652 !         ij=(j-1)*iim+i 
    653 !         phi(ij,l) = phi(ij,l)+ phi(ij,l-1) 
    654 !       ENDDO 
    655 !     ENDDO 
    656 !   ENDDO 
    657 !! --> IMPLICIT FLUSH on phi 
    658        
    659 !!! Compute exner function and geopotential 
    660  
    661 ! geopot=phis for first layer 
    662 !$OMP DO SCHEDULE(STATIC)  
    663     DO j=jj_begin-1,jj_end+1 
    664        DO i=ii_begin-1,ii_end+1 
    665           ij=(j-1)*iim+i 
    666           geopot(ij,1) = phis(ij) 
    667        ENDDO 
    668     ENDDO 
    669 ! for other layers 
    670     DO l = 1,llm 
    671 !$OMP DO SCHEDULE(STATIC)  
    672        DO j=jj_begin-1,jj_end+1 
    673           DO i=ii_begin-1,ii_end+1 
    674              ij=(j-1)*iim+i 
    675              p_ik = ptop + mass_ak(l) + mass_bk(l)*ps(ij) ! FIXME : leave ps for the moment ; change ps to Ms later 
    676              !         p_ik = ptop + g*(mass_ak(l)+ mass_bk(l)*ps(i,j)) 
    677              exner_ik = cpp * (p_ik/preff) ** kappa 
    678              pk(ij,l) = exner_ik 
    679              ! specific volume v = kappa*theta*pi/p = dphi/g/rhodz 
    680              geopot(ij,l+1) = geopot(ij,l) + (g*kappa)*rhodz(ij,l)*theta(ij,l)*exner_ik/p_ik  
    681           ENDDO 
    682        ENDDO 
    683     ENDDO 
    684  
    685560!!! Compute bernouilli term = Kinetic Energy + geopotential 
    686561   DO l=ll_begin,ll_end 
    … …  
    702577 
    703578  
    704 !!! gradients of Bernoulli and Exner functions 
     579!!! Add gradients of Bernoulli and Exner functions to du 
    705580 
    706581   DO l=ll_begin,ll_end 
    … …  
    729604  ENDDO 
    730605 
    731     
     606  CALL trace_end("compute_caldyn_horiz") 
     607 
     608END SUBROUTINE compute_caldyn_horiz 
     609 
     610SUBROUTINE compute_caldyn_vert(u,theta,rhodz,divm, wflux,wwuu, dps,dtheta_rhodz,du) 
     611  USE icosa 
     612  USE disvert_mod 
     613  USE exner_mod 
     614  USE trace 
     615  USE omp_para 
     616  IMPLICIT NONE 
     617    REAL(rstd),INTENT(IN)  :: u(iim*3*jjm,llm) 
     618    REAL(rstd),INTENT(IN)  :: theta(iim*jjm,llm) 
     619    REAL(rstd),INTENT(IN)  :: rhodz(iim*jjm,llm) 
     620    REAL(rstd),INTENT(INOUT)  :: divm(iim*jjm,llm)  ! mass flux divergence 
     621 
     622    REAL(rstd),INTENT(OUT) :: wflux(iim*jjm,llm+1) ! vertical mass flux (kg/m2/s) 
     623    REAL(rstd),INTENT(OUT) :: wwuu(iim*3*jjm,llm+1) 
     624    REAL(rstd),INTENT(OUT) :: du(iim*3*jjm,llm) 
     625    REAL(rstd),INTENT(OUT) :: dtheta_rhodz(iim*jjm,llm) 
     626    REAL(rstd),INTENT(OUT) :: dps(iim*jjm) 
     627 
     628! temporary variable     
     629    INTEGER :: i,j,ij,l 
     630    REAL(rstd) :: p_ik, exner_ik 
     631 
     632  CALL trace_start("compute_caldyn_vert") 
     633 
     634!$OMP BARRIER    
     635!!! cumulate mass flux divergence from top to bottom 
     636  DO  l = llm-1, 1, -1 
     637    IF (caldyn_conserv==energy) CALL test_message(req_qu)  
     638!$OMP DO SCHEDULE(STATIC)  
     639    DO j=jj_begin,jj_end 
     640      DO i=ii_begin,ii_end 
     641        ij=(j-1)*iim+i 
     642        divm(ij,l) = divm(ij,l) + divm(ij,l+1) 
     643      ENDDO 
     644    ENDDO 
     645  ENDDO 
     646 
     647! IMPLICIT FLUSH on divm 
     648!!!!!!!!!!!!!!!!!!!!!!!!! 
     649 
     650! compute dps 
     651  IF (omp_first) THEN 
     652    DO j=jj_begin,jj_end 
     653      DO i=ii_begin,ii_end 
     654        ij=(j-1)*iim+i 
     655        ! dps/dt = -int(div flux)dz 
     656        dps(ij)=-divm(ij,1) * g  
     657      ENDDO 
     658    ENDDO 
     659  ENDIF 
     660   
     661!!! Compute vertical mass flux (l=1,llm+1 done by caldyn_BC) 
     662  DO l=ll_beginp1,ll_end 
     663    IF (caldyn_conserv==energy) CALL test_message(req_qu)  
     664    DO j=jj_begin,jj_end 
     665      DO i=ii_begin,ii_end 
     666        ij=(j-1)*iim+i 
     667        ! w = int(z,ztop,div(flux)dz) + B(eta)dps/dt 
     668        ! => w>0 for upward transport 
     669        wflux( ij, l ) = divm( ij, l ) - bp(l) * divm( ij, 1 ) 
     670      ENDDO 
     671    ENDDO 
     672  ENDDO 
     673 
    732674  DO l=ll_begin,ll_endm1 
    733675    DO j=jj_begin,jj_end 
    … …  
    748690  ENDDO 
    749691   
    750   IF (omp_first) THEN 
    751     DO j=jj_begin,jj_end 
    752       DO i=ii_begin,ii_end 
    753         ij=(j-1)*iim+i 
    754         wwuu(ij+u_right,1)=0    
    755         wwuu(ij+u_lup,1)=0    
    756         wwuu(ij+u_ldown,1)=0 
    757       ENDDO 
    758     ENDDO 
    759   ENDIF    
    760  
    761   IF (omp_last) THEN 
    762     DO j=jj_begin,jj_end 
    763       DO i=ii_begin,ii_end 
    764         ij=(j-1)*iim+i 
    765         wwuu(ij+u_right,llm+1)=0    
    766         wwuu(ij+u_lup,llm+1)=0    
    767         wwuu(ij+u_ldown,llm+1)=0 
    768       ENDDO 
    769     ENDDO 
    770   ENDIF    
    771    
     692! Compute vertical transport 
    772693  DO l=ll_beginp1,ll_end 
    773694    DO j=jj_begin,jj_end 
    … …  
    784705 !$OMP BARRIER 
    785706 
     707! Add vertical transport to du 
    786708  DO l=ll_begin,ll_end 
    787709    DO j=jj_begin,jj_end 
    … …  
    795717  ENDDO       
    796718   
    797   CALL trace_end("compute_caldyn") 
    798  
    799   END SUBROUTINE compute_caldyn 
     719  CALL trace_end("compute_caldyn_vert") 
     720 
     721  END SUBROUTINE compute_caldyn_vert 
    800722 
    801723!-------------------------------- Diagnostics ---------------------------- 

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

    r154 r157  
    199199    LOGICAL :: fluxt_zero(ndomain) ! set to .TRUE. to start accumulating fluxes in time 
    200200    LOGICAL, PARAMETER :: check=.FALSE. 
     201 
     202    CALL caldyn_BC(f_phis, f_wflux) 
    201203 
    202204!$OMP BARRIER 
Note: See TracChangeset for help on using the changeset viewer.