Context Navigation

limit_read_mod.F90 @ 5209

Last change on this file since 5209 was 5084, checked in by Laurent Fairhead, 4 months ago
Reverting to r4065. Updating fortran standard broke too much stuff. Will do it by smaller chunks AB, LF
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 17.8 KB

Rev	Line
[996]	1	!
[3435]	2	! $Id: limit_read_mod.F90 5084 2024-07-19 16:40:44Z fairhead $
[996]	3	!
	4	MODULE limit_read_mod
	5	!
	6	! This module reads the fichier "limit.nc" containing fields for surface forcing.
	7	!
	8	! Module subroutines :
	9	! limit_read_frac : call limit_read_tot and return the fractions
	10	! limit_read_rug_alb : return rugosity and albedo, if coupled ocean call limit_read_tot first
	11	! limit_read_sst : return sea ice temperature
	12	! limit_read_tot : read limit.nc and store the fields in local modules variables
	13	!
	14	IMPLICIT NONE
	15
	16	REAL, ALLOCATABLE, DIMENSION(:,:), SAVE, PRIVATE :: pctsrf
[1001]	17	!$OMP THREADPRIVATE(pctsrf)
[996]	18	REAL, ALLOCATABLE, DIMENSION(:), SAVE, PRIVATE :: rugos
[1001]	19	!$OMP THREADPRIVATE(rugos)
[996]	20	REAL, ALLOCATABLE, DIMENSION(:), SAVE, PRIVATE :: albedo
[1001]	21	!$OMP THREADPRIVATE(albedo)
[996]	22	REAL, ALLOCATABLE, DIMENSION(:), SAVE, PRIVATE :: sst
[1001]	23	!$OMP THREADPRIVATE(sst)
	24	LOGICAL,SAVE :: read_continents=.FALSE.
	25	!$OMP THREADPRIVATE(read_continents)
[996]	26
	27	CONTAINS
	28	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	29	!!
	30	!! Public subroutines :
	31	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	32
[3435]	33
	34	SUBROUTINE init_limit_read(first_day)
	35	USE mod_grid_phy_lmdz
	36	USE surface_data
	37	USE mod_phys_lmdz_para
[4619]	38	USE lmdz_xios
	39
[3435]	40	IMPLICIT NONE
	41	INTEGER, INTENT(IN) :: first_day
	42
	43
	44	IF ( type_ocean /= 'couple') THEN
	45	IF (grid_type==unstructured) THEN
	46	IF (is_omp_master) CALL xios_set_file_attr("limit_read",enabled=.TRUE.,record_offset=first_day)
	47	ENDIF
	48	ENDIF
	49
	50	END SUBROUTINE init_limit_read
	51
[996]	52	SUBROUTINE limit_read_frac(itime, dtime, jour, pctsrf_new, is_modified)
	53	!
	54	! This subroutine is called from "change_srf_frac" for case of
	55	! ocean=force or from ocean_slab_frac for ocean=slab.
	56	! The fraction for all sub-surfaces at actual time step is returned.
	57
	58	USE dimphy
[1785]	59	USE indice_sol_mod
[996]	60
	61	! Input arguments
	62	!****************************************************************************************
	63	INTEGER, INTENT(IN) :: itime ! time step
	64	INTEGER, INTENT(IN) :: jour ! current day
	65	REAL , INTENT(IN) :: dtime ! length of time step
	66
	67	! Output arguments
	68	!****************************************************************************************
	69	REAL, DIMENSION(klon,nbsrf), INTENT(OUT) :: pctsrf_new ! sub surface fractions
	70	LOGICAL, INTENT(OUT) :: is_modified ! true if pctsrf is modified at this time step
	71
	72	! End declaration
	73	!****************************************************************************************
	74
	75	! 1) Read file limit.nc
	76	CALL limit_read_tot(itime, dtime, jour, is_modified)
	77
	78	! 2) Return the fraction read in limit_read_tot
	79	pctsrf_new(:,:) = pctsrf(:,:)
	80
	81	END SUBROUTINE limit_read_frac
	82
	83	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	84
	85	SUBROUTINE limit_read_rug_alb(itime, dtime, jour, &
	86	knon, knindex, &
	87	rugos_out, alb_out)
	88	!
	89	! This subroutine is called from surf_land_bucket.
	90	! The flag "ok_veget" must can not be true. If coupled run, "ocean=couple"
	91	! then this routine will call limit_read_tot.
	92	!
	93	USE dimphy
	94	USE surface_data
	95
	96	! Input arguments
	97	!****************************************************************************************
	98	INTEGER, INTENT(IN) :: itime ! numero du pas de temps courant
	99	INTEGER, INTENT(IN) :: jour ! jour a lire dans l'annee
	100	REAL , INTENT(IN) :: dtime ! pas de temps de la physique (en s)
	101	INTEGER, INTENT(IN) :: knon ! nomber of points on compressed grid
	102	INTEGER, DIMENSION(klon), INTENT(IN) :: knindex ! grid point number for compressed grid
	103	! Output arguments
	104	!****************************************************************************************
	105	REAL, DIMENSION(klon), INTENT(OUT) :: rugos_out
	106	REAL, DIMENSION(klon), INTENT(OUT) :: alb_out
	107
	108	! Local variables
	109	!****************************************************************************************
	110	INTEGER :: i
	111	LOGICAL :: is_modified
	112	!****************************************************************************************
	113
