Changeset 17


Ignore:
Timestamp:
07/16/12 10:24:35 (12 years ago)
Author:
ymipsl
Message:

Merge advection scheme from sarvesh in standard version

YM

Location:
codes/icosagcm/trunk/src
Files:
9 added
10 edited

Legend:

Unmodified
Added
Removed

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

    r15 r17  
    44  USE transfert_mod 
    55 
     6  PRIVATE  
    67  TYPE(t_field),POINTER :: f_out(:) 
    78  REAL(rstd),POINTER :: out(:,:) 
    … …  
    1213  
    1314  INTEGER :: itau_out 
    14    
     15 
     16  PUBLIC init_caldyn, caldyn 
     17 
    1518CONTAINS 
    1619   
    … …  
    1922  IMPLICIT NONE 
    2023    REAL(rstd),INTENT(IN) :: dt 
    21     REAL(rstd) :: write_period 
     24    INTEGER :: write_period 
     25 
    2226    CALL allocate_caldyn 
    23  
    2427    CALL getin('write_period',write_period) 
    2528     
    2629    itau_out=INT(write_period/dt) 
     30     
     31    CALL allocate_caldyn 
    2732   
    2833  END SUBROUTINE init_caldyn 
    … …  
    3439  USE metric 
    3540  IMPLICIT NONE 
    36     INTEGER :: ind,i,j 
    3741 
    3842    CALL allocate_field(f_out,field_t,type_real,llm)  

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

    r12 r17  
    11MODULE disvert_mod 
    22  USE prec 
    3   REAL(rstd), SAVE, ALLOCATABLE :: ap(:) 
    4   REAL(rstd), SAVE, ALLOCATABLE :: bp(:) 
    5   REAL(rstd), SAVE, ALLOCATABLE :: presnivs(:) 
     3  REAL(rstd), SAVE, POINTER :: ap(:) 
     4  REAL(rstd), SAVE, POINTER :: bp(:) 
     5  REAL(rstd), SAVE, POINTER :: presnivs(:) 
    66 
    77CONTAINS 
    88 
    99  SUBROUTINE init_disvert 
    10   USE grid_param 
     10  USE disvert_std_mod, ONLY: ap_std=>ap, bp_std=>bp, presnivs_std=>presnivs, init_disvert_std=>init_disvert 
     11  USE disvert_ncar_mod, ONLY: ap_ncar=>ap, bp_ncar=>bp, presnivs_ncar=>presnivs, init_disvert_ncar=>init_disvert 
     12  USE ioipsl 
    1113  IMPLICIT NONE 
    12    
    13     ALLOCATE(ap(llm+1)) 
    14     ALLOCATE(bp(llm+1)) 
    15     ALLOCATE(presnivs(llm)) 
     14    CHARACTER(LEN=255) :: disvert_type = 'std' 
     15     
     16    CALL getin("disvert",disvert_type) 
     17     
     18    SELECT CASE (TRIM(disvert_type)) 
     19      CASE('std') 
     20     
     21        CALL init_disvert_std 
     22        ap=>ap_std 
     23        bp=>bp_std 
     24        presnivs=>presnivs_std 
     25       
     26      CASE ('ncar') 
    1627 
    17     CALL disvert(ap,bp,presnivs)     
     28        CALL init_disvert_ncar 
     29        ap=>ap_ncar 
     30        bp=>bp_ncar 
     31        presnivs=>presnivs_ncar 
     32         
     33      CASE default 
     34        PRINT*,'Bad selector for variable disvert : <', TRIM(disvert_type),"> options are <std>, <ncar>"  
     35        STOP 
     36         
     37    END SELECT 
    1838 
    1939  END SUBROUTINE init_disvert   
    20  
    21  
    22   SUBROUTINE disvert(ap,bp,presnivs) 
    23   USE earth_const 
    24   USE math_const  
    25   USE grid_param 
    26   IMPLICIT NONE 
    27   REAL(rstd),INTENT(OUT) :: ap(:) 
    28   REAL(rstd),INTENT(OUT) :: bp(:) 
    29   REAL(rstd),INTENT(OUT) :: presnivs(:) 
    30    
    31   REAL(rstd) :: dsig(llm) 
    32   REAL(rstd) :: sig(llm+1) 
    33   REAL(rstd) :: snorm 
    34   INTEGER :: l 
    35    
    36     snorm  = 0. 
    37     DO l = 1, llm 
    38       dsig(l) = 1.0 + 7.0 * SIN( Pi*(l-0.5)/(llm+1) )**2 
    39       snorm = snorm + dsig(l) 
    40     ENDDO     
    41      
    42     DO l = 1, llm 
    43       dsig(l) = dsig(l)/snorm 
    44     ENDDO 
    45  
    46     sig(llm+1) = 0. 
    47     DO l = llm, 1, -1 
    48       sig(l) = sig(l+1) + dsig(l) 
    49     ENDDO 
    50  
    51     bp(llm+1) =   0. 
    52     DO l = 1, llm 
    53       bp(l) = EXP( 1. -1./( sig(l)*sig(l)) ) 
    54       ap(l) = pa * ( sig(l) - bp(l) ) 
    55     ENDDO 
    56     bp(1)=1. 
    57     ap(1)=0. 
    58     ap(llm+1) = pa * ( sig(llm+1) - bp(llm+1) ) 
    59      
    60     PRINT*,'ap',ap 
    61     PRINT*,'bp',bp 
    62      
    63     PRINT*, 'Niveaux de pressions approximatifs aux centres des' 
    64     PRINT*, 'couches calcules pour une pression de surface =', preff 
    65     PRINT*, 'et altitudes equivalentes pour une hauteur d echelle de' 
    66     PRINT*, '8km' 
    67      
    68     DO l = 1, llm 
    69       presnivs(l) = 0.5 *( ap(l)+bp(l)*preff + ap(l+1)+bp(l+1)*preff ) 
    70    
    71       PRINT*, 'PRESNIVS(',l,')=',presnivs(l),'  Z ~ ',log(preff/presnivs(l))*8.,       & 
    72               ' DZ ~ ',8.*log((ap(l)+bp(l)*preff)/ max(ap(l+1)+bp(l+1)*preff,1.e-10)) 
    73     ENDDO 
    74    
    75   END SUBROUTINE disvert 
    7640   
    7741END MODULE disvert_mod 

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

    r12 r17  
    11MODULE etat0_mod 
    22 
    3   USE etat0_academic_mod 
    4   USE etat0_jablonowsky06_mod 
    5   USE etat0_williamson_mod 
     3 
     4CONTAINS 
     5   
     6  SUBROUTINE etat0(f_ps,f_phis,f_theta_rhodz,f_u, f_q) 
     7  USE field_mod 
     8  USE domain_mod 
     9  USE domain_mod 
     10  USE dimensions 
     11  USE grid_param 
     12  USE geometry 
     13  USE etat0_jablonowsky06_mod, ONLY : etat0_jablonowsky06=>etat0 
     14  USE etat0_academic_mod, ONLY : etat0_academic=>etat0   
     15  USE etat0_ncar_mod, ONLY : etat0_ncar=>etat0   
     16  USE ioipsl 
     17  IMPLICIT NONE 
     18    TYPE(t_field),POINTER :: f_ps(:) 
     19    TYPE(t_field),POINTER :: f_phis(:) 
     20    TYPE(t_field),POINTER :: f_theta_rhodz(:) 
     21    TYPE(t_field),POINTER :: f_u(:) 
     22    TYPE(t_field),POINTER :: f_q(:) 
     23 
     24    CHARACTER(len=255) :: etat0_type 
     25     
     26    etat0_type='jablonowsky06' 
     27    CALL getin("etat0",etat0_type) 
     28     
     29    SELECT CASE (TRIM(etat0_type)) 
     30      CASE ('jablonowsky06') 
     31        CALL etat0_jablonowsky06(f_ps,f_phis,f_theta_rhodz,f_u, f_q) 
     32      CASE ('academic') 
     33        CALL etat0_academic(f_ps,f_phis,f_theta_rhodz,f_u, f_q) 
     34      CASE ('ncar') 
     35        CALL etat0_ncar(f_ps,f_phis,f_theta_rhodz,f_u, f_q) 
     36      CASE DEFAULT 
     37        PRINT*, 'Bad selector for varaiable etat0 <',etat0_type,'> options are <jablonowsky06>, <academic>, <ncar> ' 
     38        STOP 
     39    END SELECT  
     40 
     41  END SUBROUTINE etat0   
     42          
    643END MODULE etat0_mod 

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

    r13 r17  
    1919    TYPE(t_field),POINTER :: f_theta_rhodz(:) 
    2020    TYPE(t_field),POINTER :: f_u(:) 
     21    TYPE(t_field),POINTER :: f_q(:) 
    2122    TYPE(t_field),POINTER :: f_Ki(:) 
    2223    TYPE(t_field),POINTER :: f_temp(:) 
    … …  
    3435    CALL allocate_field(f_temp,field_t,type_real) 
    3536     
    36     CALL etat0_academic(f_ps,f_phis,f_theta_rhodz,f_u) 
     37    CALL etat0(f_ps,f_phis,f_theta_rhodz,f_u, f_q) 
    3738 
    3839    CALL kinetic(f_u,f_Ki) 
    … …  
    4647    
    4748     
    48   SUBROUTINE etat0_academic(f_ps,f_phis,f_theta_rhodz,f_u) 
     49  SUBROUTINE etat0(f_ps,f_phis,f_theta_rhodz,f_u, f_q) 
    4950  USE field_mod 
    5051  USE domain_mod 
    … …  
    5859    TYPE(t_field),POINTER :: f_theta_rhodz(:) 
    5960    TYPE(t_field),POINTER :: f_u(:) 
     61    TYPE(t_field),POINTER :: f_q(:) 
    6062   
    6163    REAL(rstd),POINTER :: ps(:) 
    … …  
    7577    ENDDO 
    7678 
    77   END SUBROUTINE etat0_academic  
     79  END SUBROUTINE etat0 
    7880   
    7981  SUBROUTINE compute_etat0_academic(ps, phis, theta_rhodz, u) 

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

    r15 r17  
    1212  REAL(rstd),PARAMETER :: Gamma=0.005 
    1313  REAL(rstd),PARAMETER :: up0=1 
    14   PUBLIC  test_etat0_jablonowsky06, etat0_jablonowsky06, compute_etat0_jablonowsky06 
     14  PUBLIC  test_etat0_jablonowsky06, etat0, compute_etat0_jablonowsky06 
    1515CONTAINS 
    1616   
    … …  
    3939    TYPE(t_field),POINTER :: f_phi(:) 
    4040    TYPE(t_field),POINTER :: f_vort(:) 
     41    TYPE(t_field),POINTER :: f_q(:) 
    4142   
    4243    REAL(rstd),POINTER :: Ki(:,:) 
    … …  
    5758    CALL allocate_field(f_vort,field_z,type_real,llm) 
    5859     
    59     CALL etat0_jablonowsky06(f_ps,f_phis,f_theta_rhodz,f_u) 
     60    CALL etat0(f_ps,f_phis,f_theta_rhodz,f_u,f_q) 
    6061 
    6162    CALL kinetic(f_u,f_Ki) 
    … …  
    7778    
    7879     
    79   SUBROUTINE etat0_jablonowsky06(f_ps,f_phis,f_theta_rhodz,f_u) 
     80  SUBROUTINE etat0(f_ps,f_phis,f_theta_rhodz,f_u, f_q) 
    8081  USE field_mod 
    8182  USE domain_mod 
    … …  
    8990    TYPE(t_field),POINTER :: f_theta_rhodz(:) 
    9091    TYPE(t_field),POINTER :: f_u(:) 
     92    TYPE(t_field),POINTER :: f_q(:) 
    9193   
    9294    REAL(rstd),POINTER :: ps(:) 
    … …  
    9496    REAL(rstd),POINTER :: theta_rhodz(:,:) 
    9597    REAL(rstd),POINTER :: u(:,:) 
     98    REAL(rstd),POINTER :: q(:,:,:) 
    9699    INTEGER :: ind 
    97100     
    … …  
    103106      theta_rhodz=f_theta_rhodz(ind) 
    104107      u=f_u(ind) 
     108      q=f_q(ind) 
     109      q=0 
    105110      CALL compute_etat0_jablonowsky06(ps, phis, theta_rhodz, u) 
    106111    ENDDO 
    107112 
    108   END SUBROUTINE etat0_jablonowsky06  
     113  END SUBROUTINE etat0 
    109114   
    110115  SUBROUTINE compute_etat0_jablonowsky06(ps, phis, theta_rhodz, u) 

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

    r15 r17  
    115115  USE transfert_mod 
    116116  USE vector 
     117  USE ioipsl 
    117118  IMPLICIT NONE 
    118     INTEGER,PARAMETER :: nb_it=3000 
     119    INTEGER :: nb_it=0 
    119120    TYPE(t_domain),POINTER :: d 
    120121    INTEGER :: ind,it,i,j,n,k 
    … …  
    124125    REAL(rstd) :: sum 
    125126    LOGICAL    :: check 
     127     
     128     
     129    CALL getin('optim_it',nb_it) 
    126130     
    127131    DO ind=1,ndomain 

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

    r15 r17  
    44  INTEGER,PARAMETER  :: nb_face=10 
    55  INTEGER  :: llm=19 
     6  INTEGER  :: nqtot 
    67 
    78CONTAINS 
    … …  
    1112  IMPLICIT NONE 
    1213    CALL getin('nbp',iim_glo) 
    13     jjm_glo=40 
     14    jjm_glo=iim_glo 
    1415    CALL getin('llm',llm) 
     16     
     17    nqtot=1 
     18    CALL getin('nqtot',nqtot) 
    1519     
    1620  END SUBROUTINE  init_grid_param 

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

    r15 r17  
    33  USE transfert_mod 
    44  USE etat0_mod 
    5    
    6   INTEGER,PARAMETER :: euler=1, leapfrog=2, leapfrog_matsuno=3, adam_bashforth=4 
    75   
    86CONTAINS 
    … …  
    1816  USE dissip_gcm_mod 
    1917  USE ioipsl 
    20   USE caldyn_gcm_mod 
     18  USE caldyn_mod 
    2119  USE theta2theta_rhodz_mod 
    2220  USE etat0_mod 
     21  USE guided_mod 
     22  USE advect_tracer_mod 
     23   
    2324  IMPLICIT NONE 
    2425  TYPE(t_field),POINTER :: f_phis(:) 
    2526  TYPE(t_field),POINTER :: f_theta(:) 
     27  TYPE(t_field),POINTER :: f_q(:) 
    2628  TYPE(t_field),POINTER :: f_dtheta(:) 
    2729  TYPE(t_field),POINTER :: f_ps(:),f_psm1(:), f_psm2(:) 
    … …  
    3335 
    3436  REAL(rstd),POINTER :: phis(:) 
     37  REAL(rstd),POINTER :: q(:,:,:) 
    3538  REAL(rstd),POINTER :: ps(:) ,psm1(:), psm2(:) 
    3639  REAL(rstd),POINTER :: u(:,:) , um1(:,:), um2(:,:) 
    … …  
    4245  INTEGER :: ind 
    4346  INTEGER :: it,i,j,n 
    44   INTEGER :: scheme 
     47  CHARACTER(len=255) :: scheme 
    4548  INTEGER :: matsuno_period 
    4649  INTEGER :: itaumax 
    47    
     50  INTEGER :: write_period   
     51  INTEGER :: itau_out 
    4852 
    4953  dt=90. 
    … …  
    5256  itaumax=100 
    5357  CALL getin('itaumax',itaumax) 
    54    
    55   scheme=leapfrog_matsuno 
     58 
     59  write_period=0 
     60  CALL getin('write_period',write_period) 
     61  itau_out=INT(write_period/dt) 
     62   
     63  scheme='adam_bashforth' 
    5664  CALL getin('scheme',scheme) 
    5765   
    5866  matsuno_period=5 
    5967  CALL getin('matsuno_period',matsuno_period) 
    60   IF (scheme==leapfrog) matsuno_period=itaumax+1 
     68  IF (TRIM(scheme)=='leapfrog') matsuno_period=itaumax+1 
    6169 
    6270  CALL allocate_field(f_phis,field_t,type_real) 
    … …  
    7987  CALL allocate_field(f_dtheta,field_t,type_real,llm) 
    8088 
     89  CALL allocate_field(f_q,field_t,type_real,llm,nqtot) 
     90 
    8191  CALL allocate_field(f_theta_rhodz,field_t,type_real,llm) 
    8292  CALL allocate_field(f_theta_rhodzm1,field_t,type_real,llm) 
    … …  
    8898  CALL init_dissip(dt) 
    8999  CALL init_caldyn(dt) 
    90    
    91 !  CALL etat0_academic(f_ps,f_phis,f_theta_rhodz,f_u) 
    92   CALL etat0_jablonowsky06(f_ps,f_phis,f_theta_rhodz,f_u) 
    93 !  CALL test_etat0_jablonowsky06 
     100  CALL init_guided(dt) 
     101  CALL init_advect_tracer(dt) 
     102   
     103  CALL etat0(f_ps,f_phis,f_theta_rhodz,f_u, f_q) 
    94104   
    95105  DO it=0,itaumax 
    96106    PRINT *,"It No :",It 
    97107 
     108    CALL guided(it,f_ps,f_theta_rhodz,f_u,f_q) 
    98109    CALL caldyn(it,f_phis,f_ps,f_theta_rhodz,f_u, f_dps, f_dtheta_rhodz, f_du) 
    99  
    100     IF (scheme==Euler) THEN 
    101       CALL  euler_scheme 
    102     ELSE IF (scheme==leapfrog) THEN 
    103       CALL leapfrog_scheme 
    104     ELSE IF (scheme==leapfrog_matsuno) THEN 
    105       CALL  leapfrog_matsuno_scheme 
    106     ELSE IF (scheme==adam_bashforth) THEN 
    107       CALL dissip(f_u,f_du,f_ps,f_theta_rhodz,f_dtheta_rhodz) 
    108       CALL adam_bashforth_scheme 
     110    CALL advect_tracer(f_ps,f_u,f_q) 
     111     
     112    SELECT CASE (TRIM(scheme)) 
     113      CASE('euler') 
     114        CALL  euler_scheme 
     115 
     116      CASE ('leapfrog') 
     117        CALL leapfrog_scheme 
     118 
     119      CASE ('leapfrog_matsuno') 
     120        CALL  leapfrog_matsuno_scheme 
     121 
     122      CASE ('adam_bashforth') 
     123        CALL dissip(f_u,f_du,f_ps,f_theta_rhodz,f_dtheta_rhodz) 
     124        CALL adam_bashforth_scheme 
     125 
     126      CASE default 
     127        PRINT*,'Bad selector for variable scheme : <', TRIM(scheme),"> options are <euler>, <leapfrog>, <leapfrog_matsuno>, <adam_bashforth>"  
     128        STOP 
     129         
     130    END SELECT 
     131 
     132     
     133    IF ( itau_out>0 .AND. MOD(it,itau_out)==0) THEN 
     134      CALL writefield("q",f_q) 
     135      CALL writefield("ps",f_ps) 
    109136    ENDIF 
    110137 

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

    r15 r17  
    22  USE genmod 
    33  USE transfert_mod 
    4   USE etat0_mod 
     4  USE etat0_williamson_mod 
    55   
    66  INTEGER,PARAMETER :: euler=1, leapfrog=2, leapfrog_matsuno=3, adam_bashforth=4 

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

    r15 r17  
    33implicit none 
    44 
     5  TYPE ncvar 
     6    INTEGER :: size 
     7    INTEGER,POINTER :: nc_id(:) 
     8  END TYPE ncvar 
     9 
    510  integer, parameter :: MaxWriteField = 1000 
    611  integer, dimension(MaxWriteField),save :: FieldId 
    7   integer, dimension(MaxWriteField),save :: FieldVarId 
     12  TYPE(ncvar), dimension(MaxWriteField),save :: FieldVarId 
    813  integer, dimension(MaxWriteField),save :: FieldIndex 
    914  character(len=255), dimension(MaxWriteField) ::  FieldName  
    … …  
    8792      TYPE(t_domain),POINTER :: d 
    8893      INTEGER :: Index 
    89       INTEGER :: ind,i,j,k,n,ncell 
     94      INTEGER :: ind,i,j,k,n,ncell,q 
    9095      INTEGER :: iie,jje,iin,jjn 
    9196      INTEGER :: status 
    … …  
    149154            ENDDO 
    150155          ENDDO 
    151           status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index),Field_val2d,start=(/ ncell,FieldIndex(Index) /),count=(/n,1 /)) 
     156          status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val2d,start=(/ ncell,FieldIndex(Index) /),count=(/n,1 /)) 
    152157          DEALLOCATE(field_val2d) 
    153158        ELSE IF (field(ind)%ndim==3) THEN 
    … …  
    163168            ENDDO 
    164169          ENDDO 
    165            status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index),Field_val3d,start=(/ ncell,1,FieldIndex(Index) /),   & 
     170           status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val3d,start=(/ ncell,1,FieldIndex(Index) /),   & 
    166171                               count=(/n,size(field(1)%rval3d,2),1 /)) 
    167172          DEALLOCATE(field_val3d) 
    168173        ELSE IF (field(1)%ndim==4) THEN 
    169           ALLOCATE(Field_val4d(n,size(field(ind)%rval4d,2),size(field(ind)%rval4d,3))) 
    170           n=0 
    171           DO j=d%jj_begin-halo_size,d%jj_end+halo_size 
    172             DO i=d%ii_begin-halo_size,d%ii_end+halo_size 
    173               k=d%iim*(j-1)+i 
    174               IF (d%own(i,j) .OR. single) THEN 
    175                 n=n+1 
    176                 Field_val4d(n,:,:)=field(ind)%rval4d(k,:,:) 
    177               ENDIF 
    178             ENDDO 
    179           ENDDO 
    180  
    181           status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index),Field_val4d,start=(/ ncell,1,1,FieldIndex(Index) /),   & 
    182                               count=(/n,size(field(1)%rval4d,2),size(field(1)%rval4d,3),1 /)) 
    183           DEALLOCATE(field_val4d) 
     174 
     175          DO q=1,FieldVarId(index)%size 
     176           
     177            ALLOCATE(Field_val3d(n,size(field(ind)%rval4d,2))) 
     178            n=0 
     179            DO j=d%jj_begin-halo_size,d%jj_end+halo_size 
     180              DO i=d%ii_begin-halo_size,d%ii_end+halo_size 
     181                k=d%iim*(j-1)+i 
     182                IF (d%own(i,j) .OR. single) THEN 
     183                  n=n+1 
     184                  Field_val3d(n,:)=field(ind)%rval4d(k,:,q) 
     185                ENDIF 
     186              ENDDO 
     187            ENDDO 
     188            status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(q),Field_val3d,start=(/ ncell,1,FieldIndex(Index) /),   & 
     189                               count=(/n,size(field(1)%rval4d,2),1 /)) 
     190            DEALLOCATE(field_val3d) 
     191          ENDDO          
    184192        ENDIF 
    185193         
    … …  
    230238          ENDDO 
    231239 
    232           status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index),Field_val2d,start=(/ ncell,FieldIndex(Index) /),count=(/n,1 /)) 
     240          status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val2d,start=(/ ncell,FieldIndex(Index) /),count=(/n,1 /)) 
    233241          DEALLOCATE(field_val2d) 
    234242 
    … …  
    251259            ENDDO 
    252260          ENDDO 
    253            status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index),Field_val3d,start=(/ ncell,1,FieldIndex(Index) /),   & 
     261           status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val3d,start=(/ ncell,1,FieldIndex(Index) /),   & 
    254262                               count=(/n,size(field(1)%rval3d,2),1 /)) 
    255263          DEALLOCATE(field_val3d) 
    256264        ELSE IF (field(1)%ndim==4) THEN 
    257           ALLOCATE(Field_val4d(n,size(field(ind)%rval4d,2),size(field(ind)%rval4d,3))) 
    258           n=0 
    259           DO j=jj_begin+1,jj_end 
    260             DO i=ii_begin,ii_end-1 
    261               n=n+1 
    262               k=iim*(j-1)+i 
    263               Field_val4d(n,:,:)=field(ind)%rval4d(k+z_down,:,:) 
    264             ENDDO 
    265           ENDDO 
    266  
    267           DO j=jj_begin,jj_end-1 
    268             DO i=ii_begin+1,ii_end 
    269               n=n+1 
    270               k=iim*(j-1)+i 
    271               Field_val4d(n,:,:)=field(ind)%rval4d(k+z_up,:,:) 
    272             ENDDO 
    273           ENDDO 
    274  
    275           status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index),Field_val4d,start=(/ ncell,1,1,FieldIndex(Index) /),   & 
    276                               count=(/n,size(field(1)%rval4d,2),size(field(1)%rval4d,3),1 /)) 
    277           DEALLOCATE(field_val4d) 
     265 
     266          DO q=1,FieldVarId(index)%size 
     267            ALLOCATE(Field_val3d(n,size(field(ind)%rval4d,2))) 
     268            n=0 
     269            DO j=jj_begin+1,jj_end 
     270              DO i=ii_begin,ii_end-1 
     271                n=n+1 
     272                k=iim*(j-1)+i 
     273                Field_val3d(n,:)=field(ind)%rval4d(k+z_down,:,q) 
     274              ENDDO 
     275            ENDDO 
     276 
     277            DO j=jj_begin,jj_end-1 
     278              DO i=ii_begin+1,ii_end 
     279                n=n+1 
     280                k=iim*(j-1)+i 
     281                Field_val3d(n,:)=field(ind)%rval4d(k+z_up,:,q) 
     282              ENDDO 
     283            ENDDO 
     284 
     285            status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(q),Field_val3d,start=(/ ncell,1,1,FieldIndex(Index) /),   & 
     286                                count=(/n,size(field(1)%rval4d,2),1 /)) 
     287            DEALLOCATE(field_val3d) 
     288          ENDDO 
    278289        ENDIF 
    279290         
    … …  
    307318      INTEGER :: dim3id,dim4id 
    308319      INTEGER :: status 
    309       INTEGER :: ind,i,j,k,n 
     320      INTEGER :: ind,i,j,k,n,q 
    310321      INTEGER :: iie,jje,iin,jjn 
    311322      INTEGER :: ind_b,ind_e 
    … …  
    351362        status = NF90_DEF_DIM(ncid,'nvert',nvert,nvertid) 
    352363 
    353         IF (Field(ind_b)%ndim==3)  THEN 
     364        IF (Field(ind_b)%ndim==2)  THEN 
     365          FieldVarId(NbField)%size=1 
     366          ALLOCATE(FieldVarId(NbField)%nc_id(1)) 
     367        ELSE IF (Field(ind_b)%ndim==3)  THEN 
     368          FieldVarId(NbField)%size=1 
     369          ALLOCATE(FieldVarId(NbField)%nc_id(1)) 
    354370          status = NF90_DEF_DIM(ncid,'Z',size(field(ind_b)%rval3d,2),dim3id) 
    355371        ELSE IF (Field(1)%ndim==4) THEN 
     372          FieldVarId(NbField)%size=size(field(ind_b)%rval4d,3) 
     373          ALLOCATE(FieldVarId(NbField)%nc_id(FieldVarId(NbField)%size)) 
    356374          status = NF90_DEF_DIM(ncid,'Z',size(field(ind_b)%rval4d,2),dim3id) 
    357           status = NF90_DEF_DIM(ncid,'Q',size(field(ind_b)%rval4d,3),dim4id) 
     375!          status = NF90_DEF_DIM(ncid,'Q',size(field(ind_b)%rval4d,3),dim4id) 
    358376        ENDIF 
    359377       
    … …  
    372390 
    373391        IF (Field(ind_b)%ndim==2) THEN 
    374           status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),NF90_DOUBLE,(/ ncellId,timeId /),FieldVarId(NbField)) 
     392          status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),NF90_DOUBLE,(/ ncellId,timeId /),FieldVarId(NbField)%nc_id(1)) 
     393          status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon lat") 
    375394        ELSE IF (Field(ind_b)%ndim==3) THEN 
    376           status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),NF90_DOUBLE,(/ ncellId,dim3id,timeId /),FieldVarId(NbField)) 
     395          status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),NF90_DOUBLE,(/ ncellId,dim3id,timeId /),FieldVarId(NbField)%nc_id(1)) 
     396          status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon lat") 
    377397        ELSE IF (Field(ind_b)%ndim==4) THEN 
    378           status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),NF90_DOUBLE,(/ ncellId,dim3id,dim4id,timeId /),FieldVarId(NbField)) 
     398          DO i=1,FieldVarId(NbField)%size 
     399            status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name))//int2str(i),NF90_DOUBLE,(/ ncellId,dim3id,timeId /),FieldVarId(NbField)%nc_id(i)) 
     400            status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(i),"coordinates","lon lat") 
     401          ENDDO         
    379402        ENDIF  
    380403   
    381         status = NF90_PUT_ATT(ncid,FieldVarId(NbField),"coordinates","lon lat") 
    382          
     404       
    383405        status = NF90_ENDDEF(ncid)       
    384406 
    … …  
    445467        status = NF90_DEF_DIM(ncid,'nvert',nvert,nvertid) 
    446468 
    447         IF (Field(ind_b)%ndim==3)  THEN 
     469        IF (Field(ind_b)%ndim==2)  THEN 
     470          FieldVarId(NbField)%size=1 
     471          ALLOCATE(FieldVarId(NbField)%nc_id(1)) 
     472        ELSE IF (Field(ind_b)%ndim==3)  THEN 
     473          FieldVarId(NbField)%size=1 
     474          ALLOCATE(FieldVarId(NbField)%nc_id(1)) 
    448475          status = NF90_DEF_DIM(ncid,'Z',size(field(ind_b)%rval3d,2),dim3id) 
    449476        ELSE IF (Field(1)%ndim==4) THEN 
     477          FieldVarId(NbField)%size=size(field(ind_b)%rval4d,3) 
     478          ALLOCATE(FieldVarId(NbField)%nc_id(FieldVarId(NbField)%size)) 
    450479          status = NF90_DEF_DIM(ncid,'Z',size(field(ind_b)%rval4d,2),dim3id) 
    451           status = NF90_DEF_DIM(ncid,'Q',size(field(ind_b)%rval4d,3),dim4id) 
    452480        ENDIF 
     481 
     482 
    453483       
    454484        status = NF90_DEF_DIM(ncid,'time_counter',NF90_UNLIMITED,timeId) 
    … …  
    465495        status = NF90_DEF_VAR(ncid,'bounds_lat',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_latId) 
    466496 
     497 
    467498        IF (Field(ind_b)%ndim==2) THEN 
    468           status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),NF90_DOUBLE,(/ ncellId,timeId /),FieldVarId(NbField)) 
     499          status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),NF90_DOUBLE,(/ ncellId,timeId /),FieldVarId(NbField)%nc_id(1)) 
     500          status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon lat") 
    469501        ELSE IF (Field(ind_b)%ndim==3) THEN 
    470           status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),NF90_DOUBLE,(/ ncellId,dim3id,timeId /),FieldVarId(NbField)) 
     502          status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),NF90_DOUBLE,(/ ncellId,dim3id,timeId /),FieldVarId(NbField)%nc_id(1)) 
     503          status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon lat") 
    471504        ELSE IF (Field(ind_b)%ndim==4) THEN 
    472           status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),NF90_DOUBLE,(/ ncellId,dim3id,dim4id,timeId /),FieldVarId(NbField)) 
     505          DO q=1,FieldVarId(NbField)%size 
     506            status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)//int2str(q)),NF90_DOUBLE,(/ ncellId,dim3id,timeId /),FieldVarId(NbField)%nc_id(q)) 
     507            status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(q),"coordinates","lon lat") 
     508          ENDDO         
    473509        ENDIF  
    474    
    475         status = NF90_PUT_ATT(ncid,FieldVarId(NbField),"coordinates","lon lat") 
    476510         
    477511        status = NF90_ENDDEF(ncid)       
Note: See TracChangeset for help on using the changeset viewer.