Context Navigation

dynetat0.F @ 1713

Last change on this file since 1713 was 1422, checked in by milmd, 10 years ago
In GENERIC, MARS and COMMON models replace some include files by modules (usefull for decoupling physics with dynamics).
File size: 15.3 KB

Rev	Line
[1416]	1	SUBROUTINE dynetat0(fichnom,vcov,ucov,
[999]	2	. teta,q,masse,ps,phis,time0)
[38]	3
	4	use netcdf
[1416]	5	use infotrac, only: tname,nqtot
[1130]	6	use control_mod, only: timestart
[1422]	7	USE comvert_mod, ONLY: pa,preff
	8	USE comconst_mod, ONLY: im,jm,lllm,daysec,dtvr,
	9	. rad,omeg,r,cpp,kappa,g
	10	USE logic_mod, ONLY: fxyhypb,ysinus
	11	USE serre_mod, ONLY: clon,clat,grossismx,grossismy,dzoomx,dzoomy,
	12	. taux,tauy
	13	USE temps_mod, ONLY: day_ini,hour_ini
	14	USE ener_mod, ONLY: etot0,ptot0,ztot0,stot0,ang0
[38]	15
	16	IMPLICIT NONE
	17
	18	c=======================================================================
	19	c
	20	c Auteur: P. Le Van / L.Fairhead
	21	c -------
	22	c
	23	c objet:
	24	c ------
	25	c
	26	c Lecture de l'etat initial
	27	c
	28	c Modifs: Oct.2008 read in tracers by name. Ehouarn Millour
	29	c Aug.2010 use NetCDF90 to load variables (enables using
	30	c r4 or r8 restarts independently of having compiled
	31	c the GCM in r4 or r8)
[999]	32	c June 2013 TN : Possibility to read files with a time axis.
	33	c NOW defrun_new MUST BE CALLED BEFORE dynetat0.
[38]	34	c=======================================================================
	35	c-----------------------------------------------------------------------
	36	c Declarations:
	37	c -------------
	38
	39	#include "dimensions.h"
	40	#include "paramet.h"
	41	#include "comgeom2.h"
[1036]	42	!#include "advtrac.h"
[1130]	43	!#include "control.h"
[38]	44
	45	c Arguments:
	46	c ----------
	47
	48	CHARACTER() fichnom
	49	REAL vcov(iip1,jjm,llm),ucov(iip1,jjp1,llm),teta(iip1,jjp1,llm)
[1416]	50	REAL q(iip1,jjp1,llm,nqtot)
[38]	51	REAL masse(iip1,jjp1,llm)
	52	REAL ps(iip1,jjp1),phis(iip1,jjp1)
	53
[999]	54	REAL time0 ! fraction of day the fields correspond to
[38]	55
	56	c Variables
	57	c
	58	REAL traceur(iip1,jjp1,llm) ! to temporarily store a tracer
	59	INTEGER length,iq,i,j,l
	60	PARAMETER (length = 100)
	61	REAL tab_cntrl(length) ! tableau des parametres du run
	62	INTEGER ierr, nid, nvarid, nqold
	63	CHARACTER str33,yes1
[999]	64
	65	REAL,ALLOCATABLE :: time(:) ! times stored in start
	66	INTEGER timelen ! number of times stored in the file
	67	INTEGER indextime ! index of selected time
	68	!REAL hour_ini ! fraction of day of stored date. Equivalent of day_ini, but 0=<hour_ini<1
[38]	69
[999]	70	INTEGER edges(4),corner(4)
	71
[38]	72	c-----------------------------------------------------------------------
	73
	74	c Ouverture NetCDF du fichier etat initial
	75
	76	ierr=nf90_open(fichnom,NF90_NOWRITE,nid)
	77	IF (ierr.NE.nf90_noerr) THEN
	78	write(6,*)' Pb d''ouverture du fichier ',trim(fichnom)
	79	write(,)trim(nf90_strerror(ierr))
	80	CALL ABORT
	81	ENDIF
	82
	83	c
	84	ierr=nf90_inq_varid(nid,"controle",nvarid)
	85	IF (ierr .NE. nf90_noerr) THEN
	86	PRINT*, "dynetat0: Le champ <controle> est absent"
	87	write(,)trim(nf90_strerror(ierr))
	88	CALL abort
	89	ENDIF
	90	ierr=nf90_get_var(nid,nvarid,tab_cntrl)
	91	IF (ierr .NE. nf90_noerr) THEN
	92	PRINT*, "dynetat0: Lecture echoue pour <controle>"
	93	write(,)trim(nf90_strerror(ierr))
	94	CALL abort
	95	ENDIF
	96
	97	im = tab_cntrl(1)
	98	jm = tab_cntrl(2)
	99	lllm = tab_cntrl(3)
	100	day_ini = tab_cntrl(4)
	101	rad = tab_cntrl(5)
	102	omeg = tab_cntrl(6)
	103	g = tab_cntrl(7)
	104	cpp = tab_cntrl(8)
	105	kappa = tab_cntrl(9)
	106	daysec = tab_cntrl(10)
	107	dtvr = tab_cntrl(11)
	108	etot0 = tab_cntrl(12)
	109	ptot0 = tab_cntrl(13)
	110	ztot0 = tab_cntrl(14)
	111	stot0 = tab_cntrl(15)
	112	ang0 = tab_cntrl(16)
	113	pa = tab_cntrl(17)
	114	preff = tab_cntrl(18)