[2209]	114	IF (type_ocean == 'couple'.OR. &
	115	(type_ocean == 'slab' .AND. version_ocean == 'sicINT')) THEN
[996]	116	! limit.nc has not yet been read. Do it now!
	117	CALL limit_read_tot(itime, dtime, jour, is_modified)
	118	END IF
	119
	120	DO i=1,knon
	121	rugos_out(i) = rugos(knindex(i))
	122	alb_out(i) = albedo(knindex(i))
	123	END DO
	124
	125	END SUBROUTINE limit_read_rug_alb
	126
	127	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	128
	129	SUBROUTINE limit_read_sst(knon, knindex, sst_out)
	130	!
	131	! This subroutine returns the sea surface temperature already read from limit.nc.
	132	!
	133	USE dimphy, ONLY : klon
	134
	135	INTEGER, INTENT(IN) :: knon ! nomber of points on compressed grid
	136	INTEGER, DIMENSION(klon), INTENT(IN) :: knindex ! grid point number for compressed grid
	137	REAL, DIMENSION(klon), INTENT(OUT) :: sst_out
	138
	139	INTEGER :: i
	140
	141	DO i = 1, knon
	142	sst_out(i) = sst(knindex(i))
	143	END DO
	144
	145	END SUBROUTINE limit_read_sst
	146
	147	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	148	!!
	149	!! Private subroutine :
	150	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	151
	152	SUBROUTINE limit_read_tot(itime, dtime, jour, is_modified)
	153	!
	154	! Read everything needed from limit.nc
	155	!
	156	! 0) Initialize
	157	! 1) Open the file limit.nc, if it is time
	158	! 2) Read fraction, if not type_ocean=couple
	159	! 3) Read sea surface temperature, if not type_ocean=couple
	160	! 4) Read albedo and rugosity for land surface, only in case of no vegetation model
	161	! 5) Close file and distribuate variables to all processus
	162
	163	USE dimphy
	164	USE mod_grid_phy_lmdz
	165	USE mod_phys_lmdz_para
	166	USE surface_data, ONLY : type_ocean, ok_veget
[5084]	167	USE netcdf
[1785]	168	USE indice_sol_mod
[2665]	169	USE phys_cal_mod, ONLY : calend, year_len
	170	USE print_control_mod, ONLY: lunout, prt_level
[4619]	171	USE lmdz_xios, ONLY: xios_recv_field, using_xios
[3435]	172
[996]	173	IMPLICIT NONE
	174
	175	! In- and ouput arguments
	176	!****************************************************************************************
	177	INTEGER, INTENT(IN) :: itime ! numero du pas de temps courant
	178	INTEGER, INTENT(IN) :: jour ! jour a lire dans l'annee
	179	REAL , INTENT(IN) :: dtime ! pas de temps de la physique (en s)
	180
	181	LOGICAL, INTENT(OUT) :: is_modified ! true if pctsrf is modified at this time step
	182
	183	! Locals variables with attribute SAVE
	184	!****************************************************************************************
	185	! frequence de lecture des conditions limites (en pas de physique)
[1001]	186	INTEGER,SAVE :: lmt_pas
	187	!$OMP THREADPRIVATE(lmt_pas)
[996]	188	LOGICAL, SAVE :: first_call=.TRUE.
[1577]	189	!$OMP THREADPRIVATE(first_call)
	190	INTEGER, SAVE :: jour_lu = -1
	191	!$OMP THREADPRIVATE(jour_lu)
[996]	192	! Locals variables
	193	!****************************************************************************************
