Index: LMDZ5/branches/testing/libf/bibio/wxios.F90 =================================================================== --- LMDZ5/branches/testing/libf/bibio/wxios.F90 (revision 1864) +++ LMDZ5/branches/testing/libf/bibio/wxios.F90 (revision 1864) @@ -0,0 +1,523 @@ +! $Id: wxios.F90 $ +#ifdef CPP_XIOS +MODULE wxios + USE xios + USE iaxis + USE iaxis_attr + USE icontext_attr + USE idate + USE idomain_attr + USE ifield_attr + USE ifile_attr + USE ixml_tree + + !Variables disponibles pendant toute l'execution du programme: + + INTEGER, SAVE :: g_comm + CHARACTER(len=100), SAVE :: g_ctx_name + TYPE(xios_context), SAVE :: g_ctx +!$OMP THREADPRIVATE(g_comm,g_cts_name,g_ctx) + LOGICAL, SAVE :: g_flag_xml = .FALSE. + CHARACTER(len=100) :: g_field_name = "nofield" +!$OMP THREADPRIVATE(g_flag_xml,g_field_name) + + + CONTAINS + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! str + i => str_i !!!!!!!!!!!!!!!!!!!! + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + SUBROUTINE concat(str, str2, str_str2) + CHARACTER(len=*), INTENT(IN) :: str, str2 + CHARACTER(len=20), INTENT(OUT) :: str_str2 + + + str_str2 = TRIM(ADJUSTL(str//"_"//TRIM(ADJUSTL(str2)))) + !IF (prt_level >= 10) WRITE(lunout,*) "Xios. ",str,"+",str2,"=",str_str2 + END SUBROUTINE concat + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! 36day => 36d etc !!!!!!!!!!!!!!!!!!!! + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + SUBROUTINE reformadate(odate, ndate) + CHARACTER(len=*), INTENT(IN) :: odate + CHARACTER(len=100), INTENT(OUT) :: ndate + + INTEGER :: i = 0 + !!!!!!!!!!!!!!!!!! + ! Pour XIOS: + ! year : y + ! month : mo + ! day : d + ! hour : h + ! minute : mi + ! second : s + !!!!!!!!!!!!!!!!!! + + i = INDEX(odate, "day") + IF (i > 0) THEN + ndate = odate(1:i-1)//"d" + END IF + + i = INDEX(odate, "hr") + IF (i > 0) THEN + ndate = odate(1:i-1)//"h" + END IF + + i = INDEX(odate, "mth") + IF (i > 0) THEN + ndate = odate(1:i-1)//"mo" + END IF + + !IF (prt_level >= 10) WRITE(lunout,*) "Xios. ", odate, " => ", ndate + END SUBROUTINE reformadate + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! ave(X) => average etc !!!!!!!!!!!!!!! + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + CHARACTER(len=7) FUNCTION reformaop(op) + CHARACTER(len=*), INTENT(IN) :: op + + INTEGER :: i = 0 + reformaop = "average" + + IF (op.EQ."inst(X)") THEN + reformaop = "instant" + END IF + + IF (op.EQ."once") THEN + reformaop = "once" + END IF + + IF (op.EQ."t_max(X)") THEN + reformaop = "maximum" + END IF + + IF (op.EQ."t_min(X)") THEN + reformaop = "minimum" + END IF + + !IF (prt_level >= 10) WRITE(lunout,*) "Xios. ", op, " => ", reformaop + END FUNCTION reformaop + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! Routine d'initialisation !!!!!!!!!!!!! + ! A lancer juste après mpi_init !!!!!!!!!!!!! + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + SUBROUTINE wxios_init(xios_ctx_name, locom, outcom) + IMPLICIT NONE + INCLUDE 'iniprint.h' + + CHARACTER(len=*), INTENT(IN) :: xios_ctx_name + INTEGER, INTENT(IN), OPTIONAL :: locom + INTEGER, INTENT(OUT), OPTIONAL :: outcom + + + TYPE(xios_context) :: xios_ctx + INTEGER :: xios_comm + + IF (prt_level >= 10) WRITE(lunout,*) "Xios. Initialization" + + + + IF (PRESENT(locom)) THEN + IF (prt_level >= 10) WRITE(lunout,*) "Xios. ctx=",xios_ctx_name,"local_comm=",locom,",return_comm=",xios_comm + CALL xios_initialize(xios_ctx_name, local_comm = locom, return_comm = xios_comm ) + ELSE + IF (prt_level >= 10) WRITE(lunout,*) "Xios. ctx=",xios_ctx_name,"return_comm=",outcom + CALL xios_initialize(xios_ctx_name, return_comm = xios_comm ) + END IF + + IF (PRESENT(outcom)) THEN + outcom = xios_comm + END IF + + !Enregistrement des variables globales: + g_comm = xios_comm + g_ctx_name = xios_ctx_name + + CALL wxios_context_init() + + END SUBROUTINE wxios_init + + SUBROUTINE wxios_context_init() + IMPLICIT NONE + INCLUDE 'iniprint.h' + + TYPE(xios_context) :: xios_ctx + + !Initialisation du contexte: + CALL xios_context_initialize(g_ctx_name, g_comm) + CALL xios_get_handle(g_ctx_name, xios_ctx) !Récupération + CALL xios_set_current_context(xios_ctx) !Activation + g_ctx = xios_ctx + + IF (prt_level >= 10) WRITE(lunout,*) "Xios. Current context is ", g_ctx_name + + !Une première analyse des héritages: + CALL xios_solve_inheritance() + END SUBROUTINE wxios_context_init + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! Routine de paramétrisation !!!!!!!!!!!!!!!!!! + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + SUBROUTINE wxios_set_cal(pasdetemps, calendrier, annee, mois, jour, heure) + IMPLICIT NONE + INCLUDE 'iniprint.h' + + !Paramètres: + CHARACTER(len=*), INTENT(IN) :: calendrier + INTEGER, INTENT(IN) :: annee, mois, jour + REAL, INTENT(IN) :: pasdetemps, heure + + !Variables: + CHARACTER(len=80) :: abort_message + CHARACTER(len=19) :: date + INTEGER :: njour = 1 + + !Variables pour xios: + TYPE(xios_time) :: mdtime + !REAL(kind = 8) :: year=0, month=0, day=0, hour=0, minute=0, second=0 + + mdtime = xios_time(0, 0, 0, 0, 0, pasdetemps) + + !Réglage du calendrier: + SELECT CASE (calendrier) + CASE('earth_360d') + CALL xios_set_context_attr_hdl(g_ctx, calendar_type= "D360") + IF (prt_level >= 10) WRITE(lunout,*) 'Xios. Calendrier terrestre a 360 jours/an' + CASE('earth_365d') + CALL xios_set_context_attr_hdl(g_ctx, calendar_type= "NoLeap") + IF (prt_level >= 10) WRITE(lunout,*) 'Xios. Calendrier terrestre a 365 jours/an' + CASE('earth_366d') + CALL xios_set_context_attr_hdl(g_ctx, calendar_type= "Gregorian") + IF (prt_level >= 10) WRITE(lunout,*) 'Xios. Calendrier gregorien' + CASE DEFAULT + abort_message = 'Xios. Mauvais choix de calendrier' + CALL abort_gcm('Gcm:Xios',abort_message,1) + END SELECT + + !Formatage de la date de départ: + WRITE(date, "(i4.4,'-',i2.2,'-',i2.2,' 00:00:00')") annee, mois, jour + + IF (prt_level >= 10) WRITE(lunout,*) "Xios. Initial time: ", date + + CALL xios_set_context_attr_hdl(g_ctx, start_date= date) + + !Et enfin,le pas de temps: + CALL xios_set_timestep(mdtime) + IF (prt_level >= 10) WRITE(lunout,*) "Xios. ts=",mdtime + END SUBROUTINE wxios_set_cal + + SUBROUTINE wxios_set_timestep(ts) + REAL, INTENT(IN) :: ts + TYPE(xios_time) :: mdtime + + mdtime = xios_time(0, 0, 0, 0, 0, ts) + + CALL xios_set_timestep(mdtime) + END SUBROUTINE wxios_set_timestep + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! Pour initialiser un domaine !!!!!!!!!!!!!!!!!!!! + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + SUBROUTINE wxios_domain_param(dom_id, is_sequential, ni, nj, ni_glo, nj_glo, & + ibegin, iend, ii_begin, ii_end, jbegin, jend, & + data_ni, data_ibegin, data_iend, & + io_lat, io_lon) + + IMPLICIT NONE + INCLUDE 'iniprint.h' + + CHARACTER (len=*), INTENT(IN) :: dom_id + LOGICAL, INTENT(IN) :: is_sequential + INTEGER, INTENT(IN) :: ni, nj, ni_glo, nj_glo, ibegin, iend, ii_begin, ii_end, jbegin, jend + INTEGER, INTENT(IN) :: data_ni, data_ibegin, data_iend + REAL, DIMENSION(:) :: io_lat, io_lon + + + TYPE(xios_domain) :: dom + LOGICAL :: boool + + !Masque pour les problèmes de recouvrement MPI: + LOGICAL :: mask(ni,nj) + + !On récupère le handle: + CALL xios_get_domain_handle(dom_id, dom) + + IF (prt_level >= 10) WRITE(lunout,*) "Xios. ni:",ni," ni_glo:", ni_glo, " nj:", nj, " nj_glo:", nj_glo + IF (prt_level >= 10) WRITE(lunout,*) "Xios. ibegin:",ibegin," iend:", iend, " jbegin:", jbegin, " jend:", jend + IF (prt_level >= 10) WRITE(lunout,*) "Xios. Size lon:", SIZE(io_lon(ibegin:iend)), " lat:", SIZE(io_lat(jbegin:jend)) + + !On parametrise le domaine: + CALL xios_set_domain_attr_hdl(dom, ni_glo=ni_glo, ibegin=ibegin, ni=ni) + CALL xios_set_domain_attr_hdl(dom, nj_glo=nj_glo, jbegin=jbegin, nj=nj, data_dim=2) + CALL xios_set_domain_attr_hdl(dom, lonvalue=io_lon(ibegin:iend), latvalue=io_lat(jbegin:jend)) + + IF (.NOT.is_sequential) THEN + mask(:,:)=.TRUE. + mask(1:ii_begin-1,1) = .FALSE. + mask(ii_end+1:ni,nj) = .FALSE. + IF (prt_level >= 10) WRITE(lunout,*) "Xios. mask" + !CALL xios_set_domain_attr_hdl(dom, mask=mask) + END IF + + CALL xios_is_defined_domain_attr_hdl(dom,ni_glo=boool) + !Vérification: + IF (xios_is_valid_domain(dom_id)) THEN + IF (prt_level >= 10) WRITE(lunout,*) "Xios. Domain initialized: ", dom_id, boool + ELSE + IF (prt_level >= 10) WRITE(lunout,*) "Xios. Invalid domain: ", dom_id + END IF + END SUBROUTINE wxios_domain_param + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! Pour déclarer un axe vertical !!!!!!!!!!!!!!! + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + SUBROUTINE wxios_add_vaxis(axisgroup_id, axis_file, axis_size, axis_value) + IMPLICIT NONE + INCLUDE 'iniprint.h' + + CHARACTER (len=*), INTENT(IN) :: axisgroup_id, axis_file + INTEGER, INTENT(IN) :: axis_size + REAL, DIMENSION(axis_size), INTENT(IN) :: axis_value + + TYPE(xios_axisgroup) :: axgroup + TYPE(xios_axis) :: ax + CHARACTER(len=20) :: axis_id + + + !Préparation du nom de l'axe: + CALL concat(axisgroup_id, axis_file, axis_id) + + !On récupère le groupe d'axes qui va bien: + CALL xios_get_axisgroup_handle(axisgroup_id, axgroup) + + !On ajoute l'axe correspondant à ce fichier: + CALL xios_add_axis(axgroup, ax, TRIM(ADJUSTL(axis_id))) + + !Et on le parametrise: + CALL xios_set_axis_attr_hdl(ax, size=axis_size, value=axis_value) + + !Vérification: + IF (xios_is_valid_axis(TRIM(ADJUSTL(axis_id)))) THEN + IF (prt_level >= 10) WRITE(lunout,*) "Xios. Axis created: ", TRIM(ADJUSTL(axis_id)) + ELSE + WRITE(*,*) "Xios. Invalid axis: ", TRIM(ADJUSTL(axis_id)) + END IF + + END SUBROUTINE wxios_add_vaxis + + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! Pour déclarer un fichier !!!!!!!!!!!!!!!!!!! + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + SUBROUTINE wxios_add_file(fname, ffreq, flvl) + IMPLICIT NONE + INCLUDE 'iniprint.h' + + CHARACTER(len=*), INTENT(IN) :: fname + CHARACTER(len=*), INTENT(IN) :: ffreq + INTEGER, INTENT(IN) :: flvl + + TYPE(xios_file) :: x_file + TYPE(xios_filegroup) :: x_fg + CHARACTER(len=100) :: nffreq + + !On regarde si le fichier n'est pas défini par XML: + IF (.NOT.xios_is_valid_file(fname)) THEN + !On créé le noeud: + CALL xios_get_filegroup_handle("defile", x_fg) + CALL xios_add_file(x_fg, x_file, fname) + + !On reformate la fréquence: + CALL reformadate(ffreq, nffreq) + + !On configure: + CALL xios_set_file_attr_hdl(x_file, name="X"//fname,& + output_freq=TRIM(ADJUSTL(nffreq)), output_level=flvl, enabled=.TRUE.) + + IF (xios_is_valid_file("X"//fname)) THEN + IF (prt_level >= 10) WRITE(lunout,*) "Xios. New file: ", "X"//fname + IF (prt_level >= 10) WRITE(lunout,*) "Xios. output_freq=",TRIM(ADJUSTL(nffreq)),"; output_lvl=",flvl + ELSE + WRITE(*,*) "Xios. Error, invalid file: ", "X"//fname + WRITE(*,*) "Xios. output_freq=",TRIM(ADJUSTL(nffreq)),"; output_lvl=",flvl + END IF + ELSE + IF (prt_level >= 10) WRITE(lunout,*) "Xios. Fichier ", fname, " défini par XML." + CALL xios_set_file_attr(fname, enabled=.TRUE.) + END IF + END SUBROUTINE wxios_add_file + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! Pour créer un champ !!!!!!!!!!!!!!!!!!!! + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + SUBROUTINE wxios_add_field(fieldname, fieldgroup, fieldlongname, fieldunit) + USE netcdf + + IMPLICIT NONE + INCLUDE 'iniprint.h' + + CHARACTER(len=*), INTENT(IN) :: fieldname + TYPE(xios_fieldgroup), INTENT(IN) :: fieldgroup + CHARACTER(len=*), INTENT(IN) :: fieldlongname + CHARACTER(len=*), INTENT(IN) :: fieldunit + + TYPE(xios_field) :: field + CHARACTER(len=10) :: newunit + REAL(KIND=8) :: def + + !La valeur par défaut des champs non définis: + def = nf90_fill_real + + IF (fieldunit .EQ. " ") THEN + newunit = "-" + ELSE + newunit = fieldunit + ENDIF + + !On ajoute le champ: + CALL xios_add_field(fieldgroup, field, fieldname) + !IF (prt_level >= 10) WRITE(lunout,*) "Xios. ",fieldname,fieldgroup, fieldlongname, fieldunit + + !On rentre ses paramètres: + CALL xios_set_field_attr_hdl(field, standard_name=fieldlongname, unit=newunit, default_value=def) + IF (prt_level >= 10) WRITE(lunout,*) "Xios. Champ ", fieldname, "cree:" + IF (prt_level >= 10) WRITE(lunout,*) "Xios. long_name=",fieldlongname,"; unit=",newunit,"; default_value=",nf90_fill_real + + END SUBROUTINE wxios_add_field + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! Pour déclarer un champ !!!!!!!!!!!!!!!!! + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + SUBROUTINE wxios_add_field_to_file(fieldname, fdim, fid, fname, fieldlongname, fieldunit, field_level, op) + IMPLICIT NONE + INCLUDE 'iniprint.h' + + CHARACTER(len=*), INTENT(IN) :: fieldname + INTEGER, INTENT(IN) :: fdim, fid + CHARACTER(len=*), INTENT(IN) :: fname + CHARACTER(len=*), INTENT(IN) :: fieldlongname + CHARACTER(len=*), INTENT(IN) :: fieldunit + INTEGER, INTENT(IN) :: field_level + CHARACTER(len=*), INTENT(IN) :: op + + CHARACTER(len=20) :: axis_id + CHARACTER(len=100) :: operation + TYPE(xios_file) :: f + TYPE(xios_field) :: field + TYPE(xios_fieldgroup) :: fieldgroup + LOGICAL :: bool=.FALSE. + INTEGER :: lvl =0 + + + !Préparation du nom de l'axe: + CALL concat("presnivs", fname, axis_id) + + !on prépare le nom de l'opération: + operation = reformaop(op) + + + + !On selectionne le bon groupe de champs: + IF (fdim.EQ.2) THEN + CALL xios_get_fieldgroup_handle("fields_2D", fieldgroup) + ELSE + CALL xios_get_fieldgroup_handle("fields_3D", fieldgroup) + ENDIF + + !On regarde si le champ à déjà été créé ou non: + IF (xios_is_valid_field(fieldname) .AND. .NOT. g_field_name == fieldname) THEN + !Si ce champ existe via XML (ie, dès le premier passage, ie g_field_name != fieldname) alors rien d'autre à faire + IF (prt_level >= 10) WRITE(lunout,*) "Xios. Champ ", fieldname, "existe via XML" + g_flag_xml = .TRUE. + g_field_name = fieldname + + ELSE IF (.NOT. g_field_name == fieldname) THEN + !Si premier pssage et champ indéfini, alors on le créé + + IF (prt_level >= 10) WRITE(lunout,*) "Xios. Champ ", fieldname, "nexiste pas" + + !On le créé: + CALL wxios_add_field(fieldname, fieldgroup, fieldlongname, fieldunit) + IF (xios_is_valid_field(fieldname)) THEN + IF (prt_level >= 10) WRITE(lunout,*) "Xios. Champ ", fieldname, "cree" + ENDIF + + g_flag_xml = .FALSE. + g_field_name = fieldname + + END IF + + IF (.NOT. g_flag_xml) THEN + !Champ existe déjà, mais pas XML, alors on l'ajoute + !On ajoute le champ: + CALL xios_get_file_handle(fname, f) + CALL xios_add_fieldtofile(f, field) + + + !L'operation, sa frequence: + CALL xios_set_field_attr_hdl(field, field_ref=fieldname, operation=TRIM(ADJUSTL(operation)), freq_op="1ts", prec=4) + + + !On rentre ses paramètres: + CALL xios_set_field_attr_hdl(field, level=field_level, enabled=.TRUE.) + + IF (fdim.EQ.2) THEN + !Si c'est un champ 2D: + IF (prt_level >= 10) WRITE(lunout,*) "Xios. Champ 2D ", fieldname, " de ", "X"//fname ," configure:" + IF (prt_level >= 10) WRITE(lunout,*) "Xios. op=", TRIM(ADJUSTL(operation)) + IF (prt_level >= 10) WRITE(lunout,*) "Xios. freq_op=1ts","; lvl=",field_level + ELSE + !Si 3D : + !On ajoute l'axe vertical qui va bien: + CALL xios_set_field_attr_hdl(field, axis_ref=TRIM(ADJUSTL(axis_id))) + + IF (prt_level >= 10) WRITE(lunout,*) "Xios. Champ 3D ", fieldname, " de ", "X"//fname, "configure:" + IF (prt_level >= 10) WRITE(lunout,*) "Xios. freq_op=1ts","; lvl=",field_level + IF (prt_level >= 10) WRITE(lunout,*) "Xios. axe=",TRIM(ADJUSTL(axis_id)) + END IF + + ELSE + !Sinon on se contente de l'activer: + CALL xios_set_field_attr(fieldname, enabled=.TRUE.) + ENDIF + + END SUBROUTINE wxios_add_field_to_file + + SUBROUTINE wxios_update_calendar(ito) + INTEGER, INTENT(IN) :: ito + CALL xios_update_calendar(ito) + END SUBROUTINE wxios_update_calendar + + SUBROUTINE wxios_write_2D(fieldname, fdata) + CHARACTER(len=*), INTENT(IN) :: fieldname + REAL, DIMENSION(:,:), INTENT(IN) :: fdata + + CALL xios_send_field(fieldname, fdata) + END SUBROUTINE wxios_write_2D + + SUBROUTINE wxios_write_3D(fieldname, fdata) + CHARACTER(len=*), INTENT(IN) :: fieldname + REAL, DIMENSION(:,:,:), INTENT(IN) :: fdata + + CALL xios_send_field(fieldname, fdata) + END SUBROUTINE wxios_write_3D + + SUBROUTINE wxios_closedef() + CALL xios_close_context_definition() + CALL xios_update_calendar(0) + END SUBROUTINE wxios_closedef + + SUBROUTINE wxios_close() + CALL xios_context_finalize() + CALL xios_finalize() + END SUBROUTINE wxios_close +END MODULE wxios +#endif