! $Id: limit_read_mod.F90 3435 2019-01-22 15:21:59Z fairhead $ MODULE limit_read_mod ! This module reads the fichier "limit.nc" containing fields for surface forcing. ! Module subroutines : ! limit_read_frac : CALL limit_read_tot and return the fractions ! limit_read_rug_alb : return rugosity and albedo, if coupled ocean CALL limit_read_tot first ! limit_read_sst : return sea ice temperature ! limit_read_tot : read limit.nc and store the fields in local modules variables IMPLICIT NONE REAL, ALLOCATABLE, DIMENSION(:,:), SAVE, PRIVATE :: pctsrf !$OMP THREADPRIVATE(pctsrf) REAL, ALLOCATABLE, DIMENSION(:), SAVE, PRIVATE :: rugos !$OMP THREADPRIVATE(rugos) REAL, ALLOCATABLE, DIMENSION(:), SAVE, PRIVATE :: albedo !$OMP THREADPRIVATE(albedo) REAL, ALLOCATABLE, DIMENSION(:), SAVE, PRIVATE :: sst !$OMP THREADPRIVATE(sst) #ifdef ISO REAL, ALLOCATABLE, DIMENSION(:), SAVE, PRIVATE :: tuoce !$OMP THREADPRIVATE(tuoce) #endif LOGICAL,SAVE :: read_continents=.FALSE. !$OMP THREADPRIVATE(read_continents) CONTAINS !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! !! Public subroutines : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE init_limit_read(first_day) USE lmdz_grid_phy USE surface_data USE lmdz_phys_para USE lmdz_xios IMPLICIT NONE INTEGER, INTENT(IN) :: first_day IF ( type_ocean /= 'couple') THEN IF (grid_type==unstructured) THEN IF (is_omp_master) CALL xios_set_file_attr("limit_read",enabled=.TRUE.,record_offset=first_day) ENDIF ENDIF END SUBROUTINE init_limit_read SUBROUTINE limit_read_frac(itime, dtime, jour, pctsrf_new, is_modified) ! This SUBROUTINE is called from "change_srf_frac" for case of ! ocean=force or from ocean_slab_frac for ocean=slab. ! The fraction for all sub-surfaces at actual time step is returned. USE dimphy USE indice_sol_mod ! Input arguments !**************************************************************************************** INTEGER, INTENT(IN) :: itime ! time step INTEGER, INTENT(IN) :: jour ! current day REAL , INTENT(IN) :: dtime ! length of time step ! Output arguments !**************************************************************************************** REAL, DIMENSION(klon,nbsrf), INTENT(OUT) :: pctsrf_new ! sub surface fractions LOGICAL, INTENT(OUT) :: is_modified ! true if pctsrf is modified at this time step ! End declaration !**************************************************************************************** ! 1) Read file limit.nc CALL limit_read_tot(itime, dtime, jour, is_modified) ! 2) Return the fraction read in limit_read_tot pctsrf_new(:,:) = pctsrf(:,:) END SUBROUTINE limit_read_frac !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE limit_read_rug_alb(itime, dtime, jour, & knon, knindex, & rugos_out, alb_out) ! This SUBROUTINE is called from surf_land_bucket. ! The flag "ok_veget" must can not be true. If coupled run, "ocean=couple" ! then this routine will CALL limit_read_tot. USE dimphy USE surface_data #ifdef ISO USE isotopes_mod, ONLY: P_veg #endif ! Input arguments !**************************************************************************************** INTEGER, INTENT(IN) :: itime ! numero du pas de temps courant INTEGER, INTENT(IN) :: jour ! jour a lire dans l'annee REAL , INTENT(IN) :: dtime ! pas de temps de la physique (en s) INTEGER, INTENT(IN) :: knon ! nomber of points on compressed grid INTEGER, DIMENSION(klon), INTENT(IN) :: knindex ! grid point number for compressed grid ! Output arguments !**************************************************************************************** REAL, DIMENSION(klon), INTENT(OUT) :: rugos_out REAL, DIMENSION(klon), INTENT(OUT) :: alb_out ! Local variables !**************************************************************************************** INTEGER :: i LOGICAL :: is_modified !**************************************************************************************** IF (type_ocean == 'couple'.OR. & (type_ocean == 'slab' .AND. version_ocean == 'sicINT')) THEN ! limit.nc has not yet been read. Do it now! CALL limit_read_tot(itime, dtime, jour, is_modified) END IF DO i=1,knon rugos_out(i) = rugos(knindex(i)) alb_out(i) = albedo(knindex(i)) END DO END SUBROUTINE limit_read_rug_alb !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE limit_read_sst(knon, knindex, sst_out & #ifdef ISO ,Roce,rlat & #endif ) ! This SUBROUTINE returns the sea surface temperature already read from limit.nc. USE dimphy, ONLY: klon #ifdef ISO USE infotrac_phy, ONLY: niso USE isotopes_mod, ONLY: tcorr,toce,modif_sst, & deltaTtest,sstlatcrit,deltaTtestpoles,dsstlatcrit, & iso_HTO,ok_prod_nucl_tritium #ifdef ISOVERIF USE isotopes_verif_mod, ONLY: iso_verif_egalite_vect2D,iso_verif_positif, & iso_verif_positif_nostop #endif #endif INTEGER, INTENT(IN) :: knon ! nomber of points on compressed grid INTEGER, DIMENSION(klon), INTENT(IN) :: knindex ! grid point number for compressed grid REAL, DIMENSION(klon), INTENT(OUT) :: sst_out #ifdef ISO REAL, DIMENSION(klon) :: tuoce_out ! sortie tritium surface ocean #endif INTEGER :: i #ifdef ISO REAL, INTENT(OUT), DIMENSION(niso,klon) :: Roce INTEGER :: ixt REAL, INTENT(IN),DIMENSION(klon) :: rlat REAL lat_locale #endif !#ifdef ISOVERIF ! integer iso_verif_positif_nostop !#endif DO i = 1, knon sst_out(i) = sst(knindex(i)) END DO #ifdef ISO IF (iso_HTO.gt.0) THEN IF (ok_prod_nucl_tritium) then ! si on active la production nucleaire de tritium DO i = 1, knon tuoce_out(i)=tuoce(knindex(i)) END DO endif endif #endif #ifdef ISO IF (modif_sst.ge.1) THEN do i = 1, knon lat_locale=rlat(knindex(i)) ! test: modification uniforme de la sst IF (modif_sst.EQ.1) THEN sst_out(i)= sst_out(i)+deltaTtest elseif (modif_sst.EQ.2) then !if (modif_sst.EQ.1) THEN ! pattern parabolique en dehors des tropiques (sstlatcrit) IF (abs(lat_locale).gt.sstlatcrit) THEN sst_out(i)= sst_out(i)+deltaTtestpoles & *(lat_locale**2-sstlatcrit**2) & /(90.0**2-sstlatcrit**2) endif !if (abs(lat_locale).gt.abs(sstlatcrit)) THEN ELSE IF (modif_sst.EQ.3) THEN IF (abs(lat_locale).gt.abs(sstlatcrit)) THEN IF (abs(lat_locale).gt.sstlatcrit+dsstlatcrit) THEN sst_out(i)= sst_out(i)+deltaTtestpoles else sst_out(i)= sst_out(i)+deltaTtestpoles & *(abs(lat_locale)-sstlatcrit)/dsstlatcrit endif endif endif !if (modif_sst.EQ.1) THEN enddo !do i = 1, knon endif !if (modif_sst.ge.1) THEN #endif #ifdef ISOVERIF do i=1,knon,20 CALL iso_verif_positif(sst_out(i)-100.0,'limit_read 4323') enddo #endif #ifdef ISO !* lecture de Roce ! 1) première possibilité: valeur fixe à SMOW DO i = 1, knon do ixt=1,niso Roce(ixt,i)=tcorr(ixt)*toce(ixt) enddo !do ixt=1,niso enddo !DO i = 1, knon ! 2) deuxième possibilité: lecture de la carte ! A FAIRE ! lecture pour le tritium IF ((iso_HTO.gt.0).AND.(ok_prod_nucl_tritium)) THEN ! lecture de la carte tritium ocean surface Roce(iso_HTO,i)=tcorr(iso_HTO)*tuoce_out(i)*1.E-18*2. endif #endif #ifdef ISOVERIF do i=1,knon IF (iso_verif_positif_nostop(370.0-sst_out(i), & 'limit_read 368').EQ.1) THEN WRITE(*,*) 'i,knindex,sst_out=',i,knindex,sst_out(i) stop endif enddo #endif END SUBROUTINE limit_read_sst !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! !! Private SUBROUTINE : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE limit_read_tot(itime, dtime, jour, is_modified) ! Read everything needed from limit.nc ! 0) Initialize ! 1) Open the file limit.nc, if it is time ! 2) Read fraction, if not type_ocean=couple ! 3) Read sea surface temperature, if not type_ocean=couple ! 4) Read albedo and rugosity for land surface, ONLY in case of no vegetation model ! 5) Close file and distribuate variables to all processus USE dimphy USE lmdz_grid_phy USE lmdz_phys_para USE surface_data, ONLY: type_ocean, ok_veget USE netcdf, ONLY:nf90_noerr,nf90_close,nf90_get_var,nf90_inq_varid,nf90_nowrite,& nf90_inq_dimid,nf90_inquire_dimension,nf90_open,nf90_get_att,nf90_inquire USE indice_sol_mod #ifdef ISO USE isotopes_mod, ONLY: iso_HTO,ok_prod_nucl_tritium #ifdef ISOVERIF USE isotopes_verif_mod, ONLY: iso_verif_positif_nostop #endif #endif USE phys_cal_mod, ONLY: calend, year_len USE lmdz_print_control, ONLY: lunout, prt_level USE lmdz_XIOS, ONLY: xios_recv_field IMPLICIT NONE ! In- and ouput arguments !**************************************************************************************** INTEGER, INTENT(IN) :: itime ! numero du pas de temps courant INTEGER, INTENT(IN) :: jour ! jour a lire dans l'annee REAL , INTENT(IN) :: dtime ! pas de temps de la physique (en s) LOGICAL, INTENT(OUT) :: is_modified ! true if pctsrf is modified at this time step ! Locals variables with attribute SAVE !**************************************************************************************** ! frequence de lecture des conditions limites (en pas de physique) INTEGER,SAVE :: lmt_pas !$OMP THREADPRIVATE(lmt_pas) LOGICAL, SAVE :: first_call=.TRUE. !$OMP THREADPRIVATE(first_call) INTEGER, SAVE :: jour_lu = -1 !$OMP THREADPRIVATE(jour_lu) ! Locals variables !**************************************************************************************** INTEGER :: nid, nvarid, ndimid, nn INTEGER :: ii, ierr INTEGER, DIMENSION(2) :: start, epais REAL, DIMENSION(klon_glo,nbsrf) :: pct_glo ! fraction at global grid REAL, DIMENSION(klon_glo) :: sst_glo ! sea-surface temperature at global grid REAL, DIMENSION(klon_glo) :: rug_glo ! rugosity at global grid REAL, DIMENSION(klon_glo) :: alb_glo ! albedo at global grid REAL, DIMENSION(klon_mpi,nbsrf) :: pct_mpi ! fraction at global grid REAL, DIMENSION(klon_mpi) :: sst_mpi ! sea-surface temperature at global grid REAL, DIMENSION(klon_mpi) :: rug_mpi ! rugosity at global grid REAL, DIMENSION(klon_mpi) :: alb_mpi ! albedo at global grid CHARACTER(len=20) :: modname='limit_read_mod' CHARACTER(LEN=99) :: abort_message, calendar, str #ifdef ISO REAL, DIMENSION(klon_glo) :: tuoce_glo ! sea-surface tritium et global grid #endif ! End declaration !**************************************************************************************** !**************************************************************************************** ! 0) Initialization !**************************************************************************************** IF (first_call) THEN first_call=.FALSE. ! calculate number of time steps for one day lmt_pas = NINT(86400./dtime * 1.0) ! Allocate module save variables IF ( type_ocean /= 'couple' ) THEN ALLOCATE(pctsrf(klon,nbsrf), sst(klon), stat=ierr) IF (ierr /= 0) CALL abort_physic(modname, 'PB in allocating pctsrf and sst',1) END IF IF ( .NOT. ok_veget ) THEN ALLOCATE(rugos(klon), albedo(klon), stat=ierr) IF (ierr /= 0) CALL abort_physic(modname, 'PB in allocating rugos and albedo',1) END IF !$OMP MASTER ! Only master thread IF (is_mpi_root) THEN ! Only master processus ierr = nf90_open ('limit.nc', nf90_nowrite, nid) IF (ierr /= nf90_noerr) CALL abort_physic(modname,& 'Pb d''ouverture du fichier de conditions aux limites',1) !--- WARNING IF CALENDAR IS KNOWN AND DOES NOT MATCH THE ONE OF LMDZ ierr=nf90_inq_varid(nid, 'TEMPS', nvarid) ierr=nf90_get_att(nid, nvarid, 'calendar', calendar) IF(ierr==nf90_noerr.AND.calendar/=calend.AND.prt_level>=1) THEN WRITE(lunout,*)'BEWARE: gcm and limit.nc calendars differ: ' WRITE(lunout,*)' '//TRIM(calend)//' for gcm' WRITE(lunout,*)' '//TRIM(calendar)//' for limit.nc file' END IF !--- ERROR IF FILE RECORDS NUMBER IS NOT EQUAL TO EXPECTED NUMBER OF DAYS IF (grid_type==unstructured) THEN ierr=nf90_inq_dimid(nid,"time_year",ndimid) ELSE ierr=nf90_inquire(nid, UnlimitedDimID=ndimid) ENDIF ierr=nf90_inquire_dimension(nid, ndimid, len=nn) WRITE(abort_message,'(a,2(i0,a))')'limit.nc records number (',nn,') does no'//& 't match year length (',year_len,')' IF(nn/=year_len) CALL abort_physic(modname,abort_message,1) !--- ERROR IF FILES AND LMDZ HORIZONTAL RESOLUTIONS DO NOT MATCH IF (grid_type==unstructured) THEN ierr=nf90_inq_dimid(nid, 'cell', ndimid) ELSE ierr=nf90_inq_dimid(nid, 'points_physiques', ndimid) ENDIF ierr=nf90_inquire_dimension(nid, ndimid, len=nn) WRITE(abort_message,'(a,2(i0,a))')'limit.nc horizontal number of cells (',nn, & ') does not match LMDZ klon_glo (',klon_glo,')' IF(nn/=klon_glo) CALL abort_physic(modname,abort_message,1) ierr = nf90_close(nid) IF (ierr /= nf90_noerr) CALL abort_physic(modname,'Pb when closing file', 1) END IF ! is_mpi_root !$OMP END MASTER !$OMP BARRIER END IF !**************************************************************************************** ! 1) Open the file limit.nc if it is the right moment to read, once a day. ! The file is read only by the master thread of the master mpi process(is_mpi_root) ! Check by the way if the number of records is correct. !**************************************************************************************** is_modified = .FALSE. !ym IF (MOD(itime-1, lmt_pas) == 0 .OR. jour_lu /= jour ) THEN ! time to read ! not REALLY PERIODIC IF (MOD(itime-1, lmt_pas) == 0 .OR. (jour_lu /= jour .AND. grid_type /= unstructured)) THEN ! time to read ! IF (MOD(itime-1, lmt_pas) == 0) THEN ! time to read jour_lu = jour is_modified = .TRUE. IF (grid_type==unstructured) THEN IF ( type_ocean /= 'couple') THEN IF (is_omp_master) CALL xios_recv_field("foce_limin",pct_mpi(:,is_oce)) IF (is_omp_master) CALL xios_recv_field("fsic_limin",pct_mpi(:,is_sic)) ! IF (read_continents .OR. itime == 1) THEN IF (is_omp_master) CALL xios_recv_field("fter_limin",pct_mpi(:,is_ter)) IF (is_omp_master) CALL xios_recv_field("flic_limin",pct_mpi(:,is_lic)) ! ENDIF ENDIF! type_ocean /= couple IF ( type_ocean /= 'couple') THEN IF (is_omp_master) CALL xios_recv_field("sst_limin",sst_mpi) ENDIF IF (.NOT. ok_veget) THEN IF (is_omp_master) CALL xios_recv_field("alb_limin",alb_mpi) IF (is_omp_master) CALL xios_recv_field("rug_limin",rug_mpi) ENDIF IF ( type_ocean /= 'couple') THEN CALL Scatter_omp(sst_mpi,sst) CALL Scatter_omp(pct_mpi(:,is_oce),pctsrf(:,is_oce)) CALL Scatter_omp(pct_mpi(:,is_sic),pctsrf(:,is_sic)) ! IF (read_continents .OR. itime == 1) THEN CALL Scatter_omp(pct_mpi(:,is_ter),pctsrf(:,is_ter)) CALL Scatter_omp(pct_mpi(:,is_lic),pctsrf(:,is_lic)) ! END IF END IF IF (.NOT. ok_veget) THEN CALL Scatter_omp(alb_mpi, albedo) CALL Scatter_omp(rug_mpi, rugos) END IF ELSE ! grid_type==regular !$OMP MASTER ! Only master thread IF (is_mpi_root) THEN ! Only master processus! ierr = nf90_open ('limit.nc', nf90_nowrite, nid) IF (ierr /= nf90_noerr) CALL abort_physic(modname,& 'Pb d''ouverture du fichier de conditions aux limites',1) ! La tranche de donnees a lire: start(1) = 1 start(2) = jour epais(1) = klon_glo epais(2) = 1 !**************************************************************************************** ! 2) Read fraction if not type_ocean=couple !**************************************************************************************** IF ( type_ocean /= 'couple') THEN ! Ocean fraction ierr = nf90_inq_varid(nid, 'FOCE', nvarid) IF (ierr /= nf90_noerr) CALL abort_physic(modname, 'Le champ est absent',1) ierr = nf90_get_var(nid,nvarid,pct_glo(:,is_oce),start,epais) IF (ierr /= nf90_noerr) CALL abort_physic(modname,'Lecture echouee pour ' ,1) ! Sea-ice fraction ierr = nf90_inq_varid(nid, 'FSIC', nvarid) IF (ierr /= nf90_noerr) CALL abort_physic(modname,'Le champ est absent',1) ierr = nf90_get_var(nid,nvarid,pct_glo(:,is_sic),start,epais) IF (ierr /= nf90_noerr) CALL abort_physic(modname,'Lecture echouee pour ' ,1) ! Read land and continentals fraction only if asked for IF (read_continents .OR. itime == 1) THEN ! Land fraction ierr = nf90_inq_varid(nid, 'FTER', nvarid) IF (ierr /= nf90_noerr) CALL abort_physic(modname,'Le champ est absent',1) ierr = nf90_get_var(nid,nvarid,pct_glo(:,is_ter),start,epais) IF (ierr /= nf90_noerr) CALL abort_physic(modname,'Lecture echouee pour ',1) ! Continentale ice fraction ierr = nf90_inq_varid(nid, 'FLIC', nvarid) IF (ierr /= nf90_noerr) CALL abort_physic(modname,'Le champ est absent',1) ierr = nf90_get_var(nid,nvarid,pct_glo(:,is_lic),start,epais) IF (ierr /= nf90_noerr) CALL abort_physic(modname,'Lecture echouee pour ',1) END IF END IF ! type_ocean /= couple !**************************************************************************************** ! 3) Read sea-surface temperature, if not coupled ocean !**************************************************************************************** IF ( type_ocean /= 'couple') THEN ierr = nf90_inq_varid(nid, 'SST', nvarid) IF (ierr /= nf90_noerr) CALL abort_physic(modname,'Le champ est absent',1) ierr = nf90_get_var(nid,nvarid,sst_glo,start,epais) IF (ierr /= nf90_noerr) CALL abort_physic(modname,'Lecture echouee pour ',1) #ifdef ISO IF ((iso_HTO.gt.0).AND.(ok_prod_nucl_tritium)) THEN ierr = nf90_inq_varid(nid, 'TUOCE', nvarid) IF (ierr /= nf90_noerr) CALL abort_gcm(modname,'Le champ est absent',1) ierr = nf90_get_var(nid,nvarid,tuoce_glo,start,epais) IF (ierr /= nf90_noerr) CALL abort_gcm(modname,'Lecture echouee pour ',1) END IF #ifdef ISOVERIF do ii=1,klon_glo IF (iso_verif_positif_nostop(370.0-sst_glo(ii), & 'limit_read 384').EQ.1) THEN WRITE(*,*) 'ii,sst_glo=',ii,sst_glo(ii) WRITE(*,*) 'jour,start,epais=',jour,start,epais stop endif enddo #endif #endif END IF !**************************************************************************************** ! 4) Read albedo and rugosity for land surface, ONLY in case of no vegetation model !**************************************************************************************** IF (.NOT. ok_veget) THEN ! Read albedo ierr = nf90_inq_varid(nid, 'ALB', nvarid) IF (ierr /= nf90_noerr) CALL abort_physic(modname,'Le champ est absent',1) ierr = nf90_get_var(nid,nvarid,alb_glo,start,epais) IF (ierr /= nf90_noerr) CALL abort_physic(modname,'Lecture echouee pour ',1) ! Read rugosity ierr = nf90_inq_varid(nid, 'RUG', nvarid) IF (ierr /= nf90_noerr) CALL abort_physic(modname,'Le champ est absent',1) ierr = nf90_get_var(nid,nvarid,rug_glo,start,epais) IF (ierr /= nf90_noerr) CALL abort_physic(modname,'Lecture echouee pour ',1) END IF !**************************************************************************************** ! 5) Close file and distribuate variables to all processus !**************************************************************************************** ierr = nf90_close(nid) IF (ierr /= nf90_noerr) CALL abort_physic(modname,'Pb when closing file', 1) ENDIF ! is_mpi_root !$OMP END MASTER !$OMP BARRIER IF ( type_ocean /= 'couple') THEN CALL Scatter(sst_glo,sst) CALL Scatter(pct_glo(:,is_oce),pctsrf(:,is_oce)) CALL Scatter(pct_glo(:,is_sic),pctsrf(:,is_sic)) IF (read_continents .OR. itime == 1) THEN CALL Scatter(pct_glo(:,is_ter),pctsrf(:,is_ter)) CALL Scatter(pct_glo(:,is_lic),pctsrf(:,is_lic)) END IF #ifdef ISO IF ((iso_HTO.gt.0).AND.(ok_prod_nucl_tritium)) THEN CALL Scatter(tuoce_glo,tuoce) END IF #endif END IF IF (.NOT. ok_veget) THEN CALL Scatter(alb_glo, albedo) CALL Scatter(rug_glo, rugos) END IF ENDIF ! Grid type ENDIF ! time to read END SUBROUTINE limit_read_tot END MODULE limit_read_mod