source: trunk/LMDZ.MARS/libf/phymars/phyetat0_mod.F90 @ 3249

Last change on this file since 3249 was 3183, checked in by jbclement, 2 years ago

Mars PCM:
Some changes to prepare the introduction of slopes in 1D: in particular, the subroutine "getslopes.F90" and "param_slope.F90" are now inside the module "slope_mod.F90" + Few small cleanings throughout the code.
JBC
File size: 29.3 KB
Line 
1module phyetat0_mod
2
3implicit none
4  real,save :: tab_cntrl_mod(100)
5
6!$OMP THREADPRIVATE(tab_cntrl_mod)
7
8contains
9
10subroutine phyetat0 (fichnom,tab0,Lmodif,nsoil,ngrid,nlay,nq,nqsoil, &
11                     day_ini,time0,tsurf,tsoil,albedo,emis,q2,qsurf,qsoil, &
12                     tauscaling,totcloudfrac,wstar,watercap,perennial_co2ice, &
13                     def_slope,def_slope_mean,subslope_dist)
14
15  use tracer_mod, only: noms ! tracer names
16  use surfdat_h, only: phisfi, albedodat, z0, z0_default,&
17                       zmea, zstd, zsig, zgam, zthe, hmons, summit, base,&
18                       watercaptag
19  use iostart, only: nid_start, open_startphy, close_startphy, &
20                     get_field, get_var, inquire_field, &
21                     inquire_dimension, inquire_dimension_length
22  use nonoro_gwd_ran_mod, only: du_nonoro_gwd, dv_nonoro_gwd
23  use nonoro_gwd_mix_mod, only: du_eddymix_gwd, dv_eddymix_gwd, de_eddymix_rto, &
24                                df_eddymix_flx !dr_depflux_gwd
25
26  use compute_dtau_mod, only: dtau
27  use dust_rad_adjust_mod, only: dust_rad_adjust_prev,dust_rad_adjust_next
28  use dust_param_mod, only: dustscaling_mode
29  USE ioipsl_getin_p_mod, ONLY : getin_p
30  use comsoil_h, only: flux_geo
31  USE comslope_mod, ONLY: nslope, major_slope
32  USE paleoclimate_mod, ONLY: paleoclimate, h2o_ice_depth, lag_co2_ice, d_coef
33  USE comcstfi_h, only: pi
34  use geometry_mod, only: latitude
35  implicit none
36 
37  include "callkeys.h"
38!======================================================================
39! Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
40!  Adaptation � Mars : Yann Wanherdrick
41! Objet: Lecture de l etat initial pour la physique
42! Modifs: Aug.2010 EM : use NetCDF90 to load variables (enables using
43!                      r4 or r8 restarts independently of having compiled
44!                      the GCM in r4 or r8)
45!         June 2013 TN : Possibility to read files with a time axis
46!         November 2013 EM : Enabeling parallel, using iostart module
47!         March 2020 AD: Enabling initialization of physics without startfi
48!                        flag: startphy_file
49!======================================================================
50  INTEGER nbsrf !Mars nbsrf a 1 au lieu de 4
51  PARAMETER (nbsrf=1) ! nombre de sous-fractions pour une maille
52!======================================================================
53!  Arguments:
54!  ---------
55!  inputs:
56!  logical,intent(in) :: startphy_file ! .true. if reading start file
57  character*(*),intent(in) :: fichnom ! "startfi.nc" file
58  integer,intent(in) :: tab0
59  integer,intent(in) :: Lmodif
60  integer,intent(in) :: nsoil ! # of soil layers
61  integer,intent(in) :: ngrid ! # of atmospheric columns
62  integer,intent(in) :: nlay ! # of atmospheric layers
63  integer,intent(in) :: nq
64  integer,intent(in) :: nqsoil ! # of tracers in the soil
65  integer :: day_ini
66  real :: time0
67
68!  outputs:
69  real,intent(out) :: tsurf(ngrid,nslope) ! surface temperature
70  real,intent(out) :: tsoil(ngrid,nsoil,nslope) ! soil temperature
71  real,intent(out) :: albedo(ngrid,2,nslope) ! surface albedo
72  real,intent(out) :: emis(ngrid,nslope) ! surface emissivity
73  real,intent(out) :: q2(ngrid,nlay+1) !
74  real,intent(out) :: qsurf(ngrid,nq,nslope) ! tracers on surface
75  real,intent(out) :: qsoil(ngrid,nsoil,nqsoil,nslope) ! tracers in the subsurface
76  real,intent(out) :: tauscaling(ngrid) ! dust conversion factor
77  real,intent(out) :: totcloudfrac(ngrid) ! total cloud fraction
78  real,intent(out) :: wstar(ngrid) ! Max vertical velocity in thermals (m/s)
79  real,intent(out) :: watercap(ngrid,nslope) ! h2o_ice_cover
80  real,intent(out) :: perennial_co2ice(ngrid,nslope) ! perennial co2 ice(kg/m^2)
81  real,intent(out) :: def_slope(nslope+1) !boundaries for bining of the slopes
82  real,intent(out) :: def_slope_mean(nslope)
83  real,intent(out) :: subslope_dist(ngrid,nslope) !undermesh statistics
84!======================================================================
85!  Local variables:
86
87      real surffield(ngrid) ! to temporarily store a surface field
88      real xmin,xmax ! to display min and max of a field
89!
90      INTEGER ig,iq,lmax,islope
91      INTEGER nid, nvarid
92      INTEGER ierr, i, nsrf
93!      integer isoil
94!      INTEGER length
95!      PARAMETER (length=100)
96      CHARACTER*7 str7
97      CHARACTER*2 str2
98      CHARACTER*1 yes
99!
100      REAL p_rad,p_omeg,p_g,p_mugaz,p_daysec
101      INTEGER nqold
102
103! flag which identifies if 'startfi.nc' file is using old names (qsurf01,...)
104      logical :: oldtracernames=.false.
105      integer :: count
106      character(len=30) :: txt ! to store some text
107
108! specific for time
109      REAL,ALLOCATABLE :: time(:) ! times stored in start
110      INTEGER timelen ! number of times stored in the file
111      INTEGER indextime ! index of selected time
112
113      INTEGER :: edges(3),corner(3)
114      LOGICAL :: found
115
116      REAL :: timestart ! to pick which initial state to start from
117      REAL :: surfemis  ! constant emissivity when no startfi
118      REAL :: surfalbedo  ! constant albedo when no startfi
119
120      REAL :: watercaptag_tmp(ngrid)
121
122! Sub-grid scale slopes
123      LOGICAL :: startphy_slope !to be retrocompatible and add the nslope dimension
124      REAL, ALLOCATABLE :: default_def_slope(:)
125      REAL :: sum_dist
126      REAL :: current_max !var to find max distrib slope
127! Variables for CO2 index
128      INTEGER :: igcm_co2_tmp
129
130      CHARACTER(len=5) :: modname="phyetat0"
131
132write(*,*) "phyetat0: startphy_file", startphy_file
133
134if (startphy_file) then
135   ! open physics initial state file:
136   call open_startphy(fichnom)
137   ! possibility to modify tab_cntrl in tabfi
138   write(*,*)
139   write(*,*) 'TABFI in phyeta0: Lmodif=',Lmodif," tab0=",tab0
140   call tabfi (nid_start,Lmodif,tab0,day_ini,lmax,p_rad, &
141               p_omeg,p_g,p_mugaz,p_daysec,time0)
142else ! "academic" initialization of planetary parameters via tabfi
143   call tabfi (0,0,0,day_ini,lmax,p_rad, &
144               p_omeg,p_g,p_mugaz,p_daysec,time0)
145endif ! of if (startphy_file)
146
147if (startphy_file) then
148   call get_var("def_slope",def_slope,found)
149   if (.not. found) then
150       startphy_slope=.false.
151       write(*,*)'slope_settings: Problem while reading <def_slope>'
152       write(*,*)'Checking that nslope=1'
153       if (nslope == 1) then
154           write(*,*)'We take default def_slope and subslope dist'
155           allocate(default_def_slope(nslope + 1))
156           default_def_slope(1) = -90.
157           default_def_slope(2) = 90.
158           subslope_dist = 1.
159           def_slope = default_def_slope
160       else
161           write(*,*)'Variable def_slope is not present in the start and'
162           write(*,*)'you are trying to run with nslope!=1'
163           write(*,*)'This is not possible, you have to go through newstart'
164       endif
165   else
166       startphy_slope=.true.
167       call get_field("subslope_dist",subslope_dist,found,indextime)
168       if (.not. found) then
169           write(*,*)'slope_settings: Problem while reading <subslope_dist>'
170           write(*,*)'We have to abort.'
171           write(*,*)'Please check that nslope is coherent, you can modify this parameter with newstart'
172           call abort_physic(modname,"phyetat0: Failed loading <subslope_dist>",1)
173       endif
174   endif
175else ! (no startphy_file)
176    if (nslope == 1) then
177        allocate(default_def_slope(2))
178        default_def_slope = 0.
179        subslope_dist = 1.
180        def_slope = default_def_slope
181    else
182        write(*,*)'Without startfi, lets first run with nslope=1'
183        call abort_physic(modname,"phyetat0: No startfi and nslope!=1",1)
184    endif
185endif   
186
187do islope = 1,nslope
188    def_slope_mean(islope) = (def_slope(islope) + def_slope(islope + 1))/2.
189enddo
190 
191DO ig = 1,ngrid
192    sum_dist = 0.
193    DO islope = 1,nslope
194        sum_dist = sum_dist + subslope_dist(ig,islope)
195    ENDDO
196    DO islope = 1,nslope
197        subslope_dist(ig,islope) = subslope_dist(ig,islope)/sum_dist
198    ENDDO
199ENDDO
200
201!Now determine the major subslope, ie. the maximal distribution
202
203DO ig=1,ngrid
204  major_slope(ig)=1
205  current_max=subslope_dist(ig,1)
206  DO islope=2,nslope
207    if(subslope_dist(ig,islope).GT.current_max) then
208      major_slope(ig)=islope
209      current_max=subslope_dist(ig,islope)
210    ENDIF
211  ENDDO
212ENDDO
213
214if (startphy_file) then
215   ! Load surface geopotential:
216   call get_field("phisfi",phisfi,found)
217   if (.not.found) then
218     call abort_physic(modname, &
219                "phyetat0: Failed loading <phisfi>",1)
220   endif
221else
222  phisfi(:)=0.
223endif ! of if (startphy_file)
224write(*,*) "phyetat0: surface geopotential <phisfi> range:", &
225               minval(phisfi), maxval(phisfi)
226
227
228if (startphy_file) then
229   ! Load bare ground albedo:
230   call get_field("albedodat",albedodat,found)
231   if (.not.found) then
232     call abort_physic(modname, &
233                "phyetat0: Failed loading <albedodat>",1)
234   endif
235else ! If no startfi file, use parameter surfalbedo in def file
236  surfalbedo=0.1
237  call getin_p("surfalbedo_without_startfi",surfalbedo)
238  print*,"surfalbedo_without_startfi",surfalbedo
239  albedodat(:)=surfalbedo
240endif ! of if (startphy_file)
241write(*,*) "phyetat0: Bare ground albedo <albedodat> range:", &
242             minval(albedodat), maxval(albedodat)
243
244! ZMEA
245if (startphy_file) then
246   call get_field("ZMEA",zmea,found)
247   if (.not.found) then
248     call abort_physic(modname, &
249                "phyetat0: Failed loading <ZMEA>",1)
250   endif
251else
252  zmea(:)=0.
253endif ! of if (startphy_file)
254write(*,*) "phyetat0: <ZMEA> range:", &
255            minval(zmea), maxval(zmea)
256
257
258! ZSTD
259if (startphy_file) then
260   call get_field("ZSTD",zstd,found)
261   if (.not.found) then
262     call abort_physic(modname, &
263                "phyetat0: Failed loading <ZSTD>",1)
264   endif
265else
266  zstd(:)=0.
267endif ! of if (startphy_file)
268write(*,*) "phyetat0: <ZSTD> range:", &
269            minval(zstd), maxval(zstd)
270
271
272! ZSIG
273if (startphy_file) then
274   call get_field("ZSIG",zsig,found)
275   if (.not.found) then
276     call abort_physic(modname, &
277                "phyetat0: Failed loading <ZSIG>",1)
278   endif
279else
280  zsig(:)=0.
281endif ! of if (startphy_file)
282write(*,*) "phyetat0: <ZSIG> range:", &
283            minval(zsig), maxval(zsig)
284
285
286! ZGAM
287if (startphy_file) then
288   call get_field("ZGAM",zgam,found)
289   if (.not.found) then
290     call abort_physic(modname, &
291                "phyetat0: Failed loading <ZGAM>",1)
292   endif
293else
294  zgam(:)=0.
295endif ! of if (startphy_file)
296write(*,*) "phyetat0: <ZGAM> range:", &
297            minval(zgam), maxval(zgam)
298
299
300! ZTHE
301if (startphy_file) then
302   call get_field("ZTHE",zthe,found)
303   if (.not.found) then
304     call abort_physic(modname, &
305                "phyetat0: Failed loading <ZTHE>",1)
306   endif
307else
308  zthe(:)=0.
309endif ! of if (startphy_file)
310write(*,*) "phyetat0: <ZTHE> range:", &
311             minval(zthe), maxval(zthe)
312
313! hmons
314if (startphy_file) then
315   call get_field("hmons",hmons,found)
316   if (.not.found) then
317     write(*,*) "WARNING: phyetat0: Failed loading <hmons>"
318     if (rdstorm) then
319     call abort_physic(modname, &
320                "phyetat0: Failed loading <hmons>",1)
321     else
322       write(*,*) "will continue anyway..."
323       write(*,*) "because you may not need it."
324       hmons(:)=0.
325     end if ! if (rdstorm)
326   else
327     do ig=1,ngrid
328       if (hmons(ig) .eq. -999999.)  hmons(ig)=0.
329     enddo
330   endif ! (.not.found)
331else
332   hmons(:)=0.
333endif ! if (startphy_file)
334write(*,*) "phyetat0: <hmons> range:", &
335            minval(hmons), maxval(hmons)
336
337
338! summit
339if (startphy_file) then
340   call get_field("summit",summit,found)
341   if (.not.found) then
342     write(*,*) "WARNING: phyetat0: Failed loading <summit>"
343     if (rdstorm) then
344     call abort_physic(modname, &
345                "phyetat0: Failed loading <summit>",1)
346     else
347       write(*,*) "will continue anyway..."
348       write(*,*) "because you may not need it."
349       summit(:)=0.
350     end if
351   else
352     do ig=1,ngrid
353       if (summit(ig) .eq. -999999.)  summit(ig)=0.
354     enddo
355   endif ! if (.not.found)
356else
357   summit(:)=0. 
358endif ! if (startphy_file)
359write(*,*) "phyetat0: <summit> range:", &
360            minval(summit), maxval(summit)
361
362
363! base
364if (startphy_file) then
365   call get_field("base",base,found)
366   if (.not.found) then
367     write(*,*) "WARNING: phyetat0: Failed loading <base>"
368     if (rdstorm) then
369     call abort_physic(modname, &
370                "phyetat0: Failed loading <base>",1)
371     else
372       write(*,*) "will continue anyway..."
373       write(*,*) "because you may not need it."
374       base(:)=0.
375     end if
376   else
377     do ig=1,ngrid
378       if (base(ig) .eq. -999999.)  base(ig)=0.
379     enddo
380   endif ! if(.not.found)
381else
382   base(:)=0.
383endif ! if (startphy_file)
384write(*,*) "phyetat0: <base> range:", &
385            minval(base), maxval(base)
386
387! Time axis
388! obtain timestart from run.def
389timestart=-9999 ! default value
390call getin_p("timestart",timestart)
391if (startphy_file) then
392   found=inquire_dimension("Time")
393   if (.not.found) then
394     indextime = 1
395     write(*,*) "phyetat0: No time axis found in "//trim(fichnom)
396   else
397     write(*,*) "phyetat0: Time axis found in "//trim(fichnom)
398     timelen=inquire_dimension_length("Time")
399     allocate(time(timelen))
400     ! load "Time" array:
401     call get_var("Time",time,found)
402     if (.not.found) then
403     call abort_physic(modname, &
404                "phyetat0: Failed loading <Time>",1)
405     endif
406     ! seclect the desired time index
407     IF (timestart .lt. 0) THEN  ! default: we use the last time value
408       indextime = timelen
409     ELSE  ! else we look for the desired value in the time axis
410       indextime = 0
411       DO i=1,timelen
412         IF (abs(time(i) - timestart) .lt. 0.01) THEN
413           indextime = i
414           EXIT
415         ENDIF
416       ENDDO
417       IF (indextime .eq. 0) THEN
418         PRINT*, "Time", timestart," is not in "//trim(fichnom)//"!!"
419         PRINT*, "Stored times are:"
420         DO i=1,timelen
421            PRINT*, time(i)
422         ENDDO
423         call abort_physic(modname,"phyetat0: Time error",1)
424       ENDIF
425     ENDIF ! of IF (timestart .lt. 0)
426     ! In startfi the absolute date is day_ini + time0 + time
427     ! For now on, in the GCM physics, it is day_ini + time0
428     time0 = time(indextime) + time0
429     day_ini = day_ini + INT(time0)
430     time0 = time0 - INT(time0)   
431     PRINT*, "phyetat0: Selected time ",time(indextime), &
432             " at index ",indextime
433     DEALLOCATE(time)
434   endif ! of if Time not found in file
435
436   call ini_tab_controle_dyn_xios(day_ini)
437
438else
439  indextime = 1
440endif ! if (startphy_file)
441
442! Dust conversion factor
443if (startphy_file) then
444   call get_field("tauscaling",tauscaling,found,indextime)
445   if (.not.found) then
446     write(*,*) "phyetat0: <tauscaling> not in file"
447     tauscaling(:) = 1
448   endif
449else
450   tauscaling(:) = 1
451endif ! if (startphy_file)
452write(*,*) "phyetat0: dust conversion factor <tauscaling> range:", &
453            minval(tauscaling), maxval(tauscaling)
454
455! dust_rad_adjust_* for radiative rescaling of dust
456if (dustscaling_mode==2) then
457 if (startphy_file) then
458   call get_field("dust_rad_adjust_prev",dust_rad_adjust_prev,found,indextime)
459   if (.not.found) then
460     write(*,*) "phyetat0: <dust_rad_adjust_prev> not in file; set to 1"
461     dust_rad_adjust_prev(:) = 1
462   endif
463   call get_field("dust_rad_adjust_next",dust_rad_adjust_next,found,indextime)
464   if (.not.found) then
465     write(*,*) "phyetat0: <dust_rad_adjust_next> not in file; set to 1"
466     dust_rad_adjust_next(:) = 1
467   endif
468 else
469   dust_rad_adjust_prev(:)= 0
470   dust_rad_adjust_next(:)= 0
471 endif ! if (startphy_file)
472 write(*,*) "phyetat0: radiative scaling coeff <dust_rad_adjust_prev> range:", &
473            minval(dust_rad_adjust_prev), maxval(dust_rad_adjust_prev)
474 write(*,*) "phyetat0: radiative scaling coeff <dust_rad_adjust_next> range:", &
475            minval(dust_rad_adjust_next), maxval(dust_rad_adjust_next)
476endif ! of if (dustscaling_mode==2)
477
478! dtau: opacity difference between GCM and dust scenario
479if (startphy_file) then
480   call get_field("dtau",dtau,found,indextime)
481   if (.not.found) then
482     write(*,*) "phyetat0: <dtau> not in file; set to zero"
483     dtau(:) = 0
484   endif
485else
486   dtau(:)= 0
487endif ! if (startphy_file)
488write(*,*) "phyetat0: opacity diff wrt scenario <dtau> range:", &
489            minval(dtau), maxval(dtau)
490
491
492! Sub-grid cloud fraction
493if (startphy_file) then
494   call get_field("totcloudfrac",totcloudfrac,found,indextime)
495   if (.not.found) then
496     write(*,*) "phyetat0: <totcloudfrac> not in file WARNING put to 1"
497     totcloudfrac(:) = 1.0 !valeur par defaut (CLFfixval par defaut)
498   endif
499else
500   totcloudfrac(:)=1.0
501endif ! if (startphy_file)
502write(*,*) "phyetat0: total cloud fraction <totcloudfrac> range:", &
503            minval(totcloudfrac), maxval(totcloudfrac)
504
505
506! Max vertical velocity in thermals
507if (startphy_file) then
508   call get_field("wstar",wstar,found,indextime)
509   if (.not.found) then
510     write(*,*) "phyetat0: <wstar> not in file! Set to zero"
511     wstar(:)=0
512   endif
513else
514   wstar(:)=0
515endif ! if (startphy_file)
516write(*,*) "phyetat0: Max vertical velocity in thermals <wstar> range:", &
517            minval(wstar),maxval(wstar)
518
519
520! Surface temperature :
521if (startphy_file) then !tsurf
522   call get_field("tsurf",tsurf,found,indextime)
523   if (.not.found) then
524     call abort_physic(modname, &
525                "phyetat0: Failed loading <tsurf>",1)
526   endif
527else
528  tsurf(:,:)=0. ! will be updated afterwards in physiq !
529endif ! of if (startphy_file)
530write(*,*) "phyetat0: Surface temperature <tsurf> range:", &
531            minval(tsurf), maxval(tsurf)
532
533! Surface albedo
534if (startphy_file) then
535   call get_field("albedo",albedo(:,1,:),found,indextime)
536   if (.not.found) then
537     write(*,*) "phyetat0: Failed loading <albedo>"
538     do islope=1,nslope
539       albedo(:,1,islope)=albedodat(:)
540     enddo
541   endif
542else
543   do islope=1,nslope
544     albedo(:,1,islope)=albedodat(:)
545   enddo
546endif ! of if (startphy_file)
547write(*,*) "phyetat0: Surface albedo <albedo> range:", &
548            minval(albedo(:,1,:)), maxval(albedo(:,1,:))
549albedo(:,2,:)=albedo(:,1,:)
550
551! Surface emissivity
552if (startphy_file) then
553   call get_field("emis",emis,found,indextime)
554   if (.not.found) then
555     call abort_physic(modname, &
556                "phyetat0: Failed loading <emis>",1)
557   endif
558else
559  ! If no startfi file, use parameter surfemis in def file
560  surfemis=0.95
561  call getin_p("surfemis_without_startfi",surfemis)
562  print*,"surfemis_without_startfi",surfemis
563  emis(:,:)=surfemis
564endif ! of if (startphy_file)
565write(*,*) "phyetat0: Surface emissivity <emis> range:", &
566            minval(emis), maxval(emis)
567
568
569! surface roughness length (NB: z0 is a common in surfdat_h)
570if (startphy_file) then
571   call get_field("z0",z0,found)
572   if (.not.found) then
573     write(*,*) "phyetat0: Failed loading <z0>"
574     write(*,*) 'will use constant value of z0_default:',z0_default
575     z0(:)=z0_default
576   endif
577else
578   z0(:)=z0_default
579endif ! of if (startphy_file)
580write(*,*) "phyetat0: Surface roughness <z0> range:", &
581            minval(z0), maxval(z0)
582
583
584! pbl wind variance
585if (startphy_file) then
586   call get_field("q2",q2,found,indextime)
587   if (.not.found) then
588     call abort_physic(modname, &
589                "phyetat0: Failed loading <q2>",1)
590   endif
591else
592  q2(:,:)=0.
593endif ! of if (startphy_file)
594write(*,*) "phyetat0: PBL wind variance <q2> range:", &
595            minval(q2), maxval(q2)
596
597! Non-orographic gravity waves
598if (startphy_file) then
599   call get_field("du_nonoro_gwd",du_nonoro_gwd,found,indextime)
600   if (.not.found) then
601      write(*,*) "phyetat0: <du_nonoro_gwd> not in file"
602      du_nonoro_gwd(:,:)=0.
603   endif
604else
605du_nonoro_gwd(:,:)=0.
606endif ! if (startphy_file)
607write(*,*) "phyetat0: Memory of zonal wind tendency due to non-orographic GW"
608write(*,*) " <du_nonoro_gwd> range:", &
609             minval(du_nonoro_gwd), maxval(du_nonoro_gwd)
610
611if (startphy_file) then
612   call get_field("dv_nonoro_gwd",dv_nonoro_gwd,found,indextime)
613   if (.not.found) then
614      write(*,*) "phyetat0: <dv_nonoro_gwd> not in file"
615      dv_nonoro_gwd(:,:)=0.
616   endif
617else ! ! if (startphy_file)
618dv_nonoro_gwd(:,:)=0.
619endif ! if (startphy_file)
620write(*,*) "phyetat0: Memory of meridional wind tendency due to non-orographic GW"
621write(*,*) " <dv_nonoro_gwd> range:", &
622             minval(dv_nonoro_gwd), maxval(dv_nonoro_gwd)
623if (startphy_file) then
624   call get_field("du_eddymix_gwd",du_eddymix_gwd,found,indextime)
625   if (.not.found) then
626      write(*,*) "phyetat0: <du_eddymix_gwd> not in file"
627      du_eddymix_gwd(:,:)=0.
628   endif
629else
630du_eddymix_gwd(:,:)=0.
631endif ! if (startphy_file)
632write(*,*) "phyetat0: Memory of zonal wind tendency due to non-orographic GW mixing"
633write(*,*) " <du_eddymix_gwd> range:", &
634             minval(du_eddymix_gwd), maxval(du_eddymix_gwd)
635 
636
637if (startphy_file) then
638   call get_field("dv_nonoro_gwd",dv_nonoro_gwd,found,indextime)
639   if (.not.found) then
640      write(*,*) "phyetat0: <dv_nonoro_gwd> not in file"
641      dv_nonoro_gwd(:,:)=0.
642   endif
643else ! ! if (startphy_file)
644dv_nonoro_gwd(:,:)=0.
645endif ! if (startphy_file)
646write(*,*) "phyetat0: Memory of meridional wind tendency due to non-orographic GW"
647write(*,*) " <dv_nonoro_gwd> range:", &
648             minval(dv_nonoro_gwd), maxval(dv_nonoro_gwd)
649
650if (startphy_file) then
651   call get_field("dv_eddymix_gwd",dv_eddymix_gwd,found,indextime)
652   if (.not.found) then
653      write(*,*) "phyetat0: <dv_eddymix_gwd> not in file"
654      dv_eddymix_gwd(:,:)=0.
655   endif
656else ! ! if (startphy_file)
657dv_eddymix_gwd(:,:)=0.
658endif ! if (startphy_file)
659write(*,*) "phyetat0: Memory of meridional wind tendency due to non-orographic GW mixing"
660write(*,*) " <dv_eddymix_gwd> range:", &
661             minval(dv_eddymix_gwd), maxval(dv_eddymix_gwd)
662
663!if (startphy_file) then
664!   call get_field("dr_depflux_gwd",dr_depflux_gwd,found,indextime)
665!   if (.not.found) then
666!      write(*,*) "phyetat0: <dr_depflux_gwd> not in file"
667!      dr_depflux_gwd(:,:,:)=0.
668!   endif
669!else ! ! if (startphy_file)
670!!dr_depflux_gwd(:,:,:)=0.
671!endif ! if (startphy_file)
672!write(*,*) "phyetat0: Memory of meridional wind tendency due to non-orographic GW mixing"
673!write(*,*) " <dr_depflux_gwd> range:", &
674!             minval(dr_depflux_gwd), maxval(dr_depflux_gwd)
675
676if (startphy_file) then
677   call get_field("de_eddymix_rto",de_eddymix_rto,found,indextime)
678   if (.not.found) then
679      write(*,*) "phyetat0: <de_eddymix_rto> not in file"
680      de_eddymix_rto(:,:)=0.
681   endif
682else ! ! if (startphy_file)
683de_eddymix_rto(:,:)=0.
684endif ! if (startphy_file)
685write(*,*) "phyetat0: Memory of meridional wind tendency due to non-orographic GW mixing"
686write(*,*) " <de_eddymix_rto> range:", &
687             minval(de_eddymix_rto), maxval(de_eddymix_rto)
688
689if (startphy_file) then
690   call get_field("df_eddymix_flx ",df_eddymix_flx ,found,indextime)
691   if (.not.found) then
692      write(*,*) "phyetat0: <df_eddymix_flx > not in file"
693      df_eddymix_flx (:,:)=0.
694   endif
695else ! ! if (startphy_file)
696df_eddymix_flx (:,:)=0.
697endif ! if (startphy_file)
698write(*,*) "phyetat0: Memory of meridional wind tendency due to non-orographic GW mixing"
699write(*,*) " <df_eddymix_flx > range:", &
700             minval(df_eddymix_flx ), maxval(df_eddymix_flx )
701
702
703
704! tracer on surface
705if (nq.ge.1) then
706  do iq=1,nq
707    txt=noms(iq)
708    if (txt.eq."h2o_vap") then
709      ! There is no surface tracer for h2o_vap;
710      ! "h2o_ice" should be loaded instead
711      txt="h2o_ice"
712      write(*,*) 'phyetat0: loading surface tracer', &
713                           ' h2o_ice instead of h2o_vap'
714      write(*,*) 'iq = ', iq
715    endif
716
717    if (hdo) then
718      if (txt.eq."hdo_vap") then
719        txt="hdo_ice"
720        write(*,*) 'phyetat0: loading surface tracer', &
721                           ' hdo_ice instead of hdo_vap'
722      endif
723    endif !hdo
724
725    if (startphy_file) then
726      if (txt.eq."co2") then
727        ! We first check if co2ice exist in the startfi file (old way)
728        ! CO2 ice cover
729        call get_field("co2ice",qsurf(:,iq,:),found,indextime)
730        ! If not, we load the corresponding tracer in qsurf
731        if (.not.found) then
732          call get_field(txt,qsurf(:,iq,:),found,indextime)
733          if (.not.found) then
734            call abort_physic(modname, &
735                "phyetat0: Failed loading co2ice. There is neither the variable co2ice nor qsurf!",1)
736          endif
737        endif
738      else ! (not the co2 tracer)
739        call get_field(txt,qsurf(:,iq,:),found,indextime)
740        if (.not.found) then
741          write(*,*) "phyetat0: Failed loading <",trim(txt),">"
742          write(*,*) "         ",trim(txt)," is set to zero"
743          qsurf(:,iq,:)=0.
744        endif
745      endif !endif co2
746    else !(not startphy_file)
747      qsurf(:,iq,:)=0. ! co2ice is set to 0
748    endif ! of if (startphy_file)
749    write(*,*) "phyetat0: Surface tracer <",trim(txt),"> range:", &
750                 minval(qsurf(:,iq,:)), maxval(qsurf(:,iq,:))
751  enddo ! of do iq=1,nq
752endif ! of if (nq.ge.1)
753
754if (startphy_file) then
755   call get_field("watercap",watercap,found,indextime)
756   if (.not.found) then
757     write(*,*) "phyetat0: Failed loading <watercap> : ", &
758                          "<watercap> is set to zero"
759     watercap(:,:)=0
760     write(*,*) 'Now transfer negative surface water ice to', &
761                ' watercap !'
762     if (nq.ge.1) then
763      do iq=1,nq
764       txt=noms(iq)
765       if (txt.eq."h2o_ice") then
766         do ig=1,ngrid
767          do islope=1,nslope
768            if (qsurf(ig,iq,islope).lt.0.0) then
769              watercap(ig,islope) = qsurf(ig,iq,islope)
770              qsurf(ig,iq,islope) = 0.0
771            end if
772          enddo
773         end do
774       endif
775
776       if (txt.eq."hdo_ice") then
777         do ig=1,ngrid
778          do islope=1,nslope
779            if (qsurf(ig,iq,islope).lt.0.0) then
780              qsurf(ig,iq,islope) = 0.0
781            end if
782          enddo
783         end do
784       endif
785
786      enddo
787     endif ! of if (nq.ge.1)
788   endif ! of if (.not.found)
789else
790   watercap(:,:)=0
791endif ! of if (startphy_file)
792write(*,*) "phyetat0: Surface water ice <watercap> range:", &
793            minval(watercap), maxval(watercap)
794
795if (startphy_file) then
796  ! Call to soil_settings, in order to read soil temperatures,
797  ! as well as thermal inertia and volumetric heat capacity
798  call soil_settings(nid_start,ngrid,nsoil,nqsoil,tsurf,tsoil,qsoil,indextime)
799else
800    flux_geo(:,:) = 0.
801endif ! of if (startphy_file)
802
803if (startphy_file) then
804   call get_field("watercaptag",watercaptag_tmp,found,indextime)
805   if (.not.found) then
806     write(*,*) "phyetat0: Failed loading <watercaptag> : ", &
807                          "<watercaptag> is set as defined by icelocationmode in surfini.F"
808     watercaptag(:)=.false.
809   else
810     do ig=1,ngrid
811       if(watercaptag_tmp(ig).lt.0.5) then
812          watercaptag(ig)=.false.
813       else
814          watercaptag(ig)=.true.
815       endif
816     enddo
817   endif
818endif !startphy_file
819
820if (paleoclimate) then
821  do iq=1,nq
822   txt=noms(iq)
823   if (txt.eq."co2") igcm_co2_tmp = iq
824  enddo
825
826  if (startphy_file) then
827! Depth of H2O lag
828   call get_field("h2o_ice_depth",h2o_ice_depth,found,indextime)
829   if (.not.found) then
830     write(*,*) "phyetat0: Failed loading <h2o_ice_depth> : ", &
831                          "<h2o_ice_depth> is set as -1 (no subsurface ice)"
832     h2o_ice_depth(:,:) = -1.
833   endif
834
835! Diffusion coeficent
836   call get_field("d_coef",d_coef,found,indextime)
837   if (.not.found) then
838     write(*,*) "phyetat0: Failed loading <d_coef> : ", &
839                          "<d_coef> is set as 4e-4 (defualt value)"
840     d_coef(:,:) = 4e-4
841   endif
842
843! Depth of CO2 lag
844   call get_field("lag_co2_ice",lag_co2_ice,found,indextime)
845   if (.not.found) then
846     write(*,*) "phyetat0: Failed loading <lag_co2_ice> : ", &
847                          "<lag_co2_ice> is set as -1 (no subsurface ice)"
848     lag_co2_ice(:,:) = -1.
849   endif
850
851! Perennial CO2 ice
852    perennial_co2ice = 0.
853    call get_field("perennial_co2ice",perennial_co2ice,found,indextime)
854    if (.not. found) then
855        write(*,*) "phyetat0: Failed loading <perennial_co2ice> : ", &
856                   "<perennial_co2ice> is set as 10m at the South Pole"
857        if (abs(latitude(ngrid) - (-pi/2.)) < 1.e-5) then
858            perennial_co2ice(ngrid,:) = 10*1.6e3 ! 10m which is convert to kg/m^2
859            qsurf(ngrid,igcm_co2_tmp,:) = qsurf(ngrid - 1,igcm_co2_tmp,:) + perennial_co2ice(ngrid,:) ! perennial ice + frost
860        endif
861    endif ! not found
862  else ! no startfiphyle
863    h2o_ice_depth = -1.
864    lag_co2_ice = -1.
865    d_coef = 4.e-4
866    perennial_co2ice = 0.
867    if (abs(latitude(ngrid) - (-pi/2.)) < 1.e-5) then
868        perennial_co2ice(ngrid,:) = 10*1.6e3 ! 10m which is convert to kg/m^2
869        qsurf(ngrid,igcm_co2_tmp,:) = qsurf(ngrid - 1,igcm_co2_tmp,:) + perennial_co2ice(ngrid,:) ! perennial ice + frost
870    endif
871  endif !startphy_file
872else
873   h2o_ice_depth = -1.
874   lag_co2_ice = -1.
875   d_coef = 4.e-4
876   perennial_co2ice = 0.
877endif !paleoclimate
878
879! close file:
880if (startphy_file) call close_startphy
881
882end subroutine phyetat0
883
884
885subroutine ini_tab_controle_dyn_xios(idayref)
886
887      USE comcstfi_h, only: g, mugaz, omeg, rad, rcp
888      USE time_phylmdz_mod, ONLY: hour_ini, daysec, dtphys
889      USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat, nbp_lev
890      IMPLICIT NONE
891
892
893      INTEGER*4,intent(in) :: idayref  ! date (initial date for this run)
894
895
896      INTEGER length,l
897      parameter (length = 100)
898      REAL tab_cntrl(length) ! run parameters are stored in this array
899
900      DO l=1,length
901         tab_cntrl(l)=0.
902      ENDDO
903      tab_cntrl(1)  = real(nbp_lon)
904      tab_cntrl(2)  = real(nbp_lat-1)
905      tab_cntrl(3)  = real(nbp_lev)
906      tab_cntrl(4)  = real(idayref)
907      tab_cntrl(5)  = rad
908      tab_cntrl(6)  = omeg
909      tab_cntrl(7)  = g
910      tab_cntrl(8)  = mugaz
911      tab_cntrl(9)  = rcp
912      tab_cntrl(10) = daysec
913      tab_cntrl(11) = dtphys
914
915      tab_cntrl(27) = hour_ini
916
917      tab_cntrl_mod = tab_cntrl
918
919end subroutine ini_tab_controle_dyn_xios
920
921end module phyetat0_mod
Note: See TracBrowser for help on using the repository browser.