[999]	115
	116	hour_ini = tab_cntrl(29)
	117
[38]	118	c
	119	clon = tab_cntrl(19)
	120	clat = tab_cntrl(20)
	121	grossismx = tab_cntrl(21)
	122	grossismy = tab_cntrl(22)
	123	c
	124	IF ( tab_cntrl(23).EQ.1. ) THEN
	125	fxyhypb = . TRUE .
	126	dzoomx = tab_cntrl(24)
	127	dzoomy = tab_cntrl(25)
	128	taux = tab_cntrl(27)
	129	tauy = tab_cntrl(28)
	130	ELSE
	131	fxyhypb = . FALSE .
	132	ysinus = . FALSE .
	133	IF( tab_cntrl(26).EQ.1. ) ysinus = . TRUE.
	134	ENDIF
[999]	135
[38]	136	c .................................................................
	137	c
	138	c
	139	PRINT*,'dynetat0: rad,omeg,g,cpp,kappa',rad,omeg,g,cpp,kappa
	140
	141	IF( im.ne.iim ) THEN
	142	PRINT 1,im,iim
	143	STOP
	144	ELSE IF( jm.ne.jjm ) THEN
	145	PRINT 2,jm,jjm
	146	STOP
	147	ELSE IF( lllm.ne.llm ) THEN
	148	PRINT 3,lllm,llm
	149	STOP
	150	ENDIF
	151
	152	ierr=nf90_inq_varid(nid,"rlonu",nvarid)
	153	IF (ierr .NE. nf90_noerr) THEN
	154	PRINT*, "dynetat0: Le champ <rlonu> est absent"
	155	write(,)trim(nf90_strerror(ierr))
	156	CALL abort
	157	ENDIF
	158	ierr=nf90_get_var(nid,nvarid,rlonu)
	159	IF (ierr .NE. nf90_noerr) THEN
	160	PRINT*, "dynetat0: Lecture echouee pour <rlonu>"
	161	write(,)trim(nf90_strerror(ierr))
	162	CALL abort
	163	ENDIF
	164
	165	ierr=nf90_inq_varid(nid,"rlatu",nvarid)
	166	IF (ierr .NE. nf90_noerr) THEN
	167	PRINT*, "dynetat0: Le champ <rlatu> est absent"
	168	write(,)trim(nf90_strerror(ierr))
	169	CALL abort
	170	ENDIF
	171	ierr=nf90_get_var(nid,nvarid,rlatu)
	172	IF (ierr .NE. nf90_noerr) THEN
	173	PRINT*, "dynetat0: Lecture echouee pour <rlatu>"
	174	write(,)trim(nf90_strerror(ierr))
	175	CALL abort
	176	ENDIF
	177
	178	ierr=nf90_inq_varid(nid,"rlonv",nvarid)
	179	IF (ierr .NE. nf90_noerr) THEN
	180	PRINT*, "dynetat0: Le champ <rlonv> est absent"
	181	write(,)trim(nf90_strerror(ierr))
	182	CALL abort
	183	ENDIF
	184	ierr=nf90_get_var(nid,nvarid,rlonv)
	185	IF (ierr .NE. nf90_noerr) THEN
	186	PRINT*, "dynetat0: Lecture echouee pour <rlonv>"
	187	write(,)trim(nf90_strerror(ierr))
	188	CALL abort
	189	ENDIF
	190
	191	ierr=nf90_inq_varid(nid,"rlatv",nvarid)
	192	IF (ierr .NE. nf90_noerr) THEN
	193	PRINT*, "dynetat0: Le champ <rlatv> est absent"
	194	write(,)trim(nf90_strerror(ierr))
	195	CALL abort
	196	ENDIF
	197	ierr=nf90_get_var(nid,nvarid,rlatv)
	198	IF (ierr .NE. nf90_noerr) THEN
	199	PRINT*, "dynetat0: Lecture echouee pour rlatv"
	200	write(,)trim(nf90_strerror(ierr))
	201	CALL abort
	202	ENDIF
	203
	204	ierr=nf90_inq_varid(nid,"cu",nvarid)
	205	IF (ierr .NE. nf90_noerr) THEN
	206	PRINT*, "dynetat0: Le champ <cu> est absent"
	207	write(,)trim(nf90_strerror(ierr))
	208	CALL abort
	209	ENDIF
	210	ierr=nf90_get_var(nid,nvarid,cu)
	211	IF (ierr .NE. nf90_noerr) THEN
	212	PRINT*, "dynetat0: Lecture echouee pour <cu>"
	213	write(,)trim(nf90_strerror(ierr))
	214	CALL abort
	215	ENDIF
	216
	217	ierr=nf90_inq_varid(nid,"cv",nvarid)
	218	IF (ierr .NE. nf90_noerr) THEN
	219	PRINT*, "dynetat0: Le champ <cv> est absent"
	220	write(,)trim(nf90_strerror(ierr))
	221	CALL abort
	222	ENDIF
	223	ierr=nf90_get_var(nid,nvarid,cv)
	224	IF (ierr .NE. nf90_noerr) THEN
	225	PRINT*, "dynetat0: Lecture echouee pour <cv>"
	226	write(,)trim(nf90_strerror(ierr))
	227	CALL abort
	228	ENDIF
	229
	230	ierr=nf90_inq_varid(nid,"aire",nvarid)
	231	IF (ierr .NE. nf90_noerr) THEN
	232	PRINT*, "dynetat0: Le champ <aire> est absent"
	233	write(,)trim(nf90_strerror(ierr))
	234	CALL abort
	235	ENDIF
	236	ierr=nf90_get_var(nid,nvarid,aire)
	237	IF (ierr .NE. nf90_noerr) THEN
	238	PRINT*, "dynetat0: Lecture echouee pour <aire>"
	239	write(,)trim(nf90_strerror(ierr))
	240	CALL abort
	241	ENDIF
	242
	243	ierr=nf90_inq_varid(nid,"phisinit",nvarid)
	244	IF (ierr .NE. nf90_noerr) THEN
	245	PRINT*, "dynetat0: Le champ <phisinit> est absent"
	246	write(,)trim(nf90_strerror(ierr))
	247	CALL abort
	248	ENDIF
	249	ierr=nf90_get_var(nid,nvarid,phis)
	250	IF (ierr .NE. nf90_noerr) THEN
	251	PRINT*, "dynetat0: Lecture echouee pour <phisinit>"
	252	write(,)trim(nf90_strerror(ierr))
	253	CALL abort
	254	ENDIF
	255
