Changeset 50


Ignore:
Timestamp:
07/30/12 08:46:28 (12 years ago)
Author:
dubos
Message:

Temporary fix to problem with Exner function
Now also outputs velocity lat-lon, geopotential, pressure (at interfaces), Exner pressure, potential temp
Tested : test case DCMIP-3-1 with nbp=20

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

Legend:

Unmodified
Added
Removed

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

    r32 r50  
    22  USE icosa 
    33 
    4   PRIVATE  
    5   TYPE(t_field),POINTER :: f_out(:) 
    6   REAL(rstd),POINTER :: out(:,:) 
     4  INTEGER :: itau_out 
    75  TYPE(t_field),POINTER :: f_out_u(:) 
    86  REAL(rstd),POINTER :: out_u(:,:) 
    9   TYPE(t_field),POINTER :: f_out_z(:) 
    10   REAL(rstd),POINTER :: out_z(:,:) 
    11   
    12   INTEGER :: itau_out 
    13  
    14   PUBLIC init_caldyn, caldyn 
     7 
     8  TYPE(t_field),POINTER :: f_buf_i(:), f_buf_ulon(:), f_buf_ulat(:), f_buf_u3d(:) 
     9  TYPE(t_field),POINTER :: f_buf_v(:), f_buf_s(:), f_buf_p(:) 
     10 
     11  PUBLIC init_caldyn, caldyn, write_output_fields 
    1512 
    1613CONTAINS 
    1714   
    1815  SUBROUTINE init_caldyn(dt) 
    19   USE icosa 
    20   IMPLICIT NONE 
     16    USE icosa 
     17    IMPLICIT NONE 
    2118    REAL(rstd),INTENT(IN) :: dt 
    2219    INTEGER :: write_period 
    23  
    24     CALL allocate_caldyn 
     20     
     21    write_period=0 
    2522    CALL getin('write_period',write_period) 
    2623    write_period=write_period/scale_factor 
    27      
    2824    itau_out=INT(write_period/dt) 
    29      
     25   
    3026    CALL allocate_caldyn 
    3127   
    … …  
    3632  IMPLICIT NONE 
    3733 
    38     CALL allocate_field(f_out,field_t,type_real,llm)  
    39     CALL allocate_field(f_out_u,field_u,type_real,llm)  
    40     CALL allocate_field(f_out_z,field_z,type_real,llm)  
    41      
     34  CALL allocate_field(f_out_u,field_u,type_real,llm)  
     35   
     36  CALL allocate_field(f_buf_i,field_t,type_real,llm) 
     37  CALL allocate_field(f_buf_p,field_t,type_real,llm+1) 
     38  CALL allocate_field(f_buf_u3d,field_t,type_real,3,llm)  ! 3D vel at cell centers 
     39  CALL allocate_field(f_buf_ulon,field_t,type_real,llm) 
     40  CALL allocate_field(f_buf_ulat,field_t,type_real,llm) 
     41  CALL allocate_field(f_buf_v,field_z,type_real,llm) 
     42  CALL allocate_field(f_buf_s,field_t,type_real) 
     43 
    4244  END SUBROUTINE allocate_caldyn 
    4345   
    … …  
    4547  IMPLICIT NONE 
    4648    INTEGER,INTENT(IN) :: ind 
    47     out=f_out(ind)  
     49!    out=f_out(ind)  
    4850    out_u=f_out_u(ind)  
    49     out_z=f_out_z(ind)  
     51!    out_z=f_out_z(ind)  
    5052       
    5153  END SUBROUTINE swap_caldyn 
    … …  
    8890 
    8991  END SUBROUTINE check_mass_conservation 
    90    
    91    
    9292   
    9393  SUBROUTINE caldyn(it,f_phis, f_ps, f_theta_rhodz, f_u, f_dps, f_dtheta_rhodz, f_du) 
    … …  
    159159    ENDDO 
    160160 
    161     CALL transfert_request(f_out_u,req_e1) 
     161!    CALL transfert_request(f_out_u,req_e1) 
    162162!    CALL transfert_request(f_out_u,req_e1)  
    163163 
    … …  
    165165 
    166166    IF (mod(it,itau_out)==0 ) THEN 
    167       CALL writefield("ps",f_ps) 
    168       CALL writefield("dps",f_dps) 
    169 !      CALL writefield("theta_rhodz",f_theta_rhodz) 
    170 !    CALL kinetic(f_u,f_out) 
    171 !      CALL writefield("Ki",f_out) 
    172 !      CALL writefield("dtheta_rhodz",f_dtheta_rhodz) 
    173       CALL vorticity(f_u,f_out_z) 
    174       CALL writefield("vort",f_out_z) 
    175 !      CALL writefield("theta",f_out) 
    176       CALL theta_rhodz2temperature(f_ps,f_theta_rhodz,f_out) ; 
    177       CALL writefield("T",f_out) 
    178      
    179 !      CALL writefield("out",f_out) 
    180 !      DO ind=1,ndomain 
    181 !        CALL writefield("Ki",f_out,ind) 
    182 !        CALL writefield("vort",f_out_z,ind) 
    183 !        CALL writefield("dps",f_dps,ind) 
    184 !      ENDDO 
    185  
    186     ENDIF 
     167       PRINT *,'CALL write_output_fields' 
     168       CALL write_output_fields(f_ps, f_phis, f_dps, f_u, f_theta_rhodz, & 
     169            f_buf_i, f_buf_v, f_buf_u3d, f_buf_ulon, f_buf_ulat, f_buf_s, f_buf_p) 
     170    END IF 
     171 
     172 
    187173!    CALL check_mass_conservation(f_ps,f_dps) 
    188174 
    … …  
    193179  USE icosa 
    194180  USE disvert_mod 
     181  USE exner_mod 
    195182  IMPLICIT NONE 
    196183    REAL(rstd),INTENT(IN) :: phis(iim*jjm) 
    … …  
    297284    ENDDO 
    298285       
    299  
    300  
    301   
    302  !!! Compute Exnher function 
    303    
    304  !! Compute Alpha and Beta 
    305    
    306  ! for llm layer 
    307 !$OMP DO 
    308    DO j=jj_begin-1,jj_end+1 
    309      DO i=ii_begin-1,ii_end+1 
    310        ij=(j-1)*iim+i 
    311        alpha(ij,llm) = 0. 
    312        beta (ij,llm) = 1./ (1+ 2*kappa) 
    313      ENDDO 
    314    ENDDO 
    315     
    316  ! for other layer    
    317    DO l = llm-1 , 2 , -1 
    318 !$OMP DO 
    319      DO j=jj_begin-1,jj_end+1 
    320        DO i=ii_begin-1,ii_end+1 
    321           ij=(j-1)*iim+i 
    322           delta = p(ij,l)* (1+2*kappa) + p(ij,l+1)* ( beta(ij,l+1)- (1+2*kappa) ) 
    323           alpha(ij,l)  = - p(ij,l+1) / delta * alpha(ij,l+1) 
    324           beta (ij,l)  =   p(ij,l  ) / delta    
    325        ENDDO 
    326      ENDDO 
    327    ENDDO 
    328     
    329  !! Compute pk 
    330   
    331  ! for first layer 
    332 !$OMP DO 
    333    DO j=jj_begin-1,jj_end+1 
    334      DO i=ii_begin-1,ii_end+1 
    335         ij=(j-1)*iim+i 
    336         pks(ij) = cpp * ( ps(ij)/preff ) ** kappa 
    337         pk(ij,1) = ( p(ij,1)*pks(ij) - 0.5*alpha(ij,2)*p(ij,2) )     /               & 
    338                    (  p(ij,1)* (1.+kappa) + 0.5*( beta(ij,2)-(1.+2*kappa) )* p(ij,2) ) 
    339      ENDDO 
    340    ENDDO 
    341     
    342  ! for other layers 
    343       
    344    DO l = 2, llm 
    345 !$OMP DO 
    346      DO j=jj_begin-1,jj_end+1 
    347        DO i=ii_begin-1,ii_end+1 
    348          ij=(j-1)*iim+i 
    349          pk(ij,l) = alpha(ij,l) + beta(ij,l) * pk(ij,l-1) 
    350        ENDDO 
    351      ENDDO 
    352    ENDDO 
    353     
    354    
     286 !!! Compute Exner function 
     287    CALL compute_exner(ps,p,pks,pk,1)  
     288 
    355289!!! Compute mass 
    356290   DO l = 1, llm 
    … …  
    484418   
    485419 
    486   DO l = 1, llm 
    487 !$OMP DO 
    488    DO j=jj_begin,jj_end 
    489     DO i=ii_begin,ii_end 
    490       ij=(j-1)*iim+i 
    491       out(ij,l)=theta(ij,l)-288 
    492     ENDDO 
    493   ENDDO 
    494  ENDDO 
     420!  DO l = 1, llm 
     421!!$OMP DO 
     422!   DO j=jj_begin,jj_end 
     423!    DO i=ii_begin,ii_end 
     424!      ij=(j-1)*iim+i 
     425!      out(ij,l)=theta(ij,l)-288 
     426!    ENDDO 
     427!  ENDDO 
     428! ENDDO 
    495429 
    496430 
    … …  
    680614    ENDDO 
    681615  ENDDO 
     616 
    682617!!! contribution due to vertical advection 
    683618  
    … …  
    743678  END SUBROUTINE compute_caldyn 
    744679   
    745    
     680  SUBROUTINE write_output_fields(f_ps, f_phis, f_dps, f_u, f_theta_rhodz, & 
     681       f_buf_i, f_buf_v, f_buf_i3, f_buf1_i, f_buf2_i, f_buf_s, f_buf_p) 
     682    USE icosa 
     683    USE vorticity_mod 
     684    USE theta2theta_rhodz_mod 
     685    USE pression_mod 
     686    USE write_field 
     687    TYPE(t_field),POINTER :: f_ps(:), f_phis(:), f_u(:), f_theta_rhodz(:), f_dps(:), & 
     688         f_buf_i(:), f_buf_v(:), f_buf_i3(:), f_buf1_i(:), f_buf2_i(:), f_buf_s(:), f_buf_p(:) 
     689     
     690    CALL writefield("dps",f_dps) 
     691    CALL vorticity(f_u,f_buf_v) 
     692    CALL writefield("vort",f_buf_v) 
     693     
     694    ! Temperature 
     695    CALL theta_rhodz2temperature(f_ps,f_theta_rhodz,f_buf_i) ; 
     696    CALL writefield("T",f_buf_i) 
     697     
     698    ! velocity components 
     699    CALL un2ulonlat(f_u, f_buf_i3, f_buf1_i, f_buf2_i) 
     700    CALL writefield("ulon",f_buf1_i) 
     701    CALL writefield("ulat",f_buf2_i) 
     702     
     703    ! geopotential 
     704    CALL thetarhodz2geopot(f_ps,f_phis,f_theta_rhodz, f_buf_s,f_buf_p,f_buf1_i,f_buf2_i,f_buf_i) 
     705    CALL writefield("p",f_buf_p) 
     706    CALL writefield("phi",f_buf_i) 
     707    CALL writefield("theta",f_buf1_i) ! potential temperature 
     708    CALL writefield("pk",f_buf2_i) ! Exner pressure 
     709   
     710  END SUBROUTINE write_output_fields 
     711   
     712  SUBROUTINE thetarhodz2geopot(f_ps,f_phis,f_theta_rhodz, f_pks,f_p,f_theta,f_pk,f_phi)  
     713    USE field_mod 
     714    USE pression_mod 
     715    USE exner_mod 
     716    USE geopotential_mod 
     717    USE theta2theta_rhodz_mod 
     718    TYPE(t_field), POINTER :: f_ps(:), f_phis(:), f_theta_rhodz(:), &  ! IN 
     719         f_pks(:), f_p(:), f_theta(:), f_pk(:), f_phi(:)               ! OUT 
     720    REAL(rstd),POINTER :: pk(:,:), p(:,:), theta(:,:), theta_rhodz(:,:), & 
     721         phi(:,:), phis(:), ps(:), pks(:) 
     722    INTEGER :: ind 
     723 
     724    DO ind=1,ndomain 
     725       CALL swap_dimensions(ind) 
     726       CALL swap_geometry(ind) 
     727       ps = f_ps(ind) 
     728       p  = f_p(ind) 
     729       CALL compute_pression(ps,p,0) 
     730       pk = f_pk(ind) 
     731       pks = f_pks(ind) 
     732       CALL compute_exner(ps,p,pks,pk,0) 
     733       theta_rhodz = f_theta_rhodz(ind) 
     734       theta = f_theta(ind) 
     735       CALL compute_theta_rhodz2theta(ps, theta_rhodz,theta,0) 
     736       phis = f_phis(ind) 
     737       phi = f_phi(ind) 
     738       CALL compute_geopotential(phis,pks,pk,theta,phi,0) 
     739    END DO 
     740 
     741  END SUBROUTINE thetarhodz2geopot 
     742   
     743  SUBROUTINE un2ulonlat(f_u, f_u3d, f_ulon, f_ulat) 
     744    USE field_mod 
     745    USE wind_mod 
     746    TYPE(t_field), POINTER :: f_u(:), & ! IN  : normal velocity components on edges 
     747         f_u3d(:), f_ulon(:), f_ulat(:) ! OUT : velocity reconstructed at hexagons 
     748    REAL(rstd),POINTER :: u(:,:), u3d(:,:,:), ulon(:,:), ulat(:,:) 
     749    INTEGER :: ind 
     750    DO ind=1,ndomain 
     751       CALL swap_dimensions(ind) 
     752       CALL swap_geometry(ind) 
     753       u=f_u(ind) 
     754       u3d=f_u3d(ind) 
     755       CALL compute_wind_centered(u,u3d) 
     756       ulon=f_ulon(ind) 
     757       ulat=f_ulat(ind) 
     758       CALL compute_wind_centered_lonlat_compound(u3d, ulon, ulat) 
     759    END DO 
     760  END SUBROUTINE un2ulonlat 
    746761  
    747762END MODULE caldyn_gcm_mod 

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

    r19 r50  
    4343    REAL(rstd) :: delta  
    4444     
     45    IF(.FALSE.) THEN          ! LMD-Z style calculation of Exner pressure 
     46       !! Compute Alpha and Beta 
     47        
     48       ! for llm layer 
     49       !$OMP DO 
     50       DO j=jj_begin-offset,jj_end+offset 
     51          DO i=ii_begin-offset,ii_end+offset 
     52             ij=(j-1)*iim+i 
     53             alpha(ij,llm) = 0. 
     54             beta (ij,llm) = 1./ (1+ 2*kappa) 
     55          ENDDO 
     56       ENDDO 
     57        
     58       ! for other layer    
     59       DO l = llm-1 , 2 , -1 
     60          !$OMP DO 
     61          DO j=jj_begin-offset,jj_end+offset 
     62             DO i=ii_begin-offset,ii_end+offset 
     63                ij=(j-1)*iim+i 
     64                delta = p(ij,l)* (1+2*kappa) + p(ij,l+1)* ( beta(ij,l+1)- (1+2*kappa) ) 
     65                alpha(ij,l)  = - p(ij,l+1) / delta * alpha(ij,l+1) 
     66                beta (ij,l)  =   p(ij,l  ) / delta    
     67             ENDDO 
     68          ENDDO 
     69       ENDDO 
     70        
     71       !! Compute pk 
     72        
     73       ! for first layer 
     74       !$OMP DO 
     75       DO j=jj_begin-offset,jj_end+offset 
     76          DO i=ii_begin-offset,ii_end+offset 
     77             ij=(j-1)*iim+i 
     78             pks(ij) = cpp * ( ps(ij)/preff ) ** kappa 
     79             pk(ij,1) = ( p(ij,1)*pks(ij) - 0.5*alpha(ij,2)*p(ij,2) )     /    & 
     80                  (  p(ij,1)* (1.+kappa) + 0.5*( beta(ij,2)-(1.+2*kappa) )* p(ij,2) ) 
     81          ENDDO 
     82       ENDDO 
     83        
     84       ! for other layers 
     85        
     86       DO l = 2, llm 
     87          !$OMP DO 
     88          DO j=jj_begin-offset,jj_end+offset 
     89             DO i=ii_begin-offset,ii_end+offset 
     90                ij=(j-1)*iim+i 
     91                pk(ij,l) = alpha(ij,l) + beta(ij,l) * pk(ij,l-1) 
     92             ENDDO 
     93          ENDDO 
     94       ENDDO 
     95        
     96    ELSE ! Simple calculation of Exner pressure based on centered average 
     97       ! surface : pks 
     98       !$OMP DO 
     99       DO j=jj_begin-offset,jj_end+offset 
     100          DO i=ii_begin-offset,ii_end+offset 
     101             ij=(j-1)*iim+i 
     102             pks(ij) = cpp * ( ps(ij)/preff ) ** kappa 
     103          ENDDO 
     104       ENDDO 
    45105 
    46    !! Compute Alpha and Beta 
    47    
    48    ! for llm layer 
    49      DO j=jj_begin-offset,jj_end+offset 
    50        DO i=ii_begin-offset,ii_end+offset 
    51          ij=(j-1)*iim+i 
    52          alpha(ij,llm) = 0. 
    53          beta (ij,llm) = 1./ (1+ 2*kappa) 
     106       ! 3D : pk 
     107       DO l = 1, llm 
     108          !$OMP DO 
     109          DO j=jj_begin-offset,jj_end+offset 
     110             DO i=ii_begin-offset,ii_end+offset 
     111                ij=(j-1)*iim+i 
     112                pk(ij,l) = cpp * ((.5/preff)*(p(ij,l)+p(ij,l+1))) ** kappa 
     113             ENDDO 
     114          ENDDO 
    54115       ENDDO 
    55      ENDDO 
    56     
    57    ! for other layer    
    58      DO l = llm-1 , 2 , -1 
    59        DO j=jj_begin-offset,jj_end+offset 
    60          DO i=ii_begin-offset,ii_end+offset 
    61             ij=(j-1)*iim+i 
    62             delta = p(ij,l)* (1+2*kappa) + p(ij,l+1)* ( beta(ij,l+1)- (1+2*kappa) ) 
    63             alpha(ij,l)  = - p(ij,l+1) / delta * alpha(ij,l+1) 
    64             beta (ij,l)  =   p(ij,l  ) / delta    
    65          ENDDO 
    66        ENDDO 
    67      ENDDO      
    68        
    69    !! Compute pk 
    70   
    71    ! for first layer 
    72      DO j=jj_begin-offset,jj_end+offset 
    73        DO i=ii_begin-offset,ii_end+offset 
    74           ij=(j-1)*iim+i 
    75           pks(ij) = cpp * ( ps(ij)/preff ) ** kappa 
    76           pk(ij,1) = ( p(ij,1)*pks(ij) - 0.5*alpha(ij,2)*p(ij,2) )     /               & 
    77                      (  p(ij,1)* (1.+kappa) + 0.5*( beta(ij,2)-(1.+2*kappa) )* p(ij,2) ) 
    78        ENDDO 
    79      ENDDO 
    80     
    81    ! for other layers 
    82       
    83      DO l = 2, llm 
    84        DO j=jj_begin-offset,jj_end+offset 
    85          DO i=ii_begin-offset,ii_end+offset 
    86            ij=(j-1)*iim+i 
    87            pk(ij,l) = alpha(ij,l) + beta(ij,l) * pk(ij,l-1) 
    88          ENDDO 
    89        ENDDO 
    90      ENDDO 
     116    END IF 
    91117     
    92118  END SUBROUTINE compute_exner 
    93  
     119   
    94120END MODULE exner_mod 

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

    r46 r50  
    11MODULE timeloop_gcm_mod 
    2  
    32   
    43CONTAINS 
    … …  
    87  USE dissip_gcm_mod 
    98  USE caldyn_mod 
    10   USE theta2theta_rhodz_mod 
    119  USE etat0_mod 
    1210  USE guided_mod 
    … …  
    3331  REAL(rstd),POINTER :: du(:,:), dum1(:,:), dum2(:,:) 
    3432  REAL(rstd),POINTER :: dtheta_rhodz(:,:),dtheta_rhodzm1(:,:),dtheta_rhodzm2(:,:) 
     33 
    3534  REAL(rstd) :: dt, run_length 
    3635  INTEGER :: ind 
    … …  
    3938  INTEGER :: matsuno_period 
    4039  INTEGER :: itaumax 
    41   INTEGER :: write_period   
    42   INTEGER :: itau_out 
    4340 
    4441  dt=90. 
    … …  
    5451  dt=dt/scale_factor 
    5552 
    56   write_period=0 
    57   CALL getin('write_period',write_period) 
    58   write_period=write_period/scale_factor 
    59    
    60   itau_out=INT(write_period/dt) 
    61    
    6253  scheme='adam_bashforth' 
    6354  CALL getin('scheme',scheme) 
    … …  
    10192   
    10293  CALL etat0(f_ps,f_phis,f_theta_rhodz,f_u, f_q) 
    103    
     94 
    10495  DO it=0,itaumax 
    10596    PRINT *,"It No :",It,"   t :",dt*It 
    … …  
    242233        dpsm1(:)=dps(:) ; dum1(:,:)=du(:,:) ; dtheta_rhodzm1(:,:)=dtheta_rhodz(:,:) 
    243234 
    244       ENDDO 
    245       
     235      ENDDO      
    246236       
    247237    END SUBROUTINE adam_bashforth_scheme 
    248238 
    249   END SUBROUTINE timeloop 
    250      
    251      
     239  END SUBROUTINE timeloop     
    252240   
    253241END MODULE timeloop_gcm_mod 
Note: See TracChangeset for help on using the changeset viewer.