Changeset 162


Ignore:
Timestamp:
06/27/13 18:37:27 (11 years ago)
Author:
dubos
Message:

Lagrangian vertical coordinate tested with test4.1 - 60 MPI procs

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

Legend:

Unmodified
Added
Removed

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

    r159 r162  
    3232   
    3333  SUBROUTINE caldyn(write_out,f_phis, f_ps, f_mass, f_theta_rhodz, f_u, f_q, & 
    34        f_hflux, f_wflux, f_dps, f_dtheta_rhodz, f_du) 
     34       f_hflux, f_wflux, f_dps, f_dmass, f_dtheta_rhodz, f_du) 
    3535  USE icosa 
    3636  USE caldyn_gcm_mod, ONLY : caldyn_gcm=>caldyn 
    … …  
    4747  TYPE(t_field),POINTER :: f_wflux(:) 
    4848  TYPE(t_field),POINTER :: f_dps(:) 
     49  TYPE(t_field),POINTER :: f_dmass(:) 
    4950  TYPE(t_field),POINTER :: f_dtheta_rhodz(:) 
    5051  TYPE(t_field),POINTER :: f_du(:) 
    … …  
    5354      CASE('gcm') 
    5455        CALL caldyn_gcm(write_out,f_phis, f_ps, f_mass, f_theta_rhodz, f_u, f_q, & 
    55              f_hflux, f_wflux, f_dps, f_dtheta_rhodz, f_du) 
     56             f_hflux, f_wflux, f_dps, f_dmass, f_dtheta_rhodz, f_du) 
    5657      CASE('adv') 
    5758        CALL caldyn_adv(write_out,f_phis, f_ps, f_theta_rhodz, f_u, f_q, & 

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

    r159 r162  
    55 
    66  INTEGER, PARAMETER :: energy=1, enstrophy=2 
    7   TYPE(t_field),POINTER :: f_out_u(:), f_qu(:) 
     7  TYPE(t_field),POINTER :: f_out_u(:), f_qu(:), f_qv(:) 
    88  REAL(rstd),POINTER :: out_u(:,:), p(:,:), qu(:,:) 
    99 
    … …  
    1515  TYPE(t_field),POINTER :: f_pk(:) 
    1616  TYPE(t_field),POINTER :: f_geopot(:) 
    17   TYPE(t_field),POINTER :: f_divm(:) 
    1817  TYPE(t_field),POINTER :: f_wwuu(:) 
    1918    
    2019  INTEGER :: caldyn_conserv 
    2120   
    22   TYPE(t_message) :: req_ps, req_theta_rhodz, req_u, req_qu 
    23  
    24   PUBLIC init_caldyn, caldyn_BC, caldyn, write_output_fields,req_ps, & 
    25        caldyn_eta, eta_mass, eta_lag 
     21  TYPE(t_message) :: req_ps, req_mass, req_theta_rhodz, req_u, req_qu 
     22 
     23  PUBLIC init_caldyn, caldyn_BC, caldyn, write_output_fields, & 
     24       req_ps, req_mass, caldyn_eta, eta_mass, eta_lag 
    2625 
    2726CONTAINS 
    … …  
    5857    CALL allocate_field(f_out_u,field_u,type_real,llm)  
    5958    CALL allocate_field(f_qu,field_u,type_real,llm)  
     59    CALL allocate_field(f_qv,field_z,type_real,llm)  
    6060   
    6161    CALL allocate_field(f_buf_i,   field_t,type_real,llm) 
    … …  
    7070    CALL allocate_field(f_pk,    field_t,type_real,llm,  name='pk') 
    7171    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 
    7372    CALL allocate_field(f_wwuu,  field_u,type_real,llm+1,name='wwuu') 
    7473 
    … …  
    121120    
    122121  SUBROUTINE caldyn(write_out,f_phis, f_ps, f_mass, f_theta_rhodz, f_u, f_q, & 
    123        f_hflux, f_wflux, f_dps, f_dtheta_rhodz, f_du) 
     122       f_hflux, f_wflux, f_dps, f_dmass, f_dtheta_rhodz, f_du) 
    124123    USE icosa 
     124    USE disvert_mod 
    125125    USE vorticity_mod 
    126126    USE kinetic_mod 
    … …  
    139139    TYPE(t_field),POINTER :: f_hflux(:), f_wflux(:) 
    140140    TYPE(t_field),POINTER :: f_dps(:) 
     141    TYPE(t_field),POINTER :: f_dmass(:) 
    141142    TYPE(t_field),POINTER :: f_dtheta_rhodz(:) 
    142143    TYPE(t_field),POINTER :: f_du(:) 
    … …  
    146147    REAL(rstd),POINTER :: u(:,:), du(:,:), hflux(:,:), wflux(:,:) 
    147148    REAL(rstd),POINTER :: qu(:,:) 
     149    REAL(rstd),POINTER :: qv(:,:) 
    148150 
    149151! temporary shared variable 
    … …  
    151153    REAL(rstd),POINTER  :: pk(:,:) 
    152154    REAL(rstd),POINTER  :: geopot(:,:) 
    153     REAL(rstd),POINTER  :: divm(:,:)  
     155    REAL(rstd),POINTER  :: convm(:,:)  
    154156    REAL(rstd),POINTER  :: wwuu(:,:) 
    155157         
    … …  
    158160!$OMP THREADPRIVATE(first) 
    159161     
    160     IF (first) THEN 
     162    ! MPI messages need to be sent at first call to caldyn 
     163    ! This is needed only once : the next ones will be sent by timeloop 
     164    IF (first) THEN  
    161165      first=.FALSE. 
    162       CALL init_message(f_ps,req_i1,req_ps) 
     166      IF(caldyn_eta==eta_mass) THEN 
     167         CALL init_message(f_ps,req_i1,req_ps) 
     168      ELSE 
     169         CALL init_message(f_mass,req_i1,req_mass) 
     170      END IF 
    163171      CALL init_message(f_theta_rhodz,req_i1,req_theta_rhodz) 
    164172      CALL init_message(f_u,req_e1_vect,req_u) 
    165173      CALL init_message(f_qu,req_e1_scal,req_qu) 
    166       CALL send_message(f_ps,req_ps)  
     174      IF(caldyn_eta==eta_mass) THEN 
     175         CALL send_message(f_ps,req_ps)  
     176      ELSE 
     177         CALL send_message(f_mass,req_mass)  
     178      END IF 
    167179    ENDIF 
    168180     
    … …  
    183195          theta = f_theta(ind) 
    184196          qu=f_qu(ind) 
    185           CALL compute_pvort(ps,u,theta_rhodz, mass,theta,qu) 
     197          qv=f_qv(ind) 
     198          CALL compute_pvort(ps,u,theta_rhodz, mass,theta,qu,qv) 
    186199       ENDDO 
    187200 
    … …  
    201214          CALL compute_geopot(ps,mass,theta, pk,geopot) 
    202215          hflux=f_hflux(ind) 
    203           divm = f_divm(ind) 
     216          convm = f_dmass(ind) 
    204217          dtheta_rhodz=f_dtheta_rhodz(ind) 
    205218          du=f_du(ind) 
    206           CALL compute_caldyn_horiz(u,mass,qu,theta,pk,geopot, hflux,divm,dtheta_rhodz,du) 
    207           wflux=f_wflux(ind) 
    208           wwuu=f_wwuu(ind) 
    209           dps=f_dps(ind) 
    210           CALL compute_caldyn_vert(u,theta,mass,divm, wflux,wwuu, dps, dtheta_rhodz, du) 
     219          CALL compute_caldyn_horiz(u,mass,qu,theta,pk,geopot, hflux,convm,dtheta_rhodz,du) 
     220          IF(caldyn_eta==eta_mass) THEN 
     221             wflux=f_wflux(ind) 
     222             wwuu=f_wwuu(ind) 
     223             dps=f_dps(ind) 
     224             CALL compute_caldyn_vert(u,theta,mass,convm, wflux,wwuu, dps, dtheta_rhodz, du) 
     225          END IF 
    211226       ENDDO        
    212227        
    … …  
    221236          theta = f_theta(ind) 
    222237          qu=f_qu(ind) 
    223           CALL compute_pvort(ps,u,theta_rhodz, mass,theta,qu) 
     238          CALL compute_pvort(ps,u,theta_rhodz, mass,theta,qu,qv) 
    224239          pk = f_pk(ind) 
    225240          geopot = f_geopot(ind)   
    226241          CALL compute_geopot(ps,mass,theta, pk,geopot) 
    227242          hflux=f_hflux(ind) 
    228           divm = f_divm(ind) 
     243          convm = f_dmass(ind) 
    229244          dtheta_rhodz=f_dtheta_rhodz(ind) 
    230245          du=f_du(ind) 
    231           CALL compute_caldyn_horiz(u,mass,qu,theta,pk,geopot, hflux,divm,dtheta_rhodz,du) 
    232           wflux=f_wflux(ind) 
    233           wwuu=f_wwuu(ind) 
    234           dps=f_dps(ind) 
    235           CALL compute_caldyn_vert(u,theta,mass,divm, wflux,wwuu, dps, dtheta_rhodz, du) 
     246          CALL compute_caldyn_horiz(u,mass,qu,theta,pk,geopot, hflux,convm,dtheta_rhodz,du) 
     247          IF(caldyn_eta==eta_mass) THEN 
     248             wflux=f_wflux(ind) 
     249             wwuu=f_wwuu(ind) 
     250             dps=f_dps(ind) 
     251             CALL compute_caldyn_vert(u,theta,mass,convm, wflux,wwuu, dps, dtheta_rhodz, du) 
     252          END IF 
    236253       ENDDO 
    237254        
    … …  
    252269       CALL writefield("dps",f_dps) 
    253270       CALL writefield("mass",f_mass) 
     271       CALL writefield("dmass",f_dmass) 
    254272       CALL writefield("vort",f_qu) 
    255273       CALL writefield("theta",f_theta) 
     274       CALL writefield("exner",f_pk) 
     275       CALL writefield("pv",f_qv) 
    256276 
    257277    END IF 
    … …  
    263283END SUBROUTINE caldyn 
    264284     
    265 SUBROUTINE compute_pvort(ps,u,theta_rhodz, rhodz,theta,qu) 
     285SUBROUTINE compute_pvort(ps,u,theta_rhodz, rhodz,theta,qu,qv) 
    266286  USE icosa 
    267287  USE disvert_mod 
    … …  
    273293  REAL(rstd),INTENT(IN)  :: ps(iim*jjm) 
    274294  REAL(rstd),INTENT(IN)  :: theta_rhodz(iim*jjm,llm) 
    275   REAL(rstd),INTENT(OUT) :: rhodz(iim*jjm,llm) 
     295  REAL(rstd),INTENT(INOUT) :: rhodz(iim*jjm,llm) 
    276296  REAL(rstd),INTENT(OUT) :: theta(iim*jjm,llm) 
    277297  REAL(rstd),INTENT(OUT) :: qu(iim*3*jjm,llm) 
     298  REAL(rstd),INTENT(OUT) :: qv(iim*2*jjm,llm) 
    278299   
    279300  INTEGER :: i,j,ij,l 
    280   REAL(rstd) :: etav,hv 
    281   REAL(rstd) :: qv(2*iim*jjm,llm)     ! potential velocity   
     301  REAL(rstd) :: etav,hv, m 
     302!  REAL(rstd) :: qv(2*iim*jjm,llm)     ! potential velocity   
    282303   
    283304 
    284305  CALL trace_start("compute_pvort")   
    285306 
    286   CALL wait_message(req_ps)   
     307  IF(caldyn_eta==eta_mass) THEN 
     308     CALL wait_message(req_ps)   
     309  ELSE 
     310     CALL wait_message(req_mass) 
     311  END IF 
    287312  CALL wait_message(req_theta_rhodz)  
    288313 
    289 !!! Compute mass & theta 
    290    DO l = ll_begin,ll_end 
    291      CALL test_message(req_u)  
    292      DO j=jj_begin-1,jj_end+1 
    293         DO i=ii_begin-1,ii_end+1 
    294            ij=(j-1)*iim+i 
    295            rhodz(ij,l) = ( mass_dak(l)+ps(ij)*mass_dbk(l) )/g 
    296            theta(ij,l) = theta_rhodz(ij,l)/rhodz(ij,l) 
     314  IF(caldyn_eta==eta_mass) THEN ! Compute mass & theta 
     315     DO l = ll_begin,ll_end 
     316        CALL test_message(req_u)  
     317        DO j=jj_begin-1,jj_end+1 
     318           DO i=ii_begin-1,ii_end+1 
     319              ij=(j-1)*iim+i 
     320              m = ( ap(l) - ap(l+1) + (bp(l)-bp(l+1))*ps(ij) )/g  
     321!              m = ( mass_dak(l)+ps(ij)*mass_dbk(l) )/g 
     322              rhodz(ij,l) = m 
     323              theta(ij,l) = theta_rhodz(ij,l)/rhodz(ij,l) 
     324           ENDDO 
    297325        ENDDO 
    298326     ENDDO 
    299   ENDDO 
     327  ELSE ! Compute only theta 
     328     DO l = ll_begin,ll_end 
     329        CALL test_message(req_u)  
     330        DO j=jj_begin-1,jj_end+1 
     331           DO i=ii_begin-1,ii_end+1 
     332              ij=(j-1)*iim+i 
     333              theta(ij,l) = theta_rhodz(ij,l)/rhodz(ij,l) 
     334           ENDDO 
     335        ENDDO 
     336     ENDDO 
     337  END IF 
    300338 
    301339  CALL wait_message(req_u)    
    … …  
    377415                pk(ij,l) = exner_ik 
    378416                ! 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  
     417                geopot(ij,l+1) = geopot(ij,l) + (g*kappa)*rhodz(ij,l)*theta(ij,l)*exner_ik/p_ik 
    380418             ENDDO 
    381419          ENDDO 
    … …  
    413451          END DO 
    414452       END DO 
    415  
    416453       DO l = 1,llm 
    417454          !$OMP DO SCHEDULE(STATIC)  
    … …  
    434471  END SUBROUTINE compute_geopot 
    435472 
    436   SUBROUTINE compute_caldyn_horiz(u,rhodz,qu,theta,pk,geopot, hflux,divm, dtheta_rhodz, du) 
     473  SUBROUTINE compute_caldyn_horiz(u,rhodz,qu,theta,pk,geopot, hflux,convm, dtheta_rhodz, du) 
    437474  USE icosa 
    438475  USE disvert_mod 
    … …  
    449486 
    450487    REAL(rstd),INTENT(OUT) :: hflux(iim*3*jjm,llm) ! hflux in kg/s 
    451     REAL(rstd),INTENT(OUT) :: divm(iim*jjm,llm)  ! mass flux divergence 
     488    REAL(rstd),INTENT(OUT) :: convm(iim*jjm,llm)  ! mass flux convergence 
    452489    REAL(rstd),INTENT(OUT) :: dtheta_rhodz(iim*jjm,llm) 
    453490    REAL(rstd),INTENT(OUT) :: du(iim*3*jjm,llm) 
    … …  
    483520      DO i=ii_begin,ii_end 
    484521        ij=(j-1)*iim+i   
    485         ! divm = +div(mass flux), sign convention as in Ringler et al. 2012, eq. 21 
    486         divm(ij,l)= 1./Ai(ij)*(ne_right*hflux(ij+u_right,l)   +  & 
     522        ! convm = -div(mass flux), sign convention as in Ringler et al. 2012, eq. 21 
     523        convm(ij,l)= -1./Ai(ij)*(ne_right*hflux(ij+u_right,l)   +  & 
    487524                                ne_rup*hflux(ij+u_rup,l)       +  &   
    488525                                ne_lup*hflux(ij+u_lup,l)       +  &   
    … …  
    659696END SUBROUTINE compute_caldyn_horiz 
    660697 
    661 SUBROUTINE compute_caldyn_vert(u,theta,rhodz,divm, wflux,wwuu, dps,dtheta_rhodz,du) 
     698SUBROUTINE compute_caldyn_vert(u,theta,rhodz,convm, wflux,wwuu, dps,dtheta_rhodz,du) 
    662699  USE icosa 
    663700  USE disvert_mod 
    … …  
    669706    REAL(rstd),INTENT(IN)  :: theta(iim*jjm,llm) 
    670707    REAL(rstd),INTENT(IN)  :: rhodz(iim*jjm,llm) 
    671     REAL(rstd),INTENT(INOUT)  :: divm(iim*jjm,llm)  ! mass flux divergence 
     708    REAL(rstd),INTENT(INOUT)  :: convm(iim*jjm,llm)  ! mass flux convergence 
    672709 
    673710    REAL(rstd),INTENT(OUT) :: wflux(iim*jjm,llm+1) ! vertical mass flux (kg/m2/s) 
    … …  
    684721 
    685722!$OMP BARRIER    
    686 !!! cumulate mass flux divergence from top to bottom 
     723!!! cumulate mass flux convergence from top to bottom 
    687724  DO  l = llm-1, 1, -1 
    688725    IF (caldyn_conserv==energy) CALL test_message(req_qu)  
    … …  
    691728      DO i=ii_begin,ii_end 
    692729        ij=(j-1)*iim+i 
    693         divm(ij,l) = divm(ij,l) + divm(ij,l+1) 
     730        convm(ij,l) = convm(ij,l) + convm(ij,l+1) 
    694731      ENDDO 
    695732    ENDDO 
    696733  ENDDO 
    697734 
    698 ! IMPLICIT FLUSH on divm 
     735! IMPLICIT FLUSH on convm 
    699736!!!!!!!!!!!!!!!!!!!!!!!!! 
    700737 
    … …  
    705742        ij=(j-1)*iim+i 
    706743        ! dps/dt = -int(div flux)dz 
    707         dps(ij)=-divm(ij,1) * g  
     744        dps(ij) = convm(ij,1) * g  
    708745      ENDDO 
    709746    ENDDO 
    … …  
    718755        ! w = int(z,ztop,div(flux)dz) + B(eta)dps/dt 
    719756        ! => w>0 for upward transport 
    720         wflux( ij, l ) = divm( ij, l ) - bp(l) * divm( ij, 1 ) 
     757        wflux( ij, l ) = bp(l) * convm( ij, 1 ) - convm( ij, l )  
    721758      ENDDO 
    722759    ENDDO 

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

    r159 r162  
    99  INTEGER  :: itau_sync=10 
    1010 
    11   TYPE(t_message) :: req_ps0, req_theta_rhodz0, req_u0, req_q0 
     11  TYPE(t_message) :: req_ps0, req_mass0, req_theta_rhodz0, req_u0, req_q0 
    1212 
    1313  TYPE(t_field),POINTER :: f_q(:) 
    … …  
    130130    CALL allocate_field(f_hfluxt,field_u,type_real,llm)   ! mass "fluxes" accumulated in time 
    131131    CALL allocate_field(f_wflux,field_t,type_real,llm+1)  ! vertical mass fluxes 
     132    CALL allocate_field(f_dmass,field_t,type_real,llm, name='dmass') 
    132133 
    133134    IF(caldyn_eta == eta_mass) THEN ! eta = mass coordinate (default) 
    … …  
    136137       CALL allocate_field(f_wfluxt,field_t,type_real,llm+1,name='wfluxt') 
    137138       ! the following are unused but must point to something 
    138        f_massm1 => f_mass 
    139        f_dmass => f_mass 
     139!       f_massm1 => f_mass 
    140140    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) 
     141       CALL allocate_field(f_massm1,field_t,type_real,llm, name='massm1') 
    144142       ! the following are unused but must point to something 
    145143       f_wfluxt => f_wflux 
    … …  
    147145       f_psm1 => f_phis 
    148146    END IF 
    149  
    150147 
    151148    def='runge_kutta' 
    … …  
    197194 
    198195    CALL init_message(f_ps,req_i0,req_ps0) 
     196    CALL init_message(f_mass,req_i0,req_mass0) 
    199197    CALL init_message(f_theta_rhodz,req_i0,req_theta_rhodz0) 
    200198    CALL init_message(f_u,req_e0_vect,req_u0) 
    … …  
    208206  USE disvert_mod 
    209207  USE caldyn_mod 
    210   USE caldyn_gcm_mod, ONLY : req_ps 
     208  USE caldyn_gcm_mod, ONLY : req_ps, req_mass 
    211209  USE etat0_mod 
    212210  USE guided_mod 
    … …  
    237235     CALL swap_dimensions(ind) 
    238236     CALL swap_geometry(ind) 
    239      rhodz=f_rhodz(ind); ps=f_ps(ind) 
    240      CALL compute_rhodz(.TRUE., ps, rhodz) ! save rhodz for transport scheme before dynamics update ps 
     237     rhodz=f_rhodz(ind); mass=f_mass(ind); ps=f_ps(ind) 
     238     IF(caldyn_eta==eta_mass) THEN 
     239        CALL compute_rhodz(.TRUE., ps, rhodz) ! save rhodz for transport scheme before dynamics update ps 
     240     ELSE 
     241        rhodz(:,:)=mass(:,:) 
     242     END IF 
    241243  END DO 
    242244  fluxt_zero=.TRUE. 
    … …  
    245247    IF (MOD(it,itau_sync)==0) THEN 
    246248      CALL send_message(f_ps,req_ps0) 
     249      CALL send_message(f_mass,req_mass0) 
    247250      CALL send_message(f_theta_rhodz,req_theta_rhodz0)  
    248251      CALL send_message(f_u,req_u0) 
    249252      CALL send_message(f_q,req_q0)  
    250253      CALL wait_message(req_ps0) 
     254      CALL wait_message(req_mass0) 
    251255      CALL wait_message(req_theta_rhodz0)  
    252256      CALL wait_message(req_u0) 
    … …  
    266270       CALL caldyn((stage==1) .AND. (MOD(it,itau_out)==0), & 
    267271            f_phis,f_ps,f_mass,f_theta_rhodz,f_u, f_q, & 
    268             f_hflux, f_wflux, f_dps, f_dtheta_rhodz, f_du) 
     272            f_hflux, f_wflux, f_dps, f_dmass, f_dtheta_rhodz, f_du) 
    269273       SELECT CASE (scheme) 
    270274       CASE(euler) 
    … …  
    326330       CALL swap_dimensions(ind) 
    327331       CALL swap_geometry(ind) 
    328        IF(with_dps) THEN 
    329          ps=f_ps(ind) ; dps=f_dps(ind) ;  
    330  
    331          IF (omp_first) THEN 
    332            DO j=jj_begin,jj_end 
    333              DO i=ii_begin,ii_end 
    334                ij=(j-1)*iim+i 
    335                ps(ij)=ps(ij)+dt*dps(ij) 
    336              ENDDO 
    337            ENDDO 
    338          ENDIF 
    339           
    340          hflux=f_hflux(ind);     hfluxt=f_hfluxt(ind) 
    341          wflux=f_wflux(ind);     wfluxt=f_wfluxt(ind) 
    342          CALL accumulate_fluxes(hflux,wflux,hfluxt,wfluxt,dt,fluxt_zero(ind)) 
    343        END IF 
     332 
     333       IF(with_dps) THEN ! update ps/mass 
     334          IF(caldyn_eta==eta_mass) THEN ! update ps 
     335             ps=f_ps(ind) ; dps=f_dps(ind) ;               
     336             IF (omp_first) THEN 
     337                DO j=jj_begin,jj_end 
     338                   DO i=ii_begin,ii_end 
     339                      ij=(j-1)*iim+i 
     340                      ps(ij)=ps(ij)+dt*dps(ij) 
     341                   ENDDO 
     342                ENDDO 
     343             ENDIF  
     344          ELSE ! update mass 
     345             mass=f_mass(ind) ; dmass=f_dmass(ind) ;               
     346             DO l=1,llm 
     347                DO j=jj_begin,jj_end 
     348                   DO i=ii_begin,ii_end 
     349                      ij=(j-1)*iim+i 
     350                      mass(ij,l)=mass(ij,l)+dt*dmass(ij,l) 
     351                   ENDDO 
     352                ENDDO 
     353             END DO 
     354          END IF 
     355 
     356          hflux=f_hflux(ind);     hfluxt=f_hfluxt(ind) 
     357          wflux=f_wflux(ind);     wfluxt=f_wfluxt(ind) 
     358          CALL accumulate_fluxes(hflux,wflux,hfluxt,wfluxt,dt,fluxt_zero(ind)) 
     359       END IF ! update ps/mass 
    344360        
    345361       u=f_u(ind) ; theta_rhodz=f_theta_rhodz(ind) 
    … …  
    377393      IF(caldyn_eta==eta_mass) THEN 
    378394         IF (omp_first) THEN 
     395 
    379396            DO ind=1,ndomain 
    380397               CALL swap_dimensions(ind) 
    381398               CALL swap_geometry(ind) 
    382                ps=f_ps(ind)    
    383                psm1=f_psm1(ind)  
    384                dps=f_dps(ind)  
     399               ps=f_ps(ind) ; psm1=f_psm1(ind) ; dps=f_dps(ind)  
    385400                
    386401               IF (stage==1) THEN ! first stage : save model state in XXm1 
    … …  
    402417            ENDDO 
    403418         ENDIF 
    404     
    405419         CALL send_message(f_ps,req_ps) 
     420       
     421      ELSE ! Lagrangian coordinate, deal with mass 
     422         DO ind=1,ndomain 
     423            CALL swap_dimensions(ind) 
     424            CALL swap_geometry(ind) 
     425            mass=f_mass(ind); dmass=f_dmass(ind); massm1=f_massm1(ind) 
     426 
     427            IF (stage==1) THEN ! first stage : save model state in XXm1 
     428               DO l=ll_begin,ll_end 
     429                  DO j=jj_begin,jj_end 
     430                     DO i=ii_begin,ii_end 
     431                        ij=(j-1)*iim+i 
     432                        massm1(ij,l)=mass(ij,l) 
     433                     ENDDO 
     434                  ENDDO 
     435               ENDDO 
     436            END IF 
     437 
     438            ! updates are of the form x1 := x0 + tau*f(x1) 
     439            DO l=ll_begin,ll_end 
     440               DO j=jj_begin,jj_end 
     441                  DO i=ii_begin,ii_end 
     442                     ij=(j-1)*iim+i 
     443                     mass(ij,l)=massm1(ij,l)+tau*dmass(ij,l) 
     444                  ENDDO 
     445               ENDDO 
     446            ENDDO 
     447         END DO 
     448         CALL send_message(f_mass,req_mass) 
     449 
    406450      END IF 
    407451 
    … …  
    410454         CALL swap_dimensions(ind) 
    411455         CALL swap_geometry(ind) 
    412          ps=f_ps(ind)   ; u=f_u(ind)   ; theta_rhodz=f_theta_rhodz(ind) 
    413          psm1=f_psm1(ind) ; um1=f_um1(ind) ; theta_rhodzm1=f_theta_rhodzm1(ind) 
    414          dps=f_dps(ind) ; du=f_du(ind) ; dtheta_rhodz=f_dtheta_rhodz(ind) 
     456         u=f_u(ind)      ; du=f_du(ind)      ; um1=f_um1(ind)  
     457         theta_rhodz=f_theta_rhodz(ind) 
     458         theta_rhodzm1=f_theta_rhodzm1(ind) 
     459         dtheta_rhodz=f_dtheta_rhodz(ind) 
    415460          
    416461         IF (stage==1) THEN ! first stage : save model state in XXm1 
    417              
    418462           DO l=ll_begin,ll_end 
    419463             DO j=jj_begin,jj_end 
    … …  
    427471             ENDDO 
    428472           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 
    440  
    441          END IF 
    442          ! updates are of the form x1 := x0 + tau*f(x1) 
    443           
     473         END IF         
     474 
    444475         DO l=ll_begin,ll_end 
    445476           DO j=jj_begin,jj_end 
    … …  
    453484           ENDDO 
    454485         ENDDO 
    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           
     486 
    466487         IF(stage==nb_stage) THEN ! accumulate mass fluxes at last stage 
    467488            hflux=f_hflux(ind);     hfluxt=f_hfluxt(ind) 
    … …  
    562583          ENDDO 
    563584       END IF 
     585 
    564586    ELSE 
    565587 
    … …  
    590612  END SUBROUTINE accumulate_fluxes 
    591613   
     614  FUNCTION maxval_i(p) 
     615    USE icosa 
     616    IMPLICIT NONE 
     617    REAL(rstd), DIMENSION(iim*jjm) :: p 
     618    REAL(rstd) :: maxval_i 
     619    INTEGER :: j, ij 
     620     
     621    maxval_i=p((jj_begin-1)*iim+ii_begin) 
     622     
     623    DO j=jj_begin-1,jj_end+1 
     624       ij=(j-1)*iim 
     625       maxval_i = MAX(maxval_i, MAXVAL(p(ij+ii_begin:ij+ii_end))) 
     626    END DO 
     627  END FUNCTION maxval_i 
     628 
     629  FUNCTION maxval_ik(p) 
     630    USE icosa 
     631    IMPLICIT NONE 
     632    REAL(rstd) :: p(iim*jjm, llm) 
     633    REAL(rstd) :: maxval_ik(llm) 
     634    INTEGER :: l,j, ij 
     635     
     636    DO l=1,llm 
     637       maxval_ik(l)=p((jj_begin-1)*iim+ii_begin,l) 
     638       DO j=jj_begin-1,jj_end+1 
     639          ij=(j-1)*iim 
     640          maxval_ik(l) = MAX(maxval_ik(l), MAXVAL(p(ij+ii_begin:ij+ii_end,l))) 
     641       END DO 
     642    END DO 
     643  END FUNCTION maxval_ik 
     644 
    592645END MODULE timeloop_gcm_mod 
Note: See TracChangeset for help on using the changeset viewer.