[999]	256
	257	! read time axis
	258	ierr = nf90_inq_dimid(nid,"Time",nvarid)
	259	ierr = nf90_inquire_dimension(nid,nvarid,len=timelen)
	260
	261	ALLOCATE(time(timelen))
	262
[38]	263	ierr=nf90_inq_varid(nid,"Time",nvarid)
	264	IF (ierr .NE. nf90_noerr) THEN
	265	ierr=nf90_inq_varid(nid,"temps", nvarid)
	266	IF (ierr .NE. nf90_noerr) THEN
	267	PRINT*, "dynetat0: <Time> or <temps> absent"
	268	write(,)trim(nf90_strerror(ierr))
	269	CALL abort
	270	ENDIF
	271	ENDIF
	272	ierr=nf90_get_var(nid,nvarid,time)
	273	IF (ierr .NE. nf90_noerr) THEN
	274	PRINT*, "dynetat0: Lecture echouee <Time>/<temps>"
	275	write(,)trim(nf90_strerror(ierr))
	276	CALL abort
	277	ENDIF
[999]	278
	279	! select desired time
	280	IF (timestart .lt. 0) THEN ! default: we use the last time value
	281	indextime = timelen
	282	ELSE ! else we look for the desired value in the time axis
	283	indextime = 0
	284	DO i=1,timelen
	285	IF (abs(time(i) - timestart) .lt. 0.01) THEN
	286	indextime = i
	287	EXIT
	288	ENDIF
	289	ENDDO
	290	IF (indextime .eq. 0) THEN
	291	PRINT*, "Time", timestart," is not in "//fichnom//"!!"
	292	PRINT*, "Stored times are:"
	293	DO i=1,timelen
	294	PRINT*, time(i)
	295	ENDDO
	296	CALL abort
	297	ENDIF
	298	ENDIF
	299
	300	! In start the absolute date is day_ini + hour_ini + time
	301	! For now on, in the GCMdynamics, it is day_ini + time0
	302	time0 = time(indextime) + hour_ini
	303	day_ini = day_ini + INT(time0)
	304	time0 = time0 - INT(time0) ! time0 devient le nouveau hour_ini
	305	hour_ini = time0
	306
	307	PRINT*, "dynetat0: Selected time ",time(indextime),
	308	. " at index ",indextime
	309
	310	DEALLOCATE(time)