[2665]	194	INTEGER :: nid, nvarid, ndimid, nn
[996]	195	INTEGER :: ii, ierr
	196	INTEGER, DIMENSION(2) :: start, epais
	197	REAL, DIMENSION(klon_glo,nbsrf) :: pct_glo ! fraction at global grid
	198	REAL, DIMENSION(klon_glo) :: sst_glo ! sea-surface temperature at global grid
	199	REAL, DIMENSION(klon_glo) :: rug_glo ! rugosity at global grid
	200	REAL, DIMENSION(klon_glo) :: alb_glo ! albedo at global grid
[3435]	201
	202	REAL, DIMENSION(klon_mpi,nbsrf) :: pct_mpi ! fraction at global grid
	203	REAL, DIMENSION(klon_mpi) :: sst_mpi ! sea-surface temperature at global grid
	204	REAL, DIMENSION(klon_mpi) :: rug_mpi ! rugosity at global grid
	205	REAL, DIMENSION(klon_mpi) :: alb_mpi ! albedo at global grid
	206
[996]	207	CHARACTER(len=20) :: modname='limit_read_mod'
[2665]	208	CHARACTER(LEN=99) :: abort_message, calendar, str
[996]	209
	210	! End declaration
	211	!****************************************************************************************
	212
	213	!****************************************************************************************
	214	! 0) Initialization
	215	!
	216	!****************************************************************************************
	217	IF (first_call) THEN
[2770]	218	first_call=.FALSE.
[996]	219	! calculate number of time steps for one day
	220	lmt_pas = NINT(86400./dtime * 1.0)
	221
	222	! Allocate module save variables
	223	IF ( type_ocean /= 'couple' ) THEN
	224	ALLOCATE(pctsrf(klon,nbsrf), sst(klon), stat=ierr)
[2311]	225	IF (ierr /= 0) CALL abort_physic(modname, 'PB in allocating pctsrf and sst',1)
[996]	226	END IF
	227
	228	IF ( .NOT. ok_veget ) THEN
	229	ALLOCATE(rugos(klon), albedo(klon), stat=ierr)
[2311]	230	IF (ierr /= 0) CALL abort_physic(modname, 'PB in allocating rugos and albedo',1)
[996]	231	END IF
	232
	233	!$OMP MASTER ! Only master thread
	234	IF (is_mpi_root) THEN ! Only master processus
	235	ierr = NF90_OPEN ('limit.nc', NF90_NOWRITE, nid)
[2311]	236	IF (ierr /= NF90_NOERR) CALL abort_physic(modname,&
[996]	237	'Pb d''ouverture du fichier de conditions aux limites',1)
[2665]	238
	239	!--- WARNING IF CALENDAR IS KNOWN AND DOES NOT MATCH THE ONE OF LMDZ
	240	ierr=NF90_INQ_VARID(nid, 'TEMPS', nvarid)
	241	ierr=NF90_GET_ATT(nid, nvarid, 'calendar', calendar)
	242	IF(ierr==NF90_NOERR.AND.calendar/=calend.AND.prt_level>=1) THEN
	243	WRITE(lunout,*)'BEWARE: gcm and limit.nc calendars differ: '
	244	WRITE(lunout,*)' '//TRIM(calend)//' for gcm'
	245	WRITE(lunout,*)' '//TRIM(calendar)//' for limit.nc file'
	246	END IF
	247
[3435]	248	!--- ERROR IF FILE RECORDS NUMBER IS NOT EQUAL TO EXPECTED NUMBER OF DAYS
	249	IF (grid_type==unstructured) THEN
	250	ierr=NF90_INQ_DIMID(nid,"time_year",ndimid)
	251	ELSE
	252	ierr=NF90_INQUIRE(nid, UnlimitedDimID=ndimid)
	253	ENDIF
[2665]	254	ierr=NF90_INQUIRE_DIMENSION(nid, ndimid, len=nn)
[3435]	255	WRITE(abort_message,'(a,2(i0,a))')'limit.nc records number (',nn,') does no'//&
[2768]	256	't match year length (',year_len,')'
[2665]	257	IF(nn/=year_len) CALL abort_physic(modname,abort_message,1)
	258
	259	!--- ERROR IF FILES AND LMDZ HORIZONTAL RESOLUTIONS DO NOT MATCH
[3435]	260	IF (grid_type==unstructured) THEN
	261	ierr=NF90_INQ_DIMID(nid, 'cell', ndimid)
	262	ELSE
	263	ierr=NF90_INQ_DIMID(nid, 'points_physiques', ndimid)
	264	ENDIF
