Index: LMDZ5/trunk/libf/dyn3dmem/calfis_loc.F =================================================================== --- LMDZ5/trunk/libf/dyn3dmem/calfis_loc.F (revision 1656) +++ LMDZ5/trunk/libf/dyn3dmem/calfis_loc.F (revision 1657) @@ -1,4 +1,4 @@ ! -! $Id: calfis_p.F 1299 2010-01-20 14:27:21Z fairhead $ +! $Id$ ! C @@ -108,4 +108,5 @@ #include "comvert.h" #include "comgeom2.h" +#include "iniprint.h" #ifdef CPP_MPI include 'mpif.h' @@ -180,8 +181,31 @@ REAL,SAVE,ALLOCATABLE :: flxwfi_omp(:,:) ! Flux de masse verticale sur la grille physiq +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! Introduction du splitting (FH) +! Question pour Yann : +! J'ai été surpris au début que les tableaux zufi_omp, zdufi_omp n'co soitent +! en SAVE. Je crois comprendre que c'est parce que tu voulais qu'il +! soit allocatable (plutot par exemple que de passer une dimension +! dépendant du process en argument des routines) et que, du coup, +! le SAVE évite d'avoir à refaire l'allocation à chaque appel. +! Tu confirmes ? +! J'ai suivi le même principe pour les zdufic_omp +! Mais c'est surement bien que tu controles. +! + + REAL,ALLOCATABLE,SAVE :: zdufic_omp(:,:) + REAL,ALLOCATABLE,SAVE :: zdvfic_omp(:,:) + REAL,ALLOCATABLE,SAVE :: zdtfic_omp(:,:) + REAL,ALLOCATABLE,SAVE :: zdqfic_omp(:,:,:) + REAL jH_cur_split,zdt_split + LOGICAL debut_split,lafin_split + INTEGER isplit +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + c$OMP THREADPRIVATE(zplev_omp,zplay_omp,zphi_omp,zphis_omp, c$OMP+ presnivs_omp,zufi_omp,zvfi_omp,ztfi_omp, c$OMP+ zqfi_omp,zdufi_omp,zdvfi_omp, -c$OMP+ zdtfi_omp,zdqfi_omp,zdpsrf_omp,flxwfi_omp) +c$OMP+ zdtfi_omp,zdqfi_omp,zdpsrf_omp,flxwfi_omp, +c$OMP+ zdufic_omp,zdvfic_omp,zdtfic_omp,zdqfic_omp) LOGICAL,SAVE :: first_omp=.true. @@ -235,9 +259,10 @@ debut = .TRUE. IF (ngridmx.NE.2+(jjm-1)*iim) THEN - PRINT*,'STOP dans calfis' - PRINT*,'La dimension ngridmx doit etre egale a 2 + (jjm-1)*iim' - PRINT*,' ngridmx jjm iim ' - PRINT*,ngridmx,jjm,iim - STOP + write(lunout,*) 'STOP dans calfis' + write(lunout,*) + & 'La dimension ngridmx doit etre egale a 2 + (jjm-1)*iim' + write(lunout,*) ' ngridmx jjm iim ' + write(lunout,*) ngridmx,jjm,iim + STOP ENDIF c$OMP MASTER @@ -533,4 +558,8 @@ allocate(zdtfi_omp(klon,llm)) allocate(zdqfi_omp(klon,llm,nqtot)) + allocate(zdufic_omp(klon,llm)) + allocate(zdvfic_omp(klon,llm)) + allocate(zdtfic_omp(klon,llm)) + allocate(zdqfic_omp(klon,llm,nqtot)) allocate(zdpsrf_omp(klon)) allocate(flxwfi_omp(klon,llm)) @@ -635,11 +664,29 @@ if (planet_type=="earth") then #ifdef CPP_EARTH + + +!$OMP MASTER + write(lunout,*) 'PHYSIQUE AVEC NSPLIT_PHYS=',nsplit_phys +!$OMP END MASTER + zdt_split=dtphys/nsplit_phys + zdufic_omp(:,:)=0. + zdvfic_omp(:,:)=0. + zdtfic_omp(:,:)=0. + zdqfic_omp(:,:,:)=0. + + do isplit=1,nsplit_phys + + jH_cur_split=jH_cur+(isplit-1) * dtvr / (daysec *nsplit_phys) + debut_split=debut.and.isplit==1 + lafin_split=lafin.and.isplit==nsplit_phys + + CALL physiq (klon, . llm, - . debut, - . lafin, + . debut_split, + . lafin_split, . jD_cur, - . jH_cur, - . dtphys, + . jH_cur_split, + . zdt_split, . zplev_omp, . zplay_omp, @@ -663,4 +710,22 @@ . pducov, . PVteta) + + zufi_omp(:,:)=zufi_omp(:,:)+zdufi_omp(:,:)*zdt_split + zvfi_omp(:,:)=zvfi_omp(:,:)+zdvfi_omp(:,:)*zdt_split + ztfi_omp(:,:)=ztfi_omp(:,:)+zdtfi_omp(:,:)*zdt_split + zqfi_omp(:,:,:)=zqfi_omp(:,:,:)+zdqfi_omp(:,:,:)*zdt_split + + zdufic_omp(:,:)=zdufic_omp(:,:)+zdufi_omp(:,:) + zdvfic_omp(:,:)=zdvfic_omp(:,:)+zdvfi_omp(:,:) + zdtfic_omp(:,:)=zdtfic_omp(:,:)+zdtfi_omp(:,:) + zdqfic_omp(:,:,:)=zdqfic_omp(:,:,:)+zdqfi_omp(:,:,:) + + enddo + + zdufi_omp(:,:)=zdufic_omp(:,:)/nsplit_phys + zdvfi_omp(:,:)=zdvfic_omp(:,:)/nsplit_phys + zdtfi_omp(:,:)=zdtfic_omp(:,:)/nsplit_phys + zdqfi_omp(:,:,:)=zdqfic_omp(:,:,:)/nsplit_phys + #endif endif !of if (planet_type=="earth") @@ -1116,4 +1181,6 @@ #else write(*,*) "calfis_p: for now can only work with parallel physics" + write(lunout,*) + & "calfis_p: for now can only work with parallel physics" stop #endif Index: LMDZ5/trunk/libf/dyn3dmem/conf_gcm.F =================================================================== --- LMDZ5/trunk/libf/dyn3dmem/conf_gcm.F (revision 1656) +++ LMDZ5/trunk/libf/dyn3dmem/conf_gcm.F (revision 1657) @@ -1,4 +1,4 @@ ! -! $Id: conf_gcm.F 1299 2010-01-20 14:27:21Z fairhead $ +! $Id: conf_gcm.F 1357 2010-04-14 14:03:19Z emillour $ ! c @@ -174,4 +174,12 @@ CALL getin('day_step',day_step) +!Config Key = nsplit_phys +!Config Desc = nombre d'iteration de la physique +!Config Def = 240 +!Config Help = nombre d'itration de la physique +! + nsplit_phys = 1 + CALL getin('nsplit_phys',nsplit_phys) + !Config Key = iperiod !Config Desc = periode pour le pas Matsuno @@ -594,4 +602,19 @@ CALL getin('ok_dynzon',ok_dynzon) +!Config Key = ok_dyn_ins +!Config Desc = sorties instantanees dans la dynamique +!Config Def = n +!Config Help = +!Config + ok_dyn_ins = .FALSE. + CALL getin('ok_dyn_ins',ok_dyn_ins) + +!Config Key = ok_dyn_ave +!Config Desc = sorties moyennes dans la dynamique +!Config Def = n +!Config Help = +!Config + ok_dyn_ave = .FALSE. + CALL getin('ok_dyn_ave',ok_dyn_ave) write(lunout,*)' #########################################' @@ -604,4 +627,5 @@ write(lunout,*)' day_step = ', day_step write(lunout,*)' iperiod = ', iperiod + write(lunout,*)' nsplit_phys = ', nsplit_phys write(lunout,*)' iconser = ', iconser write(lunout,*)' iecri = ', iecri @@ -633,4 +657,6 @@ write(lunout,*)' config_inca = ', config_inca write(lunout,*)' ok_dynzon = ', ok_dynzon + write(lunout,*)' ok_dyn_ins = ', ok_dyn_ins + write(lunout,*)' ok_dyn_ave = ', ok_dyn_ave RETURN @@ -765,4 +791,20 @@ ok_dynzon = .FALSE. CALL getin('ok_dynzon',ok_dynzon) + +!Config Key = ok_dyn_ins +!Config Desc = sorties instantanees dans la dynamique +!Config Def = n +!Config Help = +!Config + ok_dyn_ins = .FALSE. + CALL getin('ok_dyn_ins',ok_dyn_ins) + +!Config Key = ok_dyn_ave +!Config Desc = sorties moyennes dans la dynamique +!Config Def = n +!Config Help = +!Config + ok_dyn_ave = .FALSE. + CALL getin('ok_dyn_ave',ok_dyn_ave) !Config Key = use_filtre_fft @@ -859,4 +901,6 @@ write(lunout,*)' config_inca = ', config_inca write(lunout,*)' ok_dynzon = ', ok_dynzon + write(lunout,*)' ok_dyn_ins = ', ok_dyn_ins + write(lunout,*)' ok_dyn_ave = ', ok_dyn_ave write(lunout,*)' use_filtre_fft = ', use_filtre_fft write(lunout,*)' use_mpi_alloc = ', use_mpi_alloc Index: LMDZ5/trunk/libf/dyn3dmem/control_mod.F90 =================================================================== --- LMDZ5/trunk/libf/dyn3dmem/control_mod.F90 (revision 1656) +++ LMDZ5/trunk/libf/dyn3dmem/control_mod.F90 (revision 1657) @@ -11,10 +11,17 @@ REAL :: periodav - INTEGER :: nday,day_step,iperiod,iapp_tracvl + INTEGER :: nday,day_step,iperiod,iapp_tracvl,nsplit_phys INTEGER :: iconser,iecri,idissip,iphysiq,iecrimoy INTEGER :: dayref,anneeref, raz_date, ip_ebil_dyn - LOGICAL :: offline, output_grads_dyn + LOGICAL :: offline CHARACTER (len=4) :: config_inca - CHARACTER (len=10) :: planet_type + CHARACTER (len=10) :: planet_type ! planet type ('earth','mars',...) + LOGICAL output_grads_dyn ! output dynamics diagnostics in + ! binary grads file 'dyn.dat' (y/n) + LOGICAL ok_dynzon ! output zonal transports in dynzon.nc file + LOGICAL ok_dyn_ins ! output instantaneous values of fields + ! in the dynamics in NetCDF files dyn_hist*nc + LOGICAL ok_dyn_ave ! output averaged values of fields in the dynamics + ! in NetCDF files dyn_hist*ave.nc LOGICAL ok_dynzon ! output zonal transports in dynzon.nc file Index: LMDZ5/trunk/libf/dyn3dmem/dynetat0.F =================================================================== --- LMDZ5/trunk/libf/dyn3dmem/dynetat0.F (revision 1656) +++ LMDZ5/trunk/libf/dyn3dmem/dynetat0.F (revision 1657) @@ -1,7 +1,8 @@ ! -! $Header$ +! $Id $ ! SUBROUTINE dynetat0(fichnom,vcov,ucov, . teta,q,masse,ps,phis,time) + USE infotrac IMPLICIT NONE @@ -33,4 +34,5 @@ #include "serre.h" #include "logic.h" +#include "iniprint.h" c Arguments: @@ -52,10 +54,11 @@ c----------------------------------------------------------------------- + c Ouverture NetCDF du fichier etat initial ierr = NF_OPEN (fichnom, NF_NOWRITE,nid) IF (ierr.NE.NF_NOERR) THEN - write(6,*)' Pb d''ouverture du fichier start.nc' - write(6,*)' ierr = ', ierr + write(lunout,*)'dynetat0: Pb d''ouverture du fichier start.nc' + write(lunout,*)' ierr = ', ierr CALL ABORT ENDIF @@ -64,5 +67,5 @@ ierr = NF_INQ_VARID (nid, "controle", nvarid) IF (ierr .NE. NF_NOERR) THEN - PRINT*, "dynetat0: Le champ est absent" + write(lunout,*)"dynetat0: Le champ est absent" CALL abort ENDIF @@ -73,5 +76,5 @@ #endif IF (ierr .NE. NF_NOERR) THEN - PRINT*, "dynetat0: Lecture echoue pour " + write(lunout,*)"dynetat0: Lecture echoue pour " CALL abort ENDIF @@ -119,5 +122,6 @@ c c - PRINT*,'rad,omeg,g,cpp,kappa',rad,omeg,g,cpp,kappa + write(lunout,*)'dynetat0: rad,omeg,g,cpp,kappa', + & rad,omeg,g,cpp,kappa IF( im.ne.iim ) THEN @@ -134,5 +138,5 @@ ierr = NF_INQ_VARID (nid, "rlonu", nvarid) IF (ierr .NE. NF_NOERR) THEN - PRINT*, "dynetat0: Le champ est absent" + write(lunout,*)"dynetat0: Le champ est absent" CALL abort ENDIF @@ -143,5 +147,5 @@ #endif IF (ierr .NE. NF_NOERR) THEN - PRINT*, "dynetat0: Lecture echouee pour " + write(lunout,*)"dynetat0: Lecture echouee pour " CALL abort ENDIF @@ -149,5 +153,5 @@ ierr = NF_INQ_VARID (nid, "rlatu", nvarid) IF (ierr .NE. NF_NOERR) THEN - PRINT*, "dynetat0: Le champ est absent" + write(lunout,*)"dynetat0: Le champ est absent" CALL abort ENDIF @@ -158,5 +162,5 @@ #endif IF (ierr .NE. NF_NOERR) THEN - PRINT*, "dynetat0: Lecture echouee pour " + write(lunout,*)"dynetat0: Lecture echouee pour " CALL abort ENDIF @@ -164,5 +168,5 @@ ierr = NF_INQ_VARID (nid, "rlonv", nvarid) IF (ierr .NE. NF_NOERR) THEN - PRINT*, "dynetat0: Le champ est absent" + write(lunout,*)"dynetat0: Le champ est absent" CALL abort ENDIF @@ -173,5 +177,5 @@ #endif IF (ierr .NE. NF_NOERR) THEN - PRINT*, "dynetat0: Lecture echouee pour " + write(lunout,*)"dynetat0: Lecture echouee pour " CALL abort ENDIF @@ -179,5 +183,5 @@ ierr = NF_INQ_VARID (nid, "rlatv", nvarid) IF (ierr .NE. NF_NOERR) THEN - PRINT*, "dynetat0: Le champ est absent" + write(lunout,*)"dynetat0: Le champ est absent" CALL abort ENDIF @@ -188,5 +192,5 @@ #endif IF (ierr .NE. NF_NOERR) THEN - PRINT*, "dynetat0: Lecture echouee pour rlatv" + write(lunout,*)"dynetat0: Lecture echouee pour rlatv" CALL abort ENDIF @@ -194,5 +198,5 @@ ierr = NF_INQ_VARID (nid, "cu", nvarid) IF (ierr .NE. NF_NOERR) THEN - PRINT*, "dynetat0: Le champ est absent" + write(lunout,*)"dynetat0: Le champ est absent" CALL abort ENDIF @@ -203,5 +207,5 @@ #endif IF (ierr .NE. NF_NOERR) THEN - PRINT*, "dynetat0: Lecture echouee pour " + write(lunout,*)"dynetat0: Lecture echouee pour " CALL abort ENDIF @@ -209,5 +213,5 @@ ierr = NF_INQ_VARID (nid, "cv", nvarid) IF (ierr .NE. NF_NOERR) THEN - PRINT*, "dynetat0: Le champ est absent" + write(lunout,*)"dynetat0: Le champ est absent" CALL abort ENDIF @@ -218,5 +222,5 @@ #endif IF (ierr .NE. NF_NOERR) THEN - PRINT*, "dynetat0: Lecture echouee pour " + write(lunout,*)"dynetat0: Lecture echouee pour " CALL abort ENDIF @@ -224,5 +228,5 @@ ierr = NF_INQ_VARID (nid, "aire", nvarid) IF (ierr .NE. NF_NOERR) THEN - PRINT*, "dynetat0: Le champ est absent" + write(lunout,*)"dynetat0: Le champ est absent" CALL abort ENDIF @@ -233,5 +237,5 @@ #endif IF (ierr .NE. NF_NOERR) THEN - PRINT*, "dynetat0: Lecture echouee pour " + write(lunout,*)"dynetat0: Lecture echouee pour " CALL abort ENDIF @@ -239,5 +243,5 @@ ierr = NF_INQ_VARID (nid, "phisinit", nvarid) IF (ierr .NE. NF_NOERR) THEN - PRINT*, "dynetat0: Le champ est absent" + write(lunout,*)"dynetat0: Le champ est absent" CALL abort ENDIF @@ -248,5 +252,5 @@ #endif IF (ierr .NE. NF_NOERR) THEN - PRINT*, "dynetat0: Lecture echouee pour " + write(lunout,*)"dynetat0: Lecture echouee pour " CALL abort ENDIF @@ -254,5 +258,5 @@ ierr = NF_INQ_VARID (nid, "temps", nvarid) IF (ierr .NE. NF_NOERR) THEN - PRINT*, "dynetat0: Le champ est absent" + write(lunout,*)"dynetat0: Le champ est absent" CALL abort ENDIF @@ -263,5 +267,5 @@ #endif IF (ierr .NE. NF_NOERR) THEN - PRINT*, "dynetat0: Lecture echouee " + write(lunout,*)"dynetat0: Lecture echouee " CALL abort ENDIF @@ -269,5 +273,5 @@ ierr = NF_INQ_VARID (nid, "ucov", nvarid) IF (ierr .NE. NF_NOERR) THEN - PRINT*, "dynetat0: Le champ est absent" + write(lunout,*)"dynetat0: Le champ est absent" CALL abort ENDIF @@ -278,5 +282,5 @@ #endif IF (ierr .NE. NF_NOERR) THEN - PRINT*, "dynetat0: Lecture echouee pour " + write(lunout,*)"dynetat0: Lecture echouee pour " CALL abort ENDIF @@ -284,5 +288,5 @@ ierr = NF_INQ_VARID (nid, "vcov", nvarid) IF (ierr .NE. NF_NOERR) THEN - PRINT*, "dynetat0: Le champ est absent" + write(lunout,*)"dynetat0: Le champ est absent" CALL abort ENDIF @@ -293,5 +297,5 @@ #endif IF (ierr .NE. NF_NOERR) THEN - PRINT*, "dynetat0: Lecture echouee pour " + write(lunout,*)"dynetat0: Lecture echouee pour " CALL abort ENDIF @@ -299,5 +303,5 @@ ierr = NF_INQ_VARID (nid, "teta", nvarid) IF (ierr .NE. NF_NOERR) THEN - PRINT*, "dynetat0: Le champ est absent" + write(lunout,*)"dynetat0: Le champ est absent" CALL abort ENDIF @@ -308,14 +312,16 @@ #endif IF (ierr .NE. NF_NOERR) THEN - PRINT*, "dynetat0: Lecture echouee pour " - CALL abort - ENDIF - - + write(lunout,*)"dynetat0: Lecture echouee pour " + CALL abort + ENDIF + + + IF(nqtot.GE.1) THEN DO iq=1,nqtot ierr = NF_INQ_VARID (nid, tname(iq), nvarid) IF (ierr .NE. NF_NOERR) THEN - PRINT*, "dynetat0: Le champ <"//tname(iq)//"> est absent" - PRINT*, " Il est donc initialise a zero" + write(lunout,*)"dynetat0: Le champ <"//tname(iq)// + & "> est absent" + write(lunout,*)" Il est donc initialise a zero" q(:,:,iq)=0. ELSE @@ -326,13 +332,14 @@ #endif IF (ierr .NE. NF_NOERR) THEN - PRINT*, "dynetat0: Lecture echouee pour "//tname(iq) - CALL abort + write(lunout,*)"dynetat0: Lecture echouee pour "//tname(iq) + CALL abort ENDIF ENDIF ENDDO + ENDIF ierr = NF_INQ_VARID (nid, "masse", nvarid) IF (ierr .NE. NF_NOERR) THEN - PRINT*, "dynetat0: Le champ est absent" + write(lunout,*)"dynetat0: Le champ est absent" CALL abort ENDIF @@ -343,5 +350,5 @@ #endif IF (ierr .NE. NF_NOERR) THEN - PRINT*, "dynetat0: Lecture echouee pour " + write(lunout,*)"dynetat0: Lecture echouee pour " CALL abort ENDIF @@ -349,5 +356,5 @@ ierr = NF_INQ_VARID (nid, "ps", nvarid) IF (ierr .NE. NF_NOERR) THEN - PRINT*, "dynetat0: Le champ est absent" + write(lunout,*)"dynetat0: Le champ est absent" CALL abort ENDIF @@ -358,5 +365,5 @@ #endif IF (ierr .NE. NF_NOERR) THEN - PRINT*, "dynetat0: Lecture echouee pour " + write(lunout,*)"dynetat0: Lecture echouee pour " CALL abort ENDIF Index: LMDZ5/trunk/libf/dyn3dmem/dynetat0_loc.F =================================================================== --- LMDZ5/trunk/libf/dyn3dmem/dynetat0_loc.F (revision 1656) +++ LMDZ5/trunk/libf/dyn3dmem/dynetat0_loc.F (revision 1657) @@ -1,4 +1,4 @@ ! -! $Header$ +! $Id$ ! SUBROUTINE dynetat0_loc(fichnom,vcov,ucov, @@ -34,4 +34,5 @@ #include "serre.h" #include "logic.h" +#include "iniprint.h" c Arguments: @@ -61,6 +62,7 @@ ierr = NF_OPEN (fichnom, NF_NOWRITE,nid) IF (ierr.NE.NF_NOERR) THEN - write(6,*)' Pb d''ouverture du fichier start.nc' - write(6,*)' ierr = ', ierr + write(lunout,*) + & 'dynetat0_loc: Pb d''ouverture du fichier start.nc' + write(lunout,*)' ierr = ', ierr CALL ABORT ENDIF @@ -69,5 +71,5 @@ ierr = NF_INQ_VARID (nid, "controle", nvarid) IF (ierr .NE. NF_NOERR) THEN - PRINT*, "dynetat0: Le champ est absent" + write(lunout,*)"dynetat0_loc: Le champ est absent" CALL abort ENDIF @@ -78,5 +80,5 @@ #endif IF (ierr .NE. NF_NOERR) THEN - PRINT*, "dynetat0: Lecture echoue pour " + write(lunout,*)"dynetat0_loc: Lecture echoue pour " CALL abort ENDIF @@ -124,5 +126,6 @@ c c - PRINT*,'rad,omeg,g,cpp,kappa',rad,omeg,g,cpp,kappa + write(lunout,*)'dynetat0_loc: rad,omeg,g,cpp,kappa', + & rad,omeg,g,cpp,kappa IF( im.ne.iim ) THEN @@ -139,5 +142,5 @@ ierr = NF_INQ_VARID (nid, "rlonu", nvarid) IF (ierr .NE. NF_NOERR) THEN - PRINT*, "dynetat0: Le champ est absent" + write(lunout,*)"dynetat0_loc: Le champ est absent" CALL abort ENDIF @@ -148,5 +151,5 @@ #endif IF (ierr .NE. NF_NOERR) THEN - PRINT*, "dynetat0: Lecture echouee pour " + write(lunout,*)"dynetat0_loc: Lecture echouee pour " CALL abort ENDIF @@ -154,5 +157,5 @@ ierr = NF_INQ_VARID (nid, "rlatu", nvarid) IF (ierr .NE. NF_NOERR) THEN - PRINT*, "dynetat0: Le champ est absent" + write(lunout,*)"dynetat0_loc: Le champ est absent" CALL abort ENDIF @@ -163,5 +166,5 @@ #endif IF (ierr .NE. NF_NOERR) THEN - PRINT*, "dynetat0: Lecture echouee pour " + write(lunout,*)"dynetat0_loc: Lecture echouee pour " CALL abort ENDIF @@ -169,5 +172,5 @@ ierr = NF_INQ_VARID (nid, "rlonv", nvarid) IF (ierr .NE. NF_NOERR) THEN - PRINT*, "dynetat0: Le champ est absent" + write(lunout,*)"dynetat0_loc: Le champ est absent" CALL abort ENDIF @@ -178,5 +181,5 @@ #endif IF (ierr .NE. NF_NOERR) THEN - PRINT*, "dynetat0: Lecture echouee pour " + write(lunout,*)"dynetat0_loc: Lecture echouee pour " CALL abort ENDIF @@ -184,5 +187,5 @@ ierr = NF_INQ_VARID (nid, "rlatv", nvarid) IF (ierr .NE. NF_NOERR) THEN - PRINT*, "dynetat0: Le champ est absent" + write(lunout,*)"dynetat0_loc: Le champ est absent" CALL abort ENDIF @@ -193,5 +196,5 @@ #endif IF (ierr .NE. NF_NOERR) THEN - PRINT*, "dynetat0: Lecture echouee pour rlatv" + write(lunout,*)"dynetat0_loc: Lecture echouee pour rlatv" CALL abort ENDIF @@ -199,5 +202,5 @@ ierr = NF_INQ_VARID (nid, "cu", nvarid) IF (ierr .NE. NF_NOERR) THEN - PRINT*, "dynetat0: Le champ est absent" + write(lunout,*)"dynetat0_loc: Le champ est absent" CALL abort ENDIF @@ -208,5 +211,5 @@ #endif IF (ierr .NE. NF_NOERR) THEN - PRINT*, "dynetat0: Lecture echouee pour " + write(lunout,*)"dynetat0_loc: Lecture echouee pour " CALL abort ENDIF @@ -214,5 +217,5 @@ ierr = NF_INQ_VARID (nid, "cv", nvarid) IF (ierr .NE. NF_NOERR) THEN - PRINT*, "dynetat0: Le champ est absent" + write(lunout,*)"dynetat0_loc: Le champ est absent" CALL abort ENDIF @@ -223,5 +226,5 @@ #endif IF (ierr .NE. NF_NOERR) THEN - PRINT*, "dynetat0: Lecture echouee pour " + write(lunout,*)"dynetat0_loc: Lecture echouee pour " CALL abort ENDIF @@ -229,5 +232,5 @@ ierr = NF_INQ_VARID (nid, "aire", nvarid) IF (ierr .NE. NF_NOERR) THEN - PRINT*, "dynetat0: Le champ est absent" + write(lunout,*)"dynetat0_loc: Le champ est absent" CALL abort ENDIF @@ -238,5 +241,5 @@ #endif IF (ierr .NE. NF_NOERR) THEN - PRINT*, "dynetat0: Lecture echouee pour " + write(lunout,*)"dynetat0_loc: Lecture echouee pour " CALL abort ENDIF @@ -246,5 +249,5 @@ ierr = NF_INQ_VARID (nid, "phisinit", nvarid) IF (ierr .NE. NF_NOERR) THEN - PRINT*, "dynetat0: Le champ est absent" + write(lunout,*)"dynetat0_loc: Le champ est absent" CALL abort ENDIF @@ -255,5 +258,5 @@ #endif IF (ierr .NE. NF_NOERR) THEN - PRINT*, "dynetat0: Lecture echouee pour " + write(lunout,*)"dynetat0_loc: Lecture echouee pour " CALL abort ENDIF @@ -263,5 +266,5 @@ ierr = NF_INQ_VARID (nid, "temps", nvarid) IF (ierr .NE. NF_NOERR) THEN - PRINT*, "dynetat0: Le champ est absent" + write(lunout,*)"dynetat0_loc: Le champ est absent" CALL abort ENDIF @@ -272,5 +275,5 @@ #endif IF (ierr .NE. NF_NOERR) THEN - PRINT*, "dynetat0: Lecture echouee " + write(lunout,*)"dynetat0_loc: Lecture echouee " CALL abort ENDIF @@ -278,5 +281,5 @@ ierr = NF_INQ_VARID (nid, "ucov", nvarid) IF (ierr .NE. NF_NOERR) THEN - PRINT*, "dynetat0: Le champ est absent" + write(lunout,*)"dynetat0_loc: Le champ est absent" CALL abort ENDIF @@ -290,5 +293,5 @@ #endif IF (ierr .NE. NF_NOERR) THEN - PRINT*, "dynetat0: Lecture echouee pour " + write(lunout,*)"dynetat0_loc: Lecture echouee pour " CALL abort ENDIF @@ -300,5 +303,5 @@ ierr = NF_INQ_VARID (nid, "vcov", nvarid) IF (ierr .NE. NF_NOERR) THEN - PRINT*, "dynetat0: Le champ est absent" + write(lunout,*)"dynetat0_loc: Le champ est absent" CALL abort ENDIF @@ -309,5 +312,5 @@ #endif IF (ierr .NE. NF_NOERR) THEN - PRINT*, "dynetat0: Lecture echouee pour " + write(lunout,*)"dynetat0_loc: Lecture echouee pour " CALL abort ENDIF @@ -318,5 +321,5 @@ ierr = NF_INQ_VARID (nid, "teta", nvarid) IF (ierr .NE. NF_NOERR) THEN - PRINT*, "dynetat0: Le champ est absent" + write(lunout,*)"dynetat0_loc: Le champ est absent" CALL abort ENDIF @@ -327,5 +330,5 @@ #endif IF (ierr .NE. NF_NOERR) THEN - PRINT*, "dynetat0: Lecture echouee pour " + write(lunout,*)"dynetat0_loc: Lecture echouee pour " CALL abort ENDIF @@ -339,7 +342,8 @@ ierr = NF_INQ_VARID (nid, tname(iq), nvarid) IF (ierr .NE. NF_NOERR) THEN - PRINT*, "dynetat0: Le champ <"//tname(iq)//"> est absent" - PRINT*, " Il est donc initialise a zero" - q_glo(:,:)=0. + write(lunout,*)"dynetat0_loc: Le champ <"//tname(iq)// + & "> est absent" + write(lunout,*)" Il est donc initialise a zero" + q(:,:,iq)=0. ELSE #ifdef NC_DOUBLE @@ -349,6 +353,7 @@ #endif IF (ierr .NE. NF_NOERR) THEN - PRINT*, "dynetat0: Lecture echouee pour "//tname(iq) - CALL abort + write(lunout,*) + & "dynetat0_loc: Lecture echouee pour "//tname(iq) + CALL abort ENDIF ENDIF @@ -361,5 +366,5 @@ ierr = NF_INQ_VARID (nid, "masse", nvarid) IF (ierr .NE. NF_NOERR) THEN - PRINT*, "dynetat0: Le champ est absent" + write(lunout,*)"dynetat0_loc: Le champ est absent" CALL abort ENDIF @@ -370,5 +375,5 @@ #endif IF (ierr .NE. NF_NOERR) THEN - PRINT*, "dynetat0: Lecture echouee pour " + write(lunout,*)"dynetat0_loc: Lecture echouee pour " CALL abort ENDIF @@ -379,5 +384,5 @@ ierr = NF_INQ_VARID (nid, "ps", nvarid) IF (ierr .NE. NF_NOERR) THEN - PRINT*, "dynetat0: Le champ est absent" + write(lunout,*)"dynetat0_loc: Le champ est absent" CALL abort ENDIF @@ -388,5 +393,5 @@ #endif IF (ierr .NE. NF_NOERR) THEN - PRINT*, "dynetat0: Lecture echouee pour " + write(lunout,*)"dynetat0_loc: Lecture echouee pour " CALL abort ENDIF Index: LMDZ5/trunk/libf/dyn3dmem/dynredem_loc.F =================================================================== --- LMDZ5/trunk/libf/dyn3dmem/dynredem_loc.F (revision 1656) +++ LMDZ5/trunk/libf/dyn3dmem/dynredem_loc.F (revision 1657) @@ -27,4 +27,5 @@ #include "description.h" #include "serre.h" +#include "iniprint.h" c Arguments: @@ -64,5 +65,5 @@ if (mpi_rank==0) then - modname='dynredem0_p' + modname='dynredem0_loc' #ifdef CPP_IOIPSL @@ -132,6 +133,7 @@ ierr = NF_CREATE(fichnom, NF_CLOBBER, nid) IF (ierr.NE.NF_NOERR) THEN - WRITE(6,*)" Pb d ouverture du fichier "//fichnom - WRITE(6,*)' ierr = ', ierr + write(lunout,*)"dynredem0: Pb d ouverture du fichier " + & //trim(fichnom) + write(lunout,*)' ierr = ', ierr CALL ABORT ENDIF @@ -514,8 +516,8 @@ ierr = NF_CLOSE(nid) ! fermer le fichier - - PRINT*,'iim,jjm,llm,iday_end',iim,jjm,llm,iday_end - PRINT*,'rad,omeg,g,cpp,kappa', - , rad,omeg,g,cpp,kappa + write(lunout,*)'dynredem_loc: iim,jjm,llm,iday_end', + & iim,jjm,llm,iday_end + write(lunout,*)'dynredem_loc: rad,omeg,g,cpp,kappa', + & rad,omeg,g,cpp,kappa endif ! mpi_rank==0 @@ -540,4 +542,5 @@ #include "comgeom.h" #include "temps.h" +#include "iniprint.h" INTEGER l @@ -579,8 +582,8 @@ !$OMP MASTER if (mpi_rank==0) then - modname = 'dynredem1' + modname = 'dynredem1_loc' ierr = NF_OPEN(fichnom, NF_WRITE, nid) IF (ierr .NE. NF_NOERR) THEN - PRINT*, "Pb. d ouverture "//fichnom + write(lunout,*)"dynredem1: Pb. d ouverture "//trim(fichnom) CALL abort ENDIF @@ -591,5 +594,5 @@ ierr = NF_INQ_VARID(nid, "temps", nvarid) IF (ierr .NE. NF_NOERR) THEN - print *, NF_STRERROR(ierr) + write(lunout,*) NF_STRERROR(ierr) abort_message='Variable temps n est pas definie' CALL abort_gcm(modname,abort_message,ierr) @@ -600,5 +603,5 @@ ierr = NF_PUT_VAR1_REAL (nid,nvarid,nb,time) #endif - PRINT*, "Enregistrement pour ", nb, time + write(lunout,*) "dynredem1_loc: Enregistrement pour ", nb, time c Index: LMDZ5/trunk/libf/dyn3dmem/exner_hyb.F =================================================================== --- LMDZ5/trunk/libf/dyn3dmem/exner_hyb.F (revision 1656) +++ LMDZ5/trunk/libf/dyn3dmem/exner_hyb.F (revision 1657) @@ -1,4 +1,4 @@ ! -! $Header$ +! $Id$ ! SUBROUTINE exner_hyb ( ngrid, ps, p,alpha,beta, pks, pk, pkf ) @@ -51,5 +51,33 @@ REAL SSUM c + + if (llm.eq.1) then + ! Specific behaviour for Shallow Water (1 vertical layer) case + ! Sanity checks + if (kappa.ne.1) then + call abort_gcm("exner_hyb", + & "kappa!=1 , but running in Shallow Water mode!!",42) + endif + if (cpp.ne.r) then + call abort_gcm("exner_hyb", + & "cpp!=r , but running in Shallow Water mode!!",42) + endif + + ! Compute pks(:),pk(:),pkf(:) + + DO ij = 1, ngrid + pks(ij) = (cpp/preff) * ps(ij) + pk(ij,1) = .5*pks(ij) + ENDDO + + CALL SCOPY ( ngrid * llm, pk, 1, pkf, 1 ) + CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 ) + + ! our work is done, exit routine + return + endif ! of if (llm.eq.1) + + unpl2k = 1.+ 2.* kappa c Index: LMDZ5/trunk/libf/dyn3dmem/exner_hyb_loc.F =================================================================== --- LMDZ5/trunk/libf/dyn3dmem/exner_hyb_loc.F (revision 1656) +++ LMDZ5/trunk/libf/dyn3dmem/exner_hyb_loc.F (revision 1657) @@ -1,2 +1,5 @@ +c +c $Id$ +c SUBROUTINE exner_hyb_loc(ngrid, ps, p,alpha,beta, pks,pk,pkf) c @@ -56,4 +59,67 @@ c c$OMP BARRIER + + if (llm.eq.1) then + ! Specific behaviour for Shallow Water (1 vertical layer) case + + ! Sanity checks + if (kappa.ne.1) then + call abort_gcm("exner_hyb", + & "kappa!=1 , but running in Shallow Water mode!!",42) + endif + if (cpp.ne.r) then + call abort_gcm("exner_hyb", + & "cpp!=r , but running in Shallow Water mode!!",42) + endif + + ! Compute pks(:),pk(:),pkf(:) + ijb=ij_begin + ije=ij_end +!$OMP DO SCHEDULE(STATIC) + DO ij=ijb, ije + pks(ij)=(cpp/preff)*ps(ij) + pk(ij,1) = .5*pks(ij) + pkf(ij,1)=pk(ij,1) + ENDDO +!$OMP ENDDO + +!$OMP MASTER + if (pole_nord) then + DO ij = 1, iim + ppn(ij) = aire( ij ) * pks( ij ) + ENDDO + xpn = SSUM(iim,ppn,1) /apoln + + DO ij = 1, iip1 + pks( ij ) = xpn + pk(ij,1) = .5*pks(ij) + pkf(ij,1)=pk(ij,1) + ENDDO + endif + + if (pole_sud) then + DO ij = 1, iim + pps(ij) = aire(ij+ip1jm) * pks(ij+ip1jm ) + ENDDO + xps = SSUM(iim,pps,1) /apols + + DO ij = 1, iip1 + pks( ij+ip1jm ) = xps + pk(ij+ip1jm,1)=.5*pks(ij+ip1jm) + pkf(ij+ip1jm,1)=pk(ij+ip1jm,1) + ENDDO + endif +!$OMP END MASTER + + jjb=jj_begin + jje=jj_end + CALL filtreg_p ( pkf,jjb_u,jje_u,jjb,jje, jmp1, llm, + & 2, 1, .TRUE., 1 ) + + ! our work is done, exit routine + return + endif ! of if (llm.eq.1) + + unpl2k = 1.+ 2.* kappa c Index: LMDZ5/trunk/libf/dyn3dmem/gcm.F =================================================================== --- LMDZ5/trunk/libf/dyn3dmem/gcm.F (revision 1656) +++ LMDZ5/trunk/libf/dyn3dmem/gcm.F (revision 1657) @@ -1,4 +1,4 @@ ! -! $Id: gcm.F 1316 2010-02-22 14:51:12Z acozic $ +! $Id: gcm.F 1397 2010-06-02 12:57:39Z emillour $ ! c @@ -70,8 +70,11 @@ #include "description.h" #include "serre.h" -#include "com_io_dyn.h" +!#include "com_io_dyn.h" #include "iniprint.h" #include "tracstoke.h" +#ifdef INCA +! Only INCA needs these informations (from the Earth's physics) #include "indicesol.h" +#endif INTEGER longcles @@ -273,14 +276,21 @@ if (read_start) then ! we still need to run iniacademic to initialize some - ! constants & fields, if we run the 'newtonian' case: - if (iflag_phys.eq.2) then + ! constants & fields, if we run the 'newtonian' or 'SW' cases: + if (iflag_phys.ne.1) then CALL iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0) endif -!#ifdef CPP_IOIPSL + if (planet_type.eq."earth") then #ifdef CPP_EARTH ! Load an Earth-format start file CALL dynetat0_loc("start.nc",vcov,ucov, - . teta,q,masse,ps,phis, time_0) + & teta,q,masse,ps,phis, time_0) +#else + ! SW model also has Earth-format start files + ! (but can be used without the CPP_EARTH directive) + if (iflag_phys.eq.0) then + CALL dynetat0_loc("start.nc",vcov,ucov, + & teta,q,masse,ps,phis, time_0) + endif #endif endif ! of if (planet_type.eq."earth") @@ -326,5 +336,14 @@ C on remet le calendrier à zero si demande c - if (annee_ref .ne. anneeref .or. day_ref .ne. dayref) then + IF (raz_date == 1) THEN + annee_ref = anneeref + day_ref = dayref + day_ini = dayref + itau_dyn = 0 + itau_phy = 0 + time_0 = 0. + write(lunout,*) + . 'GCM: On reinitialise a la date lue dans gcm.def' + ELSE IF (annee_ref .ne. anneeref .or. day_ref .ne. dayref) THEN write(lunout,*) . 'GCM: Attention les dates initiales lues dans le fichier' @@ -332,22 +351,32 @@ . ' restart ne correspondent pas a celles lues dans ' write(lunout,*)' gcm.def' - write(lunout,*)' annee_ref=',annee_ref," anneeref=",anneeref - write(lunout,*)' day_ref=',day_ref," dayref=",dayref - if (raz_date .ne. 1) then - write(lunout,*) - . 'GCM: On garde les dates du fichier restart' - else - annee_ref = anneeref - day_ref = dayref - day_ini = dayref - itau_dyn = 0 - itau_phy = 0 - time_0 = 0. - write(lunout,*) - . 'GCM: On reinitialise a la date lue dans gcm.def' - endif - ELSE - raz_date = 0 - endif + write(lunout,*)' annee_ref=',annee_ref," anneeref=",anneeref + write(lunout,*)' day_ref=',day_ref," dayref=",dayref + write(lunout,*)' Pas de remise a zero' + ENDIF +c if (annee_ref .ne. anneeref .or. day_ref .ne. dayref) then +c write(lunout,*) +c . 'GCM: Attention les dates initiales lues dans le fichier' +c write(lunout,*) +c . ' restart ne correspondent pas a celles lues dans ' +c write(lunout,*)' gcm.def' +c write(lunout,*)' annee_ref=',annee_ref," anneeref=",anneeref +c write(lunout,*)' day_ref=',day_ref," dayref=",dayref +c if (raz_date .ne. 1) then +c write(lunout,*) +c . 'GCM: On garde les dates du fichier restart' +c else +c annee_ref = anneeref +c day_ref = dayref +c day_ini = dayref +c itau_dyn = 0 +c itau_phy = 0 +c time_0 = 0. +c write(lunout,*) +c . 'GCM: On reinitialise a la date lue dans gcm.def' +c endif +c ELSE +c raz_date = 0 +c endif #ifdef CPP_IOIPSL @@ -428,5 +457,5 @@ if (planet_type.eq."earth") then #ifdef CPP_EARTH - CALL iniphysiq(ngridmx,llm,daysec,day_ini,dtphys , + CALL iniphysiq(ngridmx,llm,daysec,day_ini,dtphys/nsplit_phys , , latfi,lonfi,airefi,zcufi,zcvfi,rad,g,r,cpp ) #endif @@ -473,20 +502,24 @@ #ifdef CPP_IOIPSL - if ( 1.eq.1) then time_step = zdtvr - t_ops = iecri * daysec - t_wrt = iecri * daysec -! CALL inithist_p(dynhist_file,day_ref,annee_ref,time_step, -! . t_ops, t_wrt, histid, histvid) + IF (mpi_rank==0) then + if (ok_dyn_ins) then + ! initialize output file for instantaneous outputs + ! t_ops = iecri * daysec ! do operations every t_ops + t_ops =((1.0*iecri)/day_step) * daysec + t_wrt = daysec ! iecri * daysec ! write output every t_wrt + t_wrt = daysec ! iecri * daysec ! write output every t_wrt + CALL inithist(day_ref,annee_ref,time_step, + & t_ops,t_wrt) + endif IF (ok_dyn_ave) THEN - t_ops = iperiod * time_step - t_wrt = periodav * daysec + ! initialize output file for averaged outputs + t_ops = iperiod * time_step ! do operations every t_ops + t_wrt = periodav * daysec ! write output every t_wrt CALL initdynav_loc(day_ref,annee_ref,time_step,t_ops,t_wrt) - END IF + END IF + ENDIF dtav = iperiod*dtvr/daysec - endif - - #endif ! #endif of #ifdef CPP_IOIPSL Index: LMDZ5/trunk/libf/dyn3dmem/iniacademic.F =================================================================== --- LMDZ5/trunk/libf/dyn3dmem/iniacademic.F (revision 1656) +++ LMDZ5/trunk/libf/dyn3dmem/iniacademic.F (revision 1657) @@ -1,4 +1,4 @@ ! -! $Id: iniacademic.F 1299 2010-01-20 14:27:21Z fairhead $ +! $Id: iniacademic.F 1363 2010-04-16 09:50:10Z emillour $ ! c @@ -8,4 +8,6 @@ USE filtreg_mod USE infotrac, ONLY : nqtot + USE control_mod + c%W% %G% @@ -31,5 +33,4 @@ c c======================================================================= - USE control_mod IMPLICIT NONE c----------------------------------------------------------------------- @@ -46,4 +47,5 @@ #include "temps.h" #include "iniprint.h" +#include "logic.h" c Arguments: @@ -55,5 +57,5 @@ REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants REAL teta(ip1jmp1,llm) ! temperature potentielle - REAL q(ip1jmp1,llm,nqtot) ! champs advectes + REAL q(ip1jmp1,llm,nqtot) ! champs advectes REAL ps(ip1jmp1) ! pression au sol REAL masse(ip1jmp1,llm) ! masse d'air @@ -84,5 +86,5 @@ time_0=0. day_ref=0 - annee_ref=0 + annee_ref=0 im = iim @@ -105,13 +107,31 @@ ang0 = 0. + if (llm.eq.1) then + ! specific initializations for the shallow water case + kappa=1 + endif + CALL iniconst CALL inigeom CALL inifilr - ps=0. - phis=0. + if (llm.eq.1) then + ! initialize fields for the shallow water case, if required + if (.not.read_start) then + phis(:)=0. + q(:,:,1)=1.e-10 + q(:,:,2)=1.e-15 + q(:,:,3:nqtot)=0. + CALL sw_case_williamson91_6(vcov,ucov,teta,masse,ps) + endif + endif + + if (iflag_phys.eq.2) then + ! initializations for the academic case + ps(:)=1.e5 + phis(:)=0. c--------------------------------------------------------------------- - taurappel=10.*daysec + taurappel=10.*daysec c--------------------------------------------------------------------- @@ -119,63 +139,61 @@ c -------------------------------------- - DO l=1,llm - zsig=ap(l)/preff+bp(l) - if (zsig.gt.0.3) then - lsup=l - tetarappell=1./8.*(-log(zsig)-.5) - DO j=1,jjp1 + DO l=1,llm + zsig=ap(l)/preff+bp(l) + if (zsig.gt.0.3) then + lsup=l + tetarappell=1./8.*(-log(zsig)-.5) + DO j=1,jjp1 ddsin=sin(rlatu(j))-sin(pi/20.) tetajl(j,l)=300.*(1+1./18.*(1.-3.*ddsin*ddsin)+tetarappell) - ENDDO - else + ENDDO + else c Choix isotherme au-dessus de 300 mbar - do j=1,jjp1 - tetajl(j,l)=tetajl(j,lsup)*(0.3/zsig)**kappa - enddo - endif ! of if (zsig.gt.0.3) - ENDDO ! of DO l=1,llm - - do l=1,llm - do j=1,jjp1 + do j=1,jjp1 + tetajl(j,l)=tetajl(j,lsup)*(0.3/zsig)**kappa + enddo + endif ! of if (zsig.gt.0.3) + ENDDO ! of DO l=1,llm + + do l=1,llm + do j=1,jjp1 do i=1,iip1 ij=(j-1)*iip1+i tetarappel(ij,l)=tetajl(j,l) enddo - enddo - enddo + enddo + enddo c call dump2d(jjp1,llm,tetajl,'TEQ ') - ps=1.e5 - phis=0. - CALL pression ( ip1jmp1, ap, bp, ps, p ) - CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf ) - CALL massdair(p,masse) + CALL pression ( ip1jmp1, ap, bp, ps, p ) + CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf ) + CALL massdair(p,masse) c intialisation du vent et de la temperature - teta(:,:)=tetarappel(:,:) - CALL geopot(ip1jmp1,teta,pk,pks,phis,phi) - call ugeostr(phi,ucov) - vcov=0. - q(:,:,1 )=1.e-10 - q(:,:,2 )=1.e-15 - q(:,:,3:nqtot)=0. + teta(:,:)=tetarappel(:,:) + CALL geopot(ip1jmp1,teta,pk,pks,phis,phi) + call ugeostr(phi,ucov) + vcov=0. + q(:,:,1 )=1.e-10 + q(:,:,2 )=1.e-15 + q(:,:,3:nqtot)=0. c perturbation aleatoire sur la temperature - idum = -1 - zz = ran1(idum) - idum = 0 - do l=1,llm - do ij=iip2,ip1jm + idum = -1 + zz = ran1(idum) + idum = 0 + do l=1,llm + do ij=iip2,ip1jm teta(ij,l)=teta(ij,l)*(1.+0.005*ran1(idum)) - enddo - enddo - - do l=1,llm - do ij=1,ip1jmp1,iip1 + enddo + enddo + + do l=1,llm + do ij=1,ip1jmp1,iip1 teta(ij+iim,l)=teta(ij,l) - enddo - enddo + enddo + enddo @@ -187,9 +205,10 @@ c initialisation d'un traceur sur une colonne - j=jjp1*3/4 - i=iip1/2 - ij=(j-1)*iip1+i - q(ij,:,3)=1. - + j=jjp1*3/4 + i=iip1/2 + ij=(j-1)*iip1+i + q(ij,:,3)=1. + endif ! of if (iflag_phys.eq.2) + else write(lunout,*)"iniacademic: planet types other than earth", Index: LMDZ5/trunk/libf/dyn3dmem/iniconst.F =================================================================== --- LMDZ5/trunk/libf/dyn3dmem/iniconst.F (revision 1656) +++ LMDZ5/trunk/libf/dyn3dmem/iniconst.F (revision 1657) @@ -1,8 +1,9 @@ ! -! $Id: iniconst.F 1299 2010-01-20 14:27:21Z fairhead $ +! $Id: iniconst.F 1380 2010-05-06 12:24:59Z emillour $ ! SUBROUTINE iniconst USE control_mod + IMPLICIT NONE c @@ -18,4 +19,5 @@ #include "temps.h" #include "comvert.h" +#include "iniprint.h" @@ -53,5 +55,5 @@ r = cpp * kappa - PRINT*,' R CP Kappa ', r , cpp, kappa + write(lunout,*)'iniconst: R CP Kappa ', r , cpp, kappa c c----------------------------------------------------------------------- Index: LMDZ5/trunk/libf/dyn3dmem/inter_barxy_m.F90 =================================================================== --- LMDZ5/trunk/libf/dyn3dmem/inter_barxy_m.F90 (revision 1656) +++ LMDZ5/trunk/libf/dyn3dmem/inter_barxy_m.F90 (revision 1657) @@ -118,7 +118,7 @@ IMPLICIT NONE - REAL, intent(in):: dlonid(:) - real, intent(in):: fdat(:) - real, intent(in):: rlonimod(:) + REAL, intent(in):: dlonid(:) ! dim(idatmax) + real, intent(in):: fdat(:) ! dim(idatmax) + real, intent(in):: rlonimod(:) ! dim(imodmax) real inter_barx(size(rlonimod)) @@ -176,5 +176,5 @@ DO idat = 1, idatmax - xxd(idat) = AMOD( xxd(idat) - xim0, 360. ) + xxd(idat) = MOD( xxd(idat) - xim0, 360. ) fdd(idat) = fdat (idat) ENDDO Index: LMDZ5/trunk/libf/dyn3dmem/leapfrog_loc.F =================================================================== --- LMDZ5/trunk/libf/dyn3dmem/leapfrog_loc.F (revision 1656) +++ LMDZ5/trunk/libf/dyn3dmem/leapfrog_loc.F (revision 1657) @@ -75,5 +75,5 @@ #include "description.h" #include "serre.h" -#include "com_io_dyn.h" +!#include "com_io_dyn.h" #include "iniprint.h" #include "academic.h" @@ -397,4 +397,5 @@ IF( purmats ) THEN + ! Purely Matsuno time stepping IF( MOD(itau,iconser) .EQ.0.AND. forward ) conser = .TRUE. IF( MOD(itau,idissip ).EQ.0.AND..NOT.forward ) apdiss = .TRUE. @@ -402,4 +403,5 @@ s .and. iflag_phys.EQ.1 ) apphys = .TRUE. ELSE + ! Leapfrog/Matsuno time stepping IF( MOD(itau ,iconser) .EQ. 0 ) conser = .TRUE. IF( MOD(itau+1,idissip) .EQ. 0 ) apdiss = .TRUE. @@ -407,8 +409,14 @@ END IF +! Ehouarn: for Shallow Water case (ie: 1 vertical layer), +! supress dissipation step + if (llm.eq.1) then + apdiss=.false. + endif + cym ---> Pour le moment cym apphys = .FALSE. statcl = .FALSE. - conser = .FALSE. + conser = .FALSE. ! ie: no output of control variables to stdout in // if (firstCaldyn) then @@ -1069,6 +1077,12 @@ ijb=ij_begin ije=ij_end - teta(ijb:ije,:)=teta(ijb:ije,:) - s -iphysiq*dtvr*(teta(ijb:ije,:)-tetarappel(ijb:ije,:))/taurappel +!LF teta(ijb:ije,:)=teta(ijb:ije,:) +!LF s -iphysiq*dtvr*(teta(ijb:ije,:)-tetarappel(ijb:ije,:))/taurappel +!$OMP DO SCHEDULE(STATIC,OMP_CHUNK) + do l=1,llm + teta(ijb:ije,l)=teta(ijb:ije,l) + & -iphysiq*dtvr*(teta(ijb:ije,l)-tetarappel(ijb:ije,l))/taurappel + enddo +!$OMP END DO call Register_Hallo_u(ucov,llm,0,1,1,0,Request_Physic) @@ -1077,6 +1091,7 @@ c$OMP BARRIER call WaitRequest(Request_Physic) - +c$OMP BARRIER call friction_loc(ucov,vcov,iphysiq*dtvr) +!$OMP BARRIER ENDIF ! of IF(iflag_phys.EQ.2) @@ -1312,5 +1327,5 @@ !c$OMP END MASTER !c$OMP BARRIER - END IF + END IF ! of IF(apdiss) cc$OMP END PARALLEL Index: LMDZ5/trunk/libf/dyn3dmem/sw_case_williamson91_6.F =================================================================== --- LMDZ5/trunk/libf/dyn3dmem/sw_case_williamson91_6.F (revision 1657) +++ LMDZ5/trunk/libf/dyn3dmem/sw_case_williamson91_6.F (revision 1657) @@ -0,0 +1,140 @@ +! +! $Id $ +! + SUBROUTINE sw_case_williamson91_6(vcov,ucov,teta,masse,ps) + +c======================================================================= +c +c Author: Thomas Dubos original: 26/01/2010 +c ------- +c +c Subject: +c ------ +c Realise le cas-test 6 de Williamson et al. (1991) : onde de Rossby-Haurwitz +c +c Method: +c -------- +c +c Interface: +c ---------- +c +c Input: +c ------ +c +c Output: +c ------- +c +c======================================================================= + IMPLICIT NONE +c----------------------------------------------------------------------- +c Declararations: +c --------------- + +#include "dimensions.h" +#include "paramet.h" +#include "comvert.h" +#include "comconst.h" +#include "comgeom.h" +#include "iniprint.h" + +c Arguments: +c ---------- + +c variables dynamiques + REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants + REAL teta(ip1jmp1,llm) ! temperature potentielle + REAL ps(ip1jmp1) ! pression au sol + REAL masse(ip1jmp1,llm) ! masse d'air + REAL phis(ip1jmp1) ! geopotentiel au sol + +c Local: +c ------ + + REAL p (ip1jmp1,llmp1 ) ! pression aux interfac.des couches + REAL pks(ip1jmp1) ! exner au sol + REAL pk(ip1jmp1,llm) ! exner au milieu des couches + REAL pkf(ip1jmp1,llm) ! exner filt.au milieu des couches + REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm) + + REAL :: sinth,costh,costh2, Ath,Bth,Cth, lon,dps + INTEGER i,j,ij + + REAL, PARAMETER :: rho=1 ! masse volumique de l'air (arbitraire) + REAL, PARAMETER :: K = 7.848e-6 ! K = \omega + REAL, PARAMETER :: gh0 = 9.80616 * 8e3 + INTEGER, PARAMETER :: R0=4, R1=R0+1, R2=R0+2 ! mode 4 +c NB : rad = 6371220 dans W91 (6371229 dans LMDZ) +c omeg = 7.292e-5 dans W91 (7.2722e-5 dans LMDZ) + + IF(0==0) THEN +c Williamson et al. (1991) : onde de Rossby-Haurwitz + teta = preff/rho/cpp +c geopotentiel (pression de surface) + do j=1,jjp1 + costh2 = cos(rlatu(j))**2 + Ath = (R0+1)*(costh2**2) + (2*R0*R0-R0-2)*costh2 - 2*R0*R0 + Ath = .25*(K**2)*(costh2**(R0-1))*Ath + Ath = .5*K*(2*omeg+K)*costh2 + Ath + Bth = (R1*R1+1)-R1*R1*costh2 + Bth = 2*(omeg+K)*K/(R1*R2) * (costh2**(R0/2))*Bth + Cth = R1*costh2 - R2 + Cth = .25*K*K*(costh2**R0)*Cth + do i=1,iip1 + ij=(j-1)*iip1+i + lon = rlonv(i) + dps = Ath + Bth*cos(R0*lon) + Cth*cos(2*R0*lon) + ps(ij) = rho*(gh0 + (rad**2)*dps) + enddo + enddo + write(lunout,*) 'W91 ps', MAXVAL(ps), MINVAL(ps) +c vitesse zonale ucov + do j=1,jjp1 + costh = cos(rlatu(j)) + costh2 = costh**2 + Ath = rad*K*costh + Bth = R0*(1-costh2)-costh2 + Bth = rad*K*Bth*(costh**(R0-1)) + do i=1,iip1 + ij=(j-1)*iip1+i + lon = rlonu(i) + ucov(ij,1) = (Ath + Bth*cos(R0*lon)) + enddo + enddo + write(lunout,*) 'W91 u', MAXVAL(ucov(:,1)), MINVAL(ucov(:,1)) + ucov(:,1)=ucov(:,1)*cu +c vitesse meridienne vcov + do j=1,jjm + sinth = sin(rlatv(j)) + costh = cos(rlatv(j)) + Ath = -rad*K*R0*sinth*(costh**(R0-1)) + do i=1,iip1 + ij=(j-1)*iip1+i + lon = rlonv(i) + vcov(ij,1) = Ath*sin(R0*lon) + enddo + enddo + write(lunout,*) 'W91 v', MAXVAL(vcov(:,1)), MINVAL(vcov(:,1)) + vcov(:,1)=vcov(:,1)*cv + +c ucov=0 +c vcov=0 + ELSE +c test non-tournant, onde se propageant en latitude + do j=1,jjp1 + do i=1,iip1 + ij=(j-1)*iip1+i + ps(ij) = 1e5*(1 + .1*exp(-100*(1+sin(rlatu(j)))**2) ) + enddo + enddo + +c rho = preff/(cpp*teta) + teta = .01*preff/cpp ! rho = 100 ; phi = ps/rho = 1e3 ; c=30 m/s = 2600 km/j = 23 degres / j + ucov=0. + vcov=0. + END IF + + CALL pression ( ip1jmp1, ap, bp, ps, p ) + CALL massdair(p,masse) + + END +c-----------------------------------------------------------------------