[38]	311
[999]	312
	313	! read vcov
	314	corner(1)=1
	315	corner(2)=1
	316	corner(3)=1
	317	corner(4)=indextime
	318	edges(1)=iip1
	319	edges(2)=jjm
	320	edges(3)=llm
	321	edges(4)=1
	322	ierr=nf90_inq_varid(nid,"vcov",nvarid)
[38]	323	IF (ierr .NE. nf90_noerr) THEN
[999]	324	PRINT*, "dynetat0: Le champ <vcov> est absent"
[38]	325	write(,)trim(nf90_strerror(ierr))
	326	CALL abort
	327	ENDIF
[999]	328	!ierr=nf90_get_var(nid,nvarid,vcov)
	329	ierr=nf90_get_var(nid,nvarid,vcov,corner,edges)
[38]	330	IF (ierr .NE. nf90_noerr) THEN
[999]	331	PRINT*, "dynetat0: Lecture echouee pour <vcov>"
[38]	332	write(,)trim(nf90_strerror(ierr))
	333	CALL abort
	334	ENDIF
[999]	335
	336	! read ucov
	337	corner(1)=1
	338	corner(2)=1
	339	corner(3)=1
	340	corner(4)=indextime
	341	edges(1)=iip1
	342	edges(2)=jjp1
	343	edges(3)=llm
	344	edges(4)=1
	345	ierr=nf90_inq_varid(nid,"ucov",nvarid)
[38]	346	IF (ierr .NE. nf90_noerr) THEN
[999]	347	PRINT*, "dynetat0: Le champ <ucov> est absent"
[38]	348	write(,)trim(nf90_strerror(ierr))
	349	CALL abort
	350	ENDIF
[999]	351	ierr=nf90_get_var(nid,nvarid,ucov,corner,edges)
[38]	352	IF (ierr .NE. nf90_noerr) THEN
[999]	353	PRINT*, "dynetat0: Lecture echouee pour <ucov>"
[38]	354	write(,)trim(nf90_strerror(ierr))
	355	CALL abort
	356	ENDIF
	357
[999]	358	! read teta
[38]	359	ierr=nf90_inq_varid(nid,"teta",nvarid)
	360	IF (ierr .NE. nf90_noerr) THEN
	361	PRINT*, "dynetat0: Le champ <teta> est absent"
	362	write(,)trim(nf90_strerror(ierr))
	363	CALL abort
	364	ENDIF
[999]	365	ierr=nf90_get_var(nid,nvarid,teta,corner,edges)
[38]	366	IF (ierr .NE. nf90_noerr) THEN
	367	PRINT*, "dynetat0: Lecture echouee pour <teta>"
	368	write(,)trim(nf90_strerror(ierr))
	369	CALL abort
	370	ENDIF
	371
[999]	372	! read tracers
[1416]	373	IF(nqtot.GE.1) THEN
[38]	374	write(,) 'dynetat0: loading tracers'
[1416]	375	IF(nqtot.GT.99) THEN
[38]	376	PRINT*, "Trop de traceurs"
	377	CALL abort
	378	ENDIF
[1416]	379	nqold=nqtot
	380	DO iq=1,nqtot
[38]	381	! str3(1:1)='q'
	382	! WRITE(str3(2:3),'(i2.2)') iq
	383	! ierr = NF_INQ_VARID (nid, str3, nvarid)
[1130]	384	! NB: tracers are now read in using their name ('tname' from infotrac)
	385	! write(,) " loading tracer:",trim(tname(iq))
	386	ierr=nf90_inq_varid(nid,tname(iq),nvarid)
[38]	387	IF (ierr .NE. nf90_noerr) THEN
	388	! PRINT*, "dynetat0: Le champ <"//str3//"> est absent"
[1130]	389	PRINT*, "dynetat0: Le champ <"//trim(tname(iq))//
[38]	390	& "> est absent"
	391	PRINT*, " Il est donc initialise a zero"
	392	CALL initial0(ijp1llm,q(1,1,1,iq))
	393	nqold=min(iq-1,nqold)
	394	ELSE
	395	!ierr=nf90_get_var(nid,nvarid,q(1,1,1,iq))