[2665]	265	ierr=NF90_INQUIRE_DIMENSION(nid, ndimid, len=nn)
[2768]	266	WRITE(abort_message,'(a,2(i0,a))')'limit.nc horizontal number of cells (',nn, &
	267	') does not match LMDZ klon_glo (',klon_glo,')'
[2665]	268	IF(nn/=klon_glo) CALL abort_physic(modname,abort_message,1)
	269
[2768]	270	ierr = NF90_CLOSE(nid)
	271	IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Pb when closing file', 1)
	272	END IF ! is_mpi_root
	273	!$OMP END MASTER
	274	!$OMP BARRIER
	275	END IF
	276
	277	!****************************************************************************************
	278	! 1) Open the file limit.nc if it is the right moment to read, once a day.
	279	! The file is read only by the master thread of the master mpi process(is_mpi_root)
	280	! Check by the way if the number of records is correct.
	281	!
	282	!****************************************************************************************
	283
	284	is_modified = .FALSE.
[3435]	285	!ym IF (MOD(itime-1, lmt_pas) == 0 .OR. jour_lu /= jour ) THEN ! time to read
	286	! not REALLY PERIODIC
	287	IF (MOD(itime-1, lmt_pas) == 0 .OR. (jour_lu /= jour .AND. grid_type /= unstructured)) THEN ! time to read
	288	! IF (MOD(itime-1, lmt_pas) == 0) THEN ! time to read
[2768]	289	jour_lu = jour
	290	is_modified = .TRUE.
[3435]	291
	292	IF (grid_type==unstructured) THEN
	293
	294	IF ( type_ocean /= 'couple') THEN
	295
	296	IF (is_omp_master) CALL xios_recv_field("foce_limin",pct_mpi(:,is_oce))
	297	IF (is_omp_master) CALL xios_recv_field("fsic_limin",pct_mpi(:,is_sic))
	298	! IF (read_continents .OR. itime == 1) THEN
	299	IF (is_omp_master) CALL xios_recv_field("fter_limin",pct_mpi(:,is_ter))
	300	IF (is_omp_master) CALL xios_recv_field("flic_limin",pct_mpi(:,is_lic))
	301	! ENDIF
	302	ENDIF! type_ocean /= couple
	303
	304	IF ( type_ocean /= 'couple') THEN
	305	IF (is_omp_master) CALL xios_recv_field("sst_limin",sst_mpi)
	306	ENDIF
	307
	308	IF (.NOT. ok_veget) THEN
	309	IF (is_omp_master) CALL xios_recv_field("alb_limin",alb_mpi)
	310	IF (is_omp_master) CALL xios_recv_field("rug_limin",rug_mpi)
	311	ENDIF
	312
	313	IF ( type_ocean /= 'couple') THEN
	314	CALL Scatter_omp(sst_mpi,sst)
	315	CALL Scatter_omp(pct_mpi(:,is_oce),pctsrf(:,is_oce))
	316	CALL Scatter_omp(pct_mpi(:,is_sic),pctsrf(:,is_sic))
	317	! IF (read_continents .OR. itime == 1) THEN
	318	CALL Scatter_omp(pct_mpi(:,is_ter),pctsrf(:,is_ter))
	319	CALL Scatter_omp(pct_mpi(:,is_lic),pctsrf(:,is_lic))
	320	! END IF
	321	END IF
	322
	323	IF (.NOT. ok_veget) THEN
	324	CALL Scatter_omp(alb_mpi, albedo)
	325	CALL Scatter_omp(rug_mpi, rugos)
	326	END IF
	327
	328	ELSE ! grid_type==regular
	329
[2768]	330	!$OMP MASTER ! Only master thread
[3435]	331	IF (is_mpi_root) THEN ! Only master processus!
[2768]	332
	333	ierr = NF90_OPEN ('limit.nc', NF90_NOWRITE, nid)
	334	IF (ierr /= NF90_NOERR) CALL abort_physic(modname,&
	335	'Pb d''ouverture du fichier de conditions aux limites',1)
	336
[996]	337	! La tranche de donnees a lire:
	338	start(1) = 1
	339	start(2) = jour
	340	epais(1) = klon_glo
	341	epais(2) = 1
	342
	343
	344	!****************************************************************************************
	345	! 2) Read fraction if not type_ocean=couple
	346	!
	347	!****************************************************************************************
	348
	349	IF ( type_ocean /= 'couple') THEN
	350	!
	351	! Ocean fraction
	352	ierr = NF90_INQ_VARID(nid, 'FOCE', nvarid)
[2311]	353	IF (ierr /= NF90_NOERR) CALL abort_physic(modname, 'Le champ <FOCE> est absent',1)
[996]	354
	355	ierr = NF90_GET_VAR(nid,nvarid,pct_glo(:,is_oce),start,epais)
