Changeset 362 for codes


Ignore:
Timestamp:
09/25/15 14:36:36 (9 years ago)
Author:
dubos
Message:

Introduced DEC convention for velocity - HEVI only - cleanup to follow

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

Legend:

Unmodified
Added
Removed

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

    r361 r362  
    22  USE icosa 
    33  USE transfert_mod 
     4  USE caldyn_kernels_base_mod 
    45  USE caldyn_kernels_mod 
    56  IMPLICIT NONE 
    67  PRIVATE 
    78 
    8   PUBLIC init_caldyn, caldyn_BC, caldyn, req_ps, req_mass 
    9  
     9  PUBLIC init_caldyn, caldyn_BC, caldyn 
     10   
    1011CONTAINS 
    1112   
    … …  
    1314    USE icosa 
    1415    USE observable_mod 
    15 !    USE exner_mod 
    1616    USE mpipara 
    1717    USE omp_para 

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

    r361 r362  
    22  USE icosa 
    33  USE transfert_mod 
    4   USE caldyn_kernels_mod 
     4  USE caldyn_kernels_base_mod 
     5  USE caldyn_kernels_hevi_mod 
    56  USE caldyn_gcm_mod 
    6  
    77  IMPLICIT NONE 
    88  PRIVATE 
    9  
    109  PUBLIC caldyn_hevi 
    1110 
    … …  
    7574    IF(caldyn_eta==eta_mass) THEN 
    7675       CALL send_message(f_ps,req_ps) ! COM00 
     76       CALL wait_message(req_ps) ! COM00 
    7777    ELSE 
    7878       CALL send_message(f_mass,req_mass) ! COM00 
    79     END IF 
    80  
    81     ! Block below moved from caldyn_pvort when splitting it into caldyn_theta and caldyn_pvort_only 
    82     IF(caldyn_eta==eta_mass) THEN 
    83        CALL wait_message(req_ps) ! COM00 
    84     ELSE 
    8579       CALL wait_message(req_mass) ! COM00 
    8680    END IF 
    87  
    8881    CALL send_message(f_theta_rhodz,req_theta_rhodz) ! COM01 
    8982    CALL wait_message(req_theta_rhodz) ! COM01 Moved from caldyn_pvort 
    … …  
    10093       theta = f_theta(ind) 
    10194       pk = f_pk(ind) 
    102        geopot = f_geopot(ind)   
     95       geopot = f_geopot(ind) 
    10396       CALL compute_theta(ps,theta_rhodz, mass,theta) 
    10497       CALL compute_geopot(ps,mass,theta, pk,geopot) 

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

    r361 r362  
    22  USE icosa 
    33  USE transfert_mod 
     4  USE caldyn_kernels_base_mod 
    45  IMPLICIT NONE 
    5  
    6   INTEGER, PARAMETER :: energy=1, enstrophy=2 
    7   TYPE(t_field),POINTER :: f_out_u(:), f_qu(:), f_qv(:) 
    8   REAL(rstd),SAVE,POINTER :: out_u(:,:), p(:,:), qu(:,:) 
    9   !$OMP THREADPRIVATE(out_u, p, qu) 
    10  
    11   ! temporary shared variable for caldyn 
    12   TYPE(t_field),POINTER :: f_pk(:) 
    13   TYPE(t_field),POINTER :: f_wwuu(:) 
    14   TYPE(t_field),POINTER :: f_planetvel(:) 
    15  
    16   INTEGER :: caldyn_conserv 
    17   !$OMP THREADPRIVATE(caldyn_conserv)  
    18  
    19   TYPE(t_message) :: req_ps, req_mass, req_theta_rhodz, req_u, req_qu 
    20  
     6  PRIVATE 
     7 
     8  PUBLIC :: compute_planetvel, compute_pvort, compute_geopot, & 
     9       compute_caldyn_horiz, compute_caldyn_vert 
    2110CONTAINS 
    2211 
    … …  
    4029    CALL compute_wind2D_perp_from_lonlat_compound(ulon, ulat, planetvel) 
    4130  END SUBROUTINE compute_planetvel 
    42  
    43   !********************** compute_pvort = compute_theta + compute_pvort_only ****************** 
    4431 
    4532  SUBROUTINE compute_pvort(ps,u,theta_rhodz, rhodz,theta,qu,qv) 
    … …  
    6047    INTEGER :: i,j,ij,l 
    6148    REAL(rstd) :: etav,hv, m 
    62     !  REAL(rstd) :: qv(2*iim*jjm,llm)     ! potential velocity   
    63  
    6449    CALL trace_start("compute_pvort")   
    6550 
    … …  
    130115    CALL trace_end("compute_pvort") 
    131116  END SUBROUTINE compute_pvort 
    132  
    133   SUBROUTINE compute_theta(ps,theta_rhodz, rhodz,theta) 
    134     USE icosa 
    135     USE disvert_mod, ONLY : mass_dak, mass_dbk, caldyn_eta, eta_mass 
    136     USE exner_mod 
    137     USE trace 
    138     USE omp_para 
    139     IMPLICIT NONE 
    140     REAL(rstd),INTENT(IN)  :: ps(iim*jjm) 
    141     REAL(rstd),INTENT(IN)  :: theta_rhodz(iim*jjm,llm) 
    142     REAL(rstd),INTENT(INOUT) :: rhodz(iim*jjm,llm) 
    143     REAL(rstd),INTENT(OUT) :: theta(iim*jjm,llm) 
    144     INTEGER :: ij,l 
    145     REAL(rstd) :: m 
    146     CALL trace_start("compute_theta")   
    147  
    148     IF(caldyn_eta==eta_mass) THEN ! Compute mass & theta 
    149        DO l = ll_begin,ll_end 
    150           !DIR$ SIMD 
    151           DO ij=ij_begin_ext,ij_end_ext 
    152              m = ( mass_dak(l)+ps(ij)*mass_dbk(l) )/g 
    153              rhodz(ij,l) = m 
    154              theta(ij,l) = theta_rhodz(ij,l)/rhodz(ij,l) 
    155           ENDDO 
    156        ENDDO 
    157     ELSE ! Compute only theta 
    158        DO l = ll_begin,ll_end 
    159           !DIR$ SIMD 
    160           DO ij=ij_begin_ext,ij_end_ext 
    161              theta(ij,l) = theta_rhodz(ij,l)/rhodz(ij,l) 
    162           ENDDO 
    163        ENDDO 
    164     END IF 
    165  
    166     CALL trace_end("compute_theta") 
    167   END SUBROUTINE compute_theta 
    168  
    169   SUBROUTINE compute_pvort_only(u,rhodz,qu,qv) 
    170     USE icosa 
    171     USE exner_mod 
    172     USE trace 
    173     USE omp_para 
    174     IMPLICIT NONE 
    175     REAL(rstd),INTENT(IN)  :: u(iim*3*jjm,llm) 
    176     REAL(rstd),INTENT(INOUT) :: rhodz(iim*jjm,llm) 
    177     REAL(rstd),INTENT(OUT) :: qu(iim*3*jjm,llm) 
    178     REAL(rstd),INTENT(OUT) :: qv(iim*2*jjm,llm) 
    179  
    180     INTEGER :: ij,l 
    181     REAL(rstd) :: etav,hv 
    182  
    183     CALL trace_start("compute_pvort_only")   
    184 !!! Compute shallow-water potential vorticity 
    185     DO l = ll_begin,ll_end 
    186        !DIR$ SIMD 
    187        DO ij=ij_begin_ext,ij_end_ext 
    188           etav= 1./Av(ij+z_up)*(  ne_rup        * u(ij+u_rup,l)        * de(ij+u_rup)         & 
    189                + ne_left * u(ij+t_rup+u_left,l) * de(ij+t_rup+u_left)  & 
    190                - ne_lup        * u(ij+u_lup,l)        * de(ij+u_lup) )                                
    191  
    192           hv =  Riv2(ij,vup)          * rhodz(ij,l)            & 
    193                + Riv2(ij+t_rup,vldown) * rhodz(ij+t_rup,l)     & 
    194                + Riv2(ij+t_lup,vrdown) * rhodz(ij+t_lup,l) 
    195  
    196           qv(ij+z_up,l) = ( etav+fv(ij+z_up) )/hv 
    197  
    198           etav = 1./Av(ij+z_down)*(  ne_ldown         * u(ij+u_ldown,l)          * de(ij+u_ldown)          & 
    199                + ne_right * u(ij+t_ldown+u_right,l)  * de(ij+t_ldown+u_right)  & 
    200                - ne_rdown         * u(ij+u_rdown,l)          * de(ij+u_rdown) ) 
    201  
    202           hv =  Riv2(ij,vdown)        * rhodz(ij,l)          & 
    203                + Riv2(ij+t_ldown,vrup) * rhodz(ij+t_ldown,l)  & 
    204                + Riv2(ij+t_rdown,vlup) * rhodz(ij+t_rdown,l) 
    205  
    206           qv(ij+z_down,l) =( etav+fv(ij+z_down) )/hv 
    207  
    208        ENDDO 
    209  
    210        !DIR$ SIMD 
    211        DO ij=ij_begin,ij_end 
    212           qu(ij+u_right,l) = 0.5*(qv(ij+z_rdown,l)+qv(ij+z_rup,l))   
    213           qu(ij+u_lup,l) = 0.5*(qv(ij+z_up,l)+qv(ij+z_lup,l))   
    214           qu(ij+u_ldown,l) = 0.5*(qv(ij+z_ldown,l)+qv(ij+z_down,l))   
    215        END DO 
    216  
    217     ENDDO 
    218  
    219     CALL trace_end("compute_pvort_only") 
    220  
    221   END SUBROUTINE compute_pvort_only 
    222  
    223   !**************************** Geopotential ***************************** 
    224  
    225   SUBROUTINE compute_geopot(ps,rhodz,theta, pk,geopot)  
    226     USE icosa 
    227     USE disvert_mod 
    228     USE exner_mod 
    229     USE trace 
    230     USE omp_para 
    231     IMPLICIT NONE 
    232     REAL(rstd),INTENT(INOUT) :: ps(iim*jjm) 
    233     REAL(rstd),INTENT(IN)    :: rhodz(iim*jjm,llm) 
    234     REAL(rstd),INTENT(IN)    :: theta(iim*jjm,llm)    ! potential temperature 
    235     REAL(rstd),INTENT(INOUT) :: pk(iim*jjm,llm)       ! Exner function 
    236     REAL(rstd),INTENT(INOUT) :: geopot(iim*jjm,llm+1) ! geopotential 
    237  
    238     INTEGER :: i,j,ij,l 
    239     REAL(rstd) :: p_ik, exner_ik 
    240     INTEGER    :: ij_omp_begin_ext, ij_omp_end_ext 
    241  
    242  
    243     CALL trace_start("compute_geopot") 
    244  
    245     CALL distrib_level(ij_end_ext-ij_begin_ext+1,ij_omp_begin_ext,ij_omp_end_ext) 
    246     ij_omp_begin_ext=ij_omp_begin_ext+ij_begin_ext-1 
    247     ij_omp_end_ext=ij_omp_end_ext+ij_begin_ext-1 
    248  
    249     IF(caldyn_eta==eta_mass) THEN 
    250  
    251 !!! Compute exner function and geopotential 
    252        DO l = 1,llm 
    253           !DIR$ SIMD 
    254           DO ij=ij_omp_begin_ext,ij_omp_end_ext          
    255              p_ik = ptop + mass_ak(l) + mass_bk(l)*ps(ij) ! FIXME : leave ps for the moment ; change ps to Ms later 
    256              !         p_ik = ptop + g*(mass_ak(l)+ mass_bk(l)*ps(i,j)) 
    257              exner_ik = cpp * (p_ik/preff) ** kappa 
    258              pk(ij,l) = exner_ik 
    259              ! specific volume v = kappa*theta*pi/p = dphi/g/rhodz 
    260              geopot(ij,l+1) = geopot(ij,l) + (g*kappa)*rhodz(ij,l)*theta(ij,l)*exner_ik/p_ik 
    261           ENDDO 
    262        ENDDO 
    263        !       ENDIF 
    264     ELSE  
    265        ! We are using a Lagrangian vertical coordinate 
    266        ! Pressure must be computed first top-down (temporarily stored in pk) 
    267        ! Then Exner pressure and geopotential are computed bottom-up 
    268        ! Notice that the computation below should work also when caldyn_eta=eta_mass           
    269  
    270        IF(boussinesq) THEN ! compute only geopotential : pressure pk will be computed in compute_caldyn_horiz 
    271           ! specific volume 1 = dphi/g/rhodz 
    272           !         IF (is_omp_level_master) THEN ! no openMP on vertical due to dependency 
    273           DO l = 1,llm 
    274              !DIR$ SIMD 
    275              DO ij=ij_omp_begin_ext,ij_omp_end_ext          
    276                 geopot(ij,l+1) = geopot(ij,l) + g*rhodz(ij,l) 
    277              ENDDO 
    278           ENDDO 
    279        ELSE ! non-Boussinesq, compute geopotential and Exner pressure 
    280           ! uppermost layer 
    281  
    282           !DIR$ SIMD 
    283           DO ij=ij_omp_begin_ext,ij_omp_end_ext          
    284              pk(ij,llm) = ptop + (.5*g)*rhodz(ij,llm) 
    285           END DO 
    286           ! other layers 
    287           DO l = llm-1, 1, -1 
    288              !DIR$ SIMD 
    289              DO ij=ij_omp_begin_ext,ij_omp_end_ext          
    290                 pk(ij,l) = pk(ij,l+1) + (.5*g)*(rhodz(ij,l)+rhodz(ij,l+1)) 
    291              END DO 
    292           END DO 
    293           ! surface pressure (for diagnostics) 
    294           DO ij=ij_omp_begin_ext,ij_omp_end_ext          
    295              ps(ij) = pk(ij,1) + (.5*g)*rhodz(ij,1) 
    296           END DO 
    297  
    298           ! specific volume v = kappa*theta*pi/p = dphi/g/rhodz 
    299           DO l = 1,llm 
    300              !DIR$ SIMD 
    301              DO ij=ij_omp_begin_ext,ij_omp_end_ext          
    302                 p_ik = pk(ij,l) 
    303                 exner_ik = cpp * (p_ik/preff) ** kappa 
    304                 geopot(ij,l+1) = geopot(ij,l) + (g*kappa)*rhodz(ij,l)*theta(ij,l)*exner_ik/p_ik  
    305                 pk(ij,l) = exner_ik 
    306              ENDDO 
    307           ENDDO 
    308        END IF 
    309  
    310     END IF 
    311  
    312     !ym flush geopot 
    313     !$OMP BARRIER 
    314  
    315     CALL trace_end("compute_geopot") 
    316  
    317   END SUBROUTINE compute_geopot 
    318117 
    319118  !************************* caldyn_horiz = caldyn_fast + caldyn_slow ********************** 
    … …  
    490289!!! Compute bernouilli term = Kinetic Energy + geopotential 
    491290    IF(boussinesq) THEN 
    492        ! first use hydrostatic balance with theta*rhodz to find pk (Lagrange multiplier=pressure)  
    493        ! uppermost layer 
    494        !DIR$ SIMD 
    495        DO ij=ij_begin_ext,ij_end_ext          
    496           pk(ij,llm) = ptop + (.5*g)*theta(ij,llm)*rhodz(ij,llm) 
    497        END DO 
    498        ! other layers 
    499        DO l = llm-1, 1, -1 
    500           !          !$OMP DO SCHEDULE(STATIC)  
    501           !DIR$ SIMD 
    502           DO ij=ij_begin_ext,ij_end_ext          
    503              pk(ij,l) = pk(ij,l+1) + (.5*g)*(theta(ij,l)*rhodz(ij,l)+theta(ij,l+1)*rhodz(ij,l+1)) 
    504           END DO 
    505        END DO 
    506        ! surface pressure (for diagnostics) FIXME 
    507        ! DO ij=ij_begin_ext,ij_end_ext          
    508        !   ps(ij) = pk(ij,1) + (.5*g)*theta(ij,1)*rhodz(ij,1) 
    509        ! END DO 
    510        ! now pk contains the Lagrange multiplier (pressure) 
    511  
    512291       DO l=ll_begin,ll_end 
    513292          !DIR$ SIMD 
    … …  
    572351  END SUBROUTINE compute_caldyn_horiz 
    573352 
    574   SUBROUTINE compute_caldyn_fast(tau,u,rhodz,theta,pk,geopot, du) 
    575     USE icosa 
    576     USE disvert_mod 
    577     USE exner_mod 
    578     USE trace 
    579     USE omp_para 
    580     IMPLICIT NONE 
    581     REAL(rstd), INTENT(IN) :: tau ! "solve" u-tau*du/dt = rhs 
    582     REAL(rstd),INTENT(INOUT)  :: u(iim*3*jjm,llm)    ! prognostic "velocity" 
    583     REAL(rstd),INTENT(IN)  :: rhodz(iim*jjm,llm) 
    584     REAL(rstd),INTENT(IN)  :: theta(iim*jjm,llm)  ! potential temperature 
    585     REAL(rstd),INTENT(INOUT) :: pk(iim*jjm,llm) ! Exner function 
    586     REAL(rstd),INTENT(IN)  :: geopot(iim*jjm,llm+1)    ! geopotential 
    587     REAL(rstd),INTENT(OUT) :: du(iim*3*jjm,llm) 
    588     REAL(rstd) :: berni(iim*jjm,llm)  ! Bernoulli function 
    589  
    590     INTEGER :: i,j,ij,l 
    591  
    592     CALL trace_start("compute_caldyn_fast") 
    593  
    594     !    CALL wait_message(req_theta_rhodz)  
    595      
    596     ! Compute bernouilli term 
    597     IF(boussinesq) THEN 
    598        ! first use hydrostatic balance with theta*rhodz to find pk (Lagrange multiplier=pressure)  
    599        ! uppermost layer 
    600        !DIR$ SIMD 
    601        DO ij=ij_begin_ext,ij_end_ext          
    602           pk(ij,llm) = ptop + (.5*g)*theta(ij,llm)*rhodz(ij,llm) 
    603        END DO 
    604        ! other layers 
    605        DO l = llm-1, 1, -1 
    606           !          !$OMP DO SCHEDULE(STATIC)  
    607           !DIR$ SIMD 
    608           DO ij=ij_begin_ext,ij_end_ext          
    609              pk(ij,l) = pk(ij,l+1) + (.5*g)*(theta(ij,l)*rhodz(ij,l)+theta(ij,l+1)*rhodz(ij,l+1)) 
    610           END DO 
    611        END DO 
    612        ! now pk contains the Lagrange multiplier (pressure) 
    613        DO l=ll_begin,ll_end 
    614           !DIR$ SIMD 
    615           DO ij=ij_begin,ij_end          
    616              berni(ij,l) = pk(ij,l) 
    617              ! from now on pk contains the vertically-averaged geopotential 
    618              pk(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) 
    619           ENDDO 
    620        ENDDO 
    621  
    622     ELSE ! compressible 
    623  
    624        DO l=ll_begin,ll_end 
    625           !DIR$ SIMD 
    626           DO ij=ij_begin,ij_end          
    627              berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) 
    628           ENDDO 
    629        ENDDO 
    630  
    631     END IF ! Boussinesq/compressible 
    632  
    633 !!! u:=u+tau*du, du = gradients of Bernoulli and Exner functions 
    634     DO l=ll_begin,ll_end 
    635        !DIR$ SIMD 
    636        DO ij=ij_begin,ij_end          
    637  
    638           du(ij+u_right,l) = 1/de(ij+u_right) * (             & 
    639                0.5*(theta(ij,l)+theta(ij+t_right,l))                & 
    640                *( ne_right*pk(ij,l)+ne_left*pk(ij+t_right,l))    & 
    641                + ne_right*berni(ij,l)+ne_left*berni(ij+t_right,l) ) 
    642           u(ij+u_right,l) = u(ij+u_right,l) + tau*du(ij+u_right,l) 
    643  
    644  
    645           du(ij+u_lup,l) = 1/de(ij+u_lup) * (                  & 
    646                0.5*(theta(ij,l)+theta(ij+t_lup,l))                  & 
    647                *( ne_lup*pk(ij,l)+ne_rdown*pk(ij+t_lup,l))       & 
    648                + ne_lup*berni(ij,l)+ne_rdown*berni(ij+t_lup,l) ) 
    649           u(ij+u_lup,l) = u(ij+u_lup,l) + tau*du(ij+u_lup,l) 
    650  
    651           du(ij+u_ldown,l) = 1/de(ij+u_ldown) * (             & 
    652                0.5*(theta(ij,l)+theta(ij+t_ldown,l))                & 
    653                *( ne_ldown*pk(ij,l)+ne_rup*pk(ij+t_ldown,l))     & 
    654                + ne_ldown*berni(ij,l)+ne_rup*berni(ij+t_ldown,l) ) 
    655           u(ij+u_ldown,l) = u(ij+u_ldown,l) + tau*du(ij+u_ldown,l) 
    656  
    657        ENDDO 
    658     ENDDO 
    659  
    660     CALL trace_end("compute_caldyn_fast") 
    661  
    662   END SUBROUTINE compute_caldyn_fast 
    663  
    664   SUBROUTINE compute_caldyn_slow(u,rhodz,qu,theta, hflux,convm, dtheta_rhodz, du) 
    665     USE icosa 
    666     USE disvert_mod 
    667     USE exner_mod 
    668     USE trace 
    669     USE omp_para 
    670     IMPLICIT NONE 
    671     REAL(rstd),INTENT(IN)  :: u(iim*3*jjm,llm)    ! prognostic "velocity" 
    672     REAL(rstd),INTENT(IN)  :: rhodz(iim*jjm,llm) 
    673     REAL(rstd),INTENT(IN)  :: qu(iim*3*jjm,llm) 
    674     REAL(rstd),INTENT(IN)  :: theta(iim*jjm,llm)  ! potential temperature 
    675  
    676     REAL(rstd),INTENT(OUT) :: hflux(iim*3*jjm,llm) ! hflux in kg/s 
    677     REAL(rstd),INTENT(OUT) :: convm(iim*jjm,llm)  ! mass flux convergence 
    678     REAL(rstd),INTENT(OUT) :: dtheta_rhodz(iim*jjm,llm) 
    679     REAL(rstd),INTENT(OUT) :: du(iim*3*jjm,llm) 
    680  
    681     REAL(rstd) :: cor_NT(iim*jjm,llm)  ! NT coriolis force u.(du/dPhi) 
    682     REAL(rstd) :: urel(3*iim*jjm,llm)  ! relative velocity 
    683     REAL(rstd) :: Ftheta(3*iim*jjm,llm) ! theta flux 
    684     REAL(rstd) :: berni(iim*jjm,llm)  ! Bernoulli function 
    685  
    686     INTEGER :: i,j,ij,l 
    687     REAL(rstd) :: ww,uu 
    688  
    689     CALL trace_start("compute_caldyn_slow") 
    690  
    691     !    CALL wait_message(req_theta_rhodz)  
    692  
    693     DO l = ll_begin, ll_end 
    694 !!!  Compute mass and theta fluxes 
    695        IF (caldyn_conserv==energy) CALL test_message(req_qu)  
    696        !DIR$ SIMD 
    697        DO ij=ij_begin_ext,ij_end_ext          
    698           hflux(ij+u_right,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_right,l))*u(ij+u_right,l)*le(ij+u_right) 
    699           hflux(ij+u_lup,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l))*u(ij+u_lup,l)*le(ij+u_lup) 
    700           hflux(ij+u_ldown,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l))*u(ij+u_ldown,l)*le(ij+u_ldown) 
    701  
    702           Ftheta(ij+u_right,l)=0.5*(theta(ij,l)+theta(ij+t_right,l))*hflux(ij+u_right,l) 
    703           Ftheta(ij+u_lup,l)=0.5*(theta(ij,l)+theta(ij+t_lup,l))*hflux(ij+u_lup,l) 
    704           Ftheta(ij+u_ldown,l)=0.5*(theta(ij,l)+theta(ij+t_ldown,l))*hflux(ij+u_ldown,l) 
    705        ENDDO 
    706  
    707 !!! compute horizontal divergence of fluxes 
    708        !DIR$ SIMD 
    709        DO ij=ij_begin,ij_end          
    710           ! convm = -div(mass flux), sign convention as in Ringler et al. 2012, eq. 21 
    711           convm(ij,l)= -1./Ai(ij)*(ne_right*hflux(ij+u_right,l)   +  & 
    712                ne_rup*hflux(ij+u_rup,l)       +  &   
    713                ne_lup*hflux(ij+u_lup,l)       +  &   
    714                ne_left*hflux(ij+u_left,l)     +  &   
    715                ne_ldown*hflux(ij+u_ldown,l)   +  &   
    716                ne_rdown*hflux(ij+u_rdown,l))     
    717  
    718           ! signe ? attention d (rho theta dz) 
    719           ! dtheta_rhodz =  -div(flux.theta) 
    720           dtheta_rhodz(ij,l)=-1./Ai(ij)*(ne_right*Ftheta(ij+u_right,l)    +  & 
    721                ne_rup*Ftheta(ij+u_rup,l)        +  &   
    722                ne_lup*Ftheta(ij+u_lup,l)        +  &   
    723                ne_left*Ftheta(ij+u_left,l)      +  &   
    724                ne_ldown*Ftheta(ij+u_ldown,l)    +  &   
    725                ne_rdown*Ftheta(ij+u_rdown,l)) 
    726        ENDDO 
    727  
    728     END DO 
    729  
    730 !!! Compute potential vorticity (Coriolis) contribution to du 
    731  
    732     SELECT CASE(caldyn_conserv) 
    733     CASE(energy) ! energy-conserving TRiSK 
    734  
    735        CALL wait_message(req_qu) 
    736  
    737        DO l=ll_begin,ll_end 
    738           !DIR$ SIMD 
    739           DO ij=ij_begin,ij_end          
    740  
    741              uu = wee(ij+u_right,1,1)*hflux(ij+u_rup,l)*(qu(ij+u_right,l)+qu(ij+u_rup,l))+                             & 
    742                   wee(ij+u_right,2,1)*hflux(ij+u_lup,l)*(qu(ij+u_right,l)+qu(ij+u_lup,l))+                             & 
    743                   wee(ij+u_right,3,1)*hflux(ij+u_left,l)*(qu(ij+u_right,l)+qu(ij+u_left,l))+                           & 
    744                   wee(ij+u_right,4,1)*hflux(ij+u_ldown,l)*(qu(ij+u_right,l)+qu(ij+u_ldown,l))+                         & 
    745                   wee(ij+u_right,5,1)*hflux(ij+u_rdown,l)*(qu(ij+u_right,l)+qu(ij+u_rdown,l))+                         &  
    746                   wee(ij+u_right,1,2)*hflux(ij+t_right+u_ldown,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_ldown,l))+         & 
    747                   wee(ij+u_right,2,2)*hflux(ij+t_right+u_rdown,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_rdown,l))+         & 
    748                   wee(ij+u_right,3,2)*hflux(ij+t_right+u_right,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_right,l))+         & 
    749                   wee(ij+u_right,4,2)*hflux(ij+t_right+u_rup,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_rup,l))+             & 
    750                   wee(ij+u_right,5,2)*hflux(ij+t_right+u_lup,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_lup,l)) 
    751              du(ij+u_right,l) = .5*uu/de(ij+u_right) 
    752  
    753              uu = wee(ij+u_lup,1,1)*hflux(ij+u_left,l)*(qu(ij+u_lup,l)+qu(ij+u_left,l)) +                        & 
    754                   wee(ij+u_lup,2,1)*hflux(ij+u_ldown,l)*(qu(ij+u_lup,l)+qu(ij+u_ldown,l)) +                      & 
    755                   wee(ij+u_lup,3,1)*hflux(ij+u_rdown,l)*(qu(ij+u_lup,l)+qu(ij+u_rdown,l)) +                      & 
    756                   wee(ij+u_lup,4,1)*hflux(ij+u_right,l)*(qu(ij+u_lup,l)+qu(ij+u_right,l)) +                      & 
    757                   wee(ij+u_lup,5,1)*hflux(ij+u_rup,l)*(qu(ij+u_lup,l)+qu(ij+u_rup,l)) +                          &  
    758                   wee(ij+u_lup,1,2)*hflux(ij+t_lup+u_right,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_right,l)) +          & 
    759                   wee(ij+u_lup,2,2)*hflux(ij+t_lup+u_rup,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_rup,l)) +              & 
    760                   wee(ij+u_lup,3,2)*hflux(ij+t_lup+u_lup,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_lup,l)) +              & 
    761                   wee(ij+u_lup,4,2)*hflux(ij+t_lup+u_left,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_left,l)) +            & 
    762                   wee(ij+u_lup,5,2)*hflux(ij+t_lup+u_ldown,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_ldown,l)) 
    763              du(ij+u_lup,l) = .5*uu/de(ij+u_lup) 
    764  
    765  
    766              uu = wee(ij+u_ldown,1,1)*hflux(ij+u_rdown,l)*(qu(ij+u_ldown,l)+qu(ij+u_rdown,l)) +                      & 
    767                   wee(ij+u_ldown,2,1)*hflux(ij+u_right,l)*(qu(ij+u_ldown,l)+qu(ij+u_right,l)) +                      & 
    768                   wee(ij+u_ldown,3,1)*hflux(ij+u_rup,l)*(qu(ij+u_ldown,l)+qu(ij+u_rup,l)) +                          & 
    769                   wee(ij+u_ldown,4,1)*hflux(ij+u_lup,l)*(qu(ij+u_ldown,l)+qu(ij+u_lup,l)) +                          & 
    770                   wee(ij+u_ldown,5,1)*hflux(ij+u_left,l)*(qu(ij+u_ldown,l)+qu(ij+u_left,l)) +                        &  
    771                   wee(ij+u_ldown,1,2)*hflux(ij+t_ldown+u_lup,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_lup,l)) +          & 
    772                   wee(ij+u_ldown,2,2)*hflux(ij+t_ldown+u_left,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_left,l)) +        & 
    773                   wee(ij+u_ldown,3,2)*hflux(ij+t_ldown+u_ldown,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_ldown,l)) +      & 
    774                   wee(ij+u_ldown,4,2)*hflux(ij+t_ldown+u_rdown,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_rdown,l)) +      & 
    775                   wee(ij+u_ldown,5,2)*hflux(ij+t_ldown+u_right,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_right,l)) 
    776              du(ij+u_ldown,l) = .5*uu/de(ij+u_ldown) 
    777  
    778           ENDDO 
    779        ENDDO 
    780  
    781     CASE(enstrophy) ! enstrophy-conserving TRiSK 
    782  
    783        DO l=ll_begin,ll_end 
    784           !DIR$ SIMD 
    785           DO ij=ij_begin,ij_end          
    786  
    787              uu = wee(ij+u_right,1,1)*hflux(ij+u_rup,l)+                           & 
    788                   wee(ij+u_right,2,1)*hflux(ij+u_lup,l)+                           & 
    789                   wee(ij+u_right,3,1)*hflux(ij+u_left,l)+                          & 
    790                   wee(ij+u_right,4,1)*hflux(ij+u_ldown,l)+                         & 
    791                   wee(ij+u_right,5,1)*hflux(ij+u_rdown,l)+                         &  
    792                   wee(ij+u_right,1,2)*hflux(ij+t_right+u_ldown,l)+                 & 
    793                   wee(ij+u_right,2,2)*hflux(ij+t_right+u_rdown,l)+                 & 
    794                   wee(ij+u_right,3,2)*hflux(ij+t_right+u_right,l)+                 & 
    795                   wee(ij+u_right,4,2)*hflux(ij+t_right+u_rup,l)+                   & 
    796                   wee(ij+u_right,5,2)*hflux(ij+t_right+u_lup,l) 
    797              du(ij+u_right,l) = qu(ij+u_right,l)*uu/de(ij+u_right) 
    798  
    799  
    800              uu = wee(ij+u_lup,1,1)*hflux(ij+u_left,l)+                        & 
    801                   wee(ij+u_lup,2,1)*hflux(ij+u_ldown,l)+                       & 
    802                   wee(ij+u_lup,3,1)*hflux(ij+u_rdown,l)+                       & 
    803                   wee(ij+u_lup,4,1)*hflux(ij+u_right,l)+                       & 
    804                   wee(ij+u_lup,5,1)*hflux(ij+u_rup,l)+                         &  
    805                   wee(ij+u_lup,1,2)*hflux(ij+t_lup+u_right,l)+                 & 
    806                   wee(ij+u_lup,2,2)*hflux(ij+t_lup+u_rup,l)+                   & 
    807                   wee(ij+u_lup,3,2)*hflux(ij+t_lup+u_lup,l)+                   & 
    808                   wee(ij+u_lup,4,2)*hflux(ij+t_lup+u_left,l)+                  & 
    809                   wee(ij+u_lup,5,2)*hflux(ij+t_lup+u_ldown,l) 
    810              du(ij+u_lup,l) = qu(ij+u_lup,l)*uu/de(ij+u_lup) 
    811  
    812              uu = wee(ij+u_ldown,1,1)*hflux(ij+u_rdown,l)+                      & 
    813                   wee(ij+u_ldown,2,1)*hflux(ij+u_right,l)+                      & 
    814                   wee(ij+u_ldown,3,1)*hflux(ij+u_rup,l)+                        & 
    815                   wee(ij+u_ldown,4,1)*hflux(ij+u_lup,l)+                        & 
    816                   wee(ij+u_ldown,5,1)*hflux(ij+u_left,l)+                       &  
    817                   wee(ij+u_ldown,1,2)*hflux(ij+t_ldown+u_lup,l)+                & 
    818                   wee(ij+u_ldown,2,2)*hflux(ij+t_ldown+u_left,l)+               & 
    819                   wee(ij+u_ldown,3,2)*hflux(ij+t_ldown+u_ldown,l)+              & 
    820                   wee(ij+u_ldown,4,2)*hflux(ij+t_ldown+u_rdown,l)+              & 
    821                   wee(ij+u_ldown,5,2)*hflux(ij+t_ldown+u_right,l) 
    822              du(ij+u_ldown,l) = qu(ij+u_ldown,l)*uu/de(ij+u_ldown) 
    823  
    824           ENDDO 
    825        ENDDO 
    826  
    827     CASE DEFAULT 
    828        STOP 
    829     END SELECT 
    830  
    831     ! Compute bernouilli term = Kinetic Energy 
    832     IF(boussinesq) THEN 
    833  
    834        DO l=ll_begin,ll_end 
    835           !DIR$ SIMD 
    836           DO ij=ij_begin,ij_end          
    837  
    838              berni(ij,l) = & 
    839                   1/(4*Ai(ij))*(le(ij+u_right)*de(ij+u_right)*u(ij+u_right,l)**2 +    & 
    840                   le(ij+u_rup)*de(ij+u_rup)*u(ij+u_rup,l)**2 +          & 
    841                   le(ij+u_lup)*de(ij+u_lup)*u(ij+u_lup,l)**2 +          & 
    842                   le(ij+u_left)*de(ij+u_left)*u(ij+u_left,l)**2 +       & 
    843                   le(ij+u_ldown)*de(ij+u_ldown)*u(ij+u_ldown,l)**2 +    & 
    844                   le(ij+u_rdown)*de(ij+u_rdown)*u(ij+u_rdown,l)**2 )   
    845           ENDDO 
    846        ENDDO 
    847  
    848     ELSE ! compressible 
    849  
    850        DO l=ll_begin,ll_end 
    851           !DIR$ SIMD 
    852           DO ij=ij_begin,ij_end          
    853  
    854              berni(ij,l) =  & 
    855                   + 1/(4*Ai(ij))*(le(ij+u_right)*de(ij+u_right)*u(ij+u_right,l)**2 +    & 
    856                   le(ij+u_rup)*de(ij+u_rup)*u(ij+u_rup,l)**2 +          & 
    857                   le(ij+u_lup)*de(ij+u_lup)*u(ij+u_lup,l)**2 +          & 
    858                   le(ij+u_left)*de(ij+u_left)*u(ij+u_left,l)**2 +       & 
    859                   le(ij+u_ldown)*de(ij+u_ldown)*u(ij+u_ldown,l)**2 +    & 
    860                   le(ij+u_rdown)*de(ij+u_rdown)*u(ij+u_rdown,l)**2 )   
    861           ENDDO 
    862        ENDDO 
    863  
    864     END IF ! Boussinesq/compressible 
    865  
    866 !!! Add gradients of Bernoulli and Exner functions to du 
    867     DO l=ll_begin,ll_end 
    868        !DIR$ SIMD 
    869        DO ij=ij_begin,ij_end 
    870           du(ij+u_right,l) = du(ij+u_right,l) + 1/de(ij+u_right)            & 
    871                * ( ne_right*berni(ij,l)+ne_left*berni(ij+t_right,l) ) 
    872           du(ij+u_lup,l) = du(ij+u_lup,l) +  1/de(ij+u_lup)                 & 
    873                * ( ne_lup*berni(ij,l)+ne_rdown*berni(ij+t_lup,l) ) 
    874           du(ij+u_ldown,l) = du(ij+u_ldown,l) + 1/de(ij+u_ldown)            & 
    875                * ( ne_ldown*berni(ij,l)+ne_rup*berni(ij+t_ldown,l) ) 
    876        ENDDO 
    877     ENDDO 
    878  
    879     CALL trace_end("compute_caldyn_slow") 
    880  
    881   END SUBROUTINE compute_caldyn_slow 
    882  
    883   SUBROUTINE compute_caldyn_vert(u,theta,rhodz,convm, wflux,wwuu, dps,dtheta_rhodz,du) 
    884     USE icosa 
    885     USE disvert_mod 
    886     USE exner_mod 
    887     USE trace 
    888     USE omp_para 
    889     IMPLICIT NONE 
    890     REAL(rstd),INTENT(IN)  :: u(iim*3*jjm,llm) 
    891     REAL(rstd),INTENT(IN)  :: theta(iim*jjm,llm) 
    892     REAL(rstd),INTENT(IN)  :: rhodz(iim*jjm,llm) 
    893     REAL(rstd),INTENT(INOUT)  :: convm(iim*jjm,llm)  ! mass flux convergence 
    894  
    895     REAL(rstd),INTENT(INOUT) :: wflux(iim*jjm,llm+1) ! vertical mass flux (kg/m2/s) 
    896     REAL(rstd),INTENT(INOUT) :: wwuu(iim*3*jjm,llm+1) 
    897     REAL(rstd),INTENT(INOUT) :: du(iim*3*jjm,llm) 
    898     REAL(rstd),INTENT(INOUT) :: dtheta_rhodz(iim*jjm,llm) 
    899     REAL(rstd),INTENT(OUT) :: dps(iim*jjm) 
    900  
    901     ! temporary variable     
    902     INTEGER :: i,j,ij,l 
    903     REAL(rstd) :: p_ik, exner_ik 
    904     INTEGER    :: ij_omp_begin, ij_omp_end 
    905  
    906  
    907     CALL trace_start("compute_caldyn_vert") 
    908  
    909     CALL distrib_level(ij_end-ij_begin+1,ij_omp_begin,ij_omp_end) 
    910     ij_omp_begin=ij_omp_begin+ij_begin-1 
    911     ij_omp_end=ij_omp_end+ij_begin-1 
    912  
    913     !    REAL(rstd) :: wwuu(iim*3*jjm,llm+1) ! tmp var, don't know why but gain 30% on the whole code in opemp 
    914     ! need to be understood 
    915  
    916     !    wwuu=wwuu_out 
    917     CALL trace_start("compute_caldyn_vert") 
    918  
    919     !$OMP BARRIER    
    920 !!! cumulate mass flux convergence from top to bottom 
    921     !  IF (is_omp_level_master) THEN 
    922     DO  l = llm-1, 1, -1 
    923        !    IF (caldyn_conserv==energy) CALL test_message(req_qu)  
    924  
    925 !!$OMP DO SCHEDULE(STATIC)  
    926        !DIR$ SIMD 
    927        DO ij=ij_omp_begin,ij_omp_end          
    928           convm(ij,l) = convm(ij,l) + convm(ij,l+1) 
    929        ENDDO 
    930     ENDDO 
    931     !  ENDIF 
    932  
    933     !$OMP BARRIER 
    934     ! FLUSH on convm 
    935 !!!!!!!!!!!!!!!!!!!!!!!!! 
    936  
    937     ! compute dps 
    938     IF (is_omp_first_level) THEN 
    939        !DIR$ SIMD 
    940        DO ij=ij_begin,ij_end          
    941           ! dps/dt = -int(div flux)dz 
    942           dps(ij) = convm(ij,1) * g  
    943        ENDDO 
    944     ENDIF 
    945  
    946 !!! Compute vertical mass flux (l=1,llm+1 done by caldyn_BC) 
    947     DO l=ll_beginp1,ll_end 
    948        !    IF (caldyn_conserv==energy) CALL test_message(req_qu)  
    949        !DIR$ SIMD 
    950        DO ij=ij_begin,ij_end          
    951           ! w = int(z,ztop,div(flux)dz) + B(eta)dps/dt 
    952           ! => w>0 for upward transport 
    953           wflux( ij, l ) = bp(l) * convm( ij, 1 ) - convm( ij, l )  
    954        ENDDO 
    955     ENDDO 
    956  
    957     !--> flush wflux   
    958     !$OMP BARRIER 
    959  
    960     DO l=ll_begin,ll_endm1 
    961        !DIR$ SIMD 
    962        DO ij=ij_begin,ij_end          
    963           dtheta_rhodz(ij, l ) = dtheta_rhodz(ij, l )  -   0.5 * (  wflux(ij,l+1) * (theta(ij,l) + theta(ij,l+1))) 
    964        ENDDO 
    965     ENDDO 
    966  
    967     DO l=ll_beginp1,ll_end 
    968        !DIR$ SIMD 
    969        DO ij=ij_begin,ij_end          
    970           dtheta_rhodz(ij, l ) = dtheta_rhodz(ij, l )  +   0.5 * ( wflux(ij,l  ) * (theta(ij,l-1) + theta(ij,l) ) ) 
    971        ENDDO 
    972     ENDDO 
    973  
    974  
    975     ! Compute vertical transport 
    976     DO l=ll_beginp1,ll_end 
    977        !DIR$ SIMD 
    978        DO ij=ij_begin,ij_end          
    979           wwuu(ij+u_right,l) = 0.5*( wflux(ij,l) + wflux(ij+t_right,l)) * (u(ij+u_right,l) - u(ij+u_right,l-1)) 
    980           wwuu(ij+u_lup,l) = 0.5* ( wflux(ij,l) + wflux(ij+t_lup,l)) * (u(ij+u_lup,l) - u(ij+u_lup,l-1)) 
    981           wwuu(ij+u_ldown,l) = 0.5*( wflux(ij,l) + wflux(ij+t_ldown,l)) * (u(ij+u_ldown,l) - u(ij+u_ldown,l-1)) 
    982        ENDDO 
    983     ENDDO 
    984  
    985     !--> flush wwuu   
    986     !$OMP BARRIER 
    987  
    988     ! Add vertical transport to du 
    989     DO l=ll_begin,ll_end 
    990        !DIR$ SIMD 
    991        DO ij=ij_begin,ij_end          
    992           du(ij+u_right, l )   = du(ij+u_right,l)  - (wwuu(ij+u_right,l+1)+ wwuu(ij+u_right,l)) / (rhodz(ij,l)+rhodz(ij+t_right,l)) 
    993           du(ij+u_lup, l )     = du(ij+u_lup,l)    - (wwuu(ij+u_lup,l+1)  + wwuu(ij+u_lup,l))   / (rhodz(ij,l)+rhodz(ij+t_lup,l)) 
    994           du(ij+u_ldown, l )   = du(ij+u_ldown,l)  - (wwuu(ij+u_ldown,l+1)+ wwuu(ij+u_ldown,l)) / (rhodz(ij,l)+rhodz(ij+t_ldown,l)) 
    995        ENDDO 
    996     ENDDO 
    997  
    998     !  DO l=ll_beginp1,ll_end 
    999     !!DIR$ SIMD 
    1000     !      DO ij=ij_begin,ij_end          
    1001     !        wwuu_out(ij+u_right,l) = wwuu(ij+u_right,l) 
    1002     !        wwuu_out(ij+u_lup,l)   = wwuu(ij+u_lup,l) 
    1003     !        wwuu_out(ij+u_ldown,l) = wwuu(ij+u_ldown,l) 
    1004     !     ENDDO 
    1005     !   ENDDO 
    1006  
    1007     CALL trace_end("compute_caldyn_vert") 
    1008  
    1009   END SUBROUTINE compute_caldyn_vert 
    1010  
    1011353END MODULE caldyn_kernels_mod 

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

    r360 r362  
    1919  !$OMP THREADPRIVATE(nb_stage, matsuno_period, scheme, scheme_family) 
    2020 
    21  
    22   PUBLIC :: euler_scheme, accumulate_fluxes 
     21  PUBLIC :: euler_scheme, accumulate_fluxes, legacy_to_DEC, DEC_to_legacy 
    2322CONTAINS 
    2423 
    … …  
    143142  END SUBROUTINE accumulate_fluxes 
    144143 
     144  SUBROUTINE legacy_to_DEC(f_ps, f_u) 
     145    USE icosa 
     146    USE disvert_mod 
     147    USE omp_para 
     148    USE trace 
     149    TYPE(t_field),POINTER :: f_ps(:), f_u(:) 
     150    REAL(rstd), POINTER :: ps(:), u(:,:) 
     151    INTEGER :: ind,ij,l 
     152    CALL trace_start("legacy_to_DEC") 
     153 
     154    DO ind=1,ndomain 
     155       IF (.NOT. assigned_domain(ind)) CYCLE 
     156       CALL swap_dimensions(ind) 
     157       CALL swap_geometry(ind) 
     158 
     159       IF(caldyn_eta==eta_mass .AND. is_omp_first_level) THEN ! update ps 
     160          ps=f_ps(ind) 
     161          !$SIMD 
     162          DO ij=ij_begin,ij_end 
     163             ps(ij)=(ps(ij)-ptop)/g ! convert ps to column-integrated mass 
     164          ENDDO 
     165       END IF 
     166 
     167       u=f_u(ind) 
     168       DO l=ll_begin,ll_end 
     169          !$SIMD 
     170          DO ij=ij_begin,ij_end 
     171             u(ij+u_right,l)=u(ij+u_right,l)*de(ij+u_right) 
     172             u(ij+u_lup,l)=u(ij+u_lup,l)*de(ij+u_lup) 
     173             u(ij+u_ldown,l)=u(ij+u_ldown,l)*de(ij+u_ldown) 
     174          ENDDO 
     175       ENDDO 
     176    ENDDO 
     177 
     178    CALL trace_end("legacy_to_DEC")   
     179  END SUBROUTINE Legacy_to_DEC 
     180 
     181  SUBROUTINE DEC_to_legacy(f_ps, f_u) 
     182    USE icosa 
     183    USE disvert_mod 
     184    USE omp_para 
     185    USE trace 
     186    TYPE(t_field),POINTER :: f_ps(:), f_u(:) 
     187    REAL(rstd), POINTER :: ps(:), u(:,:) 
     188    INTEGER :: ind,ij,l 
     189    CALL trace_start("legacy_to_DEC") 
     190 
     191    DO ind=1,ndomain 
     192       IF (.NOT. assigned_domain(ind)) CYCLE 
     193       CALL swap_dimensions(ind) 
     194       CALL swap_geometry(ind) 
     195 
     196       IF(caldyn_eta==eta_mass .AND. is_omp_first_level) THEN 
     197          ps=f_ps(ind) 
     198          !$SIMD 
     199          DO ij=ij_begin,ij_end 
     200             ps(ij)=ptop+ps(ij)*g ! convert column-integrated mass to ps 
     201          ENDDO 
     202       ENDIF 
     203        
     204       u=f_u(ind) 
     205       DO l=ll_begin,ll_end 
     206          !$SIMD 
     207          DO ij=ij_begin,ij_end 
     208             u(ij+u_right,l)=u(ij+u_right,l)/de(ij+u_right) 
     209             u(ij+u_lup,l)=u(ij+u_lup,l)/de(ij+u_lup) 
     210             u(ij+u_ldown,l)=u(ij+u_ldown,l)/de(ij+u_ldown) 
     211          ENDDO 
     212       ENDDO 
     213    ENDDO 
     214 
     215    CALL trace_end("DEC_to_legacy")   
     216  END SUBROUTINE DEC_to_legacy 
     217 
    145218END MODULE euler_scheme_mod 

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

    r356 r362  
    2222    TYPE(t_field),POINTER :: de(:) 
    2323    TYPE(t_field),POINTER :: le(:) 
     24    TYPE(t_field),POINTER :: le_de(:) ! le/de, 0. if de=0. 
    2425    TYPE(t_field),POINTER :: Riv(:) 
    2526    TYPE(t_field),POINTER :: Riv2(:) 
    … …  
    6970  REAL(rstd),POINTER :: le(:)          ! lenght of a edge 
    7071!$OMP THREADPRIVATE(le) 
     72  REAL(rstd),POINTER :: le_de(:)          ! le/de 
     73!$OMP THREADPRIVATE(le_de) 
    7174  REAL(rstd),POINTER :: Riv(:,:)       ! weight 
    7275!$OMP THREADPRIVATE(Riv) 
    … …  
    116119    CALL allocate_field(geom%de,field_u,type_real) 
    117120    CALL allocate_field(geom%le,field_u,type_real) 
     121    CALL allocate_field(geom%le_de,field_u,type_real) 
    118122    CALL allocate_field(geom%bi,field_t,type_real) 
    119123    CALL allocate_field(geom%Av,field_z,type_real) 
    … …  
    150154    de=geom%de(ind) 
    151155    le=geom%le(ind) 
     156    le_de=geom%le_de(ind) 
    152157    Av=geom%Av(ind) 
    153158    Riv=geom%Riv(ind) 
    … …  
    438443          CALL dist_cart(xyz_v(n+z_up,:), xyz_v(n+z_lup,:),le(n+u_lup)) 
    439444          CALL dist_cart(xyz_v(n+z_ldown,:), xyz_v(n+z_down,:),le(n+u_ldown)) 
    440           
     445 
     446          le_de(n+u_right)=le(n+u_right)/de(n+u_right) ! NaN possible but should be harmless 
     447          le_de(n+u_lup)  =le(n+u_lup)  /de(n+u_lup) 
     448          le_de(n+u_ldown)=le(n+u_ldown)/de(n+u_ldown) 
     449 
    441450          Ai(n)=0 
    442451          DO k=0,5 
    … …  
    528537            CALL surf_triangle(xyz_i(n,:), xyz_v(n+z_pos(k+1),:), xyz_i(n+t_pos(k+1),:),S1) 
    529538            CALL surf_triangle(xyz_i(n,:), xyz_v(n+z_pos(k+1),:), xyz_i(n+t_pos(MOD(k+1+6,6)+1),:),S2) 
     539!            Riv(n,k+1)=0.5*(S1+S2)*(radius**2) ! Definition modified for DEC 
    530540            Riv(n,k+1)=0.5*(S1+S2)/Ai(n) 
    531541            Riv2(n,k+1)=0.5*(S1+S2)/surf_v(k+1) 

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

    r361 r362  
    11MODULE hevi_scheme_mod 
    2   USE euler_scheme_mod 
    32  USE prec 
    43  USE domain_mod 
    54  USE field_mod 
     5  USE euler_scheme_mod 
     6  USE caldyn_kernels_base_mod, ONLY : DEC 
    67  IMPLICIT NONE 
    78  PRIVATE 
    … …  
    5354    REAL(rstd),POINTER :: hflux(:,:),wflux(:,:),hfluxt(:,:),wfluxt(:,:) 
    5455 
     56    IF(DEC) CALL legacy_to_DEC(f_ps, f_u) 
    5557    DO j=1,nb_stage 
    5658       CALL caldyn_hevi((j==1) .AND. (MOD(it,itau_out)==0), taujj(j), & 
    … …  
    7981       END DO 
    8082    END DO 
     83    IF(DEC) CALL DEC_to_legacy(f_ps, f_u) 
    8184  END SUBROUTINE HEVI_scheme 
    8285   
    … …  
    8790    REAL(rstd), POINTER :: y(:,:), dy(:,:) 
    8891    INTEGER :: ind 
    89     DO ind=1,ndomain 
    90        IF (.NOT. assigned_domain(ind)) CYCLE 
    91        CALL swap_dimensions(ind) 
    92        dy=f_dy(ind); y=f_y(ind) 
    93        CALL compute_update(w,y,dy) 
    94     ENDDO 
     92    IF(w /= 0.) THEN 
     93       DO ind=1,ndomain 
     94          IF (.NOT. assigned_domain(ind)) CYCLE 
     95          CALL swap_dimensions(ind) 
     96          dy=f_dy(ind); y=f_y(ind) 
     97          CALL compute_update(w,y,dy) 
     98       ENDDO 
     99    END IF 
    95100  END SUBROUTINE update 
    96101     
Note: See TracChangeset for help on using the changeset viewer.