[999]	396	ierr=nf90_get_var(nid,nvarid,traceur,corner,edges)
[38]	397	IF (ierr .NE. nf90_noerr) THEN
	398	! PRINT*, "dynetat0: Lecture echouee pour "//str3
[1130]	399	PRINT*,"dynetat0: Lecture echouee pour "//trim(tname(iq))
[38]	400	CALL abort
	401	ENDIF
	402	q(:,:,:,iq)=traceur(:,:,:)
	403	ENDIF
	404	ENDDO
[1416]	405	if ((nqold.lt.nqtot).and.(nqold.ge.1)) then
[38]	406	c case when new tracer are added in addition to old ones
	407	write(,)'tracers 1 to ', nqold,'were already present'
[1416]	408	write(,)'tracers ', nqold+1,' to ', nqtot,'are new'
[38]	409	! yes=' '
	410	! do while ((yes.ne.'y').and.(yes.ne.'n'))
	411	! write(,) 'Would you like to reindex tracer # 1 ->',nqold
[1416]	412	! write(,) 'to #',nqtot-nqold+1,'->', nqtot,' (y or n) ?'
[38]	413	! read(*,fmt='(a)') yes
	414	! end do
	415	! if (yes.eq.'y') then
	416	! write(,) 'OK, let s reindex the tracers'
	417	! do l=1,llm
	418	! do j=1,jjp1
	419	! do i=1,iip1
[1416]	420	! do iq=nqtot,nqtot-nqold+1,-1
	421	! q(i,j,l,iq)=q(i,j,l,iq-nqtot+nqold)
[38]	422	! end do
[1416]	423	! do iq=nqtot-nqold,1,-1
[38]	424	! q(i,j,l,iq)= 0.
	425	! end do
	426	! end do
	427	! end do
	428	! end do
	429	! end if
	430	end if
	431	ENDIF
	432
[999]	433	! read masse
[38]	434	ierr=nf90_inq_varid(nid,"masse",nvarid)
	435	IF (ierr .NE. nf90_noerr) THEN
	436	PRINT*, "dynetat0: Le champ <masse> est absent"
	437	write(,)trim(nf90_strerror(ierr))
	438	CALL abort
	439	ENDIF
[999]	440	ierr=nf90_get_var(nid,nvarid,masse,corner,edges)
[38]	441	IF (ierr .NE. nf90_noerr) THEN
	442	PRINT*, "dynetat0: Lecture echouee pour <masse>"
	443	write(,)trim(nf90_strerror(ierr))
	444	CALL abort
	445	ENDIF
	446
[999]	447	! read ps
	448	corner(1)=1
	449	corner(2)=1
	450	corner(3)=indextime
	451	edges(1)=iip1
	452	edges(2)=jjp1
	453	edges(3)=1
[38]	454	ierr=nf90_inq_varid(nid,"ps",nvarid)
	455	IF (ierr .NE. nf90_noerr) THEN
	456	PRINT*, "dynetat0: Le champ <ps> est absent"
	457	write(,)trim(nf90_strerror(ierr))
	458	CALL abort
	459	ENDIF
[999]	460	ierr=nf90_get_var(nid,nvarid,ps,corner,edges)
[38]	461	IF (ierr .NE. nf90_noerr) THEN
	462	PRINT*, "dynetat0: Lecture echouee pour <ps>"
	463	write(,)trim(nf90_strerror(ierr))
	464	CALL abort
	465	ENDIF
	466
	467	ierr=nf90_close(nid)
	468
[791]	469	! day_ini=day_ini+INT(time) ! obsolete stuff ; 0<time<1 anyways
	470	! time=time-INT(time)
[38]	471
	472	1 FORMAT(//10x,'la valeur de im =',i4,2x,'lue sur le fichier de dem
	473	*arrage est differente de la valeur parametree iim =',i4//)
	474	2 FORMAT(//10x,'la valeur de jm =',i4,2x,'lue sur le fichier de dem
	475	*arrage est differente de la valeur parametree jjm =',i4//)
	476	3 FORMAT(//10x,'la valeur de lmax =',i4,2x,'lue sur le fichier dema
	477	*rrage est differente de la valeur parametree llm =',i4//)
	478	4 FORMAT(//10x,'la valeur de dtrv =',i4,2x,'lue sur le fichier dema
	479	*rrage est differente de la valeur dtinteg =',i4//)
	480
	481	RETURN
	482	END

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

Download in other formats:

Original Format