[2311]	356	IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Lecture echouee pour <FOCE>' ,1)
[996]	357	!
	358	! Sea-ice fraction
	359	ierr = NF90_INQ_VARID(nid, 'FSIC', nvarid)
[2311]	360	IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Le champ <FSIC> est absent',1)
[996]	361
	362	ierr = NF90_GET_VAR(nid,nvarid,pct_glo(:,is_sic),start,epais)
[2311]	363	IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Lecture echouee pour <FSIC>' ,1)
[996]	364
	365
	366	! Read land and continentals fraction only if asked for
	367	IF (read_continents .OR. itime == 1) THEN
	368	!
	369	! Land fraction
	370	ierr = NF90_INQ_VARID(nid, 'FTER', nvarid)
[2311]	371	IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Le champ <FTER> est absent',1)
[996]	372
	373	ierr = NF90_GET_VAR(nid,nvarid,pct_glo(:,is_ter),start,epais)
[2311]	374	IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Lecture echouee pour <FTER>',1)
[996]	375	!
	376	! Continentale ice fraction
	377	ierr = NF90_INQ_VARID(nid, 'FLIC', nvarid)
[2311]	378	IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Le champ <FLIC> est absent',1)
[996]	379
	380	ierr = NF90_GET_VAR(nid,nvarid,pct_glo(:,is_lic),start,epais)
[2311]	381	IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Lecture echouee pour <FLIC>',1)
[996]	382	END IF
	383
	384	END IF ! type_ocean /= couple
	385
	386	!****************************************************************************************
	387	! 3) Read sea-surface temperature, if not coupled ocean
	388	!
	389	!****************************************************************************************
	390	IF ( type_ocean /= 'couple') THEN
	391
	392	ierr = NF90_INQ_VARID(nid, 'SST', nvarid)
[2311]	393	IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Le champ <SST> est absent',1)
[996]	394
	395	ierr = NF90_GET_VAR(nid,nvarid,sst_glo,start,epais)
[2311]	396	IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Lecture echouee pour <SST>',1)
[996]	397
	398	END IF
	399
	400	!****************************************************************************************
	401	! 4) Read albedo and rugosity for land surface, only in case of no vegetation model
	402	!
	403	!****************************************************************************************
	404
	405	IF (.NOT. ok_veget) THEN
	406	!
	407	! Read albedo
	408	ierr = NF90_INQ_VARID(nid, 'ALB', nvarid)
[2311]	409	IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Le champ <ALB> est absent',1)
[996]	410
	411	ierr = NF90_GET_VAR(nid,nvarid,alb_glo,start,epais)
[2311]	412	IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Lecture echouee pour <ALB>',1)
[996]	413	!
	414	! Read rugosity
	415	ierr = NF90_INQ_VARID(nid, 'RUG', nvarid)
[2311]	416	IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Le champ <RUG> est absent',1)
[996]	417
	418	ierr = NF90_GET_VAR(nid,nvarid,rug_glo,start,epais)
[2311]	419	IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Lecture echouee pour <RUG>',1)
[996]	420
	421	END IF
	422
	423	!****************************************************************************************
	424	! 5) Close file and distribuate variables to all processus
	425	!
	426	!****************************************************************************************
	427	ierr = NF90_CLOSE(nid)
[2311]	428	IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Pb when closing file', 1)
[996]	429	ENDIF ! is_mpi_root
	430
	431	!$OMP END MASTER
	432	!$OMP BARRIER
	433
	434	IF ( type_ocean /= 'couple') THEN
	435	CALL Scatter(sst_glo,sst)
	436	CALL Scatter(pct_glo(:,is_oce),pctsrf(:,is_oce))
	437	CALL Scatter(pct_glo(:,is_sic),pctsrf(:,is_sic))
	438	IF (read_continents .OR. itime == 1) THEN
	439	CALL Scatter(pct_glo(:,is_ter),pctsrf(:,is_ter))
	440	CALL Scatter(pct_glo(:,is_lic),pctsrf(:,is_lic))
	441	END IF
	442	END IF
	443
	444	IF (.NOT. ok_veget) THEN
	445	CALL Scatter(alb_glo, albedo)
	446	CALL Scatter(rug_glo, rugos)
	447	END IF
	448
[3435]	449	ENDIF ! Grid type
	450
[996]	451	ENDIF ! time to read
	452
	453	END SUBROUTINE limit_read_tot
	454
	455	END MODULE limit_read_